Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 16, N o 1, 2018, pp. 41 - 50 https://doi.org/10.22190/FUME171229010D © 2018 by University of Niš, Serbia | Creative Commons Licence: CC BY-NC-ND Original scientific paper SIMULATION OF FRACTURE USING A MESH-DEPENDENT FRACTURE CRITERION IN THE DISCRETE ELEMENT METHOD UDC 539.4, 519.6 Andrey Dimaki 1 , Evgeny Shilko 1 , Sergey Psakhie 1 , Valentin Popov 2 1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, Russia 2 Berlin University of Technology, Berlin, Germany Abstract. Recently, Pohrt and Popov have shown that for simulation of adhesive contacts a mesh dependent detachment criterion must be used to obtain the mesh-independent macroscopic behavior of the system. The same principle should be also applicable for the simulation of fracture processes in any method using finite discretization. In particular, in the Discrete Element Methods (DEM) the detachment criterion of particles should depend on the particle size. In the present paper, we analyze how the mesh dependent detachment criterion has to be introduced to guarantee the macroscopic invariance of mechanical behavior of a material. We find that it is possible to formulate the criterion which describes fracture both in tensile and shear experiments correctly. Key Words: Fracture, Mesh-Dependence, Discrete Element Method, Particle Size 1. INTRODUCTION Since the work of Hertz [1] it is well-known that stress distribution in a contact area is not uniform. Stress distribution in the contact between an elastic half-space and sharp- edged rigid or elastic counter-bodies of a different shape (e.g. a cylinder, frustum etc.) tends to infinity in a vicinity of the contact area boundary [2, 4-6]. In reality, stresses at the contact area boundary never achieve infinite values due to surface roughness and/or plastic deformation [3]. Nevertheless, these stresses are several times higher than in the center of the contact area. The presence of high stress concentrations near contact patch boundaries, notches etc. leads to nucleation of cracks on different scales [7, 8]. In order to describe conditions of fracture in the regions with heterogeneous stress/strain fields, non-local fracture criteria Received December 29, 2017 / Accepted February 05, 2018 Corresponding author: Andrey V. Dimaki Institute of Strength Physics and Materials Science SB RAS, 634055, Tomsk, Russia dav@ispms.tsc.ru 42 A. DIMAKI, E. SHILKO, S. PSAKHIE, V. POPOV have been developed. Perhaps one of the most popular nonlocal fracture criteria is Novozhilov fracture criterion [9] which reads 1 1 0 1 ( ) d y c x dx d    , (1) where σc is the strength of a homogeneous material under uniaxial tension without stress concentrators, σy(x) is the maximum normal stress, x axis corresponds to a “most probable” direction of crack propagation, d1 is a material constant having a dimensionality of length. Similar forms of nonlocal criterion (1) were used in [10] for composites and in [11] for notched samples. Another way of taking into account stress concentrators is to include both absolute value of stress (either maximal or average value) and its gradient in a crack processing zone: 0 ( , / ) e c f L L   , (2) where L0 is a characteristic size of material structure elements, Le is a size of stress concentration zone [12]. In the framework of this approach the size of stress concentration zone depends on a local stress gradient: Le = Le(grad(σ)). When developing a numerical model of the contact, or, more generally, of a sample with stress concentrators, it is necessary to adequately describe crack nucleation conditions since they significantly depend on stress distribution approximation accuracy that is determined, in turn, by mesh size. Correspondingly, variation of the mesh size may have influence on the values of mechanical characteristics of the simulated samples, e.g. their compressive/tensile strength, etc. The traditional way of describing heterogeneous stress/strain fields, based on mesh refinement, is not applicable to problems where “quasi-infinite” local values of stresses can occur. These problems include friction and wear, contact problems with punches having sharp edges, etc. In this case, in order to adequately describe fracture in numerical models, nonlocal fracture criteria of types (1) and (2) are widely used. Examples of their applications in the finite-element method (FEM) can be found in papers [13, 14] and in many others. Recently, Popov and Pohrt have suggested a mesh-dependent detachment criterion in the boundary element method (BEM) for a normal adhesive contact problem of linearly elastic and power-law graded elastic materials [15-17]. The main idea of the suggested approach is a modification of local ultimate stress value instead of stress field correction. The given criterion is local and it allows adequate descriptions of a normal contact between an elastic half-space and a rigid indenter in the presence of adhesive forces. Despite a wide application of non-local and mesh size dependent fracture criteria in continuum-based approaches like FEM, the size-dependent fracture criteria in DEM still remain much less studied [18]. In the paper we study the relation between tensile and shear contact strength and a discrete element size. Hereinafter the term “contact strength” means a maximal value of a force divided to the contact square that is required to detach the contacting bodies. Based on the obtained results, we suggest a mesh-dependent local fracture criterion which allows us to overcome the mentioned problem, namely to avoid dependence of the contact strength on element size. Simulation of Fracture Using a Mesh-Dependent Fracture Criterion in a Discrete Element Method 43 2. FORMULATION OF THE PROBLEM Let us consider a contact between an elastic half-space and an elastic punch, which are initially bonded (see Fig. 1). We study dependence of the contact tensile strength on element size d. A displacement with constant vertical velocity is applied to the upper layer of the punch while the lower layer of the half-space is fixed in vertical direction. Horizontal motions of elements of the top boundary of the punch and the bottom boundary of the half- space are allowed. The materials of the punch and the half-space are considered linearly elastic. We used the Movable Cellular Automaton method (MCA) [19-21] which is a DEM family representative with multi-particle distinct element formulation of the equations of elasticity and motion. Fig. 1 Schematic of the simulated sample; element size d=0.001 m The elastic modulus of the half-space was Ehalf-space = 200 GPa, the elastic modulus of the punch Epunch was varied in a certain range in order to perform a parametric study of the problem. The values of the Poisson’s ratios were νpunch = νhalf-space = 0.3. The value of the ultimate stress in von Mises fracture criterion was YMises = 200 MPa that leads to relatively small values of strains in the punch and the substrate at the moment of fracture. The strain rate was of the order of  ≈ 0.3 sec -1 that provides a quasi-static loading regime. We used a finite-size punch and a large counter-body which imitates a half-space. Evidently, the boundary conditions produce a contribution to the values of sample strength obtained during simulations. In order to diminish influence of the boundary conditions we studied the convergence of the strength value when thickness and height of a punch and a counter-body simulating a half-space increase. The values of thickness H and width W of the counter-body greater than five widths of the punch 2a0 guarantee convergence of the value of the tensile strength and, consequently, the absence of the boundaries’ influence. Punch height h>3a0 guarantees tensile strength being independent of h. The following values of the sizes of the punch and the “half-space” were used in the calculations described below: a0=0.01 m, h=0.03 m, H=0.125 m, W=0.1 m. The closed package of 44 A. DIMAKI, E. SHILKO, S. PSAKHIE, V. POPOV discrete elements was used. It should be noted that the use of a close-packed ensemble of circular discrete elements leads to the appearance of an "artificial roughness" on the boundaries of the sample (see Fig. 1). 2.1. Description of the discrete element model In the framework of MCA a solid body is considered as an ensemble of interacting particles (elements) of finite size. Contacts of interacting pairs of elements are considered to be initially bonded that simulates a consolidated solid, while the contacts between crack faces are considered as unbonded. Switching between bonded and unbounded states occurs when a given fracture or linkage criterion is satisfied. Interaction between the contacting elements (either bonded or unbonded) is realized through their common plane faces. The evolution of an ensemble of discrete elements is defined by a numerical solution of the system of Newton-Euler equations of motion: 2 2 1 1 ( ) i i N ni i i i ij ij j N i i ij j d r dv m m F F dtdt d MJ dt                , (3) where i r , i v and i  are the radius-vector, velocity vector and angular velocity of discrete element i, mi is the element mass, Ji is the moment of inertia of an equivalent disc or sphere, n ij F and ijF  are the forces of central (normal) and tangential interaction between element i and its neighbor j, and ij ij ij ijM q n F    is the torque. Equation (3) describes translational and rotational motion of the discrete elements having a finite size. Both the constituents ( n ij F and ij F  ) of the interaction force between discrete elements i and j in Eq. (3) include potential ( np ij F and p ij F  ) and dissipative/viscous ( nv ij F and v ij F  ) contributions [19]: n np nv ij ij ij p v ij ij ij F F F F F F          . (4) As follows from Eq. (4), the value of torque ijM includes both potential and dissipative constituents. The major features of the MCA method are the approximation of a homogeneous stress-strain distribution in a simply deformable element and the postulated form of the relation for the reaction force of a discrete element in response to the action of its neighboring element. In the simply deformable element approximation the state of a discrete element is determined by average stress tensor i   and average strain tensor i   [21] calculated based on the specific values of normal and tangential forces in interacting pairs of discrete elements. Tensors i   and i   are assumed to be related by an assigned constitutive law. This law determines a relationship between force and spatial interaction parameters for a pair of elements. The MCA method postulates a structural form of relations for the central and tangential interaction forces of discrete elements. We have used a generalized Hooke’s law Simulation of Fracture Using a Mesh-Dependent Fracture Criterion in a Discrete Element Method 45 for isotropic materials as constitutive relation ( ) i i     for the elastic response of a simply deformable discrete element. The developed approach allows using multi-parametric failure criteria (von Mises, Mohr–Coulomb, Drucker–Prager and others) for simulation of bond breaking between the pairs of interacting discrete elements. In this paper we used the von Mises fracture criterion in the following form: 2 2 2 21 ( ) ( ) ( ) 6 2 xx yy yy zz zz xx xy Mises Y           , (5) written for a pair of interacting discrete elements. 3. RESULTS OF SIMULATION There is a well-known analytic solution for the normal pressure distribution under an elastic flat punch contacting with an elastic half-space [3, 6]. This solution suggests the presence of peaks of normal pressure near the boundary of the contact area. Height and, more generally, shape of these peaks depend on curvature of the corners of the punch, the relation of elastic moduli of the punch and the half-space, the roughness and the coefficient of friction between the contacting bodies, etc. In any case, these peaks are much higher than the pressure on the axis of the punch. This fact leads to a curious result at the DEM simulation of detachment of a punch from a half-space. Namely, tensile strength σt of a contact between a punch and half-space becomes lower with decrease of the diameter of a discrete element in accordance with the power law: 0 0 s t t d d          , (6) where σt0 is a parameter having the dimensionality of stress, d is a discrete element size, d0 is a normalizing constant and s is the exponent. The dependencies of tensile strength on element size for different relations of elastic moduli of the punch and the half-space e = Epunch / Ehalf-space are shown in Fig. 2. Evidently, the upper limit of the element size is a punch size and the lower limit tends to zero. As one can see from Fig. 2, there is no convergence (at least, asymptotic) of the values of tensile strength on the lower limit of element size. This means the formulation of the numerical model is physically inadequate and must be improved in a way guarantying convergence of the results of simulation with decreasing of an element size. It is necessary to note the fact that the mentioned effect of element size dependence of tensile contact strength holds for different types of fracture criteria (von Mises, Drucker-Prager, detachment at a given value of normal tensile stress etc). Also this effect exists when a quadratic package of elements is used. 46 A. DIMAKI, E. SHILKO, S. PSAKHIE, V. POPOV Fig. 2 Dependencies of the contact tensile strength on the element size for the size independent local detachment criterion In order to study the influence of the direction of loading, we carried out a series of calculations in which shear loading with constant horizontal velocity was applied to the top layer of the punch, while it was fixed in vertical direction. It was found that shear strength of the contact between the punch and the half-space depends on the element size in the same manner as the tensile strength discussed above (see Fig. 3). Shear strength of the contact obeys to the power law 0 0 q shear d d          , (7) where τ0 has the dimension of stress, q is the exponent determining the slope of the logarithmical dependence of shear strength on the element size. Based on the analysis of the dependencies shown in Figs. 2 and 3 we obtained the following estimate of the exponents from Eqs. (6) and (7) for e = 1: 0.4 0.01s q   . (8) Fig. 3 Shear strength of a contact versus element size Simulation of Fracture Using a Mesh-Dependent Fracture Criterion in a Discrete Element Method 47 The estimates of s and q for different values of e are summarized in Table 1. It is evident that the highest discrepancy between the values of exponents s and q takes place for e = 10 (the stiffest punch); in other cases the difference between them does not exceed few percent. This demonstrates universality of the obtained dependencies of strength on element diameter and their applicability to construction of a mesh-dependent fracture criterion. Table 1 Values of s and q for different values of e = Epunch / Ehalf-space E Tension Shear S q 0.2 0.3 0.31 1 0.4 0.41 10 0.44 0.47 4. DISCUSSION The obtained element size dependence originates from the non-uniform stress distribution in the contact area with high values of stresses near boundaries. The peaks of stresses take place for all components of stress tensor and for equivalent (von Mises) stress which determines crack formation due to the fact that the von Mises fracture criterion was used in the performed calculations. At that, the smaller an element size the higher the peaks and the earlier fracture occurrence (see Fig. 4). Based on the hypothesis that the power-law dependence of contact strength on element size has nearly the same value of the exponent as the contact stress distribution, we have carried out the following test. In the paper [6] an approximation for normal stress distribution has been proposed: 2 2 ( ) ( , ) /( )p x FM a a x     , (9) where F – total force, applied on punch, M(λ, a) – a non-dimensional weight function, x – spatial coordinate along the contact patch measured from the punch center, λ – the exponent (see also, [23]). There is an analytic solution for λ, obtained by Rao [3] for a punch with an arbitrary angle at the corner θ:   2tan(1 ) (1 ) sin 2 sin 2(1 ) 1 cos 2(1 ) (1 ) (1 cos 2 )e                 . (10) For a cylindrical punch which is rectangular in 2D cross-section θ = π / 2 and Eq. (10) simplifies to 2 tan(1 ) sin(1 ) 1 cos(1 ) 2(1 ) 0e            . (11) It is interesting that Eqs. (10) and (11) do not contain Poisson’s ratios of a punch and a half-plane that means, in particular, that a value of λ is the same for plane stress and plane strain conditions [3]. 48 A. DIMAKI, E. SHILKO, S. PSAKHIE, V. POPOV Fig. 4 Spatial distributions of the von Mises stress on the surface of the half-space for different values of element size Having applied Eq. (9) for approximation of the von Mises stress distribution in the contact area, we obtained the estimation λ ≈ 0.39 (see Fig. 5). This value agrees well with the values of exponents s and q, which enter Eqs. (6) and (7) correspondingly. The given fact demonstrates that the character of the dependence of contact strength on element size is determined by stress distribution in a contact area. The analytic estimate of parameter λ, calculated by means of solution of Eq. (11), is λ ≈ 0.226 for e = 1. The discrepancy between the given analytic estimate and the value of λ, obtained on the basis of approximation of numerical simulation results (see. Fig. 2), has the following reason: we applied Eq. (9) for approximation of the von Mises stress distribution in the contact area although the given equation was initially obtained for description of normal stress distribution, taking no account of squeezing and bending of a material in the contact area and surroundings. Fig. 5 Normalized dependence of the von Mises stress under the punch and its analytic approximation Epunch = Ehalf-space. Simulation of Fracture Using a Mesh-Dependent Fracture Criterion in a Discrete Element Method 49 Based on the results described above, we suggest an equation for the value of ultimate stress in the von Mises failure criterion: ( ) 0 0 ( ) e Mises d Y d Y d         , (12) where λ(e) is a function of the relation e = Epunch / Ehalf-space, the value of Y0 depends on material strength, ratio e, contact size, etc. Obtaining closed-form equations for λ(e) and Y0 requires additional research. Application of Eq. (12) in fracture criterion (5) allows obtaining almost the same values of tensile strength of samples with different element size (see Figs. 6a and 6b). One can see that the loading diagrams of samples under tension are almost identical (Fig. 6b). A distinction between them is conditioned by a small difference of stiffness of samples with different element size. The deviation of tensile stress from its sample mean value does not exceed two percent (Fig. 6b). It is a well-known fact that macroscopic plasticity of brittle solids is often a consequence of microscopic failures nucleation, motion and healing [22]. This is one of the facts underlining similarity and interconnectedness of plasticity and fracture phenomena. From this point of view it is obvious to make an assumption about an existence of mesh-dependent criterion of plasticity. The latter is the topic for a future work. Fig. 6 a) The diagrams of tensile loading for samples with different element size; b) The deviation of tensile strength from its mean value vs. element size Epunch = Ehalf-space 5. CONCLUSIONS We have carried out the numerical simulation of a contact between an elastic punch and a half-space within the discrete element method. We have shown that the use of local detachment criterion without accounting for a discrete element size leads to a power-law dependence of contact strength on element size. The value of exponent in this dependence is approximately equal to 0.4±0.01 both for tensile loading and for shearing. Based on the obtained results we have proposed a mesh-dependent fracture criterion which explicitly includes the discrete element size. The suggested fracture criterion provides almost size-independent values of tensile and shear strength for a punch detachment from a half-space. Although the obtained fracture criterion is not universal, it demonstrates a 50 A. DIMAKI, E. SHILKO, S. PSAKHIE, V. POPOV perspective of development of discrete-element models for brittle and elastic-plastic materials with stress concentrators. Acknowledgements: This work was supported in parts by the Fundamental Research Program of the State Academies of Science for 2013-2020 (Russia) and by the Deutscher Akademischer Austauschdienst (DAAD). REFERENCES 1. Hertz, H., 1881, Über die Berührung fester elastischer Körper, Journal fuer die reine und angewandte Mathematik, 92, pp. 156-171. 2. Johnson, K.L., 1987, Contact Mechanics, Cambridge University Press, 452 p. 3. Shtaerman, I.Ya., The Contact Problem of the Theory of Elasticity, Moscow-Leningrad, 271 p. (in Russian). 4. Rao, A.K., 1971, Stress concentrations and singularities at interface corners, ZAMM, 51, pp. 395-406. 5. Okubo, H., 1951, On the two-dimensional problem of a semi-infinite elastic body compressed by an elastic plane, The Quarterly Journal of Mechanics and Applied Mathematics, 4(3), pp. 260-270. 6. Jordan, E.H., Urban, M.R., 1999, An approximate analytical expression for elastic stresses in flat punch problems, Wear, 236(1-2), pp. 134-143. 7. Bazant, Z., 2005, Scaling of Structural Strength: 2nd Ed., Butterworth-Heinemann, 336 p. 8. He, Z., Kotousov, A., Berto, F., Branco, R., 2016, A Brief Review of Recent Three-Dimensional Studies of Brittle Fracture, Physical Mesomechanics, 19(1), pp. 6–20. 9. Novozhilov, V.V., 1969, On a necessary and sufficient criterion for brittle strength, Journal of Applied Mathematics and Mechanics, 33(2), pp. 212-222. 10. Whitney, J.M., Nuismer, R.J., 1974, Stress fracture criteria for laminated composites containing stress concentrators, Journal of Composite Materials, 8(3), pp. 253-265. 11. Seweryn, A., 1994, Brittle fracture criterion for structures with sharp notches, Engineering Fracture Mechanics, 47(5), pp. 673–681. 12. Suknev, S.V., 2008, Formation of tensile fractures in the stress concentration zone in gypsum, Journal of Mining Science, 44(1), pp. 43-50. 13. Bobinski, J., Tejchman, J., 2005, Modelling of Concrete Behaviour with a Non-Local Continuum Damage Approach, Archives of Hydro-Engineering and Environmental Mechanics, 52(3), pp. 243-263. 14. Tovo, R., Livieri, P., Benvenuti, E., 2006, An implicit gradient type of static failure criterion for mixed mode loading, International Journal of Fracture, 141(3-4), pp. 497-511. 15. Pohrt, R., Popov, V.L., 2015, Adhesive contact simulation of elastic solids using local mesh-dependent detachment criterion in boundary elements method, Facta Universitatis-Series Mechanical Engineering, 13(1), pp. 3-10. 16. Li, Q., Popov, V.L., 2017, Boundary element method for normal non-adhesive and adhesive contacts of power-law graded elastic materials, Computational Mechanics, doi: 10.1007/s00466-017-1461-9. 17. Popov, V.L., Pohrt, R., Li, Q., 2017, Strength of adhesive contacts: Influence of contact geometry and material gradients, Friction, 5(3), pp. 308–325. 18. Tavarez, F.A., M.E. Plesha, M.E., 2007, Discrete element method for modelling solid and particulate materials, International Journal for Numerical Methods in Engineering, 70(4), pp. 379-404. 19. Psakhie, S.G., Shilko, E.V., Grigoriev, A.S., Astafurov, S.V., Dimaki, A.V., Smolin, A.Y., 2014, A mathematical model of particle–particle interaction for discrete element based modeling of deformation and fracture of heterogeneous elastic–plastic materials, Engineering Fracture Mechanics, 130, pp. 96-115. 20. Psakhie, S.G., Dimaki, A.V., Shilko, E.V., Astafurov, S.V., 2016, A coupled discrete element-finite difference approach for modeling mechanical response of fluid-saturated porous materials, International Journal for Numerical Methods in Engineering, 106(8), pp. 623-643. 21. Shilko, E.V., Psakhie, S.G., Schmauder, S., Popov, V.L., Astafurov, S.V., Smolin, \A.Yu., 2015, Overcoming the limitations of distinct element method for multiscale modelling of materials with multimodal internal structure, Computational Materials Science, 102, pp. 267-285. 22. Paterson, M.S., Wong, T.-F., 2005, Experimental rock deformation – the brittle field, Springer-Verlag, Berlin Heidelberg, Germany, 348 p. 23. Popov, V.L., Heß, M., Willert, E., 2018, Handbuch der Kontaktmechanik: Exakte Lösungen axialsymme- trischer Kontaktprobleme, Springer Vieweg, 339 p.