# Simulation of the Damage Process in Quasi-Brittle Materials by a Modified Finite Element Method Using the Consistent Embedded Discontinuity Formulation

Authors: J. Retama; A. G. Ayala
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

This paper investigates the variational finite element formulation and its numerical implementation of the damage evolution in solids, using a new discrete embedded discontinuity approach. For this purpose, the kinematically optimal symmetric (KOS) formulation, which guarantees kinematics, is consistently derived. In this formulation, rigid body motion of the parts in which the element is divided is obtained. To guarantee equilibrium at the discontinuity surfaces, the length of the discontinuity is introduced in the numerical implementation at elemental level. To illustrate and validate this approach, two examples, involving mode-I failure, are presented. Numerical results are compared with those reported from experimental tests. The presented discontinuity formulation shows a robust finite element method to simulate the damage evolution processes in quasi-brittle materials, without modifying the mesh topology when cohesive cracks propagate.


---

## Body

## 1. Introduction

The numerical simulation of the development and growth of cohesive cracks in quasi-brittle materials has posed a scientific challenge to structural engineering for decades. For this, two main paths can be distinguished: the modelling of the damage process in a smeared sense [1–3], within a continuum theory, and the modelling of cohesive cracks in a discrete form [4–6].Within the FEM context, two different types of models may be distinguished: one in which the inelastic deformation of the fracture process zone is smeared over a band of a certain width and the other in which the entire fracture process is lumped into a crack line whose behaviour is characterized by a stress-displacement law, known as cohesive zone models [1, 6–8].In the discrete modelling of fractures, the cohesive forces across a crack are modelled as a progressive degradation of the mechanical properties of the material (Figure1), and a stress-displacement relationship is used [9].Figure 1
Fictitious crack model.For the case of brittle material and considering a Rankine criterion failure, the crack begins to growth once the principal stressσI reaches the tensile strength σt0. The work energy GF needed to create a fully developed crack of unit area is as follows [1, 10]:(1)GF=∫0uσdu,where σ and u are the stress and displacement, respectively, defined on the cohesive crack (Figure 1).Recently, a different approach in which the crack crosses the element (Figure2) has been developed. This approach emerged as a mean to deal with general strain localization phenomena, of which the cohesive crack is a case (see [5, 11–17]).Figure 2
Embedded discontinuity models.All embedded discontinuity approaches have a common basic underlying; the process zone is incorporated into the variational formulation of the finite element model, by an enhancement of the relevant fields interpolation, which bears the required discontinuity only where it appears without affecting the main fields [4, 5, 18].Depending on the use of the constitutive model and the kinematics of displacements and strains, three embedded discontinuity formulations are identified: weak discontinuity, strong discontinuity, and discrete discontinuity [19].The weak discontinuity model considers that the displacement field is continuous and the strain field is discontinuous. In this model, the damage process zone is located into a bandwidth different from zero (Figure3(a)). When the bandwidth tends to zero, the strong discontinuity model is recovered, making the displacements field discontinuous, and the strains are unbounded across the discontinuity (Figure 3(b)).Figure 3
FEM mesh of a solid with an internal crack.
(a)(b)(c)The third model is the discrete discontinuity approach in which, unlike the strong discontinuity, the former considers two constitutive laws, one strain-stress for the bulk and other traction-jump for the cohesive crack zone [5, 18, 19] (Figure 3(c)).In this paper, the discrete embedded discontinuity is used to simulate damage process in quasi-brittle materials, like concrete and rocks [5, 17].


## 2. Embedded Discontinuities

The need of accurately modelling phenomena that present strong discontinuity in the fields of the state variables arises in many fields of the mechanics of materials and structures. The strong discontinuity approach (SDA) has proven effective in treating many problems in the mechanics of solids which considers discontinuous fields. Its basic idea is to introduce an enhancement in the interpolation of the relevant fields, which bears the required discontinuity only where it appears so that the model of the main field is not affected by the enhancement (Figure2).The idea is shared between methods that enrich the nodal variables (XFEM) and methods that use elements with embedded discontinuities (EED). The latter methods appear very convenient since the additional degrees of freedom enter the equations only locally and can be treated in the same way as inelastic strains [4, 5, 11–16, 20–22].There exist many formulations for EED [13, 19] that differ mainly in the numerical implementation of the method since the kinematic assumptions are basically the same as those originally proposed by Simo [23]. Furthermore, it appears that most of the formulations have been developed within the framework of constant stress elements or at least for elements with constant enhanced strains (see [24, 25]).Within the framework of the embedded discontinuity formulation, three different approaches may be identified attending to the constitutive models employed and the used kinematics of displacement and strain fields [19]. In this paper, the kinematically optimal symmetric (KOS) embedded discontinuity formulation is addressed, and attention is devoted to the kinematics of the enhancement displacement field. As explained in the paper by Jirásek [19], using the standard interpolation of the enhanced displacement field, a Ritz–Galerkin procedure produces a symmetric formulation but violates the equilibrium condition on the interface.


## 3. Kinematics of the Discontinuity

To establish the kinematics of the discontinuity interface, consider a body with domain Ω and boundary Γ (Figure4), where essential boundary conditions u¯ are imposed on Γu and natural boundary conditions t¯ are applied on Γt, such that Γ=Γu∪Γt and ∅=Γu∩Γt. Additionally, the body is crossed by a discontinuity Γd, dividing the body domain Ω in two subdomains Ω+ and Ω−. On the discontinuity boundary, the normal vector n is defined pointed towards Ω+.Figure 4
Continuum solid crossed by a discontinuityΓd.The strong discontinuity approach, addressed in this paper, can be regarded as a limit case of the weak discontinuity when the width of the localizations zonebd tends to zero, collapsing into a localized surface Γd. For this, when the strength of the material is reached, the deformation due to microcracking or other inelastic phenomena suddenly localizes into the discontinuity surface.


### 3.1. Displacement Jump

