Applications of lattice method in the simulation of crack path in heterogeneous materials

The simulation of critical and subcritical crack propagation in heterogeneous materials is not a simple problem in computational mechanics. These topics can be studied with different theoretical tools. In the crack propagation problem it is necessary to lead on the interface between the continuum and the discontinuity, and this region has different characteristics when we change the scale level point of view. In this context, this work applies a version of the lattice discrete element method (LDEM) in the study of such matters. This approach lets us to discretize the continuum with a regular tridimensional truss where the elements have an equivalent stiffness consistent with the material one wishes to model. The masses are lumped in the nodes and an uni-axial bilinear relation, inspired in the Hilleborg constitutive law, is assumed for the elements. The random characteristics of the material are introduced in the model considering the material toughness as a random field with defined statistical properties. It is important to highlight that the energy balance consistence is maintained during all the process. The spatial discretization lets us arrive to a motion equation that can be solved using an explicit scheme of integration on time. Two examples are shown in the present paper; one of them illustrates the possibilities of this method in simulating critical crack propagation in a solid mechanics problem: a simple geometry of grade material. In the second example, a simulation of subcritical crack growth is presented, when a pre-fissured quasi-brittle body is submitted to cyclic loading. In this second example, a strategy to measure crack advance in the model is proposed. Finally, obtained results and the performance of the model are discussed.


INTRODUCTION
he crack propagation in heterogeneous materials has some particular issues, as the localization phenomenon and the influence of clusters of microfissures which are some of the effects that characterize this problem.The study of this problem from a continuum mechanics point of view let advance, but it is necessary to overpass the homogenization of anisotropic damage with a non-homogeneous distribution.An alternative to the continuum approach is to simulate T damage development through a discrete element model.In this context, this paper presents a method of discrete elements formed by bars.This method consists essentially in representing the continuum through a set of uniaxial elements whose stiffness is equivalent to the solid that must be modeled, model masses are lumped in the nodes connecting the elements.The non-linear behavior is captured using a very simple bilinear constitutive law applied in each element.The area below this constitutive law is coherent with the energy balance of the system.The governing equations obtained as result of the spatial discretization are integrated in the time domain using an explicit scheme of integration (Central Difference Method).In the present paper two applications are presented.In the first one, a plate of Functionally Grade Material is studied when the variation of the critical crack path is simulated for different levels of heterogeneity in the mechanical properties.The Functionally Graded Materials (FGM) are a new generation of engineering composites characterized for having a smooth variation of mechanical/thermal, or electromagnetic properties.They are new advanced multifunctional materials, which are tailored to take advantage of their constituents, for example, in a ceramic/metal FGM, heat and corrosion resistance of ceramics works together with mechanical strength and toughness of metals.To carry out an application of FGMs, scientific knowledge of fracture and damage tolerance are important for improving their structural integrity.For this kind of problems, the Finite Element Method with cohesive interfaces is usually employed (Xu and Needleman [1]; Paulino and Zanhg [2]).In the second application, the lattice discrete element method here presented is applied in the simulation of subcritical crack growth in a heterogeneous material.In this example, the proposed strategy to measure crack propagation velocity is commented.Other approaches in the simulation of subcritical crack propagation due to fatigue can be quoted, among them the works of Xu and Yuan [3] and Unger et al [4], using the extended Finite Elements Method, and the works of Yang & Cheng [5], Yamaguchi et al [6] and Siegmund [7], that combine the Finite Elements Method with cohesive interface techniques.

DISCRETE ELEMENT MODEL DESCRIPTION
he Discrete Element Method employed in this paper was proposed by Riera [8] and is based on the representation of a solid by means of a cubic arrangement of elements able to carry only axial loads.The discrete elements representation of the orthotropic continuum was adopted to solve structural dynamic problems by means of explicit numerical integration of the equations of motion, assuming the mass lumped at the nodes.Each node has three degrees of freedom, corresponding to the nodal displacements in the three orthogonal coordinate directions.Fig. 1a,b illustrates the basic bar arrangement used in this approach.The equations that relate the properties of the elements to the elastic constants of an isotropic medium are: 18 24 3 in which E and  denote the Young's modulus and the Poisson's ratio, respectively, while l A and d A represent the areas of longitudinal and diagonal elements.The resulting equations of motion may be written in the well-known form: in which x  represents the vector of generalized nodal displacements, M is the diagonal mass matrix, C is the damping matrix, also assumed diagonal,   r F t  is the vector of internal forces acting on the nodal masses and   P t  is the vector of external forces.Obviously, if M and C are diagonal, Eq. 2 is not coupled.Then the explicit central finite differences scheme may be used to integrate Eq. 2 in the time domain.Since the nodal coordinates are updated at every time step, large displacements can be accounted for in a natural and efficient manner.

