Acta Polytechnica CTU Proceedings doi:10.14311/APP.2020.26.0039 Acta Polytechnica CTU Proceedings 26:1–6, 2020 © Czech Technical University in Prague, 2020 available online at https://ojs.cvut.cz/ojs/index.php/app MODELLING OF CRACK PROPAGATION: COMPARISON OF DISCRETE LATTICE SYSTEM AND COHESIVE ZONE MODEL Karel Mikeša, ∗, Franz Bormannb, Ondřej Rokoša, b, Ron H.J. Peerlingsb a Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic b Eindhoven University of Technology, Department of Mechanical Engineering, P.O. Box 513, 5600 MB, Eindhoven, The Netherlands ∗ corresponding author: karel.mikes@cvut.cz Abstract. Lattice models are often used to analyze materials with discrete micro-structures mainly due to their ability to accurately reflect behaviour of individual fibres or struts and capture macroscopic phenomena such as crack initiation, propagation, or branching. Due to the excessive number of discrete interactions, however, such models are often computationally expensive or even intractable for realistic problem dimensions. Simplifications therefore need to be adopted, which allow for efficient yet accurate modelling of engineering applications. For crack propagation modelling, the underlying discrete microstructure is typically replaced with an effective continuum, whereas the crack is inserted as an infinitely thin cohesive zone with a specific traction-separation law. In this work, the accuracy and efficiency of such an effective cohesive zone model is evaluated against the full lattice representation for an example of crack propagation in a three-point bending test. The variational formulation of both models is provided, and obtained results are compared for brittle and ductile behaviour of the underlying lattice in terms of force-displacement curves, crack opening diagrams, and crack length evolutions. The influence of the thickness of the process zone, which is present in the full lattice model but neglected in the effective cohesive zone model, is studied in detail. Keywords: Lattice model, damage, finite element method, cohesive zone model, three-point-bending test, crack propagation. 1. Introduction Lattice models with dissipative interactions are often used for modelling of materials with discrete micro- structures such as textile, paper, 3D printed materials, or concrete; see e.g. [1–3]. The advantage of lattice models lies in their ability to capture microstructural interactions individually, thus representing the geom- etry of the underlying micro-structure as well as the associated constitutive laws accurately and naturally. Because non-local behaviour is important for mecha- nisms occurring in the process zone at the crack tip, discrete models are suitable for problems such as dam- age localization, crack propagation, and branching. Lattice models suffer, on the other hand, from high computational and memory requirements, necessitat- ing simplifications especially in cases in which multiple evaluations are required, such as optimization, mate- rial parameter identification, or stochastic simulations. Although multiple options for building effective mod- els exist, in this contribution we focus on a continuum approximation of the underlying microstructure with an embedded cohesive zone model. Crack propagation and its cohesive properties are described in such an ap- proach by a suitable traction-separation law, whereas the bulk behavior of the continuum may often be elastic. The aim of this paper is detail comparison of crack propagation for discrete lattice model and continuous cohesive zone approach. The accuracy and perfor- mance of cohesive zone model predictions are assessed on the example of crack propagation in a three-point- bending test, as shown in Figure 1, by comparison with the exact underlying full lattice model. 2. Full lattice model An admissible configuration of a rate-independent sys- tem of interest is specified by a state variable q = (r,z), where r stores the kinematic variables (displace- ment of nodes) and z the internal variables (damage in bonds). The energetic solution q is determined using the following incremental minimization problem at time tk q(tk) ∈ arg min q̂∈Q Πk(q̂; q(tk−1)), k = 1, . . . ,nT , (1) with an initial condition q(0) = q0, where Q denotes the space of all admissible configurations, and Πk is an incremental energy defined as Πk(q; q(tk−1)) = E(tk,q) + D(z,z(tk−1)). (2) The incremental energy consists of the potential (Gibbs type) energy E and the dissipation distance D. 39 https://doi.org/10.14311/APP.2020.26.0039 https://ojs.cvut.cz/ojs/index.php/app K. Mikeš, F. Bormann, O. Rokoš, R.H.J. Peerlings Acta Polytechnica CTU Proceedings L H F ~ex ~ey ΩA Figure 1. Geometry and boundary conditions of the considered three-point-bending test with a close-up of the full lattice representation. The potential energy reads E(tk,q) = V(q) −fText(tk)r, (3) where V is the internal free energy and fext is a vector of external loads. The dissipation distance D(z2,z1) measures the minimum dissipation by a continuous transition be- tween two consecutive states z1 and z2; for further details see [4] (Section 3.2). It can be written as D(z2,z1) = 1 2 ∑ α,β∈Bα Dαβ(z2,z1), (4) where Bα denotes the set of nodes connected with node α. Dαβ is the dissipation distance for a single bond αβ, defined as follows Dαβ(z2,z1) = { Dαβ(ω2) −Dαβ(ω1) if ω2 ≥ ω1 +∞ otherwise, (5) where Dαβ(ω) reflects the amount of energy dissipated in that bond during a unidirectional damage process up to a given damage level ω. In each time step, the minimization problem of Eq. (1) is solved to obtain a local minimum that satisfies the energy balance V(q(tk)) −V(q(0)) + VarD(q, tk) = Wext(q, tk), (6) which equates the internally stored energy V(q) − V(q0) plus the dissipated energy VarD(q, tk) = k∑ i=1 D(z(ti),z(ti−1)), (7) with the work performed by the external forces Wext(q, tk) = k∑ i=1 1 2 [f(ti) + f(ti−1)]T[r(ti) −r(ti−1)]. (8) Unless stated otherwise, the same material model is considered for each lattice bond. In undamaged state, it is defined by pair potential φαβ(lαβ) = 1 2 E0A0l αβ 0 ( ε(lαβ) )2 , (9) where E0 is the Young’s modulus and A0 the cross- sectional area of the bond. l0 denotes the initial length parameter ε0 εf E0 A0 brittle 0.005 10ε0 1/ε0 1 ductile 0.005 100ε0 1/ε0 1 Table 1. Parameters of the lattice model. of the bond and ε(lαβ) = (lαβ − l0)/l0 is the bond strain. The internal energy is defined as V(r,z) = 1 2 ∑ α,β∈Bα [ (1 −ωαβ)φαβ(lαβ) ] . (10) The elastic energy of a bond φαβ is reduced by the damage variable ωαβ(ε), defined as ωαβ(ε) = { 1 −ε0 exp ( −ε−ε0 εf ) if ε ≥ ε0 0 if ε < ε0, (11) where ε0 is the limit elastic strain, and the parame- ter εf characterizes the softening branch of the asso- ciated stress-strain diagram. The particular choice of material and geometric parameters corresponding to the lattice model are specified in Table 1. For large lattice structures, the simulation must be controlled by the increment of dissipation, cf. e.g. [5], because both the (typically used) crack mouth opening displacement control (CMOD) as well as crack tip displacement opening control (CTOD) turned out to be insufficient to keep the energy balance of Eq. (6) satisfied. 3. Cohesive zone model The lattice problem as specified in Figure 1 is ap- proximated with an effective continuous finite element cohesive zone model (FE-CZ). For that purpose, the following assumptions are adopted: (1.) Only one crack is initiated; (2.) The crack trajectory is a straight line and the crack propagates vertically in the middle of the beam; (3.) All damage processes take place in the predefined crack path, whereas the rest of the domain remains elastic. 40 vol. 26/2020 Comparison of Discrete Lattice System and Cohesive Zone Model ΩAR (elastic region) L H F ΓCZ ΩAL (elastic region) Cohesive zone with traction separation law Figure 2. Geometry and boundary conditions of the finite element cohesive zone model (FE-CZ) subjected to a three-point bending test: elastic domain (grey) and cohesive zone (dashed line). 3.1. Model Description The crack is represented as one cohesive zone, Γcz, normal to the ~ex direction, which splits the considered domain Ω into two symmetric parts ΩL and ΩR: Γcz = { ~x ∈ R2 : x = xcz = 0, |y| ≤ H/2 } , ΩL = { ~x ∈ R2 : −L/2 ≤ x < xcz, |y| ≤ H/2 } , ΩR = { ~x ∈ R2 : xcz < x ≤ L/2, |y| ≤ H/2 } , Ω = ΩL ∪ ΩR, (12) where xcz = 0 denotes the horizontal position of the vertical cohesive zone, cf. Figure 2. The total free energy Ψ of the FE-CZ model consists of two contributions, i.e. Ψ = ∫ Ω\Γcz ψe dΩ + ∫ Γcz ψcz dΓcz, (13) where ψe is the elastic strain energy density ψe = 1 2 ε : C : ε, (14) corresponding to a linear elastic isotropic material under plane strain assumptions, described through a fourth-order constitutive tensor C, even though the full lattice is not isotropic. In Eq. (14), the elastic strain tensor is expressed as the symmetric part of the displacement gradient ~∇~u, i.e. ε = 1 2 ( ~∇~u + ( ~∇~u )T) . (15) The cohesive zone potential ψcz is a function of the displacement jump ~∆ = J~uK = ~uR −~uL (16) between the two sub-domains ΩR and ΩL. Due to the symmetry of the problem, the tangential component of the jump ∆t vanishes along the cohesive interface, i.e. ~uLy = ~u R y , on Γcz, (17) meaning that the cohesive zone potential ψcz becomes a function of the normal crack opening ∆n only. The constitutive parameters of both models, i.e. C and ψcz(∆n), are fitted in Section 3.2. For the solu- tion, in each loading step the resulting total potential energy of Eq. (13) is minimized with respect to the degrees of freedom associated with a discretized dis- placement field of the continuum model. In order to handle brittle behaviour and snap-back, CMOD indirect control is used. 3.2. Material Parameter Calibration Both effective models, i.e. the cohesive zone poten- tial ψcz and elasticity stiffness tensor C, are identified directly from the properties of the underlying lattice by homogenization. For that purpose, a single periodic unit cell is consid- ered, as shown in Fig. 3. Note that the cross-sectional areas of the horizontal and vertical bounding links are reduced by a factor of 0.5 to account for periodicity. l0 l0 ∆n ∆n Figure 3. The geometry of unit lattice used for calibration of elastic properties and traction-separation law. By homogenization, the effective stiffness tensor can be assembled and compared with the plane strain for- mula for Hooke’s law to obtain the elastic coefficients. In particular, the effective Young’s modulus reads Eiso = 5 6 E0A0 l0 (1 + √ 2/2), (18) whereas Poisson’s ratio ν is fixed at ν = 1/4 for the given lattice geometry; for more details see e.g. [6, Appendix A]. 41 K. Mikeš, F. Bormann, O. Rokoš, R.H.J. Peerlings Acta Polytechnica CTU Proceedings Figure 4. Normalized traction as a function of normalized opening for the two parameter sets considered, cf. Tab. 1. Figure 5. The damage (red) in bonds computed with full lattice model: evolution of the process zone one step before crack initiation (top) and final process zone and fully developed crack path in the last step of simulation (bottom). Instead of the cohesive zone potential ψcz, the cor- responding traction, T , defined as T(∆n) = dψcz(∆n) d∆n , (19) is evaluated numerically from the response of the periodic unit cell subjected to a prescribed displace- ment ∆n with the boundary conditions according to Figure 3. The computed values of tractions corre- sponding to the two different parameter sets of Tab. 1 are plotted in Figure 4. The cohesive zone poten- tial ψcz can be obtained from these tractions by nu- merical integration. But effectively, for the practical implementation, the values of ψcz are not needed be- cause the problem is solved for the force balance by the standard Newton algorithm. 4. Results In this section, the full lattice model defined in Sec- tion 2 (referred to as Full lattice), and finite element cohesive zone model defined in Section 3 (referred to as FE-CZ ), are compared on a three-point-bending test of a rectangular specimen of dimensions 2048l0×256l0, for the brittle and ductile lattice material settings (re- call Table 1 and Figure 4). For the Full lattice model, individual bonds start to damage in a relatively ductile manner, resulting in a considerable process zone, featuring multiple possible crack paths, cf. Figure 5 (top). After localization, a full (central) crack develops across the entire height of the specimen, cf. Figure 5 (bottom). The force-displacement responses are compared in Figure 6. The FE-CZ is able to exactly capture the initial elastic stiffness. However, the peak force is significantly underestimated by FE-CZ , for both ma- terials. The relative peak force error (of FE-CZ in comparison with Full lattice) is 23.3 % for the brittle and even 33.4 % for the ductile material. This is caused by the fact that for the Full lattice a large amount of inelastic strain along the bottom edge is ac- cumulated in the process zone. This complex behavior cannot be captured with one cohesive zone and oth- erwise elastic behavior naturally results in significant deviations in the FE-CZ model. The inelastic process zone in the Full lattice unloads the crack mouth (com- pared to FE-CZ ) and, therefore, higher loading force is needed to initiate the crack. To verify this, we limit the process zone in the lattice model in a further simulation (referred to as Elastic lattice). In this model, the damage can evolve only in a narrow band of width 2l0 positioned along the vertical axis of the symmetry, whereas the remaining lattice is considered as purely elastic. To this end, the critical strain parameter is modified as εi0 = { ε0 if |xi −xcz| ≤ l0 ∞ othervise (20) 42 vol. 26/2020 Comparison of Discrete Lattice System and Cohesive Zone Model Figure 6. Normalized force–displacement diagrams for a brittle (left) and a ductile (right) lattice. The results of Elastic lattice correspond to the lattice where the damage can evolve only in a narrow band of width 2l0, whereas the remaining lattice is elastic. Figure 7. Normalized opening profiles for fixed pre-peak loading force F/(E0A0ε0) = 56 (left) and for post-peak loading force F/(E0A0ε0) = 14 (right) computed for the ductile material setting. Figure 8. Normalized crack length as a function of crack mouth opening displacement (CMOD) computed for the ductile material setting. where xi −xcz is the horizontal distance between the center of the i-th bond and the symmetry axis. For all considered aspects, the FE-CZ model pro- vides very similar results to the Elastic lattice, where the influence of the process zone is removed. It reveals that the cohesive zone in FE-CZ is able to exactly cap- ture the damage properties of the underlying lattice and the error in the peak force is produced mainly by the missing process zone. This can also be observed in the difference of pre- and post-peak opening pro- files, corresponding to a fixed loading force, as shown in Figure 7 for two magnitudes of the applied force. Here we clearly see that the FE-CZ and Elastic lat- tice (both with a narrow process zone) tend to show significantly larger crack mouth openings before crack localization. Once the crack localizes and opens, how- ever, all models provide very accurate results. Also the kinematics of crack propagation, depicted in Figure 8, is accurately captured by both the Elastic lattice and FE-CZ . The FE-CZ (compared to the Full lattice model in computing time) provides, however, a substantial speed-up of the order of 50. 5. conclusion In this contribution, a cohesive-zone finite element approach to modelling of crack propagation in materi- als with discrete microstructures has been compared against the full lattice model that features phenomena observed in real materials, such as non-local behaviour, large deformations and rotations of individual struts, and distributed cracking. It has been shown that the cohesive-zone model achieves a substantial speed up, while being suffi- ciently accurate in terms of kinematic quantities such as overall displacement, crack length, or crack opening profile. On the other hand, the conjugate quantities, 43 K. Mikeš, F. Bormann, O. Rokoš, R.H.J. Peerlings Acta Polytechnica CTU Proceedings such as the peak force, are significantly underesti- mated due to the missing process zone which results in a significantly more ductile response of the con- sidered lattice system. With increasing ductility of the lattice, the accuracy further decreases. The pre- sented results have also shown that in order to correct for this deficiency, more ductile behaviour should be accounted for, e.g., by adding multiple (parallel) cohe- sive zones that cover the size of the main part of the process zone. However, the spatial positioning and orientation of such secondary cohesive zones is not known and needs to be based on prior knowledge. In the most general case, every inter-element interface can be considered as a cohesive zone, which might sig- nificantly complicate the entire procedure, introduce a dependence on the discretization (finite element size and orientation), and cause significant computational overhead. Let us note that, apart from the homogenization pro- cedures used in this paper, the effective cohesive-zone constitutive parameters can be obtained by fitting directly the overall macroscopic force–displacement curves. Although this might improve the macro- scopic accuracy (in terms of force–displacement curve) achieved by the cohesive-zone model, such a phe- nomenological approach would not directly represent the underlying microstructure. This would in turn mean that different geometries would require indepen- dent identifications and that the accuracy in kinematic quantities (such as crack opening profile) might be decreased. The cohesive zone model provides an efficient and effective tool for the description of crack propagation in discrete systems, although care must be taken to position the cohesive zone interface properly to cover not only the crack trajectories, but also phenomena important for given problem. List of symbols q State variable r Kinematic variable z Internal variable t Parametrization time nT Number of time steps Q Space of all admissible configurations Πk Incremental energy E Potential energy D Dissipation distance V Internal free energy VarD Dissipated energy fext Vector of external loads α,β Node indices Bα Nearest-neighbour nodes to α ωαβ Damage level in a bond αβ Dαβ Energy dissipated in bond αβ Wext Work performed by the external forces φ Pair potential E0 Young’s modulus of individual bonds A0 Cross-sectional area l Length of deformed bond l0 Initial lattice spacing ε,ε0 Bond strain, critical bond strain εf Damage softening parameter Ω Elastic domain Γcz Cohesive zone Ψ Total free energy ψcz Cohesive zone potential ψe Elastic strain energy density ε Elastic strain tensor σ Elastic stress tensor C Elastic stiffness tensor ~u Displacement field ∆n, ∆t Normal/tangential jump in displacements Eiso Effective Young’s modulus ν Effective Poison’s ratio F Magnitude of the point loading force u Displacement under the point loading force x,y Cartesian coordinates Acknowledgements Financial support of this work from the Grant Agency of the Czech Technical University in Prague (SGS) under project No. 19/032/OHK1/1T/11 and from the Czech Science Foundation (GAČR) under project No. 17-04150J is gratefully acknowledged. References [1] L. Beex, C. Verberne, R. Peerlings. Experimental identification of a lattice model for woven fabrics: application to electronic textile. Composites Part A: Applied Science and Manufacturing 48:82–92, 2013. doi:10.1016/j.compositesa.2012.12.014. [2] E. Bosco, R. H. Peerlings, M. G. Geers. Predicting hygro-elastic properties of paper sheets based on an idealized model of the underlying fibrous network. International Journal of Solids and Structures 56:43–52, 2015. doi:10.1016/j.ijsolstr.2014.12.006. [3] G. Cusatis, Z. P. Bažant, L. Cedolin. Confinement-shear lattice model for concrete damage in tension and compression: I. theory. Journal of Engineering Mechanics 129(12):1439–1448, 2003. doi:10.1061/(ASCE)0733-9399(2003)129:12(1439). [4] A. Mielke, T. Roubíček. Rate-Independent Systems: Theory and Application. Springer-Verlag New York, 1st edn., 2015. doi:10.1007/978-1-4939-2706-7. [5] C. V. Verhoosel, J. J. Remmers, M. A. Gutiérrez. A dissipation-based arc-length method for robust simulation of brittle and ductile failure. International Journal for Numerical Methods in Engineering 77(9):1290–1321, 2009. doi:10.1002/nme.2447. [6] K. Mikeš, M. Jirásek. Quasicontinuum method extended to irregular lattices. Computers & Structures 192:50–70, 2017. doi:10.1016/j.compstruc.2017.07.002. 44 https://doi.org/10.1016/j.compositesa.2012.12.014 https://doi.org/10.1016/j.ijsolstr.2014.12.006 https://doi.org/10.1061/(ASCE)0733-9399(2003)129:12(1439) https://doi.org/10.1007/978-1-4939-2706-7 https://doi.org/10.1002/nme.2447 https://doi.org/10.1016/j.compstruc.2017.07.002 Acta Polytechnica CTU Proceedings 26:1–6, 2020 1 Introduction 2 Full lattice model 3 Cohesive zone model 3.1 Model Description 3.2 Material Parameter Calibration 4 Results 5 conclusion List of symbols Acknowledgements References