The displacement fieldu for the body crossed by a discontinuity, shown in Figure 4, is decomposed into a continuous and discontinuous part [16, 22, 26].(2)ux=u^x+ℋΓdxux,where u^ and u are smooth and continuous functions defined on Ω, with Cartesian coordinates x, and ℋΓd is the Heaviside function centred at the discontinuity surface defined as(3)ℋΓdx=αd,∀x∈Ω+,1−αd,∀x∈Ω−.The parameterαd, 0≤αd≤1, defines the way as the jump is transmitted: for αd=0, the jump is transferred to Ω+, while for αd=1 it is transferred to Ω−.The corresponding infinitesimal strain field can be derived by taking the symmetric gradient of equation (2), and this is given as(4)ε=∇su^+ℋΓd∇su+∇ℋΓd⊗us=ε^+δΓdu⊗ns,where δΓd is the Dirac-delta distribution centred at the discontinuity Γd, satisfying(5)∫ΩδΓdϕdΩ=∫ΓdϕdΓ.The strain field, given by equation (4), is decomposed into a bounded part ε^ plus an unbounded part, irregular part, given by δΓdu⊗ns. The unbounded strains may be considered in two different ways: (1) regularizing the Dirac-delta by means of a bandwidth parameter close to zero bd≈0 (see [15, 22, 27]) or (2) with nonregularization which implies taking advantage of Dirac-delta function defined in equation (5) (see [26, 28, 29]).


### 3.2. Kinematic Interpolation

The displacement and strain fields will be constructed only from kinematic considerations like in [12, 14]. For this, consider a one-dimensional finite element crossed by the discontinuity Γd, as it is shown in Figure 5.Figure 5
One-dimensional finite element with a discontinuity.It is assumed thatΩ+ undergoes a rigid body incremental displacement u˜ with respect to Ω−. To obtain the same incremental displacement derivatives on both subdomains, the interpolations represented in Figure 5 by lines 1 and 2 may be adopted for Ω− and Ω+, respectively. Their mathematical expressions are as follows:(6)u1=Nu−01u˜for Ω− and(7)u2=Nu−01u˜+u˜for Ω+. The purpose of matrix ℋΓd=01T, stated in equations (6) and (7), is derived from Figure 5. For every point in Ω− or Ω+, the strain field is defined by(8)ε=∇su=Bu−01u˜=Bu−ℋΓdu˜,where B is the usual compatibility matrix, relating displacements and strains [30–32], and the matrix ℋΓd is defined in equation (3).For a two-dimensional case, these concepts can be properly generalized by introducing the vectoru containing the components of jump displacements associated with the discontinuity Γd:(9)u=u˜v˜.For this case, the matrixℋΓd is derived by a straightforward generalization of its counterpart in equations (6) and (7). Its form for a two-dimensional case is written as(10)ℋΓd=ℋΓd1⋮ℋΓdn,where n is the number of nodes in the finite element, and each of the submatrices ℋΓdi depends on the position of node i relative to the discontinuity given as(11)ℋΓdi=αdI,∀x∈Ω+,1−αdI,∀x∈Ω−,where I is an identity matrix with dimensions (nst×nst) and nst is the space dimension of the problem.With respect to the displacement field, this is derived from equations (6) and (7) for whole domain Ω as(12)u=Nu^+ℋΓd−NℋΓdu,and taking Nc as(13)Nc=ℋΓd−NuℋΓd,equation (12) is written in the following final form:(14)u=Nu^+Ncu.The strain kinematic is derived, in a finite element framework, by means of the displacement field given by equation (14) as(15)ε=Bu^+Bcu,where the matrix Bc is derived by applying the strain-displacement operator to equation (13), and this is given as(16)Bc=∇sNc.Equations (14) and (15) define the displacement and strain fields, respectively, corresponding to the discrete approach of the present embedded discontinuity formulation [14, 33].For the case of the constant strain triangle finite element and considering a uniform deformation along the internal interface (Figure6), equation (16) is rewritten in the form(17)Bc=−∇sφx,where φx contains the shape functions corresponding to nodes in Ω+ (see [27, 29, 33]) defined as(18)φx=∑i=1n+Ni+,where n+ denotes the number of nodes belonging to the subdomain Ω+.Figure 6
Constant strain triangle with a displacement jump.


## 3.1. Displacement Jump

The displacement fieldu for the body crossed by a discontinuity, shown in Figure 4, is decomposed into a continuous and discontinuous part [16, 22, 26].(2)ux=u^x+ℋΓdxux,where u^ and u are smooth and continuous functions defined on Ω, with Cartesian coordinates x, and ℋΓd is the Heaviside function centred at the discontinuity surface defined as(3)ℋΓdx=αd,∀x∈Ω+,1−αd,∀x∈Ω−.The parameterαd, 0≤αd≤1, defines the way as the jump is transmitted: for αd=0, the jump is transferred to Ω+, while for αd=1 it is transferred to Ω−.The corresponding infinitesimal strain field can be derived by taking the symmetric gradient of equation (2), and this is given as(4)ε=∇su^+ℋΓd∇su+∇ℋΓd⊗us=ε^+δΓdu⊗ns,where δΓd is the Dirac-delta distribution centred at the discontinuity Γd, satisfying(5)∫ΩδΓdϕdΩ=∫ΓdϕdΓ.The strain field, given by equation (4), is decomposed into a bounded part ε^ plus an unbounded part, irregular part, given by δΓdu⊗ns. The unbounded strains may be considered in two different ways: (1) regularizing the Dirac-delta by means of a bandwidth parameter close to zero bd≈0 (see [15, 22, 27]) or (2) with nonregularization which implies taking advantage of Dirac-delta function defined in equation (5) (see [26, 28, 29]).


## 3.2. Kinematic Interpolation