Non-linear constitutive model for material damage
The softening law for quasi-brittle materials, proposed by Hillerborg [9], was adopted to handle the behavior of the materials simulated here by means of the triangular constitutive relationship for the LDEM bars, presented in Fig. 1c, which allows accounting for the irreversible effects of crack nucleation and propagation.The area under the force vs. strain curve (the area of the triangle OAB) is related to the energy density necessary to fracture the element area of influence.Thus, for a given point P on the force vs. strain curve, the area of the triangle OAP quantifies the energy density dissipated by damage.

T
Once the damage energy density equals the fracture energy, the element fails and loses its load carrying capacity.On the other hand, under compression the material is assumed linearly elastic.Thus, failure in compression is induced by indirect tension.Constitutive parameters and symbols shown in Fig. 1c are defined below: the element axial force F depends on the axial strain ε.The area associated to each element was given in Eq. 1 for longitudinal and diagonal elements.An equivalent fracture area, f i A , of each element is defined in order to satisfy the condition that the energy dissipated by fracture of the continuum and by its discrete representation are equivalent.With this purpose, fracture of a cubic sample of dimensions LLL is considered.The energy dissipated by fracture of a continuum cube due to a crack parallel to one of its faces is presented in Eq. 3: in Eq. 3, Λ represents the actual fractured area, i.e., L 2 .On the other hand, the energy dissipated when a LDEM module of dimensions LLL fractures in two parts consists of the contributions of five longitudinal elements (four coincident with the module edges and one internal element) and four diagonal elements.Then, the energy dissipated by a LDEM module can be written as shown in Eq. 4 (Kosteski et al [10]).The first term between brackets, in Eq. 4, accounts for the four edge elements, the second term for the internal longitudinal element, while the third term represents the contribution of the four diagonal elements.The coefficient c A is a scaling parameter used to establish the equivalence between Γ and ΓDEM.Thus: From Eq. 5 it follows that 3 / 22 A c  .Finally, the equivalent fracture areas of the longitudinal and the diagonal elements are presented in Eq. 6.The critical failure strain ( p  ) is defined as the largest strain attained by the element before damage initiation (point A in Fig. 1c).The relationship between p  and the specific fracture energy, f G , is given in terms of Linear Elastic Fracture Mechanics as: in which f R is the so-called failure factor, which may account for the presence of an intrinsic defect of critical size a .f R may be expressed in terms of critical fracture dimension of a , as presented in Eq. 8, where Y represents the dimensionless parameter that depends on both the specimen and the crack geometry.Notice that, if the characteristic dimension of the simulated body, L, is smaller than the intrinsic critical crack size, a , the collapse will be stable.On the other hand, if L > a , an instable global collapse is expected.Also, one could write the failure factor, f R , in terms of the fragility number, s , proposed by Carpinteri [11], this dimensionless parameter measures the fragility quality of a determined structure.In the expressions presented in Eq. 9, one can see the fragility number, s , in terms of the critical stress intensity factor, IC K , and the critical stress, p  ; or in terms of the critical specific fracture energy, f G , and the critical strain, p  ; or, using Eq.7 and Eq. 8, in terms of the LDEM parameters, as shown in Eq. 10.
If two models are built with different materials and different dimensions, one can expect similar global mechanical behavior if both models have the same s value.; The element loses its load carrying capacity when the limit strain, r  , is reached (point B in Fig. 1c).This value must satisfy the condition that, upon failure of the element, the dissipated energy density equals the product of the element fracture area, f A , and the specific fracture energy, , f G divided by the element length.Hence: in which the sub index i is replaced by l or d depending on whether the element under consideration is longitudinal or diagonal.The coefficient r K is a function of the material properties and the element length i L .In order to guarantee the stability of the algorithm, the condition 1 r K  must be satisfied.In this sense, it is interesting to define the critical element length cr L (see Eq. 13).
  . Thus, for practical purposes, a single value of the critical length can be used for longitudinal and diagonal elements.Therefore, the stability condition is expressed by Eq. 16 and the limit strain by Eq. 17.
The random distribution of material parameters in the LDEM environment is modelled, defining the toughness, f G , as a random field with a Type III (Weibull) extreme value distribution.Details of the implementation can be found in Kosteski et al [12].

