Simplified Qualitative Discrete Numerical Model to Determine Cracking Pattern in Brittle Materials by Means of Finite Element Method

This paper presents the formulation, implementation, and validation of a simplified qualitative model to determine the crack path of solids considering static loads, infinitesimal strain, and plane stress condition.Thismodel is based on finite element method with a special meshing technique, where nonlinear link elements are included between the faces of the linear triangular elements. The stiffness loss of some link elements represents the crack opening. Three experimental tests of bending beams are simulated, where the cracking pattern calculatedwith the proposed numericalmodel is similar to experimental result.The advantages of the proposed model compared to discrete crack approaches with interface elements can be the implementation simplicity, the numerical stability, and the very low computational cost. The simulation with greater values of the initial stiffness of the link elements does not affect the discontinuity path and the stability of the numerical solution. The exploded mesh procedure presented in this model avoids a complex nonlinear analysis and regenerative or adaptive meshes.


Introduction
Fracture mechanics studies the cracking process of a solid subjected to progressive external load.Particularly, computational fracture mechanics allows representing the formation and propagation of cracks in solids of general geometry by means of numerical models [1].The results obtained with these models enable studying new materials and reducing cost of experimental tests [2].
The fracture process in a solid can be represented according to the description of displacement and strain field, the material constitutive model, and the numerical approximation technique in the finite element method.Likewise, the models could be divided by the following three types: the models with propagating cohesive discontinuities, the softening continuum models with partial regularization, and the regularized softening continua models [3].
Models with propagating cohesive discontinuities assume that the fracture process zone preserves a linear elastic behavior during its formation.Cohesive forces are defined between the faces of the discontinuity.These forces disappear when the gap between faces reaches a certain distance.
Consequently, the kinematics singularity vanishes and the crack opening increases in a discrete form [4]. Some authors indicate that each change of the discontinuity path needs remeshing process [5][6][7][8].Other authors propose to enhance the trial function of the finite element in order to represent a jump of the displacement field [9][10][11].The concept of partition of unity enhances the shape function of standard finite elements and ensures a uniqueness result [3].Particularly, extended finite element method (X-FEM) is based on previous concept, choosing a function that represents the discontinuity of the displacement field on the crack faces [12][13][14].