The displacement and strain fields will be constructed only from kinematic considerations like in [12, 14]. For this, consider a one-dimensional finite element crossed by the discontinuity Γd, as it is shown in Figure 5.Figure 5
One-dimensional finite element with a discontinuity.It is assumed thatΩ+ undergoes a rigid body incremental displacement u˜ with respect to Ω−. To obtain the same incremental displacement derivatives on both subdomains, the interpolations represented in Figure 5 by lines 1 and 2 may be adopted for Ω− and Ω+, respectively. Their mathematical expressions are as follows:(6)u1=Nu−01u˜for Ω− and(7)u2=Nu−01u˜+u˜for Ω+. The purpose of matrix ℋΓd=01T, stated in equations (6) and (7), is derived from Figure 5. For every point in Ω− or Ω+, the strain field is defined by(8)ε=∇su=Bu−01u˜=Bu−ℋΓdu˜,where B is the usual compatibility matrix, relating displacements and strains [30–32], and the matrix ℋΓd is defined in equation (3).For a two-dimensional case, these concepts can be properly generalized by introducing the vectoru containing the components of jump displacements associated with the discontinuity Γd:(9)u=u˜v˜.For this case, the matrixℋΓd is derived by a straightforward generalization of its counterpart in equations (6) and (7). Its form for a two-dimensional case is written as(10)ℋΓd=ℋΓd1⋮ℋΓdn,where n is the number of nodes in the finite element, and each of the submatrices ℋΓdi depends on the position of node i relative to the discontinuity given as(11)ℋΓdi=αdI,∀x∈Ω+,1−αdI,∀x∈Ω−,where I is an identity matrix with dimensions (nst×nst) and nst is the space dimension of the problem.With respect to the displacement field, this is derived from equations (6) and (7) for whole domain Ω as(12)u=Nu^+ℋΓd−NℋΓdu,and taking Nc as(13)Nc=ℋΓd−NuℋΓd,equation (12) is written in the following final form:(14)u=Nu^+Ncu.The strain kinematic is derived, in a finite element framework, by means of the displacement field given by equation (14) as(15)ε=Bu^+Bcu,where the matrix Bc is derived by applying the strain-displacement operator to equation (13), and this is given as(16)Bc=∇sNc.Equations (14) and (15) define the displacement and strain fields, respectively, corresponding to the discrete approach of the present embedded discontinuity formulation [14, 33].For the case of the constant strain triangle finite element and considering a uniform deformation along the internal interface (Figure6), equation (16) is rewritten in the form(17)Bc=−∇sφx,where φx contains the shape functions corresponding to nodes in Ω+ (see [27, 29, 33]) defined as(18)φx=∑i=1n+Ni+,where n+ denotes the number of nodes belonging to the subdomain Ω+.Figure 6
Constant strain triangle with a displacement jump.


## 4. Discrete Damage Model

For modelling the evolution of material damage where mode-I type failure is dominant, a cohesive crack model is used (see [1, 6, 7, 34–36]). Microcracking and plastic flow around a macroscopic crack tip are modelled as equivalent tractions forces on crack faces (Figure 7).Figure 7
Softening curve for a cohesive crack model.A discontinuity is introduced when the maximum principal stressσI exceeds the tensile strength of the material. This mode-I criterion is commonly used for quasi-brittle materials (see [16, 26, 29]). Since the material around the discontinuity remains elastic, two constitutive models are used: one linear elastic for the bulk and other discrete damage for the discontinuity surface.For this discrete damage model, the Helmholtz free energy densityψu,ωκ, representing the free energy per unit area [16, 26] is expressed as(19)ψu,ωκ=1−ωκψ0u,with u already defined in Section 3.2, κ is an internal scalar variable which can be considered as an equivalent jump κ=ueq, and ψ0 is the linear elastic free energy given by(20)ψ0u=12uTTelu,where Tel is the elastic constitutive tensor. Equation (20) is valid for a discontinuity with a linear elastic behaviour. For adiabatic problems, the Helmholtz free energy, given by equation (19), must satisfy the Clausius–Duhem inequality (see [2, 37, 38]):(21)D=tTu˙−ψ˙≥0.Equation (21) is a dissipative potential with D being the rate of dissipation energy and t the tractions vector across the discontinuity. Expanding the dissipative potential equation (21) and grouping similar terms,(22)D=t−∂ψ∂uTu˙+ψ0ω˙≥0.Use of standard thermodynamics arguments, in equation (21), leads to the following:(1)
Constitutive equation(23)t=∂ψ∂u=1−ωTelu.(2)
In addition, the dissipation equation is reduced to(24)D=ψ0ω˙≥0.Since the elastic free energy, defined per unit area, is positive, it may be seen from equation (24) that the rate of change of the damage parameter ω cannot be negative, i.e., monotonic damage process.


### 4.1. Tangent Constitutive Tensor

Transformation of equation (23) into its incremental form, by its differentiation with respect to time, leads to the tangent constitutive tensor expressed as follows:(25)t˙=1−ωTelu˙−ω˙Tel=1−ωTelu˙−ω˙tel,where tel is the elastic traction vector and(26)ω˙=dωdκ∂κ∂uu˙=dωdκ∂κ∂unu˙n+∂κ∂usu˙s.Substitution of equation (26) into equation (25) leads to(27)t˙=1−ωTel−dωdκtel⊗∂κ∂u⋅u˙.If unloading takes place, the rate of damage is zero, i.e.,ω=˙0, and equation (27) reduces to(28)t˙=1−ωTelu˙.


### 4.2. Yield Function

Development of a discrete damage model requires the definition of a yield function for describing the loading state at the discontinuity surface (loading or unloading/reloading). The functions which describe this process is defined as(29)f=ueq−κ≤0,where ueq is an equivalent jump obtained from the displacement jump vector u. The symbols ⋅ are the McAuley brackets denoting that only the positive part of ueq is considered. The internal scalar variable κ is a history parameter equal to the largest value of ueq, and this is given as(30)κ=κu=maxueq.


### 4.3. Loading/Unloading Conditions

Iff<0, the deformation process is reversible and the dissipation rate must be zero, which implies that ω=˙0. If f=0, the damage can grow in a way that, in a subsequent state, f remains zero ω˙>0⇒f˙=0. Both cases can be conveniently expressed by the usual Kuhn–Tucker conditions [38],(31)ω˙≥0,f≤0,ω˙f=0,together with the consistency condition(32)ω˙f˙=0.


## 4.1. Tangent Constitutive Tensor

Transformation of equation (23) into its incremental form, by its differentiation with respect to time, leads to the tangent constitutive tensor expressed as follows:(25)t˙=1−ωTelu˙−ω˙Tel=1−ωTelu˙−ω˙tel,where tel is the elastic traction vector and(26)ω˙=dωdκ∂κ∂uu˙=dωdκ∂κ∂unu˙n+∂κ∂usu˙s.Substitution of equation (26) into equation (25) leads to(27)t˙=1−ωTel−dωdκtel⊗∂κ∂u⋅u˙.If unloading takes place, the rate of damage is zero, i.e.,ω=˙0, and equation (27) reduces to(28)t˙=1−ωTelu˙.


## 4.2. Yield Function

Development of a discrete damage model requires the definition of a yield function for describing the loading state at the discontinuity surface (loading or unloading/reloading). The functions which describe this process is defined as(29)f=ueq−κ≤0,where ueq is an equivalent jump obtained from the displacement jump vector u. The symbols ⋅ are the McAuley brackets denoting that only the positive part of ueq is considered. The internal scalar variable κ is a history parameter equal to the largest value of ueq, and this is given as(30)κ=κu=maxueq.