FIRST APPLICATION: SIMULATION OF CRITICAL CRACK PROPAGATION IN FUNCTIONAL GRADED MATERIALS
he present example shows a rectangular plate of Polimetilmetacrilate (PMMA) with a central fissure.The plate dimensions are show in Fig. 2.a.This plate is submitted to constant strain velocity of 5m/s in its borders.Paulino and Zhang [2] studied the influence that the Functionally Graded Material (FGM) has in dynamic crack propagation in this same geometry using the Finite Element Method with cohesive interfaces.The cohesive elements were implemented in the region where fissure propagation could take place.The quoted authors studied the influence of three different types of FGM in the crack propagation event.In case 1 the homogeneous property was considered.In case 2 a hypothetic FGM material was proposed, where just the cohesive properties graded linearly in y axis direction.Finally, in case 3 not only the cohesive properties, but also the finite element properties were graded.

DEM Model
Only the middle of the plate was modeled due to the geometry and load symmetries of the problem.The plate was discretized with 60 modules on each side and one module in the thickness direction.The boundary conditions applied permit to represent the left border as a symmetry axis.The plane strain condition was imposed fixing the displacements in T the thickness direction; the half of the fissure was discretized using 6 cubic modules.The material properties and the main model characteristics are presented in Fig. 2b.Fig. 3 shows the constitutive law adopted for the material, considered homogeneous in the vertical direction of the plate.The critical step time is related to the time that an elastic wave takes to pass through the normal elemental bar.In the present case, for homogeneous material, the step critical time was Δtcrit =1.82 E-8 s.The main difference between Paulino and Zhang's [2] analysis and the present study is the Poisson coefficient value, in DEM model when one works with the cubic arrangement (see Fig. 1) one is limited to work with Poisson's coefficient of 0.25, if the goal is to model an elastic isotropic and homogeneous material.If one wishes to model elastic and isotropic materials with other Poisson's coefficient, it will be necessary to employ another kind of geometric arrangement.In Fig. 2d the case characteristics studied by Paulino and Zhang [2] are presented.The uniaxial constitutive law used in this work in the case 1 is also presented in Fig 2 .c.