Mathematical Problems in Engineering
Softening continuum models with partial regularization represent the fracture process zone by means of the strain localization on a finite band.Although the displacement field is defined as a continua form, there are weak discontinuities in the boundaries of the band [3].The softening in fracture zone is independent of the size finite element and is also associated with the fracture energy per area unit [15].The band of fracture zone is embedded into finite element in order to ensure objectivity with respect to the orientation of the mesh [16].
Regularized softening continua models preserve continuity of the displacement and strain fields.The fracture process zone is represented by a material band, in which the softening strain increases from the band boundary until its center [3].These models have generally high computational cost; however, they have some advantage with respect to the models of partial regulation.The influence of size element is substantially reduced with fine meshes.These models show good results with adaptive meshing technique [17].
Other classification establishes two types of numerical models in order to represent cracking in brittle materials: the smeared crack approach and the discrete crack approach [18].
The smeared crack approach considers that an infinite amount of parallel cracks, each with very small opening, are assigned to the finite elements.The constitutive material model in the element is modified, such that the tangent stiffness and the stress in normal direction of the crack are reduced while the strain increases [19].
The discrete crack approach indicates that the fracture process zone is concentrated at a surface characterized by a relationship between the traction versus displacement jump, which describes the cohesion loss of the material between the crack faces after fulfilling the failure criteria.This approach has been developed on different models.Originally, cohesive forces associated with the fracture energy between the faces of a crack are appended.The initial models as the cohesive crack model prescribe the location of crack [20,21].Subsequently, the fictitious crack model predicts the formation and propagation of the crack in anywhere of the solid [22,23].The location of the crack with respect to finite element mesh establishes two different numerical techniques.The first one states that the crack is traced on the finite element sides and the cohesive forces are obtained by means of nonlinear mixed boundary conditions [18] or through interface elements with vanished or zero thickness that connect the nodes to both sides of the discontinuity [24][25][26][27][28][29][30].The second technique states that the crack crosses the finite element [31,32].
Particularly in two-dimensional mechanical problems using discrete crack models with interface elements, some authors define two triangular finite elements with high aspect ratio in the interface [29,30], which depend on a tension damage constitutive relationship and the same kinematics as the continuum strong discontinuity approach tangent stiffness factor.Other authors state zero-thickness interface elements that couple pairs of duplicate nodes, whose behavior is defined by a traction-separation law [28].In other works [25][26][27], the fracture process in the interface surface is based on failure criteria of the three-parameter hyperbolic cracking surface [24], which relates the normal and tangential stress components to the corresponding relative displacements.These last approaches request additional procedures in order to preserve the numerical stability in the nonlinear analysis.
In composite materials, the analysis of crack nucleation and growth must be addressed taking into consideration the nonlinear behavior in its microstructure.Different phenomena such as void growth, microcracking, interfacial debonding, and other nonhomogeneities are closely related to failure mechanisms that produce macroscale failure; as a consequence, macroscopic fracture models may turn out to be inappropriate to estimate crack trajectories and the structural response of those kinds of specimens [33].Different authors have been developed approaches to tackle these problems via multiscale methods, which, by means a combination of microscale phenomena and mathematical approaches, allow describing the behavior of multiphase materials.In this sense, different works must be pointed [34][35][36][37].
This paper presents the formulation, implementation, and validation of a qualitative numerical model, which describes qualitatively the cracking pattern in a brittle homogeneous material, considering static loads, plane stress condition, and infinitesimal strain.This model is based on both the discrete crack approach and the finite element method, where the overlapping faces of the triangular elements are doubled and connected to zero-length link elements.Some sides of triangular elements are part of the discontinuity path, where the nonlinear behavior is represented by the link elements.The tangent stiffness of these elements tends to infinity during the linear elastic behavior of the material and is equal to zero when the failure criteria are fulfilled.This model avoids the remeshing process, provides the implementation into finite element code, and maintains a low computational cost.However, the structural response cannot be obtained because the cohesive law in the cracking zone was not considered.The advantage of the proposed model compared to discrete crack approaches with interface elements can be the implementation simplicity, the numerical stability, and the very low computational cost.

Formulation of the Numerical Model
The proposed model is based on continuum mechanics applied to solid with discontinuities.The latter are produced by the fracture process of the brittle material.Particularly, plane stress condition, infinitesimal strain, and static load are considered.

Government Equations in the Continuum.
A solid is subjected to body forces vector b at the domain Ω and surface forces t * at the contour Γ  with normal vector n, as shown in Figure 1(a).The prescribed displacements are indicated by the vector u * at the contour Γ  .The state stress in each material point x is expressed by the tensor .Equations (1a) to (1c) present the equilibrium and boundary conditions in the solid: Figure 1: Solid with discontinuity subjected to body and surface forces: (a) general sketch and (b) detailed traction vector inside of the discontinuity.
In Figure 1(b), a discontinuity surface Γ  with normal vector n  preserves continuity of the traction vector t  , of the following form: The weak form of the mechanical problem of a solid with discontinuities is obtained by means of virtual work principle, which establishes that the following [28,38]: The external virtual work  ext is produced by the body and surface forces b and t * and is expressed as follows: where u is the virtual displacement vector.The internal virtual work  int is generated by stress tensor  and the virtual strain tensor  = ∇  u at Ω ⊐ Γ  , of the following form: Linear elastic behavior at domain Ω outside of discontinuity contour Γ  , that is, Ω ⊐ Γ  , is assumed.Consequently the stress tensor is equal to  = D : , where D is the linear elastic constitutive tensor of the material.The cohesive virtual work  coh is produced by traction vector t  between the faces of the discontinuity Γ  and the virtual displacement jump vector ⟦u⟧; thus The cohesive loss between the faces of the discontinuity exhibits a nonlinear behavior, which is expressed by the relationship between the traction rate vector ṫ  and the rate vector of the displacement jump ⟦ u ⟧; thus where T is the tangent constitutive tensor of the tractiondisplacement jump model at the discontinuity.

