Discrete element modelling of rock cutting processes interaction with evaluation of tool wear

Discrete Element Modelling of Rock Cutting Processes Interaction with Evaluation of Tool Wear J. Rojek E. Oñate F. Zárate C.A. Recarey Monograph CIMNE Nº-87, October 2003 Discrete Element Modelling of Rock Cutting Processes Interaction with Evaluation of Tool Wear J. Rojek1, E. Oñate2, F. Zárate2 and C.A. Recarey2 1 Institute of Fundamental Technological Research (IPPT) Polish Academy of Sciences Swietokrzyska 21, 00-049 Warszawa, Poland 2 International Center for Numerical Methods in Engineering (CIMNE) Universidad Politécnica de Cataluña Campus Norte UPC, 08034 Barcelona, Spain Monograph CIMNE Nº-87, October 2003 INTERNATIONAL CENTER FOR NUMERICAL METHODS IN ENGINEERING Gran Capitán s/n, 08034 Barcelona, Spain INTERNACIONAL CENTER FOR NUMERICAL METHODS IN ENGINEERING Edificio C1, Campus Norte UPC Gran Capitán s/n 08034 Barcelona, Spain www.cimne.upc.es First edition: October 2003 DISCRETE ELEMENT MODELLING OF ROCK CUTTING PROCESSES INTERACTION WITH EVALUATION OF TOOL WEAR Monograph CIMNE M87  The authors ISBN: 84-95999-44-7 Depósito legal: B-47351-2003 Contents 1 Introduction 2 2 Modelling of rock cutting 3 2.1 Physical phenomena in rock cutting . . . . . . . . . . . . . . . . . . . . 3 2.2 Analytical models of rock cutting . . . . . . . . . . . . . . . . . . . . . 4 2.3 Basic assumptions of numerical model of rockcutting................................. 5 3 Discrete Element Method formulation 7 3.1 Equationsofmotion ............................ 7 3.2 Evaluation of contact forces . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2.1 Decomposition of the contact force . . . . . . . . . . . . . . . . 9 3.2.2 Normal contact force . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.3 Tangential frictional contact . . . . . . . . . . . . . . . . . . . . 11 3.3 Backgrounddamping............................ 12 3.4 Numericalstability ............................. 13 3.5 Contactdetection.............................. 15 3.6 Contact search algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.7 Numerical implementation – programming aspects . . . . . . . . . . . . 20 4 Micromechanical models with cohesion 24 4.1 Elastic perfectly brittle contact model . . . . . . . . . . . . . . . . . . . 24 4.2 Elasto-plastic contact model with nonlinear softening . . . . . . . . . . 26 4.3 Simplified elasto-plastic contact model with softening . . . . . . . . . . 29 4.4 Contact model with elastic damage . . . . . . . . . . . . . . . . . . . . 31 4.5 Contact model with friction, wear and heat generation . . . . . . . . . 33 4.6 Validation of micromechanical models . . . . . . . . . . . . . . . . . . . 33 4.6.1 Elasto-plastic contact model with linear softening . . . . . . . . 34 4.6.2 Elastic damage contact model with linear softening . . . . . . . 38 4.7 Summary .................................. 40 i 5 Study of macroscopic material properties 43 5.1 Simulation of uniaxial compression test . . . . . . . . . . . . . . . . . . 43 5.1.1 Numericalmodel .......................... 43 5.1.2 Numericalresults.......................... 44 6 Wear evaluation 51 6.1 Basic concepts of wear . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.2 Wear of rock cutting tools . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.3 Modellingofwear.............................. 54 6.4 Wear model with thermal effects included . . . . . . . . . . . . . . . . 55 6.5 Numerical implementation of wear evaluation algorithm . . . . . . . . . 55 6.6 Estimation of wear constants from laboratory tests . . . . . . . . . . . 57 7 Discrete Element formulation for thermal and thermo-mechanical problem 60 7.1 Heatbalanceequation ........................... 60 7.2 Thermal boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 61 7.3 Solution of thermo-mechanical problem . . . . . . . . . . . . . . . . . . 62 7.4 Material thermal properties . . . . . . . . . . . . . . . . . . . . . . . . 63 7.5 Parameters characterizing heat transfer on the surface . . . . . . . . . . 64 7.6 Benchmarks of thermal analysis using DEM . . . . . . . . . . . . . . . 65 7.6.1 Internal heat generation in parallel-sided slab . . . . . . . . . . 65 7.6.2 Semi-infinite solid subjected to a unit surface heat flux . . . . . 68 7.6.3 Transient heat conduction in an infinite parallel-sided slab . . . 71 8 Numerical simulation of rock cutting and wear evaluation 74 8.1 2D mechanical simulation of rock cutting –modelI .................................. 74 8.2 2D mechanical simulation of rock cutting (model II) . . . . . . . . . . . 75 9 Wear evaluation in multi-cycle analysis of rock cutting 82 9.1 Methodology ................................ 82 9.2 Simulation of wear of the test tooth . . . . . . . . . . . . . . . . . . . . 83 9.2.1 Description of the model . . . . . . . . . . . . . . . . . . . . . . 83 9.2.2 Results................................ 83 9.3 Wear simulation of a ripper tooth . . . . . . . . . . . . . . . . . . . . . 86 9.3.1 Modeldefinition........................... 86 9.3.2 Numericalresults.......................... 86 9.3.3 Conclusions ............................. 87 ii 10 Thermo-mechanical transient analysis of rock cutting 94 10.1 2D thermo-mechanical simulation of rockcutting................................. 94 11 Thermomechanical analysis of wear of a ripper 99 11.1Modeldefinition............................... 99 11.2 Transient thermo-mechanical solution without wear . . . . . . . . . . . 100 11.3 Quasi-stationary solution of thermal problem . . . . . . . . . . . . . . . 104 11.4 Thermomechanical simulation of rock cutting with wear evaluation . . 105 11.4.1 Model I – initial temperature distribution from transient thermal analysis ............................... 105 11.4.2 Model II – initial temperature distribution from quasi-stationary thermalanalysis........................... 107 12 Simulation of dredging 108 12.1ModelI ................................... 108 12.1.1Modeldefinition........................... 108 12.1.2 Thermo-mechanical analysis of dredging . . . . . . . . . . . . . 109 12.1.3 Thermo-mechanical analysis of dredging with wear evaluation . 114 12.2ModelII................................... 121 12.2.1Modeldefinition........................... 121 12.2.2 Thermo-mechanical analysis of dredging with wear evaluation . 121 13 3D simulation of rock cutting 127 13.1 3D simulation of interaction of a tooth with granular medium . . . . . 127 13.2 3D simulation of rock cutting (model I) . . . . . . . . . . . . . . . . . . 129 13.3 3D simulation of rock cutting (model II) . . . . . . . . . . . . . . . . . 130 13.4 3D simulation of rock cutting (model III) . . . . . . . . . . . . . . . . . 132 14 Numerical examples – granular flow 136 15 Conclusions 137 iii Abstract The document presents a numerical model of rocks and soils using spherical Discrete Elements, also called Distinct Elements. The motion of spherical elements is described by means of equations of rigid body dynamics. Explicit integration in time yields high computational efficiency. Spherical elements interact among one another with contact forces, both in normal and tangential directions. Efficient contact search scheme based on the octree structures has been implemented. Special constitutive model of contact interface taking into account cohesion forces allows us to model fracture and decohesion of materials. Numerical simulation predicts wear of rock cutting tools. The developed numerical algorithm of wear evaluation allows us us to predict evolution of the shape of the tool caused by wear. Results of numerical simulation are validated by comparison with experimental data. 1 Chapter 1 Introduction This report presents numerical modelling of rock cutting processes with tool wear evaluation. The numerical model is based on the discrete element method in the version using spherical particles. This method is widely recognized as a suitable tool to model problems characterized with strong discontinuities. The discrete element model assumes that material can be represented by an assembly of rigid particles interacting among themselves. The overall behaviour of the system is determined by the cohesive/frictional contact laws. Granular and particulate materials are characterized with inherent discrete nature. Discontinuity in rocks can occur due to progressive damage during cutting. Discrete element formulation using spherical or cylindrical particles was first proposed by Cundall and Strack [1]. Similar formulation has been developed by Rojek et al. in [2] and implemented in the explicit dynamic finite element code Simpact. 2 Chapter 2 Modelling of rock cutting 2.1 Physical phenomena in rock cutting The type of failure during rock cutting depends on the type of rock. In some rocks brittle failure with characteristic chip formation occurs (Fig. 2.1a), while others fail by ductile ploughing (Fig. 2.1b). According to Gehring [3] ductility (or brittleness) of a rock can be measured by the ratio of unconfined compressive strength to unconfined tensile strength, which he called a ductility number. Typical brittle rocks have ductility number greater than 15, while typical ductile rocks are characterised with ductility numbers lower than 9. Ductile cutting is very disadvantageous since the wear of tools is bigger than in case of brittle cutting. By adequate design of tool sometimes it is possible to change the mode of failure in cutting. a) b) Figure 2.1: Modes of failure in rock cutting: a) brittle, b) ductile Formation of a chip in brittle cutting is initiated in a crushing zone near the tooth tip (Fig. 2.2). The forces are transmitted from the crushing zone to the intact rock and microcracks are initiated. Near the crushing zone localised shear zone is developed. Shear failure further in the intact rock can bifurcate into a tensile crack. This combined 3 kn kT cn m Figure 3.3: Model of the contact interface where nis the unit vector normal to the particle surface at the contact point (this implies that it lies along the line connecting the centers of the two particles) and directed outwards from the particle 1. The contact forces Fnand FTare obtained using a constitutive model formulated for the contact between two rigid spheres (Fig. 3.3). The contact interface in our formulation is characterized by the normal and tangential stiffness knand kT, the Coulomb friction coefficient µ, and the contact damping coefficient cn. 3.2.2 Normal contact force The normal contact force Fnis decomposed to the elastic part Fne and to the damping contact force Fnd Fn=Fne +Fnd .(3.15) The damping is used to decrease oscillations of the contact forces and to dissipate kinetic energy. The elastic part of the normal contact force Fne is proportional to the normal stiffness knand the penetration of the two particle surfaces urn Fne =knurn .(3.16) The penetration is calculated as urn =d−r1−r2,(3.17) where dis the distance of the particle centres, and r1,r2their radii. In the formulation used in the present study no cohesion is allowed, so no tensile normal contact forces are allowed Fne ≤0.(3.18) If urn <0, the formula (3.16) is valid, otherwise Fne = 0. The contact damping force is assumed to be of viscous type Fnd =cnvrn (3.19) 10 urT F T F n || || | | m urT F T F n || || | | m kT a) b) Figure 3.4: Friction force vs. relative tangential displacement a) Coulomb law, b) regularized Coulomb law proportional to the normal relative velocity vrn of the centres of the two particles in contact vrn = ( ˙ u2−˙ u1)·n.(3.20) The value of damping cncan be taken as a fraction of the critical damping Ccr for the system of two rigid bodies with masses m1and m2, connected with a spring of the stiffness kn(cf. [13]) Ccr = 2sm1m2kn m1+m2 .(3.21) 3.2.3 Tangential frictional contact In the absence of cohesion (if the particles were not bonded at all or after the cohesive bond has been broken) the tangential reaction FTis brought about by friction opposing the relative motion at the contact point. The relative tangential velocity at the contact point vrT is calculated from the following relationship vrT =vr−vr·n,(3.22) vr= ( ˙ u2+ω2×rc2)−(˙ u1+ω1×rc1),(3.23) where ˙ u1,˙ u2, and ω1,ω2are the translational and rotational velocities of the particles, and rc1and rc2are the vectors connecting particle centres with contact points. The relationship between the friction force kFTkand relative tangential urT displacement for the classical Coulomb model (for a constant normal force Fn) is shown in Fig. 3.4a. This relationship would produce non physical oscillations of the friction force in the numerical solution due to possible changes of the direction of sliding velocity. To prevent this the Coulomb friction model must be regularized. A possible 11 regularization procedure involves decomposition of the tangential relative velocity into a reversible and irreversible parts, vr rT and vir rT , respectively: vrT =vr rT +vir r.(3.24) This is equivalent to formulation of the frictional contact as a problem analogous to that of elastoplasticity, which can be seen clearly from the friction force-tangential displacement relationship in Fig. 3.4b. This analogy allows us to calculate the friction force employing the radial return algorithm analogous to that used in elastoplasticity. First a trial state is calculated Ftrial T=Fold T−kTvrT ∆t , (3.25) and then the slip condition is checked φtrial =kFtrial Tk−µ|Fn|.(3.26) If φtrial ≤0, we have the case of stick contact and the friction force is assigned the trial value Fnew T=Ftrial T,(3.27) otherwise (slip contact) a return mapping is performed Fnew T=µ|Fn|Ftrial T kFtrial Tk.(3.28) 3.3 Background damping A quasi-static state of equilibrium of the assembly of particles can be achieved by application of adequate damping. Described previously contact damping is a function of the relative velocity of contacting body. It is sometimes necessary to apply damping for non-contacting particles to dissipate their energy. There are two types of such damping, referred here as background, implemented in our formulation, one of viscous type and the other of non-viscous type. In both cases damping terms Fdamp iand Tdamp i are added to equations of motion (3.1) and (3.2) mi¨ ui=Fi+Fdamp i,(3.29) Ii˙ ωi=Ti+Tdamp i.(3.30) with damping terms given by: •for viscous damping Fdamp i=−αvtmi˙ ui,(3.31) Tdamp i=−αvrIiωi,(3.32) 12 •for non-viscous damping Fdamp i=−αnvtkFik˙ ui k˙ uik,(3.33) Tdamp i=−αnvrkTikωi kωik.(3.34) where αvt,αvr,αnvt and αnvr are respective damping constants. It can be seen from Eqs. (3.31)–(3.34) that non-viscous like viscous damping is opposite to velocity, the difference consists in the evaluation of the magnitude of damping force – viscous damping is proportional to the magnitude of velocity, while non-viscous damping is proportional to the magnitude of resultant force and moment. 3.4 Numerical stability Explicit integration in time yields high computational efficiency. Therefore the method enables us to analyse large models. The known disadvantage of the explicit integration scheme is its conditional numerical stability imposing the limitation on the time step ∆t. The time step ∆tmust not be larger than a critical time step ∆tcr ∆t≤∆tcr (3.35) determined by the highest natural frequency of the system ωmax ∆tcr =2 ωmax .(3.36) If damping exists, the critical time increment is given by ∆tcr =2 ωmax µq1 + ξ2−ξ¶,(3.37) where ξis the fraction of the critical damping corresponding to the highest frequency ωmax. Exact determination of the highest frequency ωmax would require solution of the eigenvalue problem defined for the whole system of connected rigid particles. In an approximate solution procedure, eigenvalue problems can be defined separately for every rigid particle using the linearized equations of motion mi¨ ri+kiri=0,(3.38) where mi={mimimiIiIiIi}T,ri={(ux)i(uy)i(uz)i(θx)i(θy)i(θz)i}T,(3.39) 13 and kiis the stiffness matrix accounting for the contributions from all the penalty constraints active for the i-th particle. Equation (3.39) defines the vectors miand rifor a spherical particle in three-dimensional space. For a cylindrical particle in a two-dimensional model the respective vectors are defined as follows: mi={mimiIi}T,ri={(ux)i(uy)i(θz)i}T.(3.40) Equation (3.38) leads to an eigenproblem kiri=λjmiri,(3.41) where eigenvalues λj(j= 1, . . . , 6 in 3D case, and j= 1,2,3 for 2D case) are the squared frequencies of free vibrations: λj=ω2 j.(3.42) In a 3D problem, three of six frequencies ωjare translational, and the other three – rotational. In the algorithm implemented, a further simplification is assumed. The maximum frequency is estimated as the maximum of natural frequencies of mass–spring systems defined for all the particles with one translational and one rotational degree of freedom. The translational and rotational free vibrations are governed by the following equations: mi¨un+knun= 0 ,(3.43) Ii¨ θ+kθθ= 0 ,(3.44) where it is assumed that the translational motion is due to the contact interaction in the normal direction (the spring stiffness knrepresents the penalty stiffness in the normal direction), and the rotational stiffness is due to the contact stiffness (penalty) in the tangential direction. Given the tangential penalty kT, it can be shown that the rotational stiffness kθcan be obtained as kθ=kTr2,(3.45) where ris the length of the vector connecting the mass centre to the contact point. The natural frequency of the translational vibrations is given by the following equation: ωn=skn mi ,(3.46) while the rotational frequency ωθcan be obtained from the formula ωθ=skθ Ii .(3.47) 14 With the rotational inertia of a sphere I=2 5mr2(3.48) and kθgiven by Eq. (3.45), the rotational frequency can be calculated as ωθ=s5kT 2mi .(3.49) If kT=kn, the rotational frequency ωθis considerably higher than the translational frequency ωnobtained from Eq. (3.46), which results in a smaller critical time increment, cf. Eq. (3.36). To avoid the determination of a critical time step by the rotational frequencies, the rotational inertia terms are scaled adequately. The concept of scaling rotational inertia terms is commonly used for shell elements, cf. [14]. 3.5 Contact detection Changing contact pairs of elements during the analysis process must be automatically detected. The simple (“brute-force”) approach to identify interaction pairs by checking every sphere against every other sphere would be very inefficient, with the computational time proportional to n2(of order O(n2)), where nis the number of elements. A number of more effective methods has been developed to determine the particle interaction. An overview of contact detection methods developed for discrete element method can be found in [9] and [15]. In usually used contact detection schemes prior to the contact resolution objects are spatially ordered using an appropriate sorting algorithm, cf. [9]. Spatial sorting enables efficient determination of neighbouring objects, so the subsequent contact check can be limited to the pairs of objects lying close to each other. In the absence of cohesion the contact is found if the penetration between two bodies is found. In case of spheres or cylinders this is expressed by the simple condition urn ≤0,(3.50) where urn is calculated according to the formula (3.17). If cohesive bonds are active the interaction between the spheres can also occur if urn >0. The condition (3.50) is then replaced by urn ≤u+ rn,max ,(3.51) where u+ rn,max is the separation of two spheres or cylinders corresponding to the complete breakage of cohesive bonds, which must be evaluated for a given constitutive model and material properties. 15 Since the verification of the contact condition (3.50) or (3.51) is simple, we can say, that basically the effectiveness of the contact detection scheme in our case depends on the effectiveness of the sorting algorithm. Various spatial sorting algorithms are known. Most popular of them can be classified into following groups: •grid subdivision (L¨ohner [16]) •binary trees (Knuth [17], Bonet [18]) •quad-trees (in 2D) and octrees (in 3D) (Knuth [17], Samet [19]) •body based cells (Greengard [20]) •spatial heapsort (Williams [9]) Grid subdivision Figure 3.5: Grid subdivision The problem domain in this method is discretized into equal rectangular (in 2D) or hexahedral (in 3D) cells shown in Fig. 3.5 for 2D case. Cells are sometimes called bins, cf. [16]. Objects are associated with cells (bins) based on their coordinates. The performance of this method depends on the trade-off between the cell size and the number of objects per cell. The efficiency of this method is assured if the objects are evenly distributed among cells. Adaptive grid method can be employed in case if the distribution of the objects is not uniform. Binary trees An effective sorting technique for an uneven space distribution of objects is based on the binary tree structures. The main concept of this scheme is shown in Fig. 3.6. The 16 Figure 3.6: Binary trees domain with objects is divided into rectilinear cells which has maximum 2 objects. In case the number of objects exceeds 2 the cell is subdivided into two cells. The subdivision is usually done alternately along x,yand zaxes, cf. [18]. The domain with objects discretized in this way is represented by a binary tree. The structure allows us to identify easily objects lying in a particular subdomain by traversing the tree from the top downwards. The cost of this search is of order O(Nlog2N). Quadtrees and octrees Figure 3.7: Quadtree sorting Another effective sorting technique for an uneven space distribution of objects is based on quadtree (2D) and octree (3D) structures. The main concept of the quadtree sorting is shown in Fig. 3.7. In the quadtree sorting scheme a 2D domain with objects is divided into rectilinear cells (quads) with up to 4 objects assigned. In case the number of objects exceeds 4, the cell is subdivided into four cells (quads), and the objects relocated into new cells. The domain with objects discretized in this way is 17 represented by a quad-tree. The structure allows us to identify easily objects lying in a particular subdomain by traversing the quadtree from the top downwards. The cost of this search is of order O(Nlog4N). Extension of the quadtree scheme on on the octree scheme for 3D contact search is straightforward. Now a 3D domain is discretized recursively into brick cells (octants) with up to eight objects. In the number of objects in an octant exceeds eight, the octant is subdivided into eight cells, and the objects are assigned to these cells. Objects lying in a particular subdomain are found by traversing the octree from the top downwards. The cost of this search is of order O(Nlog8N). This method will be used in our search algorithm combined with body based cells. Body based cells Figure 3.8: Body based cells Another search scheme takes advantage of the assumption that object configuration evolves slowly and the new contacts can be formed only between objects that are sufficiently close at a certain stage. The list of potential contacts for each object includes the objects lying within a cell surrounding an object under consideration. The concept of this method is presented graphically in Fig. 3.8. This method will be used in our search algorithm combined with quadtree and octree sorting scheme. Spatial heapsort The basis of spatial heapsort method is sorting of the objects according to their coordinates, cf. [9]. Sorted objects are stored in a binary tree structure. Heapsort is of order O(Nlog2N) 18 3.6 Contact search algorithm In our formulation the contact search is based on the quad-tree and octree sorting coupled with the body based cell technique. Building octree or quad-tree structures at every time step can be quite expensive. Since time steps are very small, most of contacting pairs can be the same as those in the previous step. Using the information about contact pairs existing in the previous step can speed up the contact search. Therefore the algorithm of contact detection scheme consists of two stages: (i) global search, identifying the pairs of potential contacts, based on the quadtree and octree sorting, (ii) local search, verifying the list of potential contacts (typical for body based cell algorithm). Stage 1: Global search Balls are ordered in space using octree or quad-tree structures. A list of existing and potential contact pairs is created. In the list of potential contacting balls the body based cell concept is employed, each object has a circular or spherical cell assigned and neighbouring objects lying within the cell or intersecting it are included into the list of potential contacts. Global search with a quad-tree (or octree) structure being rebuilt and the list of contacting balls being updated is performed every certain number of steps. Stage 2: Local search Local search is done at each time step. Contact conditions for the ball pairs from the list of potential contacts are verified and actual contacting pairs are found. The body cell radius is assumed to be equal to the sphere/cylinder radius increased by a constant distance which is a function of a user defined parameter ctol (Fig. 3.9a). This parameter determines the length of the list of potential contacts and the frequency of global search (updating the list of potential contacts). New potential contacts are added to the list if the separation between two balls is smaller than 2√2ctol in two-dimensional problem and 2√3ctol in three-dimensional problem. Conversely if the separation in an existing contact pair gets larger than these distances the contact pair is deleted from the list. Global search interval is determined by monitoring the maximum displacement component of any ball accumulated from the time of last contact search umax = n=1, NB max i=1,nd³u(n) i´(3.52) where ndis the dimension of the problem (nd= 2 for two-dimensional problem and nd= 3 for three-dimensional problem) and NBis the number of balls. If displacement 19 ut t Rt Figure 4.2: Tangential contact force in the elastic perfectly brittle model (tensile normal contact force) s t Rn Rt Figure 4.3: Failure surface for the elastic perfectly brittle model 4.2 Elasto-plastic contact model with nonlinear softening This section summarizes the theoretical formulation of the constitutive model for geomaterials developed in the University of Padova which has been presented in more detail in the First Year Activity Report. A contact law relates the contact force acting between two spheres to their relative displacement. The model developed is based on the nonlinear elasto-plastic relationships established for the normal and tangential direction σi=ki(ui−up i), i =n, t (4.7) where: σi— normal or tangential contact force; ui— total relative displacement; up i— plastic component of the relative displacement; 26 ki— elastic parameters, which characterize contact interface stiffness n, t — normal and tangential direction; Two separate different potential functions: one for the normal and another for the tangential stress-strain relation are adopted. A general non-linear one-dimensional yield criterion has been taken in the following form: f(σ, α) = |σ|−[σY+H(α)] ,(4.8) where H(α) is the function of the plastic strain and it is expressed as a function of the softening or hardening parameter α, the so called internal hardening variable. Two different forms of H(α) can be chosen: 1. parabolic H(α) = b1α2+b2α+b3,(4.9) 2. hyperbolic (used only for tangential direction) Ht(α) = b1αd t b2+b3αd t ,(4.10) where b1,b2,b3and mare parameters defining plastic function, defined separately for tension, compression and shear action. The variation of constants defining H(α) changes significantly the shape of the yield surface and allows us to apply the model to different soils and rocks. Dependence of the tangential force on the normal force is expressed by the application of a two-dimensional failure criterion. In the proposed criterion the slip between two surfaces of a continuum is assumed to occur when the shear stress τon any plane at a point in the soil material reaches a so-called critical value. This value depends non-linearly upon the normal stress in the same plane. In the numerical implementation at each step the value of the normal force is fixed and a two-dimensional yield surface is used in order to establish the elastic range of the tangential force, as shown in Fig. 4.4. The basic concept of the strength criterion is that the shear strength of a rock material is made up of two parts varying with the normal stress: an elastic “cohesive” part and a frictional plastic part. The shear strength can be mathematically resumed as follow: τ=f(ut, up t) + f1(σn, H),(4.11) where: τ— critical tangential contact force; σn— normal contact force; 27 Figure 4.4: Scheme of application of the 2D yield surface in return mapping algorithm for one-dimensional constitutive law ut, up t— tangential total and plastic relative displacements, respectively; H— general non-linear function of plasticity. The cohesive part of the formulation is non-linear elastic regarding the shear modulus: in fact it was experimentally noticed that for some soils, such as the stiff clay, the slope of the shear-normal stresses relationship changes when the normal applied stress increases or decreases. Moreover it was deducted that there is a sort of contribute which affects the real value of the shear modulus of the material. In this way it is possible to explain the heavy dependency on the normal applied load. The plastic softening (or hardening) frictional behavior starts after that the value of τreaches the limit flow stress value, which is a function of the normal applied stress. As experimentally observed, if the normal applied stress increases, also the tangential peak stress increases. Due to the previous considerations the equation (4.11) can be more specifically rewritten as follows τ=kt(σn)(ut−up t)−[n(σn) + H1(α)] ,(4.12) where: α— internal hardening variable, n(σn) — parameter which delimits the flow range, 28 kt— shear modulus of the contact interface, and it depends on the applied normal contact force kt=m3 q||σn|−r|,(4.13) where mand rare constants which depend on the material. The parameter n(σn) in this case it is represented as function of the normal stress and is expressed trough two terms: n(σn) = η(σn+c).(4.14) This formulation of n(σn) was chosen in order to reproduce the behavior of a cohesive soil that presents the characteristic of increasing the shear strength by increasing the normal applied stress. This behavior is characteristic of a wide range of soils, from sand to some soft rock. The material parameters ηand care constant during the loading process and they respectively represent: η— parameter that gives a contribute to the frictional behavior of the material; c— is a constant which correlates the cohesive properties of the material. Flow rule The plastic strain increments for both the formulations (i.e. the normal and tangential strain) are calculated from the following flow rule: ˙up=γsign(σ),(4.15) where γ— absolute value of the slip rate (γ > 0); ˙up— derivative of the plastic relative displacement (normal or tangential); σ— contact force in the normal or tangential direction. 4.3 Simplified elasto-plastic contact model with softening A special case of the previously presented elasto-plastic contact formulation has been taken using the following simplifying assumptions: 1. Elasto-plastic contact law with linear softening is assumed for shear and normal tensile contact forces 2. Elastic linear law is assumed for compression 29 3. Contact stiffness (elastic modulus) in the tangential direction is constant and independent of normal contact force Simplified contact laws for the normal and tangential direction are shown in Figs. 4.5 and 4.6, respectively. compression tension u unloading s n tg =kan ++ tg =kan -- unmax p un p Figure 4.5: Simplified elasto-plastic contact law for the normal direction u unloading t t tg =kat utmax p ut p Figure 4.6: Simplified elasto-plastic contact law for the tangential direction The softening function H(α) given by Eq. (4.9) is now simplified to H(α) = b2α=−H0up,(4.16) where H0is the constant softening modulus (assumed positive value) and the plastic part of the relative displacement uptaken as the softening parameter α. Contact yield condition for the simplified model is shown in Fig. 4.7. Yielding of contact bonds can occur under a combined tension and shear action (if the normal contact force is tensile) or due to shear only (if the normal contact force is compressive). After contact bonds are broken due to yielding standard frictional contact can occur between the spherical elements. 30 s t Initial yield” surface Current “yield” surface Weakeningdue totensionandshear Weakeningduetoshear Figure 4.7: Contact yield condition for simplified model 4.4 Contact model with elastic damage This model can be considered as a generalization of the elastic perfectly brittle model in place of brittle failure a certain softening slope defined by softening modulus His introduced into the force-displacement relationship (Fig. 4.8). compression tension u s n kn -unmax kn H kn D Rn Figure 4.8: Normal contact force in the contact model with elastic damage The constitutive relationship for 1D elastic damage is given by: σ=kD nun= (1 −ω)knun(4.17) where: kD n— elastic damaged secant modulus, ω— scalar damage variable. 31 Scalar damage variable ωis a measure of material damage, for undamaged state ω= 0 and for damaged state 0 < ω ≤1. Scalar damage variable ωcan be written in the following form: ω=ψ(un)−1 ψ(un)(4.18) where ψ(un) is a function of total relative displacement. For linear strain softening ψ(un) is defined by ψ(un) =                      1 for un≤Rn kn k2 nun HRn+knRn−Hknun for Rn kn≤un≤Rn kn +Rn H ∞for un≥Rn kn +Rn H (4.19) where Rnis the initial tensile strength and His the softening modulus (taken as positive). Similar contact force–displacement law with damage (Fig. 4.9) can be introduced for the tangential direction. Then the bond decohesion can occur either due to tension or shear. u t t utmax kt H kt D Rt Figure 4.9: Tangential contact force in the contact model with elastic damage An optional failure criterion has been implemented in this model with failure criterion based on the tensile contact force only. This criterion is supposed to model better fracture of brittle materials, as this macroscopic failure is explained by brittle rupture of atomic bonds in tension. This microscopic mechanism explains macroscopic strain softening behaviour both under compression or tension. After contact bonds are broken due to damage (ω= 0) standard frictional contact can occur between the spherical elements. 32 4.5 Contact model with friction, wear and heat generation Contact model taking into account friction, wear and heat generation is used to model tool-rock/soil interaction. Cohesion is not included into the formulation of this model. Normal contact forces can be compressive only. Normal contact force σis calculated from the linear relationships: σ=knun(4.20) where un≤0 is the normal relative displacement (penetration of one sphere against another). Frictional contact force is evaluated using the regularized Coulomb law τ=µ|σ|(4.21) where µis the Coulomb friction coefficient. The frictional dissipation rate ˙ Dis calculated as ˙ D=τvt(4.22) where vtis the relative tangential velocity. The frictional dissipation is used to calculate wear rate ˙wand heat generated by friction Qgen. Heat generated by friction is equal Qgen =χ˙ D(4.23) where χis the part of the friction work converted to heat (χ≤1). Heat generation through frictional dissipation is assumed to be absorbed equally by the two particles in contact. Wear rate is calculated using the Archard law presented later on in Sec. 6.3. 4.6 Validation of micromechanical models Micromechanical models have been implemented in the Stampack (program developed by CIMNE) and numerical implementation has been validated by running test examples of shear, tension and compression as well as of shear combined with tension and compression. This section presents results of the tests for the simplified elasto-plastic contact model with linear softening and the contact model with elastic damage. Numerical tests for the general elasto-plastic model with non-linear softening are presented in a separate section. Numerical tests presented here allowed us to use the program to model rock specimen. Results of material model calibration are presented in an annex to WP6. Validation tests were carried out using a discrete element model of two discs (Fig. 4.10). The load was introduced under displacement control. The lower disc was fixed 33 Figure 4.10: Discrete element micromechanical model and translational motion was prescribed for the upper disc in order to obtain compression or tension at the contact interface between the discs. For the rotational motion prescribed for the upper disc the relative tangential displacement at the contact interface produced shear action at the interface. The combination of translation and rotation of the upper disc allowed us to obtain combined normal and tangential contact action. Parameters defining contact interface properties Contact interface is defined with the following parameters: kn= 2 ·1010N/m, kt= 105N, Rn= 105N, H= 1010N/m, µ= 0.839. 4.6.1 Elasto-plastic contact model with linear softening Simulation of pure tension, compression and shear The normal and tangential force were calculated for pure tension, compression and shear. Loading was applied under displacement control. Variation of the normal contact force against the normal relative displacement under monotonic loading is shown in Fig. 4.11. The peak force of 100 kN is achieved for relative displacement 0.005 mm, and then the yielding of the contact bond occurs until complete breakage at separation of 0.015 mm. Figure 4.12 shows variation of the normal contact force under tensile loading with unloading and reloading which occurs with elastic modulus. Compressive normal force under monotonic loading is shown in Fig. 4.13. Compression is treated as an elastic process. Tangential contact force vs. tangential relative displacement under shear monotonic loading is shown in Fig. 4.14. Elastic increase of the tangential force and subsequent drop due to yielding can be observed. 34 -20000 0 20000 40000 60000 80000 100000 120000 0 5e-06 1e-05 1.5e-05 2e-05 normal contact force (N) normal relative displacement (m) Figure 4.11: Normal contact force vs. normal relative displacement under monotonic tensile loading (elasto-plastic contact model with linear softening) -20000 0 20000 40000 60000 80000 100000 120000 0 5e-06 1e-05 1.5e-05 2e-05 normal contact force (N) normal relative displacement (m) Figure 4.12: Normal contact force vs. normal relative displacement under tensile loading with unloading and reloading (elasto-plastic contact model with linear softening) -600000 -500000 -400000 -300000 -200000 -100000 0 -3e-05 -2.5e-05 -2e-05 -1.5e-05 -1e-05 -5e-06 0 normal contact force (N) normal relative displacement (m) Figure 4.13: Normal contact force vs. normal relative displacement under compressive monotonic loading (elasto-plastic contact model with linear softening) 35 -20000 0 20000 40000 60000 80000 100000 0 5e-06 1e-05 1.5e-05 2e-05 normal contact force (N) normal relative displacement (m) Figure 4.26: Normal contact force vs. normal relative displacement under shear combined with tension (elastic damage contact model with linear softening) 42 Chapter 5 Study of macroscopic material properties Introduction This study is aimed to find out similarities and differences in macroscopic behaviour of a material modelled with different micromechanical models and influence of micromechanical parameters on macroscopic material behaviour. Uniaxial compression test has been chosen as a test example. 5.1 Simulation of uniaxial compression test 5.1.1 Numerical model Macroscopic response of a square material sample subjected to uniaxial compression has been studied. Figure 5.1 presents the material sample prepared for testing. Material sample 109×109 mm is represented by an assembly of randomly compacted 2100 discs of radii 1 ÷1.5 mm. The loading has been introduced under kinematic control by prescribing the motion of the right and left walls. The deformation in the y direction was free. The velocity of the displacement of the walls was 1 mm/s which was found to be sufficiently low to obtain quasi-static loading. The micromechanical parameters constant in all the numerical tests are the following: contact stiffness in the normal and tangential directions kn=kT= 20 GPa and the friction coefficient µ= 0.839. 43 Figure 5.1: Initial geometry of a material sample 5.1.2 Numerical results Numerical simulation of the compression test using elastic perfectly brittle micromechanical model with the cohesive bond strengths Rn= 0.1 MN/m RT= 1 MN/m has been found out to give macroscopic behaviour corresponding to sandstone, cf. [21]. The stress-strain relationships obtained for three different micro-mechanical models yielding similar macroscopic compressive strength are shown in Fig. 5.2. The stresses were calculated by taking the sum of the contact interactions between the walls and particles. Elastic perfectly brittle model as well as the model with elastic damage with Figure 5.2: Stress-strain relationships for three different constitutive models giving similar macroscopic compressive strength 44 linear softening were used with the same cohesive bond strengths Rn= 0.1 MN/m and RT= 1 MN/m. The softening in the latter model was defined by the softening modulus H= 2 ·1011 N/m. The parameters for the elastic-plastic micromechanical model were the following: cohesive bond strengths Rn= 0.04 MN/m and RT= 0.16 MN/m, softening modulus H= 1010 N/m. The values yielding the desired macroscopic strength were assumed here. It can be seen in Fig. 5.2 that stress-strain curves obtained with different micromechanical models are similar. It shows that similar macroscopic behaviour can be reproduced using different micromechanical models. Failure modes obtained for different micromechanical models are shown in Figs. 5.3–5.5 in the form of final deformed configurations with displacement field and distribution of broken bonds. One can see localization zones and discontinuities of the displacement fields. a) b) Figure 5.3: Failure mode obtained for the elastic perfectly brittle model, a) displacement distribution, b) broken bonds a) b) Figure 5.4: Failure mode obtained for the elastic damage model with linear softening (softening modulus H= 2 ·1011 N/m), a) displacement distribution, b) broken bonds 45 a) b) Figure 5.5: Failure mode obtained for the elastic plastic model with linear softening (softening modulus H= 1010 N/m), a) displacement distribution, b) broken bonds Comparison of Figs. 5.3-5.5 indicates that the failure modes obtained with the three models analysed are similar, strain localisation and fracture occurs along the same lines although some differences can be noted between elastic plastic model (Fig. 5.5) and the other two models (Figs. 5.3 and 5.4), broken bonds obtained with elastic-plastic model are localized in smaller zones along fracture lines. Elastic damage with different softening moduli A number of analyses have been carried out using different softening moduli. Influence of softening value on macroscopic behaviour has been studied. Stress-strain relationships for micromechanical model with elastic damage using different softening moduli are shown in Fig. 5.6. It can be seen that in a certain range of values of softening modulus the stress-strain curve does not change much, the compressive strength and Figure 5.6: Stress-strain relationships for elastic damage micromechanical model with different softening moduli 46 curves in post-critical deformation are similar for softening moduli of 2 ·1011, 2 ·1010 and 5·109N/m. Significant difference can be observed for the lowest softening modulus only. This difference can be understood by analysing the failure modes presented in Figs. 5.4 and 5.7–5.9. While the failure mode obtained for softening moduli of 2 ·1011, 2·1010 and 5·109N/m (Figs. 5.4, 5.7 and 5.8) are similar one to another, in case of the softening modulus 2 ·109N/m (Fig. 5.9) different localisation zones can be observed. a) b) Figure 5.7: Failure mode obtained for the elastic damage model with linear softening (softening modulus H= 2 ·1010 N/m), a) displacement distribution, b) broken bonds a) b) Figure 5.8: Failure mode obtained for the elastic damage model with linear softening (softening modulus H= 5 ·109N/m), a) displacement distribution, b) broken bonds 47 a) b) Figure 5.9: Failure mode obtained for the elastic damage model with linear softening (softening modulus H= 2 ·109N/m), a) displacement distribution, b) broken bonds Elastic plastic model with different softening moduli A number of analyses have been carried out using elastic plastic model with different softening moduli. Influence of softening value on macroscopic behaviour has been studied. Stress-strain relationships for different softening moduli are shown in Fig. 5.10. Figure 5.10: Stress-strain relationships for elastic-plastic micromechanical model with different softening moduli It can be seen that in this case the solution is very sensitive to softening value. Although the failure modes (Figs. 5.5, 5.11–5.13) are similar, macroscopic compressive strength increases significantly with the increase of softening modulus. In this model softening is very important parameter influencing the macroscopic strength. 48 a) b) Figure 5.11: Failure mode obtained for the elastic-plastic model with linear softening (softening modulus H= 2 ·1010 N/m), a) displacement distribution, b) broken bonds a) b) Figure 5.12: Failure mode obtained for the elastic-plastic model with linear softening (softening modulus H= 5 ·109N/m), a) displacement distribution, b) broken bonds a) b) Figure 5.13: Failure mode obtained for the elastic-plastic model with linear softening (softening modulus H= 2 ·109N/m), a) displacement distribution, b) broken bonds 49 Conclusions Numerical studies carried out with different micromechanical models show that similar macroscopic behaviour can be reproduced using different micromechanical models. Failure modes obtained using these models are quite similar. Strain localisation occurs in the same zones. Macroscopic properties were compared using stress-strain curves. Influence of strain softening on macroscopic material behaviour has been studied. It has been found out that elastic-plastic model is more sensitive to changes of strain softening. 50 Chapter 6 Wear evaluation In this section numerical model is extended on the problems of wear evaluation. Main factors influencing wear are identified. Temperature must be taken into account in wear evaluation, therefore thermal effects have been included in the numerical algorithm of wear of rock cutting. Discrete element model was extended into analysis of transient heat flow. Thermal analysis can be carried out along with the solution of mechanical problem. Thermal equations are coupled with equations of motion yielding a set of equations describing thermomechanical problem. 6.1 Basic concepts of wear Wear is the process of progressive loss of material from the surface of a solid body due to mechanical action, i.e. the contact and relative motion against a solid, liquid or gaseous counterbody. There are different wear mechanisms. These mechanisms can be classified into the following four basic groups, cf. [22]: •adhesive wear, •abrasion, •surface fatigue, •tribochemical reaction. Adhesive wear is the process of wear accompanied by formation and breaking of interfacial adhesive bonds. Adhesive wear can occur when two surfaces slide along each other (Fig. 6.1a). High local pressure between contacting asperities results in plastic deformation, adhesion and formation of local junctions. Breaking of these junctions due to shear leads to wear. Abrasion is the process of wear characterised with removal or displacement of material at a solid surface due to presence of hard particles in between or embedded in one or both of two solid surfaces in relative sliding motion. If 51 Figure 6.3: Set-up of the laboratory wear tests Table 6.2: Results of the laboratory wear tests for steels MET 71 distance (m) 100 200 cutting velocity (m/s) 0.117 0.117 normal force (N) 885 840 contact surface (mm2) 26 30 contact pessure (MPa) 34 28.2 wear depth (mm) 0.4 2 temperature (◦C) 156 146 Table 6.3: Results of the laboratory wear tests for steels MET 91 distance (m) 100 200 cutting velocity (m/s) 0.117 0.117 normal force (N) 622 1237 contact surface (mm2) 68 138 contact pessure (MPa) 9.1 9.0 wear depth (mm) 4 6 temperature (◦C) 128 213 or k=∆w pn∆sH(T) = 1.6·10−3 28.2·106·100 ·494 ·106·10 = 0.0028 (6.13) 58 An average value of wear coefficient k= 0.0025 and hardness given in Table 6.1 define the wear of steel MET 71 calculated from the Archard extended formula 6.7 and will be used in numerical simulations presented later. Using the data for steel MET 91 from Tables 6.3 and 6.1 we can calculate wear coefficient kin a similar way as above k=∆w pn∆sH(T) = 2·10−3 9.1·106·100 ·512 ·106·10 = 0.0113 (6.14) or k=∆w pn∆sH(T) = 2·10−3 9.0·106·100 ·509 ·106·10 = 0.01125 (6.15) 59 Chapter 7 Discrete Element formulation for thermal and thermo-mechanical problem 7.1 Heat balance equation Evaluation of wear requires determination of forces of cutting as well as temperature distribution. This means necessity to analyse rock cutting as a thermo-mechanical process. Temperature increases due to heat generated by friction between tool and rock. Heat is adsorbed and conducted by the tool and rock. These processes depend on thermal properties (heat capacity, thermal conductivity) of the tool and rock. Thermal phenomena during rock cutting are described by the heat balance equation. This equation can be written for a single particle in the following form: mic˙ Ti=Qi,(7.1) where the following notation has been used: c– solid heat capacity, Ti– particle temperature; Qi– heat sources or heat fluxes per single particle. Qiincludes externally supplied heat source Qext, heat generated through friction dissipation and adsorbed by the particle Qgen, heat conducted through the contact interface with another material Qcontact, heat conducted to particles of the same material Qcond and convective and radiative heat transfer between particles and environment on the free surface, Qconvection and Qradiation Qi= (Qext +Qgen)i−(Qcond +Qcontact +Qconvection +Qradiation)i(7.2) 60 It should be remembered that Qcond contains contributions from all the neighbouring particles which are in contact with the i-th particle, similarly Qgen gathers contributions from all the neighbouring particles in contact with the i-th particle. The term Qcontact is included for the particles on the contact interface, and terms Qconvection and Qradiation are evaluated for the particles on the free surface. Particle-to-particle conductive heat transfer rate Qcond is estimated as the transfer through an equivalent bar element of length dequal to the distance between the particle centres and of certain equivalent area ¯ A(function of particle radii) Qcond =hcond(Ti−Tj) = κ¯ A d(Ti−Tj) (7.3) with hcond being the heat transfer coefficient between material particles, κbeing the solid heat conductivity, and Tiand Tjbeing contacting particles temperatures. 7.2 Thermal boundary conditions Thermal boundary conditions can be specified for temperatures Tor for heat fluxes/sources Q. Thus the boundary Γ can be split into two disjoint parts, one with prescribed temperature ΓTand the other with prescribed heat flux ΓQ Γ = ΓT∪ΓQ(7.4) ΓT∩ΓQ=∅(7.5) Heat fluxes can be specified in different ways: 1. Point heat flux (representing either heat source or sink). It is represented in Eq. (7.2) by the term Qext. A special case of the prescribed heat flux is the case of insulation – with all external fluxes set to zero Qext = 0 , Qcontact = 0 , Qconvection = 0 , Qradiation = 0 (7.6) 2. Convection to the environment (on the free surface, represented by the term Qconvection) 3. Radiation to the environment (on the free surface, represented by the term Qradiation) 4. Heat transfer between contacting bodies (at the contact interface, represented by the term Qcontact). Application of the thermal boundary conditions for the free surface and contact interface requires detection of these parts of boundary in the discrete element model. The boundary itself in the model of particles changes when fractures appear in the material 61 and some particle become new surface particles. In the model of free particles (without cohesive bonds) the motion of the particles leads to change of the free boundary. In the wear analysis is carried out modification of the tool shape also leads to change of free surface. These features of the particle model make it necessary to update free surface definition, which can only be effective by employing an automatic procedure. Parts of the boundaries in contact with another body must also be detected, this however does not mean additional effort, since the detection of the contact interface is done in the solution of the mechanical problem. It is assumed that a given particle can be either on the contact interface or on the free surface. Therefore for particles with detected contact excluded from the free surface. Heat generation through frictional dissipation is calculated using the following formula Qgen =χ|Ffric vrT |,(7.7) where Ffric is the friction force, vrT is the relative tangential velocity, and χis the part of the friction work converted to heat. 7.3 Solution of thermo-mechanical problem The thermo-mechanical problem is solved by coupling the heat balance equation (7.1) with equations of motion (3.1) and (3.2). Equations (3.1), (3.2) and (7.1) describing a thermo-mechanical problem are integrated in time using a staggered scheme – solution at the n-th time step can be summarised as follows: (i) equations of motion (3.1) and (3.2) are integrated in time using a central difference scheme given by Eqs. (3.3)–(3.10); in the solution of mechanical problem frictional dissipation and heat generation is calculated according to Eq. (7.7); wear is estimated using the material hardness at a given temperature (ii) heat balance equation (7.1) is integrated in time using the explicit forward Euler scheme Tn+1 i=Tn+1 i+∆t micQn i,(7.8) Thermal problem is solved at the fixed updated geometrical configuration and the heat rate generation through frictional dissipation being supplied by the solution of the mechanical problem. Heat generation through frictional dissipation is assumed to be absorbed equally by the two particles in contact. Qgen gathers contributions from all the neighbouring particles in contact with the i-th particle. The particle temperatures obtained in the thermal step solution are in turn passed to the mechanical problem modifying material properties. 62 Explicit integration in time used in the solution of mechanical and thermal equations yields high computational efficiency. This solution scheme has been implemented in the in-house explicit dynamic code Simpact [27]. 7.4 Material thermal properties Materials considered in the model of rock cutting are steel and rock/soil. Thermomechanical analysis of rock cutting process requires knowledge of the following thermal properties of the materials such as their heat capacity and heat conductivity as well as heat transfer coefficients for the contact between rock and tool and for the convection at free surface of the tool and rock. The following thermal properties have been assumed from the literature: •steel –heat capacity c= 450 J/(kg·K), –heat conductivity κ= 60 W/(m·K), (low carbon steel 66.9 W/(m·K) –density ρ= 7830 kg/m3 –(specific heat c= 0.12 cal/(g·K)= 0.12·4.1868·1000 W/(m·K) = 502 W/(m·K) •sandstone –heat capacity c= 1970 J/(kg·K), –heat conductivity κ= 0.0125 cal/(s·cm·K) = 0.0125 ·4.1868 ·100 W/(m·K)= 5.2 W/(m·K) –density ρ= 2500 kg/m3 •sandy soil –heat capacity c = 1004 J/(kg·K), –heat conductivity κ= 0.58 W/(m·K), –density ρ= 1820 kg/m3 •marble –heat capacity c = 860 J/(kg·K), –heat conductivity κ= 3 W/(m·K), –density ρ= 2700 kg/m3 63 •tuff –heat capacity c= 950 J/(kg·K), –heat conductivity κ= 0.5÷2.5 W/(m·K) (κ= 2 W/(m·K) –density ρ= 2500 kg/m3 •sand (dry) –heat capacity c= 800 J/(kg·K), –heat conductivity κ= 0.35 W/(m·K), –density ρ= 1600 kg/m3 1cal = 4.1868 J 7.5 Parameters characterizing heat transfer on the surface Heat to be removed from the surface or the heat to be absorbed by the surface is a function of the temperature difference between the surface and surrounding media and heat transfer coefficient dependent on surrounding media (air, water) and heat transfer process (convection, radiation): •Radiation The heat transfer coefficient depending on radiation is assumed to be about 60 W/(m2K) for a metal surface temperature of 1000 ◦C. In our applications temperature is usually lower. We will consider radiation jointly with convection. •Convection Typical values of the heat transfer coefficient for convection with air or gas are typically h= 10÷100 W/(m2K). The heat transfer coefficient depending on natural convection is of the magnitude of 30 W/(m2K). Forced air-cooling may increase the convection to ≈100 W/(m2K). •Water cooling In dredging we will have heat transfer from the tooth surface to surrounding water. Depending on the flow and pressure of water, it is possible to reach values between 5000 and 50000 W/(m2K). Normal values for the water cooling are 5000 ÷25000 W/(m2K). Values above 25000 W/(m2K) are reached in special cooling equipment. 64 controlvolume V V i j(i) Figure 7.1: Definition of the control volume for density evaluation Density evaluation To evaluate the density map using the discrete element method an algorithm has been developed. For each particle a control volume is defined as it is shown in Fig. 7.1. Then an average density is defined according to the following expressions: ¯ρ(i)=ρVc−V(i) 0 Vc , V (i) 0=Vc−Vi−X j ¯ V(i) j.(7.9) In the computation of the average density associated to a particle, the intersection of the control volume with the volumes of interacting particles must be computed. An exact analytical expression is used to compute these volumes. 7.6 Benchmarks of thermal analysis using DEM Benchmarks chosen to test thermal DEM formulation are test examples for finite element thermal analysis program DOT [28]. 7.6.1 Internal heat generation in parallel-sided slab An 8-inch thick slab, infinite in extent, is initially at 0◦F. Starting at time t= 0 heat is generated internally in the slab at a uniformly distributed rate ˙q= 2000 BTU/(sec·in3). The external surfaces of the slab are maintained at 0◦F for all time. The material properties are as follows: •specific heat c= 1 BTU·in/(sec2lb·◦F), •heat conductivity κ= 16 BTU/(sec·m·◦F), •mass density ρ= 1 sec2lb/in4 The finite element model used taking advantage of the plane symmetry is shown in Fig. 7.2. Comparison between FEM and DEM solutions is shown in Fig. 7.3 for the 65 temperature distributions. A very good agreement is shown. The DEM model with temperature distributions for different time instants is shown in Fig. 7.4. 4’’ x y planeof symmetry 1’’ Figure 7.2: Internal heat generation in parallel-sided slab - finite element model 0 100 200 300 400 500 600 0 0.5 1 1.5 2 2.5 3 3.5 4 temperature x distance (inches) DEM solution, t=0.08 s DEM solution, t=0.16 s DEM solution, t=0.24 s DEM solution, t=0.32 s FEM solution, t=0.08 s FEM solution, t=0.16 s FEM solution, t=0.24 s FEM solution, t=0.32 s Figure 7.3: Internal heat generation in parallel-sided slab – temperature distribution – comparison of FEM and DEM results 66 a) t = 0.08 s b) t = 0.16 s c) t = 0.24 s d) t = 0.32 s Figure 7.4: Internal heat generation in parallel-sided slab - temperature distribution at different time instants (discrete element method results) 67 Chapter 8 Numerical simulation of rock cutting and wear evaluation 8.1 2D mechanical simulation of rock cutting – model I Figure 8.1: Initial set-up of 2D model of rock cutting – tool discretized with straight segments 2D mechanical simulation of rock cutting has been carried out using a model shown in Fig. 8.1. Material sample 109 ×109 mm is represented by an assembly of randomly compacted 2100 discs of radii 1–1.5 mm. Other model parameters are as follows: contact stiffness in the normal and tangential directions kn=kT= 20 GPa, the cohesive bond strengths Rn= 0.1 MN/m, RT= 1 MN/m, and the friction coefficient µ= 0.839. The surface of the rigid tool has been modelled with straight segments. The following parameters have been assumed for the tool-rock interface: contact stiffness modulus kn= 20 GPa, Coulomb friction coefficient µ= 0.839. 74 Cutting has been carried out with prescribed horizontal velocity of the tool 0.04 m/s. Process of cutting is shown in Fig. 8.2. In this figure we can also see the failure mode – particles with broken bonds are coloured in blue. It can be clearly seen formation of a chip during cutting typical for brittle materials. Initial wear pattern has been obtained in the analysis. The wear distribution calculated for one work cycle is shown in Fig. 8.3. Calculated wear was obtained using Eq. (6.1) with assumed parameters ¯ k= 1 and H= 1 Pa. Obtained values although not real show distribution of wear on the surface — zones with higher wear are indicated. 8.2 2D mechanical simulation of rock cutting (model II) 2D simulation of rock cutting has been carried out using the model shown in Fig. 8.4. Material sample is the same as in model I, and the rigid tool is now modelled with distinct elements. The tool is modelled with 4649 equal particles of radius 0.5 mm (average radius of particles modelling rock 1.25 mm). The following parameters have been assumed for the tool-rock interface: contact stiffness modulus kn= 20 GPa, Coulomb friction coefficient µ= 0.839. Cutting has been carried out with prescribed horizontal velocity of the tool 0.04 m/s. Analysis has been carried out without modification of tool shape. Process of cutting is shown in Fig. 8.5. Failure mode is presented in Fig. 8.5 by colouring blue particles with broken bonds. Formation of a chip in the initial phase of the cutting can also be seen in this simulation. The same model and the same process parameters have been used for the simulation of cutting with modification of tool shape due to wear. Wear has been calculated according to Eq. (6.1) using the constants ¯ k= 5 ·10−9and H= 1 Pa, which produced accelerated wear. This allowed us to test our algorithm in short simulation. Tool shape was modified by removal of particles from the surface when the accumulated wear (thickness) exceeded their diameter. Process of cutting and failure mode with formation of a chip are shown in Fig. 8.6 for different time instants. Profile of wear on the surface of the tool is shown for different stages in Fig. 8.7. Some of the particles on the cutting tip have been eliminated. This can be seen better in detail presented in Fig. 8.8. 75 a) t = 0.2 s b) t = 0.4 s c) t = 0.6 s d) t = 0.8 s e) t = 1.0 s f) t = 1.2 s Figure 8.2: Process of cutting - failure mode 76 Figure 8.3: Distribution of wear on the tool surface at the end of work cycle (model I) Figure 8.4: Initial set-up of 2D model of rock cutting – tool discretized with distinct elements 77 a) t = 0.0 s b) t = 0.1 s c) t = 0.2 s d) t = 0.3 s e) t = 0.4 s f) t = 0.5 s Figure 8.5: Process of cutting with failure mode (model II, tool shape without change) 78 a) t = 0.0 s b) t = 0.2 s c) t = 0.4 s d) t = 0.6 s e) t = 0.8 s f) t = 1.0 s Figure 8.6: Process of cutting with failure mode (model II, tool shape changed) 79 a) t = 0.2 s b) t = 0.4 s c) t = 0.6 s d) t = 0.8 s e) t = 1.0 s Figure 8.7: Process of cutting with profile of wear (model II, tool shape changed) 80 Figure 8.8: Profile of wear at time instant t = 1.0 s (model II, tool shape changed) – detail 81 Chapter 9 Wear evaluation in multi-cycle analysis of rock cutting 9.1 Methodology The objective of the present study is determination of the wear of a cutting tool due to interaction with sandstone during the excavation process. To disc assembly modelling the sandstone sample and micro-scale parameters have been previously determined. The cutting tool travels at a constant prescribed velocity of 4m/s. The wear process have been accelerated, so that lower computational time have been needed. First, we excavate the material during a brief lapse of time (0.005 s). The final configuration of this stage (Fig. 9.1) has been used as the initial configuration of the following stages. At the end of the stages we move back the cutting tool and continue excavating in the next stage. Figure 9.1: Configuration at t= 0.005 s from the beginning of the analysis 82 9.2 Simulation of wear of the test tooth 9.2.1 Description of the model The particles that form the cutting tool have been radii randomly generated in the range between 0.4 and 0.52 mm. The tooth material has been assumed steel MET 91. The constitutive parameters of the contact interface between spheres and the cutting tools are the following: •contact stiffness in the normal direction kn= 90 GPa , •Coulomb friction coefficient µ= 0.3, •wear constant kand steel hardness in the ambient temperature H(T= 20◦C) as given in Sec. 6.6. The simulation has been carried out using the constant tooth hardness. The oscillations with highest frequencies are damped out by adequate damping at the contact interface between the spheres that conform the material and the cutting tool (90% of the critical damping). We also apply a global viscous damping so that the lowest vibration modes decrease. 9.2.2 Results Tooth shapes at different stages of wear are shown in Fig. 9.2. During 0.203 seconds of excavation 40% of the tool mass has been lost. Figure 9.3 presents the tool mass loss as a fraction of the initial tool mass versus analysis time. Scaling analysis time by 7000 allows us to obtain the wear curve, which agrees with experimental curves for a ripper tooth, cf. Fig. 9.4. 83 a) t = 550 s b) t = 600 s c) t = 1300 s d) t = 1400 s Figure 9.9: Process of cutting at different stages of wear 90 a) t = 550 s b) t = 600 s c) t = 1300 s d) t = 1400 s Figure 9.10: Process of cutting — failure mode at different stages of wear 91 a) t = 550 s b) t = 600 s c) t = 1300 s d) t = 1400 s Figure 9.11: Process of cutting — wear distribution on the tool surface 92 Figure 9.12: Comparison of experimental and simulation wear curves for a ripper tooth 93 Chapter 10 Thermo-mechanical transient analysis of rock cutting 10.1 2D thermo-mechanical simulation of rock cutting Figure 10.1: Initial set-up of 2D model of rock cutting 2D thermo-mechanical simulation of rock cutting has been carried out using a model shown in Fig. 10.1. Material sample 109 ×109 mm is represented by an assembly of randomly compacted 2100 discs of radii 1–1.5 mm. Other model parameters are as follows: contact stiffness in the normal and tangential directions kn=kT= 20 GPa, the cohesive bond strengths Rn= 0.1 MN/m, RT= 1 MN/m, and the friction coefficient µ= 0.839. The rigid tool is modelled with distinct elements, 4649 equal particles of radius 0.5 mm (average radius of particles modelling rock 1.25 mm). The following parameters have been assumed for the tool-rock interface: contact stiffness modulus kn= 20 GPa, Coulomb friction coefficient µ= 0.839. The following thermal 94 properties have been assumed: •for rock heat capacity c= 3000 J/(kg·K), heat conductivity κ= 0.1 W/(m·K) •for steel heat capacity c= 450 J/(kg·K), heat conductivity κ= 60 W/(m·K) Cutting has been carried out with prescribed horizontal velocity of the tool 0.04 m/s. Process of cutting is shown in Fig. 10.2. In this figure we can also see the failure mode – particles with broken bonds are coloured in blue. It can be clearly seen formation of a chip during cutting typical for brittle materials. Wear has been calculated according to Eq. (6.9) using the constants ¯ k= 10−10,H= 100 Pa. Little wear was produced in one cycle. The objective example of the example was to test thermal analysis including heat generation and conductance. Two cases has been studied. In the first case no heat conductance was allowed, heat was assumed to be adsorbed by particles which generate heat. In the second case full thermal analysis with heat conductance by the tool and rock has been carried out. Temperature distribution at different time instants for the case without conductance is shown in Figs. 10.3 and 10.4. Temperature distribution at different time instants for the full thermal analysis is presented in Fig. 10.5. We can see increase of temperature in the zone of contact between the tool and rock. 95 a) t = 0.002 s b) t = 0.004 s c) t = 0.006 s d) t = 0.008 s Figure 10.2: Process of cutting - failure mode 96 a) t = 0.002 s b) t = 0.004 s c) t = 0.006 s d) t = 0.008 s Figure 10.3: Process of cutting – temperature distribution obtained in the solution without heat conductance a) t = 0.002 s b) t = 0.004 s c) t = 0.006 s d) t = 0.008 s Figure 10.4: Process of cutting – temperature distribution obtained in the solution without heat conductance (zoom) 97 a) t = 0.002 s b) t = 0.004 s c) t = 0.006 s d) t = 0.008 s Figure 10.5: Process of cutting – temperature distribution (full thermal analysis) 98 Chapter 11 Thermomechanical analysis of wear of a ripper 11.1 Model definition Thermo-mechanical numerical simulation of wear of a ripper tooth has been carried out Geometry of the rock sample and tooth as well as mechanical parameters have been assumed similar to those used in Sec. 9.3. Material sample (Fig. 9.5) of dimensions 109×109 mm is represented by an assembly of randomly compacted 2100 discs of radii 1 ÷1.5 mm. Mechanical model parameters are as follows: contact stiffness in the normal and tangential directions kn=kT= 90 GPa, the cohesive bond strengths Rn= 0.1 MN/m, RT= 1 MN/m, and the friction coefficient µ= 0.839. The material macroscopic behaviour corresponds to sandstone properties which has been found in previous numerical tests. The rigid tool is modelled with distinct elements, 6500 disc of radii randomly varying from 0.4 to 0.6 mm. The following parameters have been assumed for the tool-rock interface: contact stiffness modulus kn= 20 GPa, Coulomb friction coefficient µ= 0.3. The tooth material has been assumed steel MET 91. Temperature dependent steel hardness and wear coefficient given in Sec. 6.6 have been used. The thermal properties of rock (sandstone) and steel have been assumed, cf. Sec. 7.4: •for rock heat capacity c= 1970 J/(kg·K), heat conductivity κ= 5 W/(m·K) •for steel heat capacity c = 450 J/(kg·K), heat conductivity κ= 60 W/(m·K) 99 a) t = 5 s b) t = 15 s c) t = 90 s d) t = 100 s Figure 11.5: Process of cutting — tool shape evolution and temperature distribution in the tool (model II) 106 11.4.2 Model II – initial temperature distribution from quasistationary thermal analysis Wear analysis has been carried out starting with initial temperature distribution determined through quasi-stationary thermal analysis discussed in Sec. 11.3. Cutting velocity 0.4 m/s has been assumed and insulation at the tooth-crossection – corresponding initial temperature distribution is shown in Fig. 11.4d). Wear has been accelerated as accelerating wear 7000 times. Tooth shape and temperature distribution for different stages of wear are shown in Fig. 11.6. a) t = 5 s b) t = 15 s c) t = 90 s d) t = 100 s Figure 11.6: Process of cutting — tool shape evolution and temperature distribution in the tool (model I) 107 Chapter 12 Simulation of dredging 12.1 Model I 12.1.1 Model definition Numerical simulation of dredging with wear of a dredge tooth has been carried out using a 2D model shown in Fig. 12.1. Material sample is represented by an assembly of randomly compacted 92000 discs of radii 1–1.5 mm. Model parameters obtained previously for sandstone are assumed for the micromechanical model. Figure 12.1: 2D model of dredging – initial set-up The rigid tool is modelled with distinct elements, 5400 discs of equal radius of 1 mm. The following parameters have been assumed for the tool-rock interface: contact stiffness modulus kn= 50 GPa, Coulomb friction coefficient µ= 0.4. The tooth material has been assumed steel MET 91. Wear analysis has been carried out with temperature dependent hardness and wear coefficient estimated in Sec. 6.6. The 108 simulation has been carried out using the constant tooth hardness. Wear has been accelerated 7000 times by scaling the wear coefficient k. The swing velocity was assumed 0.2 m/s, and rotating velocity 1.6204 s−1, which with the distance of the tooth from the axis of rotation 0.7 m gives circumferential velocity 1.134 m/s. Mechanical process has been scaled 5 times. 12.1.2 Thermo-mechanical analysis of dredging Process of dredging has been analysed as thermo-mechanical process with subsequent cutting and cooling. During cutting the tooth was heated and then during the rotation in water the tooth was cooled. The tooth surface was cooled by water during cutting as well. The heat diffusion through the tooth material was analysed. Analysis presented in this section has been carried out with an unchanged tooth shape (zero wear). Analysis results are shown in Figs. 12.2–12.5. Failure of rock during dredging is shown in Fig. 12.2. Temperature distribution at different instants of the first tooth pass of cutting is shown in Fig. 12.3 and temperature distribution during cooling with water is presented in Fig. 12.4. Cyclic thermal load is illustrated by the curve representing change of maximum tooth temperature during subsequent cutting and cooling (Fig. 12.5). 109 a) t = 0.01 s b) t = 0.2 s c) t = 0.47 s d) t = 0.57 s Figure 12.2: Thermo-mechanical simulation of dredging – failure mode 110 a) t = 0.01 s b) t = 0.2 s c) t = 0.47 s d) t = 0.57 s Figure 12.3: Thermo-mechanical simulation of dredging – temperature distribution at different instants of the first tooth pass of cutting 111 a) t = 0.63 s b) t = 0.93 s c) t = 1.23 s d) t = 1.83 s Figure 12.4: Thermo-mechanical simulation of dredging – temperature distribution at different instants of cooling after the first tooth pass of cutting 112 20 30 40 50 60 70 80 90 0 1 2 3 4 5 6 temperature time (sec) Figure 12.5: Thermo-mechanical simulation of dredging – variation of maximum tooth temperature during three cycles of cutting and subsequent cooling 113 12.1.3 Thermo-mechanical analysis of dredging with wear evaluation Similar thermo-mechanical analysis of dredging has been carried out with extension to wear evaluation and modification of the tooth shape due to wear. Analysis results are shown in Figs. 12.6–12.12. Failure of rock during dredging is shown in Fig. 12.6. Temperature distribution at different instants of the first tooth pass of cutting is shown in Fig. 12.7 and temperature distribution during cooling with water is presented in Fig. 12.8. Cyclic thermal load is illustrated by the curve representing change of maximum tooth temperature during subsequent cutting and cooling (Fig. 12.9). Figure 12.10 shows accumulated wear on the tooth surface at different instants of the first tooth pass of cutting. Evolution of tooth shape due to wear is shown in Fig. 12.11. Wear is quantified by the curve showing the tooth mass loss in Fig. 12.12. a) t = 0.01 s b) t = 0.2 s c) t = 0.52 s d) t = 0.64 s Figure 12.6: Thermo-mechanical simulation of dredging with wear evaluation – failure mode 114 a) t = 0.01 s b) t = 0.2 s c) t = 0.52 s d) t = 0.64 s Figure 12.7: Thermo-mechanical simulation of dredging with wear evaluation – temperature distribution at different instants of the first tooth pass of cutting 115 pass is shown in Fig. 12.15 and temperature distribution during cooling with water is presented in Fig. 12.16. Figure 12.17 shows accumulated wear on the tooth surface at different instants of the first tooth pass of cutting. Change of the tooth shape due to wear after one pass of cutting is shown in Fig. 12.18. a) t = 0.01 s b) t = 0.2 s c) t = 0.4 s d) t = 0.65 s Figure 12.14: Thermo-mechanical simulation of dredging with wear evaluation – failure mode (model II) 122 a) t = 0.01 s b) t = 0.2 s c) t = 0.4 s d) t = 0.65 s Figure 12.15: Thermo-mechanical simulation of dredging with wear evaluation – temperature distribution at different instants of the first tooth pass of cutting (model II) 123 a) t = 0.68 s b) t = 0.78 s c) t = 1.48 s d) t = 2.00 s Figure 12.16: Thermo-mechanical simulation of dredging with wear evaluation – temperature distribution at different instants of cooling after the first tooth pass of cutting (model II) 124 a) t = 0.01 s b) t = 0.2 s c) t = 0.4 s d) t = 0.65 s Figure 12.17: Thermo-mechanical simulation of dredging with wear evaluation – accumulated wear on the tooth surface at different instants of the first tooth pass of cutting (model II) 125 a) before 1st pass b) after 1st pass Figure 12.18: Thermo-mechanical simulation of dredging with wear evaluation – evolution of tooth shape due to wear (model II) 126 Chapter 13 3D simulation of rock cutting 13.1 3D simulation of interaction of a tooth with granular medium 3D simulation of interaction of a tooth with granular medium has been carried out using a model shown in Fig. 13.1a. Material is represented by a collection of 11700 spherical particles of radii 9 mm placed randomly in a rectangular box (not presented in figure) and subjected to gravitational loading. The surface of the rigid tool has been modelled with triangular facets. In the contact among particles the following properties have been assumed: contact stiffness in the normal direction kn= 106N/m, friction coefficient µ= 0.1, no cohesion existed among particles. Friction between the tooth and material particles was defined by the Coulomb friction coefficient µ= 0.1, no wear calculation has been done. The tooth has been lowered, and after penetration the horizontal velocity 0.04 m/s has been prescribed. Interaction of the tooth with material particles is shown in Fig. 13.1b. 127 a) t = 0.0 s b) t = 1.1 s Figure 13.1: 3D simulation of interaction of a tooth with granular medium 128 13.2 3D simulation of rock cutting (model I) The same as above initial geometry of the tooth and particles (Fig. 13.2a) has been used to create a 3D model of rock cutting. Rock-like material properties have been obtained by introduction of cohesive bonds among contacting particles. The contact interface properties were as follows: contact stiffness in the normal and tangential directions kn=kT= 106N/m, cohesive bond strengths Rn= 500 N/m, RT= 5000 N/m, friction coefficient µ= 0.5. Friction between the tooth and material particles was defined by the Coulomb friction coefficient µ= 0.1, no wear calculation has been done. The process of cutting has been defined similarly as previously, the tooth has been lowered, and after penetration the horizontal velocity 0.04 m/s has been prescribed. Interaction of the tooth with material particles is shown in Fig. 13.2b. Formation of a chip can be seen in this figure. a) t = 0.0 s b) t = 1.1 s Figure 13.2: 3D simulation of rock cutting (model I) 129 13.3 3D simulation of rock cutting (model II) The same as above model of rock sample (Fig. 13.3a) has been used. The same contact interface properties were assumed: contact stiffness in the normal and tangential directions kn=kT= 106N/m, cohesive bond strengths Rn= 500 N/m, RT= 5000 N/m, friction coefficient µ= 0.5. The tool model in the present model has been composed of two rigid parts, one part has been modelled with rigid surfaces and the other one has been discretized with 28570 distinct elements of radii 1.5 mm (considerably smaller than the radius of rock particles, 9 mm). The small radius of particles discretizing the tool allows us to assume that the interaction of rock particles is similar like interaction with flat surface. Friction between the tooth and material particles was defined by the Coulomb friction coefficient µ= 0.1. The process of cutting has been defined similarly as previously, the tooth has been lowered, and after penetration the horizontal velocity 0.4 m/s has been prescribed. Interaction of the tooth with material particles is shown in Fig. 13.3b. Formation of a chip can be seen in this figure. a) t = 0.0 s b) t = 1.0 s Figure 13.3: 3D simulation of rock cutting (model II) Wear has been calculated according to Eq. (6.1) using the constants ¯ k= 5 ·10−9 130 and H= 1 Pa, Profile of wear of the tooth is presented in Fig. 13.4. Maximum value of wear, 0.14186 is smaller than particle diameter, 3 mm, so the tool shape was not modified. Figure 13.4: 3D simulation of rock cutting (model II) - profile of wear on the tooth at time t = 1 s 131 Acknowledgement Financial support of the European Commission through the Growth European Project GRD1-2000-25243 “Shortening Lead Times and Improving Quality by Innovative Upgrading of the Lost Foam Casting Process (FOAMCAST)” is gratefully acknowledged. 138 Bibliography [1] P.A. Cundall and O.D.L. Strack. A discrete numerical method for granular assemblies. Geotechnique, 29:47–65, 1979. [2] J. Rojek, E. Oate, F. Zarate, and J. Miquel. Modelling of rock, soil and granular materials using spherical elements. In 2nd European Conference on Computational Mechanics ECCM-2001, Cracow, 26-29 June, 2001. [3] K. Gehring. Rock testing procedures at VA’s geotechnical laboratory in Zeltweg. Technical report, Voest Alpine Zeltweg, Austria, TZU 41, 1987. [4] B.N. Whittaker, R.N. Singh, and G. Sun. Rock fracture mechanics. Amsterdam, 1987. [5] I. Evans. The force required for pointed attack picks. [6] Y. Nishimatsu. The mechanics of rock cutting. [7] C.S. Campbell. Rapid granular flows. Annual Review of Fluid Mechanics, 2:57–92, 1990. [8] Eng. Comput., 9(2), 1992. Special issue, Editor: G. Mustoe. [9] J.R. Williams and R. O’Connor. Discrete Element Simulation and the Contact Problem. Archives Comp. Meth. Engng, 6(4):279–304, 1999. [10] P.A. Cundall. Formulation of a Three Dimensional Distinct Element Model — Part I. A Scheme to Detect and Represent Contacts in a System of Many Polyhedral Blocks. Int. J. Rock Mech., Min. Sci. & Geomech. Abstr., 25(3):107–116, 1988. [11] J. Argyris. An excursion into large rotations. Comput. Meth. Appl. Mech. Eng., 32:85–155, 1982. [12] D.J. Benson and J.O. Hallquist. A simple rigid body algorithm for structural dynamics programs. Int. J. Num. Meth. Eng., 12:723–749, 1986. 139 [13] L.M. Taylor and D.S. Preece. Simulation of blasting induced rock motion. Eng. Comput., 9(2):243–252, 1992. [14] T.J.R. Hughes. The Finite Element Method. Linear Static and Dynamic Analysis. Prentice-Hall, 1987. [15] R. L¨ohner. Applied CFD techniques. An Introduction based on Finite Element Methods. Wiley, 2001. [16] R. L¨ohner and K. Morgan. An unstructured multigrid method for elliptic problems. Int. J. Num. Meth. Eng., 24:101–115, 1987. [17] D.N. Knuth. The Art of Computer Programming. Addison-Wesley, Reading, Mass., 1973. [18] J. Bonet and J. Peraire. An alternate digital tree algorithm for geometric searching and intersection problems. Int. J. Num. Meth. Eng., 31:1–17, 1991. [19] H. Samet. The quad-tree and related hierarchical data structures. Comput. Surveys, 16(2):187–285, 1984. [20] F.L. Greengard. The Rapid Evaluation of Potential Fields in Particle Systems. ACM Distinguished Dissertation 1987. ACM Distinguished Dissertation Series, 1987. [21] H. Huang. Discrete Element Modeling of Tool-Rock Interaction. PhD thesis, University of Minnesota, 1999. [22] K.H. Zum Gahr. Microstructure and wear of materials. Amsterdam, 1987. [23] P.N.W. Verhoef. Wear of rock cutting tools. Balkema, Rotterdamd, 1997. [24] J.F. Archard. Contact and rubbing of flat surfaces. J. Appl. Phys., 24(8):981–988, 1953. [25] E. Rabinowicz. Friction and wear of materials. 1995. [26] S. Stupkiewicz and Z. Mr´oz. A model of third body abrasive friction and wear in hot metal forming. [27] SIMPACT. User Manual. A finite element code for structures under dynamic and impact loadings. Report No. 236, CIMNE, Barcelona, 1997. [28] R.M. Polivka and E.L. Wilson. Finite Element Analysis of Nonlinear Heat Transfer Problems. UC SESM 76-2, 1976. 140