Obtained Results
Case 1, Homogeneous case: The homogeneous material is considered in this first case.Thus, in this case, there isn't any spatial variation in the elements constitutive law in DEM model.The Fig. 4(i) shows a fracture configuration obtained by Paulino and Zhang [2] (a), the main crack begins to branch when it reaches a length of 1.05mm, with an approximate angle of 29° in relation to the horizontal direction.In the case of the LDEM simulation (b) the branching begins when the main crack reaches a length of 1.0mm with an angle of 32° in relation to the horizontal direction.In both simulations it is possible to verify another characteristic that is seen in Fig. 4(i) (a) and (b) which is the similarity between the two final configurations, for example, it is possible to observe a secondary branching in both figures.The lack of symmetry in the final configuration of LDEM simulation characterizes an unstable propagation process.This method has been verified sensitive to disturbances that occur during the process.These disturbances could be produced by rounding errors in the calculus and due to little changes in the step integration time.Case 2, Graded Gc : In the present case only the toughness parameter c G is graded in vertical direction.The cohesive interface of the Finite Element Model used by Paulino and Zhang [2] is characterized by three parameters that are: the normal maximum tensile in the interface max n T , the critical opening displacement n and the area closed by the curve that is proportional to the toughness c G .The grade in the interface properties that was implemented to modify the curve is illustrated in Fig. 3(a).The graded material proposed in case 2, for LDEM simulation, was carried out modifying the uniaxial curve as indicated in  Case 3, Graded E and Gc : In case 3, the Elasticity Modulus E and the toughness parameter Gc are changed simultaneously in vertical direction.In Paulino and Zhang's [2] work, a graduation was proposed to the interface properties in the vertical direction, in the same manner of case 2, and the Elasticity Modulus E in the elements was modified.
In DEM approach for case 3, the constitutive elemental law was modified as shown in Fig. 3(c).In Fig. 4(iii) the comparison between Paulino and Zhang's [2] implementation (a) and LDEM implementation (b), in terms of final configurations, is shown.In the present case, the similarity between both configurations is not too clear, but both show the same general tendency.
(i) (ii) (iii) Comparison in terms of Energy Balance: Here the results in terms of Energy Balance during the whole fracture process are discussed.In Fig. 5(a) the complete balance of energy during the whole process for case 1 is shown.In Fig. 5 (b), (c), and (d), the elastic, kinetic and dissipated energy through all the damage process are shown for the three cases.
It can be verified in Fig. 5(d) that case 1 shows high level of dissipated energy due to the greater fracture area generated during the simulation.In Fig. 5(c), the abrupt change in kinetic energy in cases 2 and 3 shows a sensitive change in the crack propagation velocity during the whole process.This is in accordance with the final configuration obtained, where branching and secondary cracks appear more than in case 1.
Case 1 Case 2 Case 3 Kinetic Energy (J) t (s) 0.0 In the present case, plane strain condition was adopted; this condition is implemented in LDEM using only one cubic module of thickness and restraining of displacements in out of plane direction (dir.y).The plate was also fixed in the right boundary as indicated in Fig. 6 (left).A pre-crack was introduced in the model with the goal of favoring subcritical crack growth through the middle of the plate.It is important to highlight that the goal here is not to characterize a specific material, but to explore the capacity of LDEM in simulating the subcritical propagation of a preexistent fissure.With the considered parameters, the ratio between rupture and critical strain will be r/p = 60.
For this analysis, material toughness G f was considered as a random field characterized by means of a variation coefficient.It will also be considered that the correlation length of the random field is of the same order of the discretization level.Implementations regarding the unlinking of discretization level and correlation length of random field are presented in Puglia et al [13].The load was applied in the central nodes of the cubic cells that constitute the top and bottom surfaces of the plate, as indicated in Fig. 6 (left).In Fig. 6 (right) the graphic shows the tensile force transmitted to a normal bar element placed in the middle of the plate's high, h, and close to the pre-crack tip versus percentage of simulation time.In this model, a cubic cell length of L = 0.0075 m was adopted, the model has 100 modules in the x direction, 1 module in y direction and 41 modules in z direction.The material properties are characterized for Poisson coefficient of  = 0.25, the Young Modulus is E = 35 GPa, the density value is  = 2400 kg/m 3 and the critical strain  p = 2.18 10 -4 .The toughness random field is characterized by a mean value of(Gf) = 1155 N/m and a variation coefficient of CV(Gf)=5%.In Fig. 6 (right), it can be verified that, in a first moment, load rises according to an exponential function, until a stable magnitude is achieved, so that inertial effects are none.After that, oscillation begins and amplitude grows, also according to an exponential function, until stable loading cycle amplitude is reached.This oscillation does not induce considerable dynamic effects in the model, and after 880 cycles loading amplitude is assumed constant.The applied loads produce remote stress acting in the plate of 2.93MPa before oscillation starts, ranging later between a top value of 3.34MPa and lower value of 2.52MPa.

Obtained Results
In Fig. 7 (left), six crack configurations obtained during propagation process are ilustrated.The elements in darker shade of grey, siding the main crack, are parilly damaged.Broken elementes constitute main crack and are not shown.The applied version of LDEM constitutive law does not take residual strain into account.Constitutive laws that account for such effects are presented in Kosteski [10].The cluster of microfissures that precedes the main fissure in fatigue crack growth (characteristic in quasi-brittle materials) can be clearly seen in the partial configurations illustrated in Fig. 7 (left).It is also important to observe that crack propagation does not influence in a significant way the discrete element mesh.