Failure Criteria and Fracture
Process.The proposed model uses Rankine's failure criteria restricted by tensile stress states.These criteria establish that the failure of material takes place when the positive maximum principal stress   reaches the tensile strength of the material   .In plane stress condition, this stress is defined as follows: The fracture process in mode I appears with the failure of material, where the crack is normal to the direction of the positive maximum principal stress n  =   i +   j, where   = cos   ,   = sin   , and   is the angle between global axis and the direction of positive maximum principal stress.
The cohesion between crack faces is lost inside the cracking zone, while elastic unloading is shown at neighborhood of the cracking zone.

Implementation of the Numerical Model in Finite Element Method
The mechanical problem raised in the previous section is implemented in finite element method.The nonlinear analysis procedure, the meshing technique with potential discontinuities, the description of used type elements, and the evaluation of its tangent stiffness are presented as follows.

Nonlinear Analysis by Means of Finite Element Method.
A solid in plane stress condition can be represented with a mesh of  linear triangular finite elements connected by  link elements, where the discontinuity path may appear on the contours of the triangular elements.Each overlapping side of two triangular elements is separated, and its nodes are duplicated and connected with two link elements, as shown in Figures 2(a) and 2(b).Each link element connects two overlapping nodes; therefore, its length is null and its longitudinal axis is perpendicular to the side of associated triangular element, as shown in the virtual gap of Figure 2(c).The discontinuity path on a contour Γ  is not known a priori.Consequently, the contours of triangular elements Γ ℎ are assumed as likely discontinuity paths.The linear elastic behavior on the domain outside of possible discontinuity paths Ω ⊐ Γ ℎ is represented by linear triangular elements, while the nonlinear behavior on possible discontinuity paths Γ ℎ is represented by link elements.The latter shows displacement compatibility when there is no discontinuity and cohesive loss when there is a discontinuity.
In each triangular finite element , the transpose virtual displacement vector u   = a   N   , where N  is the shape functions matrix and a  is the nodal virtual displacement vector of the element.Likewise, the transpose virtual strain vector    = a   B   , where B  is the gradient matrix of the element  [39].
The external virtual work of the solid, expressed in (4), can be defined in a mesh of  triangular finite elements and using matrix notation; thus where the nodal virtual displacement of the mesh corresponds to a and the external force vector of mesh f ext is obtained from assembling of vector f ext  of  triangular elements, that is, The external force vector of a triangular finite element is defined as follows: The internal virtual work of the solid, expressed in (5), can be defined in a mesh of  triangular finite elements; thus where the internal force vector of the mesh f int is obtained from assembling of the vector f int  of the  triangular elements, that is, The internal force vector of a triangular element is defined of the following form: The cohesive virtual work in the discontinuity of the solid, expressed in (6), can be defined in a mesh of  link elements; thus where the cohesive force vector of the mesh f coh is obtained from assembling of the vector f coh  of the  link elements, that is, 9), (11), and ( 13) are substituted into (3), and the transpose nodal virtual displacement vector a  is factorized and canceled.The obtained equilibrium condition establishes that f int + f coh − f ext = 0.In the context of the nonlinear analysis with finite element method, the equilibrium condition is fulfilled when the residual force vector tends to zero, that is, Ψ → 0. The residual force vector Ψ is defined as follows: According to Newton-Raphson method, the residual force vector Ψ(a  + Δa  ) evaluated on nodal displacement a  plus the trial incremental displacement vector Δa  in the iteration  is equal to The derivate of the residual force vector, presented in (14), with respect to nodal displacement is The external force vector does not change with respect to nodal displacement in the iterative process; then  a f ext = 0.However, the derivative of the internal and cohesive force vectors with respect to nodal displacement are equal to the tangent stiffness matrices K int and K coh , respectively.
Nonlinear numerical solution method searches the trial incremental nodal displacement Δa  for which the residual force vector Ψ(a  + Δa  ) is equal to zero.In each iteration  of this procedure, an equations system is solved; thus The tangent stiffness matrices K int and K coh are obtained from assembling of the tangent stiffness matrices K int  and K coh  of each finite element, respectively; thus where  identifies the two-dimensional element from 1 to  and  indicates the link element from 1 to n.A linear elastic material is assumed in the domain of the triangular elements; therefore   = D    , where the constitutive matrix D  is constant with respect to the strain state and   = B  a  .The internal force vector of a triangular finite element is f int  = K int  a  , where the nodal displacement vector is a  and the stiffness matrix K int  is defined of the following form [39]: This matrix depends on thickness  of the element and Young's modulus  and Poisson's relation ] of material.
In contrast, each link element  exhibits nonlinear behavior, in which the cohesive force vector rate ḟ coh corresponds to where ȧ  is the vector rate of the displacement on the duplicate nodes of link element .Tangent stiffness matrix K coh  of this element is defined as follows: The longitudinal axis of link element is normal to the side of the neighboring triangular elements and it is defined by directional vector n  =   i +  j, where   = cos   and   = sin   .The angle between the longitudinal axis of the link element and -global axis is called   , as shown in Figure 2(c).This numerical model assumes that the crack path does not substantially depend on the cohesive law of the brittle material.Likewise, the single goal of this work is to predict the crack path, without representing the structural response of the solid.Consequently, this model proposes a simplification of the cohesive law, in which the tangent stiffness factor  coh  is equal to zero, when the cohesion is totally lost, and tends to infinity when the displacement compatibility is preserved.