## 4.3. Loading/Unloading Conditions

Iff<0, the deformation process is reversible and the dissipation rate must be zero, which implies that ω=˙0. If f=0, the damage can grow in a way that, in a subsequent state, f remains zero ω˙>0⇒f˙=0. Both cases can be conveniently expressed by the usual Kuhn–Tucker conditions [38],(31)ω˙≥0,f≤0,ω˙f=0,together with the consistency condition(32)ω˙f˙=0.


## 5. Variational and Numerical Formulation

### 5.1. Variational Formulation

For this, consider a solid divided by a discontinuityΓd in Ω− and Ω+ subdomains, as shown in Figure 8, where the normal vector n is placed on the discontinuity surface, and the relative displacement between both subdomains is defined by the jump u in the displacement field u.Figure 8
Continuum solid divided by the discontinuityΓd.The energy functional for this solid is given by(33)Πu,u=∫ΩΨε−uTbvdΩ−∫ΓtuTt¯dΓ+∫ΓdψudΓ,where ψu is an additional term which considers the elastic free energy density defined on the discontinuity derived from the field u. For a discontinuity with an elastic behaviour, this function is given by(34)ψu=∫0utTdu,and with t being the tractions vector derived from a cohesive traction-jump constitutive relation, t=Tu. Furthermore, for a discontinuity with a linear elastic behaviour, ψu is written as(35)ψu=12uTTu.Applying the first variation to the energy functional defined in equation (33), the next Euler–Lagrange equations are recovered:(36a)∇Tσ+bv=0inΩ,(36b)σTν−t=0inΓd.Equations (36a)-(36b) represent the internal equilibrium and natural boundary conditions of the body, complemented by the following equilibrium conditions across the discontinuity:(37)t=−σ+Tn+t=σ−Tn−onΓdand the kinematic condition(38)u=u¯onΓu.Equation (38) states the essential boundary conditions of the variational problem.


### 5.2. Numerical Formulation

Once the functional of the discrete approach of the embedded discontinuity formulation is defined, the discretization of the domain by the finite element method is carried out by means of the displacement and strain kinematic, developed in Section3. This is given as(39)u=Nu^+Ncu˜,(40)ε=∇su=Bu^+Bcu˜,where u^ and u˜ are the regular and not regular displacements, respectively. The matrix N contains the standard shape functions. Shape functions associated with nodes belonging to Ω+ are contained in matrix Nc.As for standard finite element method, matricesB and Bc, given in equation (40), contain derivatives of the corresponding shape functions in matrices N and Nc.Substitution of equations (39) and (40) into the energy functional defined in equation (33), and its corresponding minimization, leads to the algebraic system of equation,(41)Ku^u^Ku^u˜Ku^u˜TKu˜u˜u^u˜=fext0,with(42)Ku^u^=∫ΩBTCBdΩ,(43)Ku^u˜=∫ΩBTCBcdΩ,(44)Ku˜u˜=∫ΩBcTCBcdΩ+∫ΓdNcTTNcdΓ,and external forces vector(45)fext=∫ΩNTbvdΩ+∫ΓtNTt¯dΓ.Here, the coefficient matrix is symmetric and unlike other discontinuity formulations [26–28, 33], there are nonzeros along its leading diagonal. The inelastic response is related to the term ∫ΓtNTt¯dΓ of the submatrix Ku˜u˜.


## 5.1. Variational Formulation

For this, consider a solid divided by a discontinuityΓd in Ω− and Ω+ subdomains, as shown in Figure 8, where the normal vector n is placed on the discontinuity surface, and the relative displacement between both subdomains is defined by the jump u in the displacement field u.Figure 8
Continuum solid divided by the discontinuityΓd.The energy functional for this solid is given by(33)Πu,u=∫ΩΨε−uTbvdΩ−∫ΓtuTt¯dΓ+∫ΓdψudΓ,where ψu is an additional term which considers the elastic free energy density defined on the discontinuity derived from the field u. For a discontinuity with an elastic behaviour, this function is given by(34)ψu=∫0utTdu,and with t being the tractions vector derived from a cohesive traction-jump constitutive relation, t=Tu. Furthermore, for a discontinuity with a linear elastic behaviour, ψu is written as(35)ψu=12uTTu.Applying the first variation to the energy functional defined in equation (33), the next Euler–Lagrange equations are recovered:(36a)∇Tσ+bv=0inΩ,(36b)σTν−t=0inΓd.Equations (36a)-(36b) represent the internal equilibrium and natural boundary conditions of the body, complemented by the following equilibrium conditions across the discontinuity:(37)t=−σ+Tn+t=σ−Tn−onΓdand the kinematic condition(38)u=u¯onΓu.Equation (38) states the essential boundary conditions of the variational problem.


## 5.2. Numerical Formulation

Once the functional of the discrete approach of the embedded discontinuity formulation is defined, the discretization of the domain by the finite element method is carried out by means of the displacement and strain kinematic, developed in Section3. This is given as(39)u=Nu^+Ncu˜,(40)ε=∇su=Bu^+Bcu˜,where u^ and u˜ are the regular and not regular displacements, respectively. The matrix N contains the standard shape functions. Shape functions associated with nodes belonging to Ω+ are contained in matrix Nc.As for standard finite element method, matricesB and Bc, given in equation (40), contain derivatives of the corresponding shape functions in matrices N and Nc.Substitution of equations (39) and (40) into the energy functional defined in equation (33), and its corresponding minimization, leads to the algebraic system of equation,(41)Ku^u^Ku^u˜Ku^u˜TKu˜u˜u^u˜=fext0,with(42)Ku^u^=∫ΩBTCBdΩ,(43)Ku^u˜=∫ΩBTCBcdΩ,(44)Ku˜u˜=∫ΩBcTCBcdΩ+∫ΓdNcTTNcdΓ,and external forces vector(45)fext=∫ΩNTbvdΩ+∫ΓtNTt¯dΓ.Here, the coefficient matrix is symmetric and unlike other discontinuity formulations [26–28, 33], there are nonzeros along its leading diagonal. The inelastic response is related to the term ∫ΓtNTt¯dΓ of the submatrix Ku˜u˜.