W
In Fig. 7 (right) the graphic illustrates fluctuations of kinetic, elastic and dissipated energy during all the damaged process, where abrupt changes in such values of energy are associated to local instabilities in main crack grow, denoting the spasmodic character of simulated subcritical growth.Aiming to improve visualization of the energetic fluctuations during simulated crack growth, the dissipated energy curve (red) is multiplied by a 0.02 scaling factor and elastic energy (green) by a 0.04 scaling factor.Kinetic energy curve (blue) is shown in original scale.Observing Fig. 7 it is possible to conclude that the macro crack propagates in subcritical way from (t/tmax=0,70) untill (t/tmax=0.96)when the process becomes unstable.Within subcritical growth, in the next subsection, fatigue life is characterized according to the classic Paris law.Fatigue Life Characterization: The Paris Law, based in a potencial law, presents subcritical crack growth velocity (da/dN) as a function of the delta of stress intensity factor (K) acting at the main crack tip.For this work, was computed using the expresion presented in Gdoutos [14] which is adequate to the model geometry, boundary conditions and pre-crack configuration here adopted.The value of Kchanges (grows) during the simulation process, asK also depends on the crack length, a. Having obtained crack propagation velocity (da/dN) and stress intensity factor (K) by mean of the LDEM tool, the values can be subtituted in Paris equation [(da/dN) = CK m ] to model fatigue life during subcritical crack growth.Applying log operation in both sides of Paris expression it is possible to obtain the constants and that characterize fatigue behavior for the simulated material.The velocity (da/dN) is obtained in the LDEM model using the x coordinate of every element and the instant in time when it reaches the critical strain value  p , becoming damaged.This strategy implies in assuming that the microfissure cluster which precedes the main crack advances through the plate in the same velocity as the crack itself.Thus, (da/dN) is measured computing the time in which each bar reaches the critical strain  p , versus the x coordinate of the element central point.Also, the same procedure is realized when rupture strain r is achieved.The resulting graphic is seen in Fig. 8 (left).In Fig. 8 (left) the curve denoted by the green dots corresponds to the first elements to become damaged in a given x coordinate.Therefore, this set of points forms a curve linked with the damage propagation velocity at the crack tip and the propagation velocity of the main crack itself, according to the previous assumption made.Thus, 10 representative points were selected to carry out the curve a versus N. The cycle count begins in the normalized time t*=70, when the loading reaches permanent regime (see Fig. 6 right).Knowing the crack propagation velocity (da/dN) during the simulation process, it is possible to plot Fig. 8 (right), where the log(da/dN) versus log(K) is presented.The tendence line of this curve models the region in crack propagation known as region of stable growth or region II (see Gdoutos [14]).The applied tendence line gives the straight line equation:

CONCLUSIONS
n this work, two applications of the Lattice Discrete Element Method for studying damage development in quasi-brittle materials were presented.In the first example, unstable crack propagation in functional graded material was analyzed.In the second example, subcritical crack propagation due to fatigue damage was studied.In both examples it was possible to verify the great potential of the applied discontinuous mechanics tool in the study of damage development.

Figure 1 :
Figure 1: LDEM discretization strategy: (a) basic cubic module, (b) generation of a prismatic body, (c) Bilinear constitutive law adopted for LDEM uniaxial elements.
15)In the special case of an isotropic continuum with 0.25   , the value of the coefficients shown above are 1

Figure 2 :
Figure 2: (a) Layout and boundary condition of the rectangular plate with central fissure and imposed constant strain velocity; (b) The mechanical properties adopted for the LDEM model; (c) The constitutive uniaxial law adopted for the bars; (d) The characteristic of the analyzed cases.

Fig. 3b .Figure 3 :
Figure 3: a) The cohesive interface laws used by Paulino and Zhang [2]; The elemental constitutive law used in LDEM implementation: (a) for the graded material in case 2; (b) for the graded material in case 3.

Figure 4 :Figure 5 :
Figure 4: (i) Final configuration for the case 1, homogeneous material.Final configuration for the graded materials (ii) for the case 2, (iii) for the case 3, (in all cases (i-iii) (a) the Paulino and Zhang's [2] config.,(b) the LDEM config.
of studying the subcritical crack propagation in a rectangular plate built with quasi-brittle material, a model was developed using the Lattice Discrete Element Method (LDEM) presented before.The rectangular plate model has a pre-crack that emanates from a lateral border.The loads and boundary conditions used are illustrated inFig 6(left).

Figure 6 :
Figure 6: (Left) LDEM model used with details of how the loads and boundary conditions are applied, where b=0.75m (100 modules) and h=0.3075m (41 modules); (Right) Force transmitted to a normal bar element close to pre-crack tip versus percentage of simulation time.

Figure 7 :
Figure 7: (Left) Main crack configuration during six moments of simulation, where

Figure 8 :
Figure 8: (Left) Normalized time versus the x coordinate of the central points of elements that reached (red dots) and (blue dots) during the simulation, the curve equivalent to is denoted by the green dots; (Right) Plot versus where the straight line indicates subcritical crack growth.The velocity is expressed in µm per cycle.
log(da/dN)  0.22 log(K) -4.37, this means m = 0.22 and C = 4.26  10 -5 .Therefore, it is possible to characterize the modeled material according to Paris law as: da/dN = 4.26  10 -5 K 0.22 .Pointing out that calculation procedure to obtain the constants C and m was carried out with stress values in Pascal [Pa] and length in meters [m].As previously said, the present work focus on exploring the possibilities of the LDEM tool in simulating and characterizing fatigue damage development.Matters regarding real material modeling are topics for future works.