Generation of the Exploded Finite Element
Mesh.First, a conventional mesh of linear triangular elements and the boundary conditions are generated.Next, a new mesh of linear triangular and link elements is produced with the information of the conventional mesh.The new mesh is called exploded mesh and is made only once during the simulation.Figure 3(a) shows conventional mesh where some sides of linear triangular elements are overlapped.A virtual gap separates the overlapping sides of two triangular elements and its nodes are duplicated as shown in Figure 3(b).Each node of the conventional mesh is replaced with a set of overlapping nodes in the exploded mesh; for example, node 6 of the conventional mesh is replaced with nodes 11, 12, 13, 14, and 15.Each pair of duplicated nodes is connected with a link element; for example, link element 1 connects to nodes 11 and 12.
The triangular finite elements preserve the linear elastic behavior during the whole loading process and depend on the mechanical elastic properties of the material.Particularly, in x Figure 4: Out-scale sketch of the exploded mesh around node 6 of the conventional mesh: (a) exploded mesh in node 6 of the conventional mesh, (b) directional vector of the link elements near node 6 and comparative direction, and (c) angle between the comparative direction and the directional vector of a link element.
the first loading step, the tendency to infinity of the tangent stiffness factor of all link elements is kept, which ensures the displacement compatibility between nodes.
In the implemented numerical procedure, the numbering of triangular finite elements of the exploded and conventional meshes are the same.However, the numbering of each set of overlapping nodes of the exploded mesh is associated with the node number at the same location of the conventional mesh.

Evaluation of the Tangent Stiffness of Link Elements.
The nodal stress components   ,   , and   on each set of overlapping nodes are obtained from the average among the stress components of surrounding linear triangular elements.The principal stresses  1 and  2 and their directions n 1 and n 2 are computed as indicated in some references for plane stress condition [38,39].The positive maximum principal stress   and its direction n  are selected according to (8).The angle between the -global axis and the directional vector of the positive maximum principal stress n  is called   .
If the positive maximum principal stress   reaches the tensile strength of the material   , then the tangent stiffness in one of the link elements , . . .,  associated with this set overlapping nodes is modified.If the orientation n  of the link element has the least difference with respect to comparative direction n  , then the tangent stiffness factor  coh  of ( 21) is equal to zero.Otherwise, this factor preserves the tendency to infinity, that is, Figure 4(a) shows link elements 1 to 5 associated with the overlapping nodes 11 to 15, where the positive maximum principal stress is equal to the tensile strength of the material.
Here, the tangent stiffness factor of link element 5 is zero, because its orientation has the least difference with respect to n  , as shown Figure 4

(b).
Since there is no perfect alignment between the comparative direction n  and the link orientation n • selected because it has the least difference with respect to n  , there is an implicit error  err • in computing of crack direction.This error is defined in where the angle between -global axis and the directional vector of a link element n • is called  • and the angle between -global axis and the vector of the comparative direction n  is called   , as shown in Figure 4(c).The angle of the comparative direction   is equal to the angle of the direction principal   in the first set of overlapping nodes where   =   .In the following set of overlapping nodes, the comparative direction is defined as indicated in Section 3.4.
The side of triangular element normal to vector n  of link element with zero stiffness is part of the discontinuity path.Figure 4(a) shows link elements 1 to 5 associated with the overlapping nodes 11 to 15, where the positive maximum principal stress is equal to the tensile strength of the material.Here, the tangent stiffness factor of link element 5 is zero, because its orientation has the least difference with respect to n  , as shown in Figure 4(b).

