Continuous anisotropic damage as a twin modelling of discrete bi-dimensional fracture
CContinuous anisotropic damage as a twin modellingof discrete bi-dimensional fracture
C. Oliver-Leblond, R. Desmorat, B. Kolev
Universite´ Paris-Saclay, ENS Paris-Saclay, CNRS, LMT - Laboratoire de Me´canique et Technologie, 91190, Gif-sur-Yvette, France
Abstract
In this contribution, the use of discrete simulations to formulate an anisotropic damage model is investigated. It isproposed to use a beam-particle model to perform numerical characterization tests. Indeed, this discrete model explicitlydescribes cracking by allowing displacement discontinuities and thus capture crack induced anisotropy.Through 2D discrete simulations, the evolution of the effective elasticity tensor for various loading tests, up to failure,is obtained. The analysis of these tensors through bi-dimensional harmonic decomposition is then performed to estimatethe tensorial damage evolution. As a by-product of present work we obtain an upper bound of the distance to theorthotropic symmetry class of bi-dimensional elasticity.
Keywords:
Anisotropic Damage, Crack Density, Harmonic Decomposition, DEM, Beam-Particle, Lattice, Discretesimulation
Introduction
For most quasi-brittle materials, the initial mechanicalbehaviour of the uncracked material can be considered asisotropic. On the other hand, degradation induced bymechanical loading generally leads to anisotropy. Indeed,cracks are naturally oriented and therefore their appear-ance will not affect the material properties in an isotropicmanner. We note that the orientation of these cracks for aninitially isotropic material can be determined by the direc-tion and sign of the loading (Mazars et al., 1990; Ramtaniet al., 1992).To obtain an explicit representation of cracking – andthus of its impact on material properties – it is possibleto use discrete models of the lattice or particle type. Thefirst lattice models were historically introduced by Poisson(1828) and more recently by Hrennikoff (1941) to solveclassical elasticity problems. The elastic material is dis-cretized using 1D elements – springs or beams – whichallow for the transfer of forces between the nodes of thelattice. The development of numerical simulations has al-lowed for its extension to the study of fracture behaviourby considering a brittle behaviour for the elements form-ing this lattice (Herrmann et al., 1989). The approach wasthen applied to the quasi-brittle fracture of concrete sub-jected to tension (Schlangen and Van Mier, 1992). How-ever, these models do not allow for the representation ofcompression cracking nor of cyclic loading. Particulatemodels were proposed in 1979 to study the behaviour of
Email addresses: [email protected] (C.Oliver-Leblond), [email protected] (R.Desmorat), [email protected] (B. Kolev) granular (Cundall and Strack, 1979) joints. For these ap-plications, contact forces alone were sufficient to correctlydescribe the behaviour. A cohesive version has been pro-posed latter (Meguro and Hakuno, 1989) but does not offerthe simplicity and rapidity of the lattice models. In thepresent work, we propose to use a beam-particle modelcombining the lattice approach and the particulate ap-proach (D’Addetta et al., 2002; Delaplace, 2008; Vassauxet al., 2016). Concrete is represented via an assembly ofpolygonal particles bounded together by brittle beams.After the beams break, frictional contact forces are in-troduced between the particles. This model allows for afine and explicit description of the cracking and the asso-ciated mechanisms (initiation, propagation, closure withstiffness recovery, friction). Its application to structuralcalculations is not yet common but can be envisaged withthe implementation of high-performance calculation tech-niques. These discrete methods can be used as numericalexperimentation tools at the scale of the RepresentativeElemental Volume or at the scale of a laboratory speci-men to establish and identify the constituent equations ofa continuous model (Vassaux et al., 2015a) or simply apart of these equations (Delaplace and Desmorat, 2007).In order to obtain a representative, robust and efficientnumerical model, it is common to use macroscopic lawsthat account for the formation, propagation and coales-cence of micro-cracks by introducing an internal damagevariable. One of the first damage models for quasi-brittlematerials (Mazars, 1984) made it possible to reproducethe degradation of the concrete material with loading, viaa scalar damage variable, by considering concrete as a ho-mogeneous material at the scale of the volume element
Preprint submitted to Elsevier November 20, 2020 a r X i v : . [ phy s i c s . c l a ss - ph ] N ov f continuum mechanics. This isotropic damage model iscommonly used to study the behaviour of concrete struc-tures but its isotropic nature does not allow for it to takeinto account complex multi-axial loading (as observed byMazars and his collaborators, Ramtani et al. (1992), seealso the work of Halm and Dragon (1998)). To tacklethis issue, damage models with induced anisotropy – suchas the one proposed by Desmorat et al. (2007) and mademore robust later (Desmorat, 2016) – can be used. In thismodel, developed within the framework of the thermody-namics of irreversible processes, the damage is representedby a tensor variable of order 2. It should be noted thatthese macroscopic models are globally phenomenologicalbecause their formulation and identification are based onexperimental observations of the behaviour of quasi-brittlematerials. It is the same for their laws of evolution (Ram-tani et al., 1992; Lemaitre and Desmorat, 2005) or theirnon-local character. Their representativeness is thus lim-ited by the capacity to carry out experiments on the de-graded material.Discrete cracking models and continuous damage mod-els both have advantages and limitations, and it makessense to combine these advantages by combining the twomodelling methods. In this paper, we investigate the rel-evance to use the beam-particle model to perform virtualtests in order to formulate an anisotropic damage model.Let us recall that in (Delaplace and Desmorat, 2007) themethodology was limited to the identification of a singleparameter of the continuous model (certainly delicate tomeasure on real tests because it was the one governing theso-called ”shear-bulk” coupling of a more important effectof the damage on compressibility than on shear). More re-cent studies have proposed the formulation of continuousdamage models from a discrete analysis: by introducingscalar damage variables calculated from the macroscopicloss of stiffness in a lattice simulation (Rinaldi and Lai,2007; Rinaldi, 2013; Jivkov, 2014), or by continualizationof discrete equations (Challamel et al., 2015).In this paper, the goal is to derive a tensorial damagevariable from discrete simulations. The tools for the in-trinsic analysis of tensors, introduced by Backus (1970) inelasticity and by Leckie and Onat (1981) in damage me-chanics, mixing harmonic analysis and the notion of covari-ants (generalizing that of invariants, Olive et al. (2018b)),are used here in order to analyze the effective elasticitytensors obtained through discrete simulations, without ref-erence to a particular basis, and to achieve a general ten-sorial representation of damage. Discrete simulations willbe used to obtain the evolution of the effective elasticitytensor during the rupture of a numerical specimen underdifferent mechanical loads. The analysis of the tensorsvia the harmonic decomposition and their covariant re-construction (Olive et al., 2018a) will then allow us toestimate the evolution of the tensor damage.In section 1, definitions relating to the properties of fullysymmetrical tensors are recalled. The harmonic decom-position of a bi-dimensional elasticity tensor is then per- formed in section 2 and it is reminded that the fourth orderharmonic part of this elasticity tensor is always an har-monic square . In section 3, a covariant reconstruction of2D orthotropic elasticity tensors is used to derive an upperbound of the distance to the orthotropic symmetry class inbi-dimensional elasticity. The beam-particle model, usedin this contribution as a numerical testing tool, is pre-sented in section 4 as well as the methodology to extractthe effective elasticity tensor from the discrete simulations.A first analysis is realised in section 5 to check whether theinitial elasticity tensor, of the uncracked medium, can beconsidered as isotropic. In section 6, the analysis of theeffective elasticity tensors – and more precisely of theirharmonic part – of cracked media for various loading testup to failure is achieved. Finally, the analysis of thosetensors through harmonic decomposition/covariant recon-struction is then performed in section 7 to estimate thetensorial damage evolution.
1. Definitions
We next make use of the Euclidean structure of ℝ 𝑑 , 𝑑 = 2 ,
3, and do not make difference between covariant,contravariant or mixed tensors.
We denote by T s the totally symmetric part of a possiblynon symmetric tensor T . More precisely, if T ∈ ⊗ 𝑛 ℝ 𝑑 isof order 𝑛 , T s ( 𝑥𝑥𝑥 , . . . , 𝑥𝑥𝑥 𝑛 ) := 1 𝑛 ! ∑︁ 𝜎 ∈ S 𝑛 T ( 𝑥𝑥𝑥 𝜎 (1) , . . . , 𝑥𝑥𝑥 𝜎 ( 𝑛 ) ) , where S 𝑛 is the permutation group of 𝑛 elements.The symmetric tensor product of two tensors T and T , of respective orders 𝑛 and 𝑛 , is the symmetrizationof T ⊗ T , defining a totally symmetric tensor of order 𝑛 = 𝑛 + 𝑛 : T ⊙ T := ( T ⊗ T ) s . Contracting two subscripts 𝑖, 𝑗 of a tensor T of order 𝑛 defines a new tensor of order 𝑛 − 𝑖𝑗 T .For a totally symmetric tensor T , this operation does notdepend on a particular choice of the pair 𝑖, 𝑗 . Thus, wecan refer to this contraction just as the trace of T and wewill denote it as tr T . It is a totally symmetric tensor oforder 𝑛 −
2. Iterating the process, we definetr 𝑘 T = tr(tr( · · · (tr T ))) , which is a totally symmetric tensor of order 𝑛 − 𝑘 .Harmonic tensors H are by definition totally symmet-ric traceless tensors, i.e. such as H = H 𝑠 , tr 𝑘 H = 0.Their vector spaces are denoted ℍ 𝑛 ( ℝ 𝑑 ), with 𝑛 the orderof considered tensor, 𝑑 the dimension. In 2D ( 𝑑 = 2),dim ℍ 𝑛 ( ℝ ) = 2, while in 3D ( 𝑑 = 3), dim ℍ 𝑛 ( ℝ ) =2 𝑛 + 1. The harmonic tensors of order two ( i.e. the devi-atoric tensors) will be denoted by lowercase letter h .2 .3. Harmonic product Let H and H be two harmonic tensors of orders 𝑛 and 𝑛 respectively. The harmonic product H * H , definingan harmonic tensor of order 𝑛 = 𝑛 + 𝑛 , has been intro-duced in (Olive et al., 2018a) as the leading harmonic partof the symmetric tensor product H ⊙ H , H * H := ( H ⊙ H ) ′ ∈ ℍ 𝑛 + 𝑛 ( ℝ 𝑑 ) , 𝑑 = 2 , . In previous works, the leading harmonic part was some-times denoted ( H ⊙ H ) . As it is a generalization ofthe deviatoric part, the notation ( · ) ′ is here preferred to( · ) . It is computed at any order thanks to the harmonicdecomposition (see section 2). The harmonic product is associative and commutative , H * ( H * H ) = ( H * H ) * H , H * H = H * H . Let us particularize the harmonic product for specificcases of bi-dimensional (2D) vectors and harmonic second-order tensors (belonging thus to ℍ 𝑛 ( ℝ ) for 𝑛 = 1 , Example . For two vectors 𝑤𝑤𝑤 , 𝑤𝑤𝑤 ∈ ℍ ( ℝ ), we have 𝑤𝑤𝑤 * 𝑤𝑤𝑤 =( 𝑤𝑤𝑤 ⊙ 𝑤𝑤𝑤 ) ′ = 12 ( 𝑤𝑤𝑤 ⊗ 𝑤𝑤𝑤 + 𝑤𝑤𝑤 ⊗ 𝑤𝑤𝑤 ) −
12 ( 𝑤𝑤𝑤 · 𝑤𝑤𝑤 ) , where 𝑤𝑤𝑤 · 𝑤𝑤𝑤 = 𝑤𝑤𝑤 𝑇 𝑤𝑤𝑤 is the scalar product. Example . For two second-order harmonic (deviatoric)tensors h , h ∈ ℍ ( ℝ ), we have (Olive et al., 2018a) h * h =( h ⊙ h ) ′ = h ⊙ h −
14 tr( h h ) ⊙ . Example . The harmonic square of a 2D second orderharmonic (deviatoric) tensor h ∈ ℍ ( ℝ ) writes (Oliveet al., 2018a): h * h =( h ⊙ h ) ′ = h ⊙ h −
14 ( h : h ) ⊙ = h ⊗ h −
12 ( h : h ) J , (1.1)where J = I − ⊗ . (1.2)with 𝐼 𝑖𝑗𝑘𝑙 = ( 𝛿 𝑖𝑘 𝛿 𝑗𝑙 + 𝛿 𝑖𝑙 𝛿 𝑗𝑘 ). T (for the rotation group) The action 𝑔 ⋆ T of a rotation 𝑔 ∈ SO( 𝑑 ) on a tensor T of order 2 or 4 is( 𝑔 ⋆ T ) 𝑖𝑗 = 𝑔 𝑖𝑘 𝑔 𝑗𝑙 𝑇 𝑘𝑙 , ( 𝑔 ⋆ T ) 𝑖𝑗𝑘𝑙 = 𝑔 𝑖𝑝 𝑔 𝑗𝑞 𝑔 𝑘𝑟 𝑔 𝑙𝑠 𝑇 𝑝𝑞𝑟𝑠 . A tensor A ( T ) is said to be a covariant of tensor T forSO( 𝑑 ) if 𝑔 ⋆ A ( T ) = A ( 𝑔 ⋆ T ) , ∀ 𝑔 ∈ SO( 𝑑 ) , and dimension 𝑑 is taken next as 𝑑 = 2.The algebra of polynomial covariants of the elastic-ity tensor has been defined (and studied) in (Oliveet al., 2018b) for 𝑑 = 3 (three-dimensional case) andin (Desmorat et al., 2020), for 𝑑 = 2 (bi-dimensional case).
2. Harmonic fourth order part H of the 2D elas-ticity tensor as an harmonic square
The harmonic decomposition of tensors is a powerfulmathematical tool (Schouten, 1989; Spencer, 1970), thathas first been applied to three-dimensional elasticity ten-sors C ∈ 𝔼 la( ℝ ) by Backus (1970). Formally in 2D, it isthe equivariant decomposition C = ( 𝜇, 𝜅, d ′ , H ) ∈ 𝔼 la( ℝ ) , into two scalars (invariants) 𝜇, 𝜅 ∈ ℍ ( ℝ ) ≃ ℝ , 𝜇 be-ing the shear modulus, 𝜅 the bi-dimensional bulk mod-ulus, one harmonic (deviatoric) second order covariant d ′ = d ′ ( C ) ∈ ℍ ( ℝ ) and one harmonic fourth order co-variant H = H ( C ) ∈ ℍ ( ℝ ), such as 𝑔 ⋆ C = ( 𝜇, 𝜅, 𝑔 ⋆ d ′ , 𝑔 ⋆ H ) ∀ 𝑔 ∈ SO(2) . An explicit harmonic decomposition of C ∈ 𝔼 la( ℝ ) is: C = 2 𝜇 J + 𝜅 ⊗ + 12 ( ⊗ d ′ + d ′ ⊗ ) + H , (2.1)with J defined by (1.2). The closed form expressions ofharmonic components 𝜇 , 𝜅 , d ′ and H are gained thanks tothe definition of dilatation and Voigt second order tensors d = tr C and of v = tr C , 𝜇 = 18 (2 tr v − tr d ) ,𝜅 = 14 tr d , d ′ = d −
12 (tr d ) , H = C − 𝜇 J − 𝜅 ⊗ −
12 ( ⊗ d ′ + d ′ ⊗ ) . (2.2)In the isotropic (initial) case, one has H = 0 and d = 2 𝜅 , v = (2 𝜇 + 𝜅 ) , the Young’s modulus and Poisson’s ratio being 𝐸 = 4 𝜅𝜇𝜅 + 𝜇 , 𝜈 = 𝜅 − 𝜇𝜅 + 𝜇 . emark . In 2D the deviatoric parts of dilatation andof Voigt second order tensors are equal, v ′ = d ′ , and theterm ⊗ d ′ + d ′ ⊗ is also equal to ⊗ d ′ + d ′ ⊗ , where( a ⊗ b ) 𝑖𝑗 = ( 𝑎 𝑖𝑘 𝑏 𝑗𝑙 + 𝑎 𝑖𝑙 𝑏 𝑗𝑘 ).One can also perform the harmonic decomposition ofthe compliance tensor S = C − , and gets then d ′ ( S ) = v ′ ( S ) = (tr S ) ′ , where S = 12 𝜇 J + 14 𝜅 ⊗ + 12 ( ⊗ d ′ ( S ) + d ′ ( S ) ⊗ ) + H ( S ) , and H ( S ) is the harmonic (totally symmetric and trace-less) fourth order part of S .Recall finally that a minimal generating set (a mini-mal integrity basis) of the invariant algebra of the elas-ticity tensor in 2D, under the action of the orthogonalgroup O(2), consists in the following 5 invariants (Vianello,1997), 𝜇 = 18 (2 tr v − tr d ) ,𝜅 = 14 tr d ,𝐼 = ‖ d ′ ‖ = d ′ : d ′ ,𝐽 = ‖ H ‖ = H :: H ,𝐾 = d : H : d . (2.3) H ∈ ℍ ( ℝ ) as an harmonicsquare In 2D, it has been shown that any fourth order harmonictensor H ̸ = 0 is of the form (Desmorat and Desmorat,2015) H = 2Λ e * e , tr e = 0 , ‖ e ‖ = 1 , (2.4)where e is a unit deviatoric (second order) eigentensor as-sociated with a non zero eigenvalue Λ of the Kelvin repre-sentation of H ( i.e. such that H : e = Λ e in an orthonor-mal basis). Remark . Since H ̸ = 0 has two opposite eigenvaluesΛ + > − = − Λ + <
0, this implies that any bi-dimensional fourth order harmonic tensor is always an har-monic square (Desmorat and Desmorat, 2015, 2016). Thecorresponding explicit formula is as follows, H = h * h , h = ± √︀ + e + , (2.5)where the unit deviatoric second order tensor e = e + isthe eigentensor of H associated with the strictly positiveeigenvalue Λ + . Observe that there are two opposite har-monic square roots h of H , which are opposite to eachother. A new proof of this result, more conceptual, is pro-vided in Appendix A.
3. Covariant reconstruction of 2D orthotropic elas-ticity tensors
Bi-dimensional elasticity tensors C ∈ 𝔼 la( ℝ ) haveKelvin matrix representation[ C ] = ⎡⎣ 𝐶 𝐶 √ 𝐶 𝐶 𝐶 √ 𝐶 √ 𝐶 √ 𝐶 𝐶 ⎤⎦ . Bi-dimensional harmonic fourth order tensors H ∈ ℍ ( ℝ )have Kelvin representation[ H ] = ⎡⎣ 𝐻 − 𝐻 √ 𝐻 − 𝐻 𝐻 −√ 𝐻 √ 𝐻 −√ 𝐻 − 𝐻 ⎤⎦ . (3.1)The normal form of an orthotropic elasticity tensor cor-responds to 𝐶 = 𝐶 = 0. If C = H is moreoverharmonic, we get 𝐻 = 0.A bi-dimensional harmonic fourth order tensor H can-not be strictly orthotropic ( i.e. with exact symme-try group 𝔻 , the symmetry group of a rectangle). Anon-vanishing harmonic fourth-order tensor H has al-ways the symmetry group 𝔻 (the symmetry group of asquare) (Vannucci and Verchery, 2001; Vannucci, 2005).The covariants of H inherit its symmetry (Olive et al.,2018a,b): this implies that all the second order covariantsof H have at least the square symmetry, they are there-fore isotropic (and all the deviatoric second order covari-ants of H vanish). Combined with the fact that the har-monic product * is itself covariant, this geometrical prop-erty implies that the harmonic square root h ( i.e. such as H = h * h ) is not a covariant of H . If an elasticity tensor itself has the (exact) square sym-metry , then d ′ = 0 and H ̸ = 0 (Verchery, 1982; Vianello,1997; Vannucci, 2005), and C has no better reconstruc-tion formula by means of its covariants than its harmonicdecomposition, C = 2 𝜇 J + 𝜅 ⊗ + H . When an elasticity tensor C is orthotropic, its secondorder covariant d ′ is an eigentensor of its harmonic part H = H ( C ) (given by (2.1)), i.e. H : d ′ = Λ d ′ , (3.2)but where Λ = ± Λ + is a non-vanishing eigenvalue. Tocheck this, just consider its Kelvin normal form[ C ] = ⎡⎣ 𝐶 𝐶 𝐶 𝐶
00 0 2 𝐶 ⎤⎦ , (3.3)4hich leads to d ′ = (tr C ) ′ = 𝐶 − 𝐶 (︂ − )︂ ̸ = 0 , and, using (2.1), we get[ H ] = 𝐻 ⎡⎣ − − − ⎤⎦ . Contracting (3.2) with d ′ , altogether with the relations H : = : H = 0 , d ′ : H : d ′ = 𝐾 , we obtain finally Λ = 𝐾 𝐼 , where the invariants 𝐼 and 𝐾 are defined by (2.3). Thismeans by (2.4) that H = 2Λ ‖ d ′ ‖ d ′ * d ′ . Hence, we have the following result.
Theorem 3.1.
Any bi-dimensional orthotropic elasticitytensor C = ( 𝜇, 𝜅, d ′ , H ) ∈ 𝔼 la( ℝ ) can be reconstructed bymeans of its 4 invariants 𝜇 , 𝜅 , 𝐼 , 𝐾 , and of its (devia-toric) second order covariant d ′ , as C =2 𝜇 J + 𝜅 ⊗ + 12 ( ⊗ d ′ + d ′ ⊗ ) + H , H = 2 𝐾 𝐼 d ′ * d ′ , (3.4) where d ′ * d ′ = d ′ ⊗ d ′ −
12 ( d ′ : d ′ ) J . From theorem 3.1, one can derive an upper bound Δ forthe distance to orthotropy of a bi-dimensional elasticitytensor as defined below.
Corollary 3.2.
Let C = ( 𝜇, 𝜅, d ′ , H ) ∈ 𝔼 la( ℝ ) be a bi-dimensional elasticity tensor with no material symmetry,let 𝐼 = ‖ d ′ ‖ , 𝐽 = ‖ H ‖ and 𝐾 = d : H : d . Then, thepositive invariant Δ = ⃦⃦⃦ 𝐂 − 𝜇 𝐉 − 𝜅 ⊗ −
12 ( ⊗ 𝐝 ′ + 𝐝 ′ ⊗ ) − 𝐾 𝐼 𝐝 ′ * 𝐝 ′ ⃦⃦⃦ = √︀ 𝐽 𝐼 − 𝐾 𝐼 . is an upper bound of the distance of C to the orthotropicsymmetry class. In corollary 3.2, the fact that C has no material sym-metry implies that ‖ d ′ ‖ = 𝐼 ̸ = 0. The first equalitycorresponds to the norm of the difference between the elas-ticity tensor and its orthotropic reconstruction formula (in which necessarily d ′ ̸ = 0). The Harmonic decompositionformula (2.1) gives indeedΔ = ⃦⃦⃦⃦ H − 𝐾 𝐼 d ′ * d ′ ⃦⃦⃦⃦ . The second equality in corollary 3.2 is obtained by ex-panding the square norm, using (1.1) and the facts that J :: J = 2 and H :: J = 0, see (Desmorat and Desmorat,2015).
4. Discrete elements representative volumes
The discrete method used here to perform the virtualtesting of a quasi-brittle heterogeneous material is a beam-particle approach detailed in (Vassaux et al., 2016) andsummarized in figure 1.The representative volume is divided into an assemblyof rigid particles. The particle mesh is generated from aVoronoi tesselation of a set of randomly generated pointswithin a grid. This operation generates polygonal parti-cles.Dual of the Voronoi tesselation, the Delaunay triangu-lation associates a segment with each pair of neighbouringparticles. This segment is used as a geometric support foran elastic Euler-Bernoulli beam modelling the cohesion ofthe material. Each beam 𝑏 is parameterized by its length 𝑙 𝑏 , its section 𝐴 𝑏 , the Young modulus 𝐸 and the coefficientof inertia 𝛼 = 64 𝐼 𝑏 𝜋/𝐴 𝑏 . The first two coefficients 𝑙 𝑏 and 𝐴 𝑏 are different for each beam and imposed by the geom-etry of the mesh. The next two, 𝐸 and 𝛼 , are identicalfor all the beams and identified in order to reproduce themacroscopic elastic behaviour.In order to reproduce the fracture behaviour, a failurecriterion 𝑃 𝑝𝑞 is associated with each beam connecting twoparticles 𝑝 and 𝑞 : 𝑃 𝑝𝑞 = 𝜀 𝑝𝑞 𝜀 𝑐𝑟𝑝𝑞 + | 𝜃 𝑝 − 𝜃 𝑞 | 𝜃 𝑐𝑟𝑝𝑞 > 𝜀 𝑐𝑟𝑝𝑞 and thebreaking threshold in rotation 𝜃 𝑐𝑟𝑝𝑞 are generated for eachbeam according to a Weibull distribution. The Weibullprobability density function adopted in this study is : 𝑓 ( 𝑥 ) = 𝑘𝜆 (︁ 𝑥𝜆 )︁ 𝑘 − 𝑒 − ( 𝑥/𝜆 ) 𝑘 (4.2)with the scale factor 𝜆 and the shape factor 𝑘 . In fact,the spatial variability for the two breaking thresholds 𝜀 𝑐𝑟𝑝𝑞 , 𝜃 𝑐𝑟𝑝𝑞 , is supposed to be identical. Therefore, three parame-ters control the fracture behaviour: a shape factor 𝑘 (com-mon to both distributions), a scale factor in extension 𝜆 𝜖 𝑐𝑟 and a scale factor in rotation 𝜆 𝜃𝑐𝑟 . The combination ofthis failure criterion and of a random generation of fail-ure thresholds makes it possible to model a quasi-brittle5ehaviour in tension and compression (as proposed andshown in Vassaux et al. (2016)). These three failure pa-rameters are identified in such a way as to reproduce thenon-linear behaviour of the material.Finally, in order to be able to capture the mechanismsof crack closure and of crack sliding, frictional contact isintroduced between the particles if they overlap while theyare not connected by a beam. The detection of the softcontact between two polygonal particles, the calculationof the normal contact force as well as that of tangentialcontact via Coulomb’s law of friction follows the proposalof Tillemans and Herrmann (1995). A rewriting is howeverproposed in order to use relative displacements rather thanrelative velocities (Vassaux et al., 2015b) (figure 1e), theformer being more suitable in the quasi-static frameworkin which we place ourselves. To parametrize the friction,an angle of friction 𝜑 is introduced resulting in a coefficientof friction tan 𝜑 for Coulomb’s law. θ p θ q u q u p p q q n c t c S r u q u p p q n pq t pq a)b) c)d)e) p Figure 1: Description of the main ingredients in the beam-particlemodel (from Oliver-Leblond (2019)).
The definition of homogenized quantities on the dis-crete volume, such as stress and strain, is necessary tolink the discrete description, providing detailed informa-tion on particle movements and interaction forces, to thecontinuous description.The average Cauchy stress tensor is computed from thesymmetrization of the definition proposed by Bagi (1996): 𝜎𝜎𝜎 = 1 𝑆 𝑁 ∑︁ 𝑝 =1 f ( 𝑝 ) ⊙ x ( 𝑝 ) (4.3)where 𝑆 is the area of the discrete representative volume(the so-called Representative Volume Element or RVE in3D, the Representative Area Element or RAE in 2D), f ( 𝑝 ) is the resulting force on the particle 𝑝 and x ( 𝑝 ) is the vectorposition of the center of the particle 𝑝 . The summation is made on the 𝑁 particles constituting the boundary of theRAE.The macroscopic strain tensor is defined as the followingmean value: 𝜖𝜖𝜖 = 1 𝑆 𝑁 ∑︁ 𝑝 =1 (︂ u ( 𝑝 ) + u ( 𝑝 +1) ⊙ 𝑛𝑛𝑛 ( 𝑝,𝑝 +1) )︂ 𝑙 ( 𝑝,𝑝 +1) (4.4)where u ( 𝑝 ) is the displacement vector of the particle 𝑝 , 𝑙 ( 𝑝,𝑝 +1) is the length of the segment linking the particles 𝑝 and 𝑝 + 1 and 𝑛𝑛𝑛 ( 𝑝,𝑝 +1) the outward pointing normal vectorof this same segment.In order to obtain the effective elasticity tensor, threemeasurement loadings are performed and the correspond-ing average stress and strain vectors are extracted. Foreach loading 𝑖 , the effective elasticity tensor links the stressto the strain in the following way:^ 𝜎 ( 𝑖 ) = [ C ] ^ 𝜖 ( 𝑖 ) (4.5)where [ C ] is Kelvin matrix representation of macroscopicelasticity tensor C and^ 𝜎 ( 𝑖 ) = ⎛⎜⎝ 𝜎 ( 𝑖 )11 𝜎 ( 𝑖 )22 √ 𝜎 ( 𝑖 )12 ⎞⎟⎠ , ^ 𝜖 ( 𝑖 ) = ⎛⎜⎝ 𝜖 ( 𝑖 )11 𝜖 ( 𝑖 )22 √ 𝜖 ( 𝑖 )12 ⎞⎟⎠ . Therefore, if the measurement loadings are chosen toensure that the strain vectors are linearly independent,their determinant is non-zero and the symmetric effectiveelasticity tensor can be computed from[ C ] = (︂(︁ ^ 𝜎 (1) ^ 𝜎 (2) ^ 𝜎 (3) )︁ (︁ ^ 𝜖 (1) ^ 𝜖 (2) ^ 𝜖 (3) )︁ − )︂ 𝑠 . (4.6) Two requirements must be met to define the measure-ment loadings: ensuring the linear independence of thestrain vectors and maintaining the cracks open duringmeasurement. Several sets of measurement loadings areused and compared: • DEF WoC : Elementary strain loading without con-tact^ 𝜖 (1) = ⎛⎝ 𝜖 ⎞⎠ , ^ 𝜖 (2) = ⎛⎝ 𝜖 ⎞⎠ , ^ 𝜖 (3) = ⎛⎝ √ 𝜖 ⎞⎠ , where 𝜖 is the applied strain that must be sufficientlysmall for the (possibly cracked) RAE to remain in theelasticity domain. The following displacement fieldsare thus applied on the boundaries of the specimen: u (1) = (︂ 𝜖𝑥 )︂ , u (2) = (︂ 𝜖𝑦 )︂ , u (3) = (︂ 𝜖𝑦𝜖𝑥 )︂ , with 𝑥 and 𝑦 the positions of the center of the particleson the boundaries.6 DEF BT C : Bi-tension strain loading with contact ^ 𝜖 (1) = ⎛⎝ 𝜖 𝑏𝑡 + 𝜖𝜖 𝑏𝑡 ⎞⎠ , ^ 𝜖 (2) = ⎛⎝ 𝜖 𝑏𝑡 𝜖 𝑏𝑡 + 𝜖 ⎞⎠ , ^ 𝜖 (3) = ⎛⎝ 𝜖 𝑏𝑡 𝜖 𝑏𝑡 √ 𝜖 ⎞⎠ , where 𝜖 𝑏𝑡 is introduced to create a bi-tension loadingkeeping the cracks open. The following displacementfields are thus applied on the boundaries of the spec-imen: u (1) = (︂ ( 𝜖 𝑏𝑡 + 𝜖 ) 𝑥𝜖 𝑏𝑡 𝑦 )︂ , u (2) = (︂ 𝜖 𝑏𝑡 𝑥 ( 𝜖 𝑏𝑡 + 𝜖 ) 𝑦 )︂ , u (3) = (︂ 𝜖 𝑏𝑡 𝑥 + 𝜖𝑦𝜖𝑥 + 𝜖 𝑏𝑡 𝑦 )︂ , • SIG F WoC : Elementary stress loading with appliedforces without contact^ 𝜎 (1) = ⎛⎝ 𝜎 ⎞⎠ , ^ 𝜎 (2) = ⎛⎝ 𝜎 ⎞⎠ , ^ 𝜎 (3) = ⎛⎝ √ 𝜎 ⎞⎠ , where 𝜎 is the applied stress that must be sufficientlysmall to remain in the elasticity domain. The trac-tion vectors t = 𝜎𝜎𝜎 n are then applied on the bound-aries ( n : outer unit normal) and the central particleis completely blocked to avoid rigid motions: t (1) 𝑥𝑚𝑖𝑛 = (︂ − 𝜎 )︂ , t (1) 𝑥𝑚𝑎𝑥 = (︂ 𝜎 )︂ , t (2) 𝑦𝑚𝑖𝑛 = (︂ − 𝜎 )︂ , t (2) 𝑦𝑚𝑎𝑥 = (︂ 𝜎 )︂ , t (3) 𝑥𝑚𝑖𝑛 = (︂ − 𝜎 )︂ , t (3) 𝑥𝑚𝑎𝑥 = (︂ 𝜎 )︂ , t (3) 𝑦𝑚𝑖𝑛 = (︂ − 𝜎 )︂ , t (3) 𝑦𝑚𝑎𝑥 = (︂ 𝜎 )︂ . • SIG D WoC : Elementary stress loading with applieddisplacements without contact.The objective is to obtain the three stress states ofthe previous measurement loading by applying onlykinematic conditions. u (1) 𝑥𝑚𝑖𝑛 · 𝑒𝑒𝑒 𝑥 = 0 , u (1) 𝑥𝑚𝑎𝑥 · 𝑒𝑒𝑒 𝑥 = 𝑢, u (1) 𝑦𝑚𝑖𝑛 · 𝑒𝑒𝑒 𝑦 = 0 , u (2) 𝑦𝑚𝑖𝑛 · 𝑒𝑒𝑒 𝑦 = 0 , u (2) 𝑦𝑚𝑎𝑥 · 𝑒𝑒𝑒 𝑦 = 𝑢, u (2) 𝑥𝑚𝑖𝑛 · 𝑒𝑒𝑒 𝑥 = 0 , u (3) 𝑥𝑚𝑖𝑛 · 𝑒𝑒𝑒 𝑦 = 0 , u (3) 𝑥𝑚𝑎𝑥 · 𝑒𝑒𝑒 𝑦 = √ 𝑢, u (3) 𝑦𝑚𝑖𝑛 · 𝑒𝑒𝑒 𝑥 = 0 , u (3) 𝑦𝑚𝑎𝑥 · 𝑒𝑒𝑒 𝑥 = √ 𝑢. where u is the applied displacement that must be suf-ficiently small for the RAE to remain in the elasticitydomain. We first compare the four measurement loadings de-scribed previously on an uncracked square specimen of100 ×
100 particles (see first line of table 1). We observethat they give very close results. The dispersion can beattributed to the measurement noise. This result was ex-pected because there are enough particles in the specimento make it a Representative Area Element. Thus, the typeof boundary conditions does not influence the result, as inan infinite medium.Measurements of the effective elasticity tensors ofcracked specimens are a more important feature of thepresent study. For this purpose, two cracking tests arecarried out: a localized cracking test and a diffuse crack-ing test. In both cases, the orientation of the loads issuch that the micro-cracks are globally perpendicular tothe 𝑦 direction. For each type of cracking, we extract thecrack patterns and the associated effective elasticity ten-sors for three levels of damage 𝐷 (computed as a relativeloss of stiffness in the 𝑦 direction, Lemaitre and Chaboche(1985)): a low level of damage ( 𝐷 ≈ . 𝐷 ≈ .
4) and a high level of damage ( 𝐷 ≈ .
7) aspresented in table 1. The goal here is then to compare thefour measurement loadings in order to select one of themfor the rest of the study.Regardless of the type or amount of cracking, the twostrain loads (DEF WoC and DEF BT C) give similar re-sults. This is expected and confirms that it is sufficient todeactivate the contact within the discrete model in orderto obtain a measurement loading equivalent to a loadingkeeping the cracks open. Since the strain loading withadditional bi-tension (DEF BT C) can cause convergenceproblems, it is important to be able to validate this equiva-lence with the elementary strain loading where the contactis deactivated (DEF WoC).In the same way, we can observe that the two elementarystress loadings (SIG F WoC and SIG D WoC) give iden-tical results, up to the measurement noise, whether forcesor displacements are applied on the boundaries. Althoughthis was not the case here, stress loading with appliedforces may prove to be less robust because the blockingof the central particle, to prevent rigid body movements,can lead to convergence problems if this particle is in anarea of significant cracking.Finally, strain and stress loadings give identical resultsin the case of diffuse cracking even for a significant levelof damage. On the other hand, the results diverge for thecase of localized cracking (see figure 2). Indeed, if we pur-sue the simulation until the complete failure of the speci-men in the localized case, we notice that the 𝐶 modulusapproaches 0 for the measurement loading in stress, whichis not the case for the measurement loading in strain. Here,we encounter a limitation in the use of damage models forquasi-brittle materials such as concrete. Indeed, the dam-age is intended to represent diffuse cracking and is thusnot appropriate for cases of localized cracking. However,7 rack patterns DEF WoC DEF BT C SIG F WoC SIG D WoC Initial stateNo damage [︂ . . . . . − . . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ [︂ . . . . . . . . . ]︂ [︂ . . . . . − . . − . . ]︂ Localized crackingLow level of damage [︂ . . . . . − . . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ [︂ . . . . . − . . − . . ]︂ [︂ . . . . . − . . − . . ]︂ Localized crackingMid level of damage [︂ . . . . . . . . . ]︂ [︂ . . . . . . . . . ]︂ [︂ . . . . . . . . . ]︂ [︂ . . . . . . . . . ]︂ Localized crackingHigh level of damage [︂ . . − . . . − . − . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ Diffused crackingLow level of damage [︂ . . . . . − . . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ [︂ . . . . . . . . . ]︂ [︂ . . . . . − . . − . . ]︂ Diffused crackingMid level of damage [︂ . . . . . − . . − . . ]︂ [︂ . . − . . . − . − . − . . ]︂ [︂ . . . . . − . . − . . ]︂ [︂ . . . . . − . . − . . ]︂ Diffused crackingHigh level of damage [︂ . . . . . . . . . ]︂ not converged [︂ . . . . . . . . . ]︂ [︂ . . . . . . . . . ]︂ Table 1: Comparison of the measurement loadings used to extract the effective elasticity tensor [ 𝐂 ] (in GPa, Kelvin representation). C C C C C Figure 2: Evolution of the components of the effective elasticitytensor for the localized cracking test according to the elementarystrain loading measurement (DEF WoC in dashed lines) and theelementary stress loading measurement with applied displacements(SIG D WoC in plain lines) damage models are commonly used up to high levels ofdegradation, where the cracks are almost necessarily lo-calized in a quasi-brittle material. In order to correctlyrepresent the loss of stiffness due to this localized crackingwhile keeping the simplicity of a kinematic loading, we de-cide to choose thereafter the elementary stress loading withapplied displacements without contact (SIG D WoC).
5. Deviation from isotropy of the initial tensor
The isotropy of discrete media is generally studied froma geometric point of view (Andre´ et al., 2012). Indeed,the geometrical isotropy of the mesh is considered as agood criterion to approach the mechanical isotropy of thematerial if the elasticity parameters are homogeneous overthe domain. This geometric isotropy can be visualizedsimply by plotting the polar histogram associated with theorientation of the beams. The figure 3 presents those polarhistograms for three meshes of different densities. Onecan see that the polar histogram tends towards a circleas the density increases, which shows geometrical isotropyconvergence.
Through the harmonic decomposition of the elasticitytensors of each material samples, it is possible to quantifytheir deviation from isotropy. The evaluation of the rela-tive distance of the elasticity tensors to the isotropic classthus provides a criterion for quantifying the validity of themechanical isotropy hypothesis.Figure 4 presents the evolution with the mesh size of therelative distance to isotropy Δ 𝑖𝑠𝑜 defined as :Δ 𝑖𝑠𝑜 = ‖ C − C 𝑖𝑠𝑜 ‖‖ C ‖ (5.1) with C 𝑖𝑠𝑜 = 2 𝜇 J + 𝜅 ⊗ the isotropic projection of C .The elasticity parameters 𝜇 and 𝜅 are determined by Eq.(2.2).For each mesh density, several simulations have beenperformed (from 200 simulations for the 10 ×
10 particlesmesh to 50 simulations for a 300 ×
300 particles mesh). Asexpected, the assumption of an initial isotropic mediumis correct if the number of particles is sufficient. In Fig-ure 3, the elasticity tensors obtained for three simulationsare given along with the associated isotropic tensors andelasticity coefficients.Given the above results, it seems safe to adopt a meshdensity of 100 ×
100 particles to ensure the mechanicalisotropy of the initial, uncracked, state.
6. Multiaxial analyses of effective elasticity tensorsup to high level of damage
In this part, a 20cm × ×
100 particles, which re-sults in an average beam size 𝑙 𝑏 of 2 mm.The parameters of the discrete model, given in the ta-ble 2, are chosen so as to reproduce macroscopically thebehaviour of a mortar of Young’s modulus 𝐸 = 36 .
35 GPaand Poisson ratio 𝜈 = 0 .
22 (corresponding to 𝜅 = 23 . 𝜇 = 14 .
92 GPa), of tensile strength 𝑓 𝑡 = 3 MPaand of compressive strength 𝑓 𝑐 = 40 Mpa. Symbol Values Unit 𝑙 𝑏 𝐸 𝑏
46 GPa 𝛼 𝜆 𝜖 𝑐𝑟 − - 𝜆 𝜃 𝑐𝑟 − - 𝑘 𝜑 Table 2: Numerical parameters for the beam-particle model.
The loadings are chosen so as to obtain different crackingpaths: • Tension along the 𝑦 -axis (see figure 5); • Compression along the 𝑦 -axis with unrestrainedboundary conditions along the 𝑥 -axis (see figure 6); • Bi-tension along the 𝑥 -axis and 𝑦 -axis (see figure 7); • Simple shear (see figure 8);9 C ] = ⎡⎣ .
13 8 . − . .
42 42 .
65 0 . − .
81 0 .
48 30 . ⎤⎦ [ C 𝑖𝑠𝑜 ] = ⎡⎣ .
92 9 .
39 0 . .
39 41 .
92 0 . .
00 0 .
00 32 . ⎤⎦ 𝜅 = 25 .
65 GPa, 𝜇 = 16 .
27 Gpa 𝐸 = 38 .
82 GPa, 𝜈 = 0 . (a) 10 × [ C ] = ⎡⎣ .
72 7 . − . .
90 38 . − . − . − .
10 29 . ⎤⎦ [ C 𝑖𝑠𝑜 ] = ⎡⎣ .
16 8 .
32 0 . .
32 38 .
16 0 . .
00 0 .
00 29 . ⎤⎦ 𝜅 = 23 .
24 GPa, 𝜇 = 14 .
92 Gpa 𝐸 = 36 .
35 GPa, 𝜈 = 0 . (b) 100 × [ C ] = ⎡⎣ .
06 8 .
09 0 . .
09 38 . − . . − .
03 29 . ⎤⎦ [ C 𝑖𝑠𝑜 ] = ⎡⎣ .
82 8 .
34 0 . .
34 37 .
82 0 . .
00 0 .
00 29 . ⎤⎦ 𝜅 = 23 .
08 GPa, 𝜇 = 14 .
74 Gpa 𝐸 = 35 .
98 GPa, 𝜈 = 0 . (c) 300 × Figure 3: Polar histograms of the beams orientation for different mesh densities, elasticity tensors and associated isotropic tensors (in GPa,Kelvin representation).
Mesh density (number of particles per direction) R e l a t i v e d i s t an c e t o i s o t r op y i s o ( % ) Figure 4: Evolution of the relative distance to isotropy (5.1) withthe mesh density. • Willam loading (Willam et al., 1989) consisting on asimple tension along the 𝑦 -axis pursued up to the ten-sile stress followed by a combination of bi-tension andshear with an increase of the strain components 𝜖 𝑥𝑥 , 𝜖 𝑦𝑦 and 𝜖 𝑥𝑦 in the proportions 1.5/1.0/0.5 (see fig-ure 9). This numerical test was proposed to induce arotation of the principal stress/strain directions lead-ing to a misalignment with the material axes of or-thotropy associated with the initial crack direction.For each of them, we first observe a development of dif-fuse micro-cracks and then a localization of these micro-cracks leading to the propagation of the macro-crack.In order to follow the breakage of the specimen, the evo-lution of the ratio of broken beams is plotted in compari-son with the macroscopic response. This ratio is initiallyworth 0, when the specimen is intact, and reaches 1 at theend of loading when the specimen can no longer sustainany effort. The total number of broken beams is usuallyfar from the number of beams in the specimen.We can notice snap-back phenomena in the evolutioncurves of the ratio of broken beams. These snap-backsare related to avalanche breakage phenomena (Rinaldi andLai, 2007), classical for lattice simulations of quasi-brittlematerials, and are visible in the complete macroscopic re-sponses (see figure 5). These macroscopic responses aresmoothed on the figures 5 to 9 to get closer to the typeof data observed experimentally. However, the extractionof the effective elasticity tensors is not impacted by thispost-treatment of the curves.These loadings allow us to quantify the impact on the10 Strain yy (10 -4 ) S t r e ss yy ( M P a ) R a t i o o f b r o k en bea m s (-) Figure 5: Tensile loading: final crack pattern, macroscopic responses (complete and smoothed) and evolution of the ratio of broken beams(dashed line). -15 -10 -5 0
Strain yy (10 -4 ) -40-30-20-100 S t r e ss yy ( M P a ) R a t i o o f b r o k en bea m s (-) Figure 6: Compressive loading: final crack pattern, macroscopic response and evolution of the ratio of broken beams (dashed line).
Strain yy (10 -4 ) S t r e ss ( M P a ) R a t i o o f b r o k en bea m s (-) xxyy Figure 7: Bi-Tensile loading: final crack pattern, macroscopic response and evolution of the ratio of broken beams (dashed line).
Displacement u x (mm) F o r c e F x ( M N ) R a t i o o f b r o k en bea m s (-) Figure 8: Simple shear loading: final crack pattern, macroscopic response and evolution of the ratio of broken beams (dashed line). Strain yy (10 -4 ) -0.500.511.522.533.54 S t r e ss ( M P a ) R a t i o o f b r o k en bea m s (-) xxyyxy Figure 9: Willam loading: final crack pattern, macroscopic response and evolution of the ratio of broken beams (dashed line). effective (damaged) elasticity tensor of the presence ofmicro-cracks in different directions, which may nucleateor rotate.
For the previously defined loadings, the stiffness andcompliance tensors were extracted at different ratio of bro-ken beams. The evolution of Δ, which is an upper boundof the distance to the orthotropic class for an elasticity ten-sor as defined in corollary 3.2, is plotted in the figures 10and 11. It should be noted that it is not always possi-ble to extract the elasticity tensors for very high levels ofdamage.
Ratio of broken beams (-) E ff e c t i v e s t i ff ne ss t en s o r - ( % ) TensionCompressionBi-TensionSimple ShearWillam test
Figure 10: Upper-bound of distance to orthotropy for the effectivestiffness tensors
Even at high levels of damage, the elasticity tensors re-main close to the orthotropic symmetry class. This is eventruer in the case of the compliance tensors. It is inter-esting to note that in the case of Willam loading, whichresults in a non-orthotropic cracking pattern, the devia-tions from orthotropy remain small. The largest deviationobserved corresponds to the case of the stiffness tensor forbi-tension loading. This may be related to the coalescenceof the main two orthogonal cracks.One should note that the orthotropy of the effective elas-ticity tensor ˜ C = ˜ S − is a sufficient condition to define an Ratio of broken beams (-) E ff e c t i v e c o m p li an c e t en s o r - ( % ) TensionCompressionBi-TensionSimple ShearWillam test
Figure 11: Upper-bound of distance to orthotropy for the effectivecompliance tensors orthotropic damage tensor (see next), since the initial elas-ticity tensor C = S − is isotropic. In the bi-dimensional micromechanical approach sum-marized by Kachanov (1992), the material is consideredto be initially isotropic with an initial (undamaged) com-pliance tensor, S = 12 𝜇 J + 14 𝜅 ⊗ . (6.1)Damage is introduced using several networks of crackswith no interactions between each other. A family ofcracks 𝑝 is characterized by its orientation 𝑛𝑛𝑛 ( 𝑝 ) , which isthe normal to all the cracks in this family, a length 2 𝑙 ( 𝑝 ) and a density of micro-cracks 𝜔 ( 𝑝 ) = 𝜋𝑙 ( 𝑝 ) /𝐴 with 𝐴 thearea of the RAE.A crack density tensor is then introduced for the networkof all the families of cracks: 𝜔𝜔𝜔 = ∑︁ 𝑝 𝜔 ( 𝑝 ) 𝑛𝑛𝑛 ( 𝑝 ) ⊗ 𝑛𝑛𝑛 ( 𝑝 ) (6.2)Gibbs free enthalpy density of the cracked solid writes,12n case of open cracks (Kachanov, 1992): 𝜌𝜓 ⋆ = 12 𝜎𝜎𝜎 : ˜ S : 𝜎𝜎𝜎 = 12 𝜎𝜎𝜎 : S : 𝜎𝜎𝜎 + 1 𝐸 tr ( 𝜎𝜎𝜎 · 𝜔𝜔𝜔 · 𝜎𝜎𝜎 )= 12 𝜎𝜎𝜎 : S : 𝜎𝜎𝜎 + 12 𝐸 𝜎𝜎𝜎 : ( ⊗ 𝜔𝜔𝜔 + 𝜔𝜔𝜔 ⊗ ) : 𝜎𝜎𝜎 (6.3)with S the initial (isotropic) compliance tensor ( 𝐸 beingthe Young modulus of the isotropic uncracked solid) andwhere ˜ S is the effective compliance tensor.We recall here the harmonic decomposition of the com-pliance tensor (see remark 2.1):˜ S = ˜ S 𝑖𝑠𝑜 + 12 (︁ ⊗ d ′ (˜ S ) + d ′ (˜ S ) ⊗ )︁ + H (˜ S )= ˜ S 𝑖𝑠𝑜 + 12 (︁ ⊗ d ′ (˜ S ) + d ′ (˜ S ) ⊗ )︁ + H (˜ S ) (6.4)with d (˜ S ) = tr ˜ S (in the same way v (˜ S ) = tr ˜ S ) andwhere (see Eq. (2.2))˜ S 𝑖𝑠𝑜 = 12˜ 𝜇 J + 14˜ 𝜅 ⊗ , {︃ 𝜇 = (2 tr v (˜ S ) − tr d (˜ S )) , 𝜅 = tr d (˜ S ) . By comparing equations (6.3) and (6.4), it can be con-cluded that in the case of Kachanov micro-cracking the-ory (representing networks of non interacting open micro-cracks), one has ˜ S 𝑖𝑠𝑜 = S + 𝐸 tr 𝜔𝜔𝜔 ⊗ (with ⊗ = I = J + ⊗ ), d ′ ( S ) = 2 𝜔𝜔𝜔 ′ /𝐸 and the harmonic part ofthe compliance tensor H (˜ S ) vanishes.In order to test the validity of the property H (˜ S ) = 0,eventually when cracks interaction takes place, we plot infigures 12 to 16 the evolution of the three parts of the effec-tive compliance tensor ˜ S (computed thanks to the beam-particle method): • the isotropic part ˜ S 𝑖𝑠𝑜 , • the deviatoric dilatation part ( ⊗ d ′ (˜ S )+ d ′ (˜ S ) ⊗ ), • and the harmonic part H (˜ S ),according to the ratio of broken beams. One should re-member that the deviatoric dilatation part is always or-thotropic while the harmonic part can have the squaresymmetry or be zero.With these curves, we can distinguish four stages ofcracking: 𝑖 ) Initially, the medium is non-cracked and thereforeisotropic, which implies that the parts related to d ′ and H are null. 𝑖𝑖 ) At the beginning of loading, only the weakest beamsbreak without influence of the loading direction. Theorientation of the cracks therefore remains isotropicand here again the parts linked to d ′ and H are zero. Figure 12: Tensile loading: Evolution of the relative parts of theeffective compliance tensor with the ratio of broken beams.
Figure 13: Compressive loading: Evolution of the relative parts ofthe effective compliance tensor with the ratio of broken beams.
Figure 14: Bi-Tensile loading: Evolution of the relative parts of theeffective compliance tensor with the ratio of broken beams. Figure 15: Simple shear loading: Evolution of the relative parts ofthe effective compliance tensor with the ratio of broken beams.
Figure 16: Willam loading: Evolution of the relative parts of theeffective compliance tensor with the ratio of broken beams. 𝑖𝑖𝑖 ) After a certain level, the beams failure will be relatedto the direction of loading and the microcracks beginto orient themselves while remaining diffuse. The lossof isotropy implies that the d ′ -related part becomesnon-zero. On the other hand, the harmonic part re-mains null. This is a case of crack-induced orthotropy. 𝑖𝑣 ) Eventually, the cracks will start to interact or coa-lesce and the harmonic part becomes non-zero. Sowe deviate from the framework of Kachanov’s theory.Depending on the symmetry class of the harmonictensor, it can be induced orthotropy or no-symmetryinduced anisotropy.Those stages can be observed on the computed crackingpatterns of figures 17 to 21.We note that the results are very dependent on the load-ing studied. However, we can draw a general conclusionfrom the curves presented above: the effective compliancetensor is close to remain orthotropic even at very highlevels of cracking. Indeed, the harmonic part of the effec-tive compliance tensor remains negligible in most loadingcases. This last observation is consistent with Kachanovdilute (non interacting) micro-cracking theory, for whichthe harmonic part of the compliance tensor remains zerofor any network of diffuse and open micro-cracks. Fur-thermore, it is most often satisfied here even in cases oflocalizing and interacting cracks.When this harmonic part is not negligible, as in the casesof compression (see figure 13) or shear (see figure 15), theeffective compliance tensor does not necessarily becomesanisotropic without symmetry. Indeed, as it is observedin figure 11, ˜ S remains close to belong to the orthotropicsymmetry class (Δ <
7. Extraction of damage tensors
A natural definition of the damage variable from theeffective elasticity tensor leads to a fourth-order tensor(Chaboche, 1979; Leckie and Onat, 1980). However, theuse of a symmetric second-order tensor is common due tothe simplicity of its interpretation (Cordebois and Sidoroff,1982; Murakami, 1988).The definition of a second-order damage variable is jus-tified here since we have checked that the elasticity tensors,and thus the fourth-order damage tensors, remain close tothe orthotropic symmetric class – the symmetry class ofa generic second-order tensor– even at high levels of dam-age. Nonetheless, the definition of such a damage variableis not straightforward. After choosing between a definitionbased on the compliance tensor or the stiffness tensor, onemust select the degradations rendered by this damage vari-able, which can be related either to the bulk modulus, tothe shear modulus or possibly to a combination of both.To validate the definition of a symmetric second-ordertensor as damage variable, it must be checked that it is14 a) 60 % (b) 70 % (c) 80 % (d) 90 % (e) 100 %
Figure 17: Tensile loading: evolution of the crack pattern with the ratio of broken beams. (a) 20 % (b) 40 % (c) 60 % (d) 80 %
Figure 18: Compressive loading: evolution of the crack pattern with the ratio of broken beams. (a) 45 % (b) 55 % (c) 65 % (d) 75 % (e) 85 %
Figure 19: Bi-Tensile loading: evolution of the crack pattern with the ratio of broken beams. (a) 20 % (b) 40 % (c) 60 % (d) 80 %
Figure 20: Simple shear loading: evolution of the crack pattern with the ratio of broken beams. (a) 10 % (b) 30 % (c) 50 % (d) 70 %
Figure 21: Willam loading: evolution of the crack pattern with the ratio of broken beams.
Following M. Kachanov, a natural definition for ananisotropic damage variable is related to the compliancetensor. More precisely, we propose here to define it fromthe difference between the invert of the effective bulk mod-ulus, very sensitivitive to damage, and the invert of theinitial bulk modulus,1˜ 𝜅 − 𝜅 = tr d (˜ S − S ) = tr tr (︁ ˜ S − S )︁ . We first define a dimensionless symmetric second-ordertensor, ΩΩΩ := 𝜅 (︁ d (˜ S ) − d ( S ) )︁ = 𝜅 tr (︁ ˜ S − S )︁ , whose trace is 𝜅/ ˜ 𝜅 − ≥
0, and of eigenvalues expected tobe positive but unbounded (0 ≤ Ω 𝑖 < ∞ ). A dimensionlesssymmetric second-order damage variable is obtained as ̂︀ D := ΩΩΩ ( + ΩΩΩ) − = ( + ΩΩΩ) − ΩΩΩ . One should note that the eigenvalues of ̂︀ D remain positiveand bounded by 1 provided those of ΩΩΩ are positive.The figures 22 to 26 show the evolutions of the com-ponents of the damage tensor ̂︀ D with the ratio of brokenbeams for the studied loadings. One can see that the val-ues of the components of ̂︀ D are not bounded by 0 and1. In view of the values obtained, this is in a non sur-prising manner also the case for the eigenvalues. There-fore, this definition based solely on the invert of the initialbulk modulus is is not compatible with a thermodynamicsframework. Indeed, it does not ensure the positivity of ΩΩΩ.The most problematic cases are the compressive loading(Fig. 23) and the simple shear loading (Fig. 25), for whichthe evolution of the shear modulus might certainly not beeasily linked to the evolution of the bulk modulus. A so-lution would be to introduce this shear modulus in thedefinition of the damage variable itself. However, this so-lution would be difficult to implement because one wouldneed to introduce some parameters to combine both thebulk modulus part and the shear modulus part, parame-ters that will have to be identified later. A second – alternative and, one will see, preferred –anisotropic damage variable can be defined as related to
Figure 22: Tensile loading: Evolution of the components of the dam-age tensor ̂︀ 𝐃 with the ratio of broken beams. Figure 23: Compressive loading: Evolution of the components of thedamage tensor ̂︀ 𝐃 with the ratio of broken beams. Figure 24: Bi-Tensile loading: Evolution of the components of thedamage tensor ̂︀ 𝐃 with the ratio of broken beams. Figure 25: Simple shear loading: Evolution of the components of thedamage tensor ̂︀ 𝐃 with the ratio of broken beams. Figure 26: Willam loading: Evolution of the components of the dam-age tensor ̂︀ 𝐃 with the ratio of broken beams. the effective bulk modulus instead of its inverse, more pre-cisely as related to the difference 𝜅 − ˜ 𝜅 = 14 tr d ( C − ˜ C ) = 14 tr tr (︁ C − ˜ C )︁ We propose to define the damage variable as the dimen-sionless symmetric second order tensor: D := d ( C ) − (︁ d ( C ) − d ( ˜ C ) )︁ = (︁ d ( C ) − d ( ˜ C ) )︁ d ( C ) − as d ( C ) = tr C = 2 𝜅 is spherical due to initial isotropy,so that the previous commutativity property holds, with D = 12 𝜅 (︁ d ( C ) − d ( ˜ C ) )︁ = 12 𝜅 tr (︁ C − ˜ C )︁ (7.1)One must check that the damage tensor D thus definedhas positive eigenvalues bounded by 1 (0 ≤ 𝐷 𝑖 ≤ D with the ratio of broken beams forthe studied loadings. D D Figure 27: Tensile loading: Evolution of the components of the dam-age tensor 𝐃 with the ratio of broken beams. We can observe that, this time, the components of thesecond-order tensorial damage variable remain postive andbounded by 1. In most cases, a maximum damage closeto 1 is reached on at least one of the tensor componentsbefore the end of the loading. This may explain why it isno longer possible to extract the effective elasticity tensorseven though the ratio of broken beams has not yet reached1. In pure proportional loading we expect a strict increasein the components of the tensor damage variable. Thisis not always the case. Indeed, in the case of bi-tension(Fig. 29) and simple shear (Fig. 30), the slight decreasesof the shear component 𝐷 is mainly related to local nonproportionality (and rotation of the principal axes of ten-sor D ). The drop is clear in the Willam loading case (seefigure 31). However, the damage eigenvalues 𝐷 , 𝐷 of D always increase for this loading as illustrated in figure 32.17 D D Figure 28: Compressive loading: Evolution of the components of thedamage tensor 𝐃 with the ratio of broken beams. D D Figure 29: Bi-Tensile loading: Evolution of the components of thedamage tensor 𝐃 with the ratio of broken beams. D D Figure 30: Simple shear loading: Evolution of the components of thedamage tensor 𝐃 with the ratio of broken beams. D D Figure 31: Willam loading: Evolution of the components of the dam-age tensor 𝐃 with the ratio of broken beams. D Figure 32: Willam loading: Evolution of the eigenvalues of the dam-age tensor 𝐃 with the ratio of broken beams. . Conclusion An upper bound to the distance to orthotropy has beenobtained for bi-dimensional elasticity tensors. It natu-rally introduces a second –instead of fourth– order tensorwhich models the medium orthotropy. This has allowed usto propose and measure second order tensorial damagesvariables fully representative of the effective anisotropicdegradation due to complex cracking patterns. We haveperformed 2D discrete simulations with a beam-particlemodel of initially isotropic Representative Area Elements,allowing for strong interactions between cracks (up to theircoalescence and complete failure). We have analyzed thesesimulations in a systematic manner.The results associated with the discrete computationsconfirm the accuracy of the orthotropic approximation ofthe elasticity tensor for a 2D cracked medium, under bothproportional or non-proportional loading cases. They fullyjustify the use of a single second order damage tensor vari-able, instead of a fourth order one (Chaboche, 1979; Leckieand Onat, 1980) or of two second order damage variables(Desmorat and Desmorat, 2016), even in the strong crackinteraction case.Two definitions, in stiffness and in compliance, for asymmetric second-order damage tensor have been pro-posed and studied in this paper. They are both basedon the evolution of the bulk modulus (or of its inverse).It has been observed that the damage variable derivedfrom the stiffness tensor has the necessary properties tobe a suitable candidate for the formulation of a continuousanisotropic damage model in 2D from discrete simulations.The methodology followed does not fix a priori a par-ticular working basis, it makes accessible all the tensorial(multiaxial) components measurement of the damage vari-able are any time step, therefore of its evolution.
Appendix A. Harmonic square roots of bi-dimensional harmonic tensors
The goal of this appendix is to provide a simple proofof the following theorem, using results in (Desmorat et al.,2020).
Theorem Appendix A.1.
Any fourth-order harmonicbi-dimensional tensor H can be written either as an har-monic product h * h (an harmonic square), or as − k * k (the opposite of an harmonic square), where h and k aresecond-order bi-dimensional harmonic tensors. To prove this result, we will recall first that there is anisomorphism H ∼ = h between bi-dimensional harmonic ten-sors H of order 𝑛 and homogeneous harmonic polynomialsh of degree 𝑛 in two variables 𝑥 and 𝑦 . Rather than 𝑥 and 𝑦 , one can use the complex variables 𝑧 = 𝑥 + 𝑖𝑦 and¯ 𝑧 = 𝑥 − 𝑖𝑦 . Then, any homogeneous harmonic polyno-mial of degree 𝑛 writes as h = Re (¯ 𝑧 𝑧 𝑛 ), where 𝑧 ∈ ℂ and the harmonic product between two homogeneous har-monic polynomials h and h , of respective degree 𝑛 and 𝑛 translates intoh * h = 12 Re ( 𝑧 𝑧 𝑧 𝑛 + 𝑛 ) . For instance, the deviatoric second order tensors h = (︂ 𝑎 𝑏 𝑏 − 𝑎 )︂ and k = (︂ 𝑎 𝑏 𝑏 − 𝑎 )︂ are represented respectively by the harmonic polynomialsh = Re (¯ 𝑧 𝑧 ) , and h = Re (¯ 𝑧 𝑧 ) , where 𝑧 = 𝑎 + 𝑖𝑏 and 𝑧 = 𝑎 + 𝑖𝑏 , and an har-monic fourth order tensor H (represented by its Kelvinmatrix (3.1)) corresponds to the homogeneous harmonicpolynomial h = Re (¯ 𝑧 𝑧 ), where 𝑧 = 𝐻 + 𝑖𝐻 . Thus, the harmonic square tensorial equations H = h * h = − k * k translate, in terms of harmonic homogeneous polynomials,as h = Re (¯ 𝑧 𝑧 ) = Re (¯ 𝑧 𝑧 ) = − Re (¯ 𝑧 𝑧 ) . and the solutions are provided by roots of the algebraicequations 𝑧 = − 𝑧 = 𝑧 . (A.1) Remark
Appendix A.2 . Note that both equations H = h * h and H = − k * k have exactly two opposite solutions,when H ̸ = 0. References
Andre´, D., Iordanoff, I., Charles, J.l., Ne´auport, J., 2012. Discreteelement method to simulate continuous material by using the co-hesive beam model. Computer methods in applied mechanics andengineering 213, 113–125.Backus, G., 1970. A geometrical picture of anisotropic elastic tensors.Reviews of geophysics 8, 633–671.Bagi, K., 1996. Stress and strain in granular assemblies. Mechanicsof materials 22, 165–177.Chaboche, J.L., 1979. Le concept de contrainte effective applique´ a`l’e´lasticite´ et a` la viscoplasticite´ en pre´sence d’un endommagementanisotrope, in: Boehler, J.P. (Ed.), Colloque Int. CNRS 295, Vil-lard de Lans, Martinus Nijhoff Publishers and Editions du CNRS,1982. pp. 737–760.Challamel, N., Picandet, V., Pijaudier-Cabot, G., 2015. From dis-crete to nonlocal continuum damage mechanics: Analysis of alattice system in bending using a continualized approach. Inter-national Journal of Damage Mechanics 24, 983–1012.Cordebois, J., Sidoroff, F., 1982. Endommagement anisotrope ene´lasticite´ et plasticite´. J. Meca. Th. Appl., Special Volume , 45–65.Cundall, P.A., Strack, O.D., 1979. A discrete numerical model forgranular assemblies. geotechnique 29, 47–65. ’Addetta, G.A., Kun, F., Ramm, E., 2002. On the applicationof a discrete model to the fracture process of cohesive granularmaterials. Granular Matter 4, 77–90.Delaplace, A., 2008. Mode´lisation discre`te applique´e au comporte-ment des mate´riaux et des structures. Me´moire d’habilitation a`diriger des recherches de l’Ecole Normale Supe´rieure de Cachan .Delaplace, A., Desmorat, R., 2007. Discrete 3d model as complimen-tary numerical testing for anisotropic damage. Int J Fract 148,115–128.Desmorat, B., Desmorat, R., 2015. Tensorial polar decomposition of2d fourth-order tensors. Comptes Rendus Me´canique 343, 471–475.Desmorat, B., Desmorat, R., 2016. Second order tensorial frameworkfor 2d medium with open and closed cracks. European Journal ofMechanics-A/Solids 58, 262–277.Desmorat, B., Olive, M., Auffray, N., Desmorat, R., Kolev, B., 2020.Computation of minimal covariants bases for 2d coupled consti-tutive laws. arXiv preprint arXiv:2007.01576 .Desmorat, R., 2006. Positivity of intrinsic dissipation of a class ofnonstandard anisotropic damage models. C. R. Mecanique 334,587–592.Desmorat, R., 2016. Anisotropic damage modeling of concrete ma-terials. International Journal of Damage Mechanics 25, 818–852.Desmorat, R., Gatuingt, F., Ragueneau, F., 2007. Nonlocalanisotropic damage model and related computational aspects forquasi-brittle materials. Engineering Fracture Mechanics 74, 1539–1560.Halm, D., Dragon, A., 1998. An anisotropic model of damageand frictional sliding for brittle materials. European Journal ofMechanics-A/Solids 17, 439–460.Herrmann, H.J., Hansen, A., Roux, S., 1989. Fracture of disordered,elastic lattices in two dimensions. Physical Review B 39, 637.Hrennikoff, A., 1941. Solution of problems of elasticity by the frame-work method. J. appl. Mech. .Jivkov, A., 2014. Structure of micro-crack population and damageevolution in quasi-brittle media. Theoretical and Applied FractureMechanics 70, 1–9.Kachanov, M., 1992. Effective elastic properties of cracked solids:critical review of some basic concepts. Appl. Mech. Rev. 45, 304–335.Leckie, F., Onat, E., 1981. Tensorial nature of damage measuring in-ternal variables, in: Physical non-linearities in structural analysis.Springer, pp. 140–155.Leckie, F.A., Onat, E.T., 1980. Tensorial nature of damage mea-suring internal variables. J. Hult and J. Lemaitre eds, SpringerBerlin. chapter Physical Non-Linearities in Structural Analysis.pp. 140–155.Lemaitre, J., Chaboche, J.L., 1985. Me´canique des mate´riaux solides.Dunod, english translation 1990 ’Mechanics of Solid Materials’Cambridge University Press.Lemaitre, J., Desmorat, R., 2005. Engineering damage mechanics:ductile, creep, fatigue and brittle failures. Springer Science &Business Media.Mazars, J., 1984. Application de la me´canique de l’endommagementau comportement non line´aire et a` la rupture du be´ton de struc-ture. Ph.D. thesis. Universite´ Pierre et Marie Curie - Paris 6.Mazars, J., Berthaud, Y., Ramtani, S., 1990. The unilateral be-haviour of damaged concrete. Engineering Fracture Mechanics35, 629–635.Meguro, K., Hakuno, M., 1989. Fracture analyses of concrete struc-tures by the modified distinct element method. Doboku GakkaiRonbunshu 1989, 113–124.Murakami, S., 1988. Mechanical modeling of material damage.ASME J. Appl. Mech. 55, 280–286. doi: .Olive, M., Kolev, B., Desmorat, B., Desmorat, R., 2018a. Harmonicfactorization and reconstruction of the elasticity tensor. Journalof Elasticity 132, 67–101.Olive, M., Kolev, B., Desmorat, R., Desmorat, B., 2018b. Char-acterization of the symmetry class of an elasticity tensor usingpolynomial covariants. arXiv:1807.08996 [math.RT] .Oliver-Leblond, C., 2019. Discontinuous crack growth and tough- ening mechanisms in concrete: A numerical study based on thebeam-particle approach. Engineering Fracture Mechanics 207, 1–22.Poisson, S.D., 1828. Me´moire sur l’e´quilibre et le mouvement descorps e´lastiques. F. Didot.Ramtani, S., Berthaud, Y., Mazars, J., 1992. Orthotropic behaviorof concrete with directional aspects: modelling and experiments.Nuclear Engineering and design 133, 97–111.Rinaldi, A., 2013. Bottom-up modeling of damage in heterogeneousquasi-brittle solids. Continuum Mechanics and Thermodynamics25, 359–373.Rinaldi, A., Lai, Y.C., 2007. Statistical damage theory of 2d lat-tices: Energetics and physical foundations of damage parameter.International Journal of Plasticity 23, 1796–1825.Schlangen, E., Van Mier, J.G.M., 1992. Simple lattice model for nu-merical simulation of fracture of concrete materials and structures.Materials and Structures 25, 534–542.Schouten, J.A., 1989. Tensor analysis for physicists. Courier Corpo-ration.Spencer, A., 1970. A note on the decomposition of tensors intotraceless symmetric tensors. International Journal of EngineeringScience 8, 475–481.Tillemans, H.J., Herrmann, H.J., 1995. Simulating deformations ofgranular solids under shear. Physica A: Statistical Mechanics andits Applications 217, 261–288.Vannucci, P., 2005. Plane anisotropy by the polar method. Meccanica40, 437–454.Vannucci, P., Verchery, G., 2001. Stiffness design of laminates usingthe polar method. International Journal of Solids and Structures38, 9281–9894. doi: .Vassaux, M., Oliver-Leblond, C., Richard, B., Ragueneau, F., 2016.Beam-particle approach to model cracking and energy dissipationin concrete: Identification strategy and validation. ICement andConcrete Composites 70, 1–14.Vassaux, M., Richard, B., Ragueneau, F., Millard, A., 2015a. Reg-ularised crack behaviour effects on continuum modelling of quasi-brittle materials under cyclic loading. Engineering Fracture Me-chanics 149, 18–36.Vassaux, M., Richard, B., Ragueneau, F., Millard, A., Delaplace,A., 2015b. Lattice models applied to cyclic behavior descriptionof quasi-brittle materials: advantages of implicit integration. In-ternational Journal for Numerical and Analytical Methods in Ge-omechanics 39, 775–798.Verchery, G., 1982. Les invariants des tenseurs d’ordre 4 dutype de l’e´lasticite´, in: Mechanical behavior of anisotropicsolids/comportment Me´chanique des Solides Anisotropes.Springer, pp. 93–104.Vianello, M., 1997. An integrity basis for plane elasticity tensors.Archives of Mechanics 49, 197–208.Willam, K., Pramono, E., Sture, S., 1989. Fundamental issues ofsmeared crack models, in: Fracture of concrete and rock. Springer,pp. 142–157..Vassaux, M., Oliver-Leblond, C., Richard, B., Ragueneau, F., 2016.Beam-particle approach to model cracking and energy dissipationin concrete: Identification strategy and validation. ICement andConcrete Composites 70, 1–14.Vassaux, M., Richard, B., Ragueneau, F., Millard, A., 2015a. Reg-ularised crack behaviour effects on continuum modelling of quasi-brittle materials under cyclic loading. Engineering Fracture Me-chanics 149, 18–36.Vassaux, M., Richard, B., Ragueneau, F., Millard, A., Delaplace,A., 2015b. Lattice models applied to cyclic behavior descriptionof quasi-brittle materials: advantages of implicit integration. In-ternational Journal for Numerical and Analytical Methods in Ge-omechanics 39, 775–798.Verchery, G., 1982. Les invariants des tenseurs d’ordre 4 dutype de l’e´lasticite´, in: Mechanical behavior of anisotropicsolids/comportment Me´chanique des Solides Anisotropes.Springer, pp. 93–104.Vianello, M., 1997. An integrity basis for plane elasticity tensors.Archives of Mechanics 49, 197–208.Willam, K., Pramono, E., Sture, S., 1989. Fundamental issues ofsmeared crack models, in: Fracture of concrete and rock. Springer,pp. 142–157.