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

This paper investigates the variational ﬁnite 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 ﬁnite element method to simulate the damage evolution processes in quasi-brittle materials, without modifying the mesh topology when cohesive cracks propagate.


Introduction
e 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][2][3], within a continuum theory, and the modelling of cohesive cracks in a discrete form [4][5][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][7][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 (Figure 1), and a stress-displacement relationship is used [9].
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 . e work energy G F needed to create a fully developed crack of unit area is as follows [1,10]: 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 (Figure 2) has been developed. is approach emerged as a mean to deal with general strain localization phenomena, of which the cohesive crack is a case (see [5,[11][12][13][14][15][16][17]).
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]. e 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 (Figure 3(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)). e 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].

Embedded Discontinuities
e 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. e 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 (Figure 2). e idea is shared between methods that enrich the nodal variables (XFEM) and methods that use elements with embedded discontinuities (EED). e 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][12][13][14][15][16][20][21][22].
ere 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.

Kinematics of the Discontinuity
To establish the kinematics of the discontinuity interface, consider a body with domain Ω and boundary Γ (Figure 4), 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 Ω + . e strong discontinuity approach, addressed in this paper, can be regarded as a limit case of the weak discontinuity when the width of the localizations zone b d 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.

Displacement Jump.
e displacement field u for the body crossed by a discontinuity, shown in Figure 4, is decomposed into a continuous and discontinuous part [16,22,26].
where u and ⟦u⟧ are smooth and continuous functions defined on Ω, with Cartesian coordinates x, and H Γ d is the Heaviside function centred at the discontinuity surface defined as e 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 Ω − . e corresponding infinitesimal strain field can be derived by taking the symmetric gradient of equation (2), and this is given as where δ Γ d is the Dirac-delta distribution centred at the discontinuity Γ d , satisfying e strain field, given by equation (4), is decomposed into a bounded part ε plus an unbounded part, irregular part, given by δ Γ d (⟦u⟧ ⊗ n) s . e unbounded strains may be considered in two different ways: (1) regularizing the Diracdelta by means of a bandwidth parameter close to zero b d ≈ 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]).

Kinematic Interpolation.
e 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.
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. eir mathematical expressions are as follows: for Ω − and Ω Ω - Figure 4: Continuum solid crossed by a discontinuity Γ d .

Advances in Civil Engineering
for Ω + . e purpose of matrix H Γ d � (0 1) T , stated in equations (6) and (7), is derived from Figure 5. For every point in Ω − or Ω + , the strain field is defined by where B is the usual compatibility matrix, relating displacements and strains [30][31][32], and the matrix H Γ d is defined in equation (3). For a two-dimensional case, these concepts can be properly generalized by introducing the vector ⟦u⟧ containing the components of jump displacements associated with the discontinuity Γ d : For this case, the matrix H Γ 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 where n is the number of nodes in the finite element, and each of the submatrices H i Γ d depends on the position of node i relative to the discontinuity given as where I is an identity matrix with dimensions (n st × n st ) and n st 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 and taking N c as equation (12) is written in the following final form: e strain kinematic is derived, in a finite element framework, by means of the displacement field given by equation (14) as where the matrix B c is derived by applying the strain-displacement operator to equation (13), and this is given as 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 ( Figure 6), equation (16) is rewritten in the form where φ(x) contains the shape functions corresponding to nodes in Ω + (see [27,29,33]) defined as where n + denotes the number of nodes belonging to the subdomain Ω + .

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][35][36]). Microcracking and plastic flow around a macroscopic crack tip are modelled as equivalent tractions forces on crack faces ( Figure 7). A discontinuity is introduced when the maximum principal stress σ I exceeds the tensile strength of the material. is mode-I criterion is commonly used for quasibrittle 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 with ⟦u⟧ already defined in Section 3.2, κ is an internal scalar variable which can be considered as an equivalent jump (κ � ⟦u⟧ eq ), and ψ 0 is the linear elastic free energy given by where T el 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]) : 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, Use of standard thermodynamics arguments, in equation (21), leads to the following: (2) In addition, the dissipation equation is reduced to 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.

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: where t el is the elastic traction vector and Substitution of equation (26) into equation (25) leads to If unloading takes place, the rate of damage is zero, i.e., ω _ � 0, and equation (27) reduces to

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). e functions which describe this process is defined as where ⟦u⟧ eq is an equivalent jump obtained from the displacement jump vector ⟦u⟧. e symbols 〈·〉 are the McAuley brackets denoting that only the positive part of ⟦u⟧ eq is considered. e internal scalar variable κ is a history parameter equal to the largest value of 〈⟦u⟧ eq 〉, and this is given as

Loading/Unloading
Conditions. If f < 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 . Both cases can be conveniently expressed by the usual Kuhn-Tucker conditions [38], together with the consistency condition (32)

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

Advances in Civil Engineering
between both subdomains is defined by the jump ⟦u⟧ in the displacement field u. e energy functional for this solid is given by 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 and with t being the tractions vector derived from a cohesive traction-jump constitutive relation, t � T⟦u⟧. Furthermore, for a discontinuity with a linear elastic behaviour, ψ(⟦u⟧) is written as Applying the first variation to the energy functional defined in equation (33), the next Euler-Lagrange equations are recovered: Equations (36a)-(36b) represent the internal equilibrium and natural boundary conditions of the body, complemented by the following equilibrium conditions across the discontinuity: and the kinematic condition Equation (38) states the essential boundary conditions of the variational problem.

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 Section 3. is is given as where u and u are the regular and not regular displacements, respectively. e matrix N contains the standard shape functions. Shape functions associated with nodes belonging to Ω + are contained in matrix N c .
As for standard finite element method, matrices B and B c , given in equation (40), contain derivatives of the corresponding shape functions in matrices N and N c .
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, and external forces vector Here, the coefficient matrix is symmetric and unlike other discontinuity formulations [26][27][28]33], there are nonzeros along its leading diagonal. e inelastic response is related to the term Γ t N T tdΓ of the submatrix K uu .