Correction of the Discontinuity Path.
The crack path is considered normal to direction of positive maximum principal stress at each material point of solid.However, the discontinuity path is traced on the sides of triangular finite elements in the numerical model.In order to reduce the difference caused by the mesh alignment, the comparative direction  ()  in the set of overlapping nodes  with  ()  =   is defined as follows: where  err(−1)

𝑙
is the error angle computed in (23), which was obtained at the previous set of overlapping nodes  − 1 with  (−1)  =   .As was mentioned above, the purpose of this correction is to avoid the implicit error resulting of misalignment of the comparative directions n  with respect to mesh topography.Figure 5 shows an example of discontinuity path with and without correction in a mesh of triangular finite element.The dash line represents the expected crack path of the solid.The discontinuity path is nearer to the crack path when the correction is applied.

Application Examples
The three experimental tests developed by other authors are simulated by means of proposed numerical model.These tests correspond to notched beams with and without holes, subjected to transversal load in one or two points.The deformed shape of the finite mesh in the last loading step exhibits a path where the relative displacement between nodes is greater.This path represents the discontinuity in the solid, which is compared with the cracking pattern of the experimental test.

Example 1: Three-Point Beam with Nonconcentric Notch.
A simply supported beam of polymethylmethacrylate (PMMA) is subjected to load  at the upper center part.Figure 6 shows the geometry and load applied on the beam, where the thickness is  = 5 mm, and the notch has  of depth and separated  from the loading line on the bottom center face of the beam.The mechanical features of the material are as follows: Young's modulus  = 3100 MPa, Poisson's relation ] = 0.4, and tensile strength   = 76 MPa.Two cases with  different values of  and  are studied.In Case I,  = 10 mm and  = 60 mm are defined and three meshes with different amount of triangular elements are simulated (Table 1).Figure 7 shows the discontinuity path obtained with each mesh and the cracking pattern from the experimental test developed by Ingraffea and Grigoriu [40].The discontinuity path of mesh 1 is close to the experimental crack.However, the numerical result shows a discontinuity path almost straight because of low numbers of finite elements.
The simulation with mesh 2 captures the curvature of the experimental crack, despite the approximation being least at the beginning of the path.Likewise, the discontinuity path obtained with mesh 3 is almost equal to the experimental result, except for the zigzag line associated with the size of the orientation of the finite elements sides.In Case II,  = 15 mm and  = 50 mm are defined and a mesh of 9636 triangular elements is simulated.The notorious relative displacement between nodes is obtained by numerical model, which is similar to the cracking pattern of the experimental test [40], as shown in Figure 8(a).

Example 2: Three-Point Beam with Notch and Three Holes.
Example 2 has the same material and external geometry of the previous example but includes three internal holes which are located between the notch and the load, as shown in Figure 9.The principal stress and strain are high in the neighborhood of each hole.This produces a curvilinear cracking pattern unique [42].The depth  and the distance  are defined in two cases.Each case is simulated with three finite element meshes with different element size as indicated in Table 2.This problem has been experimentally tested by Ingraffea and Grigoriu [40] and numerically modeled by others authors [42][43][44][45][46].
A depth  = 10 mm and the distance  = 60 mm are defined in Case I.The discontinuity path of the proposed numerical model is compared with the experimental cracking pattern and with the discontinuity path computed by other authors [43][44][45][46][47].These results are shown in Figure 10, considering three different meshes of the proposed model (Table 2).The discontinuity path of meshes 1 and 2 has shown similarity with respect to experimental cracking pattern, except at its ends.The path obtained in the finest mesh or mesh 3 is very similar to the experimental result, even close to the holes.
Other approaches as the element free Galerkin method [47] or the edge-based smoothed finite element method [43] exhibit good results.Mesh 3 of proposed model, numerical approach of Nguyen-Xuan and collaborators, and the experimental test present the same crack path, in spite of the sawtooth form of the proposed model.
Figure 11 shows the notorious relative displacement between nodes of the deformed shape of the noted mesh while the applied load is increased.This represents the evolution of the discontinuity path.A depth  = 15 mm and the distance  = 50 mm are defined in Case II.The comparison between the discontinuity path of the proposed numerical model with three meshes, the numerical model of other authors, and the crack path of experimental test [40] is shown in Figure 12.The obtained numerical result of mesh 1 indicates considerable difference, in which its end does not reach the intermediate hole.Mesh 2 captures one of two curvatures of the experimental crack.The finest mesh presents the best approximation, in which the numerical discontinuity path exhibits the double curvature and the end point in the hole of the crack pattern.Furthermore, this figure shows that the work of Nguyen-Xuan and others is slightly most precise than the result of mesh 3 of the proposed model.The path of the latter is equal to the model of Ventura and collaborators.