## 6. Tractions Equilibrium

The interpolation of the displacement jump field by the finite element method of the (KOS) formulation considers the shape functionNc. This kind of interpolation is a proper way to consider that the localization band width bd is different from zero; but, when this band collapses into a line, i.e., discrete discontinuity, the matrix Nc has the meaning of the identity matrix which leads to consider tractions at the discontinuity surface in the local framework.The case, whenbd=0, may lead to a nonproper enforcement of tractions continuity which only satisfies equilibrium for cases depending in the geometry of the element. To show this, consider the element of Figure 9.Figure 9
Projection of the discontinuity to the opposite side to node 3.The discontinuity crossing the element is parallel to side1−2 with length l12; n is the normal to the discontinuity, and h is the height of the triangle. The angle β is measured from the x axis to the side 1−2, and α is the angle between the x axis and the normal to the referred side.The area of the element isA=l12h/2; hence ,2A=l12h, and from the geometry of the element and the angle β, the following relations are derived:(46)cosβ=x2−x1l12,sinβ=y2−y1l12.Furthermore, for the element shown in Figure9, the matrix Bc is defined as(47)Bc=−12Ay1−y200x2−x1x2−x1y1−y2.Rewritten equation (47) by means of relationships derived from the geometry of the element and using the angle α instead of β,(48)Bc=−1hcosα00sinαsinαcosα.It may be observed from equation (48) that it is like the normal-projection matrix, which contains the components of vector n, normal vector to the surface Γd; this is given as(49)Pn=cosα00sinαsinαcosα.The difference between matrices (48) and (48) is the term 1/h which makes that traction continuity cannot be imposed correctly. Once a discontinuity appears in the element, the equilibrium on Γd is given as(50)∫ΩBcTσu^+σu˜︸≈0dΩ+∫ΓdtdΓ=0.Hence, the equilibrium only depends of regular stressesσu^, derived from the regular displacements part u^, and tractions t developed on the discontinuity surface.After some algebraic operations, it may be shown that tractions are expressed as follows:(51a)tn=12ldy1−y2σx,(51b)ts=12ldx2−x1σy.These results show that traction continuity is not satisfied since the correct tractions aretn=σx and ts=σy. According to this, it is derived that equilibrium is satisfied only when(52a)2ld=y1−y2,(52b)2ld=x2−x1.


### 6.1. Correct Definition of the Discontinuity Length

In what follows, a consistent way for a correct enforcement of equilibrium on the discontinuity surface is given by means of geometrical relationships. The formulation is developed for a constant strain triangle with the discrete approach of the symmetric embedded discontinuity model.The first part of equation (41) gives equilibrium in the element nodes, i.e., global equilibrium, whereas the second states traction equilibrium at the discontinuity surface, elemental equilibrium, as(53)∫ΩRTBcTσu^+σu˜dΩ+∫ΓdtdΓ=0,where R is the rotation matrix, from local coordinates to global coordinates, defined as(54)R=cosθ−sinθsinθcosθ.From equation (53), it is observed that equilibrium depends of tractions t, defined on the surface Γd, and internal forces of the solitary node i (Figure 10).Figure 10
Geometry parameters for constructions of the matrixBc.According to the construction of the matrixBc, which is derived from the element geometry, the correct imposition of equilibrium on the discontinuity occurs only over a specific surface, depending on the discontinuity geometry. This is a drawback of the presented formulation which makes it dependent of the geometry of the element, unlike the nonsymmetric embedded discontinuity formulation (see [22, 29, 33, 39]) in which the equilibrium is invariant with respect to the geometry. From this condition, the localization of the discontinuity is a variable in the problem for the present KOS formulation [4, 5] (see Figure 10). To show this condition, consider the initial state for σI>σt0,(55)tn=σt0=σI,ts=0,u≈0.The state defined in equation (55) leads to the following form of equation (53):(56)tnld=Ahcosθcosασx+sinασxy+sinθsinασy+cosασxy.This last equation is rewritten by means of trigonometric functions and the next shear stress relation,(57)σxy=12tan2θσx−σy,to obtain, after some algebra,(58)tnldhA1cosθ−α=12σx−σycosθ−αcosθ+α+tan2θsinθ+α+σx+σy.The initial conditiontn=σI implies that the right-hand side of equation (58) must be equal to σI. This is satisfied if and only if the first term is equal to(59)cosθ+α+tan2θsinθ+α=1+tan22θcosθ−α.To prove this, both sides of equation (58) are squared to obtain, after some algebra work and considering trigonometric identities,(60)sin2θsin2α=sin2θsin2α,proving equality of equation (59).Once equation (59) has been proved, the right-side of equation (58) is replaced by the principal stress σI and tn=σI to obtain(61)ld=Ahcosθ−α.Equation (61) defines the surface discontinuity geometry which guarantees the equilibrium of tractions defined by equation (53).It may exhibit another way to enforce tractions equilibrium across the discontinuity surface, by projecting stresses to a surface parallel to the opposite side of the solitary nodei (Figure 10), obtaining the next discontinuity geometry definition:(62)ld=Ah.Enforcement of tractions equilibrium on the surface defined by equation (62) makes the KOS formulation strongly dependent on the mesh configuration since the discontinuity must be parallel to one side of the element.


## 6.1. Correct Definition of the Discontinuity Length