Tractions Equilibrium
e interpolation of the displacement jump field by the finite element method of the (KOS) formulation considers the shape function N c . is kind of interpolation is a proper way to consider that the localization band width b d is different from zero; but, when this band collapses into a line, i.e., discrete discontinuity, the matrix N c has the meaning of the identity matrix which leads to consider tractions at the discontinuity surface in the local framework. e case, when b d � 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. e discontinuity crossing the element is parallel to side 1 − 2 with length l 12 ; n is the normal to the discontinuity, and h is the height of the triangle. e 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. Advances in Civil Engineering e area of the element is A � l 12 h/2; hence ,2A � l 12 h, and from the geometry of the element and the angle β, the following relations are derived: Furthermore, for the element shown in Figure 9, the matrix B c is defined as Rewritten equation (47) by means of relationships derived from the geometry of the element and using the angle α instead of β, 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 e 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 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: ese results show that traction continuity is not satisfied since the correct tractions are t n � σ x and t s � σ y . According to this, it is derived that equilibrium is satisfied only when

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. e formulation is developed for a constant strain triangle with the discrete approach of the symmetric embedded discontinuity model. e 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 where R is the rotation matrix, from local coordinates to global coordinates, defined as 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).
According to the construction of the matrix B c , 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.
is 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 ,  Figure 9: Projection of the discontinuity to the opposite side to node 3. (56)

Advances in Civil Engineering
is last equation is rewritten by means of trigonometric functions and the next shear stress relation, to obtain, after some algebra, (58) e initial condition t n � σ I implies that the right-hand side of equation (58) must be equal to σ I . is is satisfied if and only if the first term is equal to cos(θ + α) + tan 2 θ sin(θ + α) � ����������������� � 1 + tan 2 2θ cos(θ − α).

(59)
To prove this, both sides of equation (58) are squared to obtain, after some algebra work and considering trigonometric identities, sin 2 θ sin 2 α � sin 2 θ sin 2 α, 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 t n � σ I to obtain 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 node i (Figure 10), obtaining the next discontinuity geometry definition: 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.

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. e discontinuity is introduced perpendicular to the direction of the maximum principal stress.

Uniaxial Tension Test.
In this example, the mechanical response of a specimen subjected only to a tension state is analysed.
is problem was studied, with numerical and experimental models, by Van Vliet [40], to get insight into the physical mechanism underlying the size effect in quasibrittle materials, such as concrete and rock. e geometry is shown in Figure 11. 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. e assumed material properties were the average of those reported for the specimen type "05C04N30"; this is given as follows: elastic modulus E � 39.8 GPa, Poisson ration ] � 0.20, tensile strength σ t0 � 2.57 MPa, and fracture energy G F � 0.1219 N/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. Figure 11 shows the reference points chosen to measure the vertical deformation; initially separated a distance L s � 0.6 D.
A structured finite element mesh, with triangle elements, is shown in Figure 12. Plane stress condition was assumed for the analysis.

Advances in Civil Engineering
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 Figure 13, the curve deformation versus load is shown. Two numerical responses are plotted for two different softening laws. e 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. ere exists a small difference between the experimental peak load and the obtained with this study. is 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 Figure 14.
e deformed mesh is shown in Figure 14(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. e 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.
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]. is fourpoint 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. Figure 15 illustrates the loading system applied to the beam, which is asymmetric, and leads to a cracking failure dominated by mode-I. e material properties as well as the geometry of the problem were taken from the work of Schlangen [49]. e geometry is shown in Figure 15 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.

Advances in Civil Engineering
With respect to the material properties, these were taken for one series of the laboratory experiments [49], for beam type SEN. e material data are elastic modulus E � 30 GPa, Poisson ratio ] � 0.20, and fracture energy G F � 0.10 N/mm. Two tensile strength were considered, σ t0 � 2.8 MPa and σ t0 � 3.0 MPa.
To compare experimental and numerical results, for the present simulation, the sliding displacement of the crack mouth (CMSD) was measured. is CMSD is the relative vertical displacement between both sides of the notch (see Figure 15).
As it may be seen from Figure 15, the supports are applied in the top of the beam. e left one has restricted just the vertical displacement, whereas that located at the right restricts vertical and horizontal displacements. e 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 Figure 16.
In zones where loads are applied, additional elements were considered to avoid a stress concentration and obtain a smooth distribution.
e 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. e 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. e crack is introduced perpendicular to the σ I direction. is 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. e numerical response, load versus CMSD, is shown in Figure 17. e 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.8 MPa with continuous thin line, and the dotted line corresponds to σ t0 � 3.0 MPa.
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 about P � 10 kN, likely when the ascendant branch ends, and from this point the response is different. is 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 of 0.1 mm 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. is 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 Figure 18, 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. e 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.
Finally, in Figure 20, 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.

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. e validity of the presented formulation is evaluated by means of two numerical examples. e 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.

Data Availability
e experimental data used to support the study have been obtained from reports of other researchers. Specifically, (a) uniaxial tension test problem: M. R. A. V. Vliet, "Size effect in tensile fracture of concrete and rock," Ph.D. esis, Delft University of Technology, Delf, the Netherlands, 2000, and (b) single-notched four-point beam: E. Schlangen, "Experimental and numerical analysis of fracture process in concrete," Ph. D. esis, Eindhoven, the Netherlands, 1993, and data are cited at relevant places as references.

Conflicts of Interest
e authors declare that they have no conflicts of interest.