Example 3: Beam with Notch and Two Loading Points.
A concrete beam is subjected to the loads 0.13 and , as shown in Figure 13.The beam has a thickness of  = 156 mm and a notch of 82 mm of depth located nearby the support.The material has a Young modulus of  = 24800 MPa, Poisson's relation of ] = 0.18, and tensile strength of   = 2.8 MPa.The geometry and the boundary conditions produce an important zone of constant shear stress [48].This test was developed by Arrea and Ingraffea [41] and was simulated with the proposed model and a mesh of 16284 triangular finite elements.
The notorious relative displacement between nodes of the deformed shape of the numerical model is shown in Figure 14(b), whereas, in Figure 14(a), both show the numerical discontinuity path and the experimental crack obtained by Arrea and Ingraffea [41].There is a qualitative similarity between the two paths, apart from inward curve of the first part of the path.
The works of Rots with smeared crack models [49] show a good approximation of this test.Particularly, the discontinuity path presents a better result near the notch with respect to proposed model, but this path has curvature less that the real crack.displacement compatibility in the overlapping nodes before the formation of the crack.
Convergence criteria of residual forces were used for these nonlinear simulations with a tolerance of 1 × 10 −3 .The maximum number of iterations per loading step was equal to 2, as shown in Table 3.

Conclusions
The main conclusions of this work are as follows.
The formulation, implementation, and validation of a discrete numerical model which predicts the cracking pattern of a solid, considering infinitesimal strain, static loads, and plane stress condition, are presented in this paper.A nonlinear analysis with the finite element method is implemented.This model assumes that the crack is produced between the sides of the linear triangular elements, connected by the link elements.The latter element type defines the activation of the fracture process in accordance with the magnitude and direction of the positive maximum principal stress.A link element with high stiffness represents the noncracking zone and a link element with null stiffness establishes the crack path.
The exploded mesh is generated from the conventional mesh one time in the simulation.This procedure separates the sides of the triangular finite elements and joins the adjacent nodes with the link elements, without modifying the original geometry.
Linear elastic behavior of the triangular elements and nonlinear behavior of the link elements in the exploded mesh are taken into consideration.Initially, the tangent stiffness of all link elements tends to infinity.In the following loading steps of the nonlinear analysis, the tangent stiffness of several link elements is reduced to zero in order to represent the crack path.The simplicity and stability of this procedure sustain a good approximation of the crack path, reduce the computational cost, and allow implementing it on a standard finite element code with few modifications.However, the structural response cannot be obtained because the cohesive law in the cracking zone was not considered.
A simple procedure of correction of the discontinuity path caused by the mesh alignment is presented.Here, the orientation of the discontinuity path is defined by the link elements with stiffness null.A link element loses the stiffness when its orientation is nearest to the comparative direction and the tensile strength is reached.The comparative direction is defined by the direction of the positive maximum principal stress.The difference between the link element orientation and the comparative direction is used in order to correct the comparative direction in the next set of overlapping nodes.
This procedure presents successful results of the discontinuity path of several examples of three-point beam with nonconcentric notch, subjected to the simultaneous action of normal and shear stresses.The topology of the numerical discontinuity path depends on the size and orientation of the finite elements.However, the finest meshes of Examples 1 and 2 show that the difference with respect to experimental test is negligible.Example 3 shows lesser accuracy because of the strong curvature of the real crack.Generally, the numerical model exhibits a good approximation of the cracking pattern, using meshes of about 400 triangular elements at height of the beam.
The discontinuity paths of Examples 2 and 3 obtained with the finest mesh of the proposed model are similar to the ones computed by other authors using different numerical models.The advantages of the proposed model compared to discrete crack approaches with interface elements can be the implementation simplicity, the numerical stability, and the very low computational cost.In the proposed model, once the conventional mesh is modified, a nonlinear problem is solved with elastic triangle elements and inelastic link elements which are normally included on finite element software.The simulation with greater values of the initial stiffness of the link elements does not affect the discontinuity path and the stability of the numerical solution.This value is limited by the precision of the computer.The exploded mesh procedure presented in this model avoids the regeneration or adaptation of the mesh during the formation of crack and consequently the computational cost is low.
The tests of structural members of brittle material reinforced with ductile material, for example, the reinforced concrete beams, exhibit several cracks.These tests could be simulated using the proposed numerical model.In future works, no models and others failure criteria of the fracture process could be implemented in order to describe the structural response.