In what follows, a consistent way for a correct enforcement of equilibrium on the discontinuity surface is given by means of geometrical relationships. The formulation is developed for a constant strain triangle with the discrete approach of the symmetric embedded discontinuity model.The first part of equation (41) gives equilibrium in the element nodes, i.e., global equilibrium, whereas the second states traction equilibrium at the discontinuity surface, elemental equilibrium, as(53)∫ΩRTBcTσu^+σu˜dΩ+∫ΓdtdΓ=0,where R is the rotation matrix, from local coordinates to global coordinates, defined as(54)R=cosθ−sinθsinθcosθ.From equation (53), it is observed that equilibrium depends of tractions t, defined on the surface Γd, and internal forces of the solitary node i (Figure 10).Figure 10
Geometry parameters for constructions of the matrixBc.According to the construction of the matrixBc, which is derived from the element geometry, the correct imposition of equilibrium on the discontinuity occurs only over a specific surface, depending on the discontinuity geometry. This is a drawback of the presented formulation which makes it dependent of the geometry of the element, unlike the nonsymmetric embedded discontinuity formulation (see [22, 29, 33, 39]) in which the equilibrium is invariant with respect to the geometry. From this condition, the localization of the discontinuity is a variable in the problem for the present KOS formulation [4, 5] (see Figure 10). To show this condition, consider the initial state for σI>σt0,(55)tn=σt0=σI,ts=0,u≈0.The state defined in equation (55) leads to the following form of equation (53):(56)tnld=Ahcosθcosασx+sinασxy+sinθsinασy+cosασxy.This last equation is rewritten by means of trigonometric functions and the next shear stress relation,(57)σxy=12tan2θσx−σy,to obtain, after some algebra,(58)tnldhA1cosθ−α=12σx−σycosθ−αcosθ+α+tan2θsinθ+α+σx+σy.The initial conditiontn=σI implies that the right-hand side of equation (58) must be equal to σI. This is satisfied if and only if the first term is equal to(59)cosθ+α+tan2θsinθ+α=1+tan22θcosθ−α.To prove this, both sides of equation (58) are squared to obtain, after some algebra work and considering trigonometric identities,(60)sin2θsin2α=sin2θsin2α,proving equality of equation (59).Once equation (59) has been proved, the right-side of equation (58) is replaced by the principal stress σI and tn=σI to obtain(61)ld=Ahcosθ−α.Equation (61) defines the surface discontinuity geometry which guarantees the equilibrium of tractions defined by equation (53).It may exhibit another way to enforce tractions equilibrium across the discontinuity surface, by projecting stresses to a surface parallel to the opposite side of the solitary nodei (Figure 10), obtaining the next discontinuity geometry definition:(62)ld=Ah.Enforcement of tractions equilibrium on the surface defined by equation (62) makes the KOS formulation strongly dependent on the mesh configuration since the discontinuity must be parallel to one side of the element.


## 7. Numerical Simulation

Two numerical examples are presented to show the numerical and theoretical consistency of the KOS formulation. In both examples, damage is simulated by considering only the mode-I of failure for quasi-brittle materials. The discontinuity is introduced perpendicular to the direction of the maximum principal stress.


### 7.1. Uniaxial Tension Test

In this example, the mechanical response of a specimen subjected only to a tension state is analysed. This problem was studied, with numerical and experimental models, by Van Vliet [40], to get insight into the physical mechanism underlying the size effect in quasi-brittle materials, such as concrete and rock.The geometry is shown in Figure11. Its characteristic dimension D is taken equal to 200 mm, radius r = 145 mm, depth equal to 100 mm, and eccentricity e=D/50 mm.Figure 11
Plate subjected to an eccentric load.The assumed material properties were the average of those reported for the specimen type “05C04N30”; this is given as follows: elastic modulusE=39.8GPa, Poisson ration ν=0.20, tensile strength σt0=2.57MPa, and fracture energy GF=0.1219N/mm.To guarantee a uniform load transfer along the top and bottom parts of the specimen, two steel plates were introduced in the laboratory experiments, allowing a free rotation. In the numerical model, this condition was reproduced by means of two layers of elements with the material properties of the steel employed in the experimental tests.In the laboratory experiments, the specimen was fully instrumented with LVDTs (linear variable displacement transducers) to obtain the deformation of its middle cross section. For the present numerical model, only the deformation along the vertical axis is considered. Figure11 shows the reference points chosen to measure the vertical deformation; initially separated a distance Ls=0.6D.A structured finite element mesh, with triangle elements, is shown in Figure12. Plane stress condition was assumed for the analysis.Figure 12
Structured mesh with triangle elements.In this example, the crack pattern was prescribed with the aim of validating the application of the KOS formulation to problems where the propagation of the crack is previously known. However, this does not invalidate its application to problems where the crack pattern is not known at initio.In Figure13, the curve deformation versus load is shown. Two numerical responses are plotted for two different softening laws. The continuous thin line corresponds to a linear softening and the continuous thick line is for an exponential softening. Experimental results, reported by Van Vliet [40], are plotted with a dotted line. It is observed that the ascendant branch is correctly reproduced by the present numerical simulation. There exists a small difference between the experimental peak load and the obtained with this study.Figure 13
Deformation versus load.This fact is due to the damage model used, which considers that mode-II occurs together with mode-I in brittle form, without energy dissipation [4, 5, 18].Regarding the rigid body motion, assumed in the present formulation of the two parts in which the solid is divided, this kinematics is observed in Figure14.Figure 14
Rigid body motion: (a) deformed mesh and (b) principal stresses.
(a)(b)The deformed mesh is shown in Figure14(a), in which it may be appreciated that element located in the middle height is completely crossed by the discontinuity, whereas the rest of elements, outside of this zone, do not undergone any damage. Similarly, in Figure 14(b), maximum principal stresses σI are plotted.The proposed embedded discontinuity formulation, in its KOS approximation, guarantees rigid body motion of the two parts in which the solid is divided by the crack.It is important to point out that, as expected from the numerical formulation and analysing equation (43), the energy dissipation only takes place on elements crossed by the discontinuity, whereas the rest unload elastically without energy dissipation.


### 7.2. Single-Notched Four-Point Beam