Figure 2 :
Figure 2: Connecting an overlapping side of two triangular elements by means of two link elements: (a) sketch of two triangular elements in real scale, (b) sketch of two triangular elements and two link elements out of scale, and (c) detailed link element and its orientation with regard to neighboring triangular elements in a virtual gap.

Figure 3 :
Figure 3: Finite element meshes: (a) conventional mesh and (b) exploded mesh out of scale.

Figure 5 :
Figure 5: Detailed discontinuity path in a mesh of triangular finite elements: (a) with correction and (b) without correction.

Figure 6 :
Figure 6: Sketch of the three points beam with nonconcentric notch.The measurements are given in millimeters.

Figure 7 :
Figure 7: Comparison between the discontinuity path in three finite element meshes and the crack path of the experimental test in Case I of Example 1.

Figure 8 :
Figure 8: Numerical and experimental crack path in Case II of Example 1: (a) comparison between the discontinuity path in numerical model and the crack path of experimental test [40] and (b) relative displacement between nodes of the deformed shape of the numerical model.

Figure 9 :
Figure 9: Sketch of the three-point beam with notch and three holes.The measurements are given in millimeters.

4. 4 .
Sensibility of the Discontinuity Path to the Initial Stiffness of Link Elements.The initial stiffness factor of link element  coh,0  tends to infinity in order to ensure the displacement compatibility between overlapping nodes.This tendency is to be expressed by means of a bounded value.This model establishes that  coh,0  is much higher than the maximum coefficient of the stiffness matrix of an equilateral triangular element with linear elastic behavior  int (max) = /( √ 3(1 − ] 2 )).A sensibility analysis of the three application examples was developed.Particularly, some simulations with different values of the initial stiffness factor of the link elements and fine meshes were carried out.Table3indicates the maximum number of iterations per loading step in column  and whether (Y) or not (N) the discontinuity path is obtained in column , for four values of  coh,0  .The simulations with initial stiffness factor between 1×10 3 and 1×10 5 times  int (max) exhibit the same discontinuity path, keeping the numerical stability of the solution.The numerical precision of the computer is overcome only when the initial stiffness factor of link element is 1 × 10 6 times  int (max) , Examples 1 and 2. Low values of initial stiffness of the link elements were not assigned to the simulations because these do not represent the Num path, proposed model mesh 1 Num path, Nguyen-Xuan et al.Num path, Ventura et al.Num path, proposed model mesh 2 Num path, Nguyen-Xuan et al.Num path, Ventura et al.Num path, proposed model mesh 3 Num path, Nguyen-Xuan et al.Num path, Ventura et al.

Figure 10 :
Figure 10: Comparison between the discontinuity path in numerical model and the crack path of experimental test in Case I of Example 2.

Figure 11 :Example 2 ,Example 2 ,Example 2 ,
Figure 11: Evolution of the discontinuity path in the numerical simulation of the finest mesh in Case I of Example 2.

Figure 12 :Figure 13 :
Figure 12: Comparison between the discontinuity path in numerical model and the crack path of experimental test in Case II of Example 2.

Figure 14 :
Figure 14: Numerical and experimental crack path in Example 3: (a) comparison between the discontinuity path in numerical model and the crack path of experimental test [41] and (b) relative displacement between nodes of the deformed shape of the numerical model.

Table 1 :
Meshes of Case I of Example 1.

Table 2 :
Finite element meshes of Example 2.