The four-point shear test, simulated by many authors, by means of continuum models [36, 41, 42], the embedded discontinuity formulation [5, 43–45], and meshless methods [46, 47], is addressed in this section.Since the pioneer experimental work of Arrea and Ingraffea [47], this type of experimental geometry has been used to study the damage process under mixed mode failure (see Schalangen [48, 49]).For quasi-brittle materials, like concrete, the fracture is locally driven by the maximum tensile stress [41]. This four-point shear problem undergoes a tension-compression state which gives rise to a curved propagation of the discontinuity, starting from the notch tip and following an inclined direction.Figure15 illustrates the loading system applied to the beam, which is asymmetric, and leads to a cracking failure dominated by mode-I. The material properties as well as the geometry of the problem were taken from the work of Schlangen [49].Figure 15
Single-notched four-point beam.The geometry is shown in Figure15 with its dimensions given in millimetres; the depth is 100 mm. At the top of the beam, and symmetrically with respect to a vertical axis, there is a notch of 5 mm width and 20 mm depth. All relevant dimensions and supports, at top and bottom, of the beam are shown in the same figure.With respect to the material properties, these were taken for one series of the laboratory experiments [49], for beam type SEN. The material data are elastic modulus E=30GPa, Poisson ratio ν=0.20, and fracture energy GF=0.10N/mm. Two tensile strength were considered, σt0=2.8MPa and σt0=3.0MPa.To compare experimental and numerical results, for the present simulation, the sliding displacement of the crack mouth (CMSD) was measured. This CMSD is the relative vertical displacement between both sides of the notch (see Figure15).As it may be seen from Figure15, the supports are applied in the top of the beam. The left one has restricted just the vertical displacement, whereas that located at the right restricts vertical and horizontal displacements.The mesh used for the discretization of the domain consists of triangle elements. Unlike the uniaxial tension test example, for this case, the mesh is nonstructures as it is observed in Figure16.Figure 16
Nonstructured mesh for the beam geometry.In zones where loads are applied, additional elements were considered to avoid a stress concentration and obtain a smooth distribution.The analysis was carried out considering a plane stress state, and a discontinuity is introduced when the principal stressσI reaches the tensile stress of the material, according to the Rankine criterion. The effect of mode-II effect was neglected.Unlike the uniaxial tension example, where crack pattern was prescribed, in this problem the crack propagation is nor a priori enforced, and its path is determined by the direction of the maximum principal stress. The crack is introduced perpendicular to theσI direction. This characteristic of introducing the crack without knowing its trajectory makes the embedded discontinuity formulation very promising in problems where the crack propagation across a continuum body is needed. Furthermore, the remeshing process associated with crack propagation, when standard finite elements are used, is eliminated.The numerical response, load versus CMSD, is shown in Figure17. The curve corresponding to the experimental test is plotted in continuous thick line. Numerical results, using an exponential softening law, are plotted for a tensile strength of σt0=2.8MPa with continuous thin line, and the dotted line corresponds to σt0=3.0MPa.Figure 17
Load versus CMSD for the four-point beam.In this figure, it is possible to appreciate that numerical results derived from the presented formulation reproduce with a great accuracy the ascendant branch of the experimental curve.Regarding the experimental results, it is observed that the numerical response coincides only until aboutP=10kN, likely when the ascendant branch ends, and from this point the response is different. This situation may be attributed in lower grade to the used material properties in the numerical simulation and second to the important influence of the modelling of the supports used in the experiments. In the simulation, different support systems were considered. In each case, the response was strongly influenced by the support system.Regarding the unloading branch, it is observed that the numerical and experimental curves are almost equal until a CMSD of0.1mm where there is a rapid descend in the numerical response attributable to the lack in the computed crack propagation, maximum principal stress criterion, when an isotropic damage model is used. This problem arises when the crack reaches the last elements located in the bottom part of the beam, leading to a hydrostatic stress state. Despite this, it is important to point out that the numerical simulation reproduces with a good accuracy the global response of the solid, especially the ductile behaviour of the material.In Figure18, the deformed mesh is shown. In this figure, it may be observed that the crack pattern is curved where the inelastic response is located, while the bulk unloads elastically. A zoom to the centre of the beam where the crack propagates is shown in Figure 19. The crack pattern which is not a priori enforced is plotted with continuous thick line. Since all the intraelement cracks are forced to pass through the centroid of the elements, the experimental crack pattern cannot be accurately reproduced. To improve this, it is necessary to enforce the crack path continuity across the elements together with a finer mesh.Figure 18
Deformed mesh for the four-point beam.Figure 19
Crack pattern across elements.Finally, in Figure20, the maximum principal stresses σI are plotted. In this figure, it is observed that the maximum values of the stresses occur in elements located in the crack pattern while those in the bulk undergo a linear elastic unloading process, giving a rigid body motion of the two parts in which the continuum is divided by the macrocrack.Figure 20
Principal stressσI, with elastic unloading of the bulk.


## 7.1. Uniaxial Tension Test

In this example, the mechanical response of a specimen subjected only to a tension state is analysed. This problem was studied, with numerical and experimental models, by Van Vliet [40], to get insight into the physical mechanism underlying the size effect in quasi-brittle materials, such as concrete and rock.The geometry is shown in Figure11. Its characteristic dimension D is taken equal to 200 mm, radius r = 145 mm, depth equal to 100 mm, and eccentricity e=D/50 mm.Figure 11
Plate subjected to an eccentric load.The assumed material properties were the average of those reported for the specimen type “05C04N30”; this is given as follows: elastic modulusE=39.8GPa, Poisson ration ν=0.20, tensile strength σt0=2.57MPa, and fracture energy GF=0.1219N/mm.To guarantee a uniform load transfer along the top and bottom parts of the specimen, two steel plates were introduced in the laboratory experiments, allowing a free rotation. In the numerical model, this condition was reproduced by means of two layers of elements with the material properties of the steel employed in the experimental tests.In the laboratory experiments, the specimen was fully instrumented with LVDTs (linear variable displacement transducers) to obtain the deformation of its middle cross section. For the present numerical model, only the deformation along the vertical axis is considered. Figure11 shows the reference points chosen to measure the vertical deformation; initially separated a distance Ls=0.6D.A structured finite element mesh, with triangle elements, is shown in Figure12. Plane stress condition was assumed for the analysis.Figure 12
Structured mesh with triangle elements.In this example, the crack pattern was prescribed with the aim of validating the application of the KOS formulation to problems where the propagation of the crack is previously known. However, this does not invalidate its application to problems where the crack pattern is not known at initio.In Figure13, the curve deformation versus load is shown. Two numerical responses are plotted for two different softening laws. The continuous thin line corresponds to a linear softening and the continuous thick line is for an exponential softening. Experimental results, reported by Van Vliet [40], are plotted with a dotted line. It is observed that the ascendant branch is correctly reproduced by the present numerical simulation. There exists a small difference between the experimental peak load and the obtained with this study.Figure 13
Deformation versus load.This fact is due to the damage model used, which considers that mode-II occurs together with mode-I in brittle form, without energy dissipation [4, 5, 18].Regarding the rigid body motion, assumed in the present formulation of the two parts in which the solid is divided, this kinematics is observed in Figure14.Figure 14
Rigid body motion: (a) deformed mesh and (b) principal stresses.
(a)(b)The deformed mesh is shown in Figure14(a), in which it may be appreciated that element located in the middle height is completely crossed by the discontinuity, whereas the rest of elements, outside of this zone, do not undergone any damage. Similarly, in Figure 14(b), maximum principal stresses σI are plotted.The proposed embedded discontinuity formulation, in its KOS approximation, guarantees rigid body motion of the two parts in which the solid is divided by the crack.It is important to point out that, as expected from the numerical formulation and analysing equation (43), the energy dissipation only takes place on elements crossed by the discontinuity, whereas the rest unload elastically without energy dissipation.


## 7.2. Single-Notched Four-Point Beam

The four-point shear test, simulated by many authors, by means of continuum models [36, 41, 42], the embedded discontinuity formulation [5, 43–45], and meshless methods [46, 47], is addressed in this section.Since the pioneer experimental work of Arrea and Ingraffea [47], this type of experimental geometry has been used to study the damage process under mixed mode failure (see Schalangen [48, 49]).For quasi-brittle materials, like concrete, the fracture is locally driven by the maximum tensile stress [41]. This four-point shear problem undergoes a tension-compression state which gives rise to a curved propagation of the discontinuity, starting from the notch tip and following an inclined direction.Figure15 illustrates the loading system applied to the beam, which is asymmetric, and leads to a cracking failure dominated by mode-I. The material properties as well as the geometry of the problem were taken from the work of Schlangen [49].Figure 15
Single-notched four-point beam.The geometry is shown in Figure15 with its dimensions given in millimetres; the depth is 100 mm. At the top of the beam, and symmetrically with respect to a vertical axis, there is a notch of 5 mm width and 20 mm depth. All relevant dimensions and supports, at top and bottom, of the beam are shown in the same figure.With respect to the material properties, these were taken for one series of the laboratory experiments [49], for beam type SEN. The material data are elastic modulus E=30GPa, Poisson ratio ν=0.20, and fracture energy GF=0.10N/mm. Two tensile strength were considered, σt0=2.8MPa and σt0=3.0MPa.To compare experimental and numerical results, for the present simulation, the sliding displacement of the crack mouth (CMSD) was measured. This CMSD is the relative vertical displacement between both sides of the notch (see Figure15).As it may be seen from Figure15, the supports are applied in the top of the beam. The left one has restricted just the vertical displacement, whereas that located at the right restricts vertical and horizontal displacements.The mesh used for the discretization of the domain consists of triangle elements. Unlike the uniaxial tension test example, for this case, the mesh is nonstructures as it is observed in Figure16.Figure 16
Nonstructured mesh for the beam geometry.In zones where loads are applied, additional elements were considered to avoid a stress concentration and obtain a smooth distribution.The analysis was carried out considering a plane stress state, and a discontinuity is introduced when the principal stressσI reaches the tensile stress of the material, according to the Rankine criterion. The effect of mode-II effect was neglected.Unlike the uniaxial tension example, where crack pattern was prescribed, in this problem the crack propagation is nor a priori enforced, and its path is determined by the direction of the maximum principal stress. The crack is introduced perpendicular to theσI direction. This characteristic of introducing the crack without knowing its trajectory makes the embedded discontinuity formulation very promising in problems where the crack propagation across a continuum body is needed. Furthermore, the remeshing process associated with crack propagation, when standard finite elements are used, is eliminated.The numerical response, load versus CMSD, is shown in Figure17. The curve corresponding to the experimental test is plotted in continuous thick line. Numerical results, using an exponential softening law, are plotted for a tensile strength of σt0=2.8MPa with continuous thin line, and the dotted line corresponds to σt0=3.0MPa.Figure 17
Load versus CMSD for the four-point beam.In this figure, it is possible to appreciate that numerical results derived from the presented formulation reproduce with a great accuracy the ascendant branch of the experimental curve.Regarding the experimental results, it is observed that the numerical response coincides only until aboutP=10kN, likely when the ascendant branch ends, and from this point the response is different. This situation may be attributed in lower grade to the used material properties in the numerical simulation and second to the important influence of the modelling of the supports used in the experiments. In the simulation, different support systems were considered. In each case, the response was strongly influenced by the support system.Regarding the unloading branch, it is observed that the numerical and experimental curves are almost equal until a CMSD of0.1mm where there is a rapid descend in the numerical response attributable to the lack in the computed crack propagation, maximum principal stress criterion, when an isotropic damage model is used. This problem arises when the crack reaches the last elements located in the bottom part of the beam, leading to a hydrostatic stress state. Despite this, it is important to point out that the numerical simulation reproduces with a good accuracy the global response of the solid, especially the ductile behaviour of the material.In Figure18, the deformed mesh is shown. In this figure, it may be observed that the crack pattern is curved where the inelastic response is located, while the bulk unloads elastically. A zoom to the centre of the beam where the crack propagates is shown in Figure 19. The crack pattern which is not a priori enforced is plotted with continuous thick line. Since all the intraelement cracks are forced to pass through the centroid of the elements, the experimental crack pattern cannot be accurately reproduced. To improve this, it is necessary to enforce the crack path continuity across the elements together with a finer mesh.Figure 18
Deformed mesh for the four-point beam.Figure 19
Crack pattern across elements.Finally, in Figure20, the maximum principal stresses σI are plotted. In this figure, it is observed that the maximum values of the stresses occur in elements located in the crack pattern while those in the bulk undergo a linear elastic unloading process, giving a rigid body motion of the two parts in which the continuum is divided by the macrocrack.Figure 20
Principal stressσI, with elastic unloading of the bulk.


## 8. Conclusions

In this paper, the kinematical optimal symmetric approach of the embedded discontinuity formulation is presented. Its variational and numerical formulations are developed in a consistent way within the finite element method framework. From this, it is shown that the field equation of the boundary value problem is derived in a natural way, and the equilibrium at the discontinuity surface is satisfied. However, when the FEM is used to approximate the solution, traction equilibrium depends on the discontinuity geometry. To guarantee the equilibrium in the numerical model, it is shown that the discontinuity length is a variable to consider in the formulation; otherwise, the equilibrium is not satisfied.The validity of the presented formulation is evaluated by means of two numerical examples. The first one is a simple tension test and the other is the classical four-point beam with a single notched at its upper face. In both examples, a model-I of failure was considered.In general, it is shown that the symmetric embedded discontinuity formulation can give robust results when the traction equilibrium across the discontinuity is evaluated in a specific surface, according to the new introduced length variable. In order to improve the numerical results, constitutive models, which consider a mixed mode of failure I-II, could be considered in the numerical implementation, guaranteeing a consistent energy dissipation.



---

*Source: 7125451-2020-11-30.xml*