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

Authors: J. Ochoa-Avendaño; D. A. Garzon-Alvarado; Dorian L. Linero; M. Cerrolaza
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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. This model 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 calculated with the proposed numerical model 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.


---

## Body

## 1. 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–8]. Other authors propose to enhance the trial function of the finite element in order to represent a jump of the displacement field [9–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–14].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–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–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–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.


## 2. 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.


### 2.1. Government Equations in the Continuum

A solid is subjected to body forces vectorb at the domain Ω and surface forces t∗ at the contour Γt with normal vector n, as shown in Figure 1(a). The prescribed displacements are indicated by the vector u∗ at the contour Γu. 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:(1a)∇·σ+b=0∀x∈Ω(1b)σ·n=t∗∀x∈Γt(1c)u=u∗∀x∈ΓuFigure 1
Solid with discontinuity subjected to body and surface forces: (a) general sketch and (b) detailed traction vector inside of the discontinuity.
(a)
(b)In Figure1(b), a discontinuity surface Γd with normal vector nd preserves continuity of the traction vector td, of the following form:(2)σ·nd=td,σ·-nd=-td∀x∈Γd.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]:(3)δWint+δWcoh-δWext=0.The external virtual work δWext is produced by the body and surface forces b and t∗ and is expressed as follows:(4)δWext=∫Ωδu·bdΩ+∫Γtδu·t∗dΓ,where δu is the virtual displacement vector. The internal virtual work δWint is generated by stress tensor σ and the virtual strain tensor δϵ=∇su at Ω⊐Γd, of the following form:(5)δWint=∫Ω⊐Γdδϵ:σdΩ.Linear elastic behavior at domain Ω outside of discontinuity contour Γd, that is, Ω⊐Γd, 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 δWcoh is produced by traction vector td between the faces of the discontinuity Γd and the virtual displacement jump vector δu; thus(6)δWcoh=∫Γdδu·tddΓ.The cohesive loss between the faces of the discontinuity exhibits a nonlinear behavior, which is expressed by the relationship between the traction rate vector t˙d and the rate vector of the displacement jump u˙; thus(7)t˙d=T·u˙,where T is the tangent constitutive tensor of the traction-displacement jump model at the discontinuity.


### 2.2. 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σp reaches the tensile strength of the material σu. In plane stress condition, this stress is defined as follows:(8)σp=max⁡σ1,σ2,0.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 stressnp=npxi+npyj, where npx=cos⁡θp, npy=sin⁡θp, and θp is the angle between x-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.


## 2.1. Government Equations in the Continuum

A solid is subjected to body forces vectorb at the domain Ω and surface forces t∗ at the contour Γt with normal vector n, as shown in Figure 1(a). The prescribed displacements are indicated by the vector u∗ at the contour Γu. 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:(1a)∇·σ+b=0∀x∈Ω(1b)σ·n=t∗∀x∈Γt(1c)u=u∗∀x∈ΓuFigure 1
Solid with discontinuity subjected to body and surface forces: (a) general sketch and (b) detailed traction vector inside of the discontinuity.
(a)
(b)In Figure1(b), a discontinuity surface Γd with normal vector nd preserves continuity of the traction vector td, of the following form:(2)σ·nd=td,σ·-nd=-td∀x∈Γd.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]:(3)δWint+δWcoh-δWext=0.The external virtual work δWext is produced by the body and surface forces b and t∗ and is expressed as follows:(4)δWext=∫Ωδu·bdΩ+∫Γtδu·t∗dΓ,where δu is the virtual displacement vector. The internal virtual work δWint is generated by stress tensor σ and the virtual strain tensor δϵ=∇su at Ω⊐Γd, of the following form:(5)δWint=∫Ω⊐Γdδϵ:σdΩ.Linear elastic behavior at domain Ω outside of discontinuity contour Γd, that is, Ω⊐Γd, 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 δWcoh is produced by traction vector td between the faces of the discontinuity Γd and the virtual displacement jump vector δu; thus(6)δWcoh=∫Γdδu·tddΓ.The cohesive loss between the faces of the discontinuity exhibits a nonlinear behavior, which is expressed by the relationship between the traction rate vector t˙d and the rate vector of the displacement jump u˙; thus(7)t˙d=T·u˙,where T is the tangent constitutive tensor of the traction-displacement jump model at the discontinuity.


## 2.2. 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σp reaches the tensile strength of the material σu. In plane stress condition, this stress is defined as follows:(8)σp=max⁡σ1,σ2,0.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 stressnp=npxi+npyj, where npx=cos⁡θp, npy=sin⁡θp, and θp is the angle between x-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.


## 3. 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.


### 3.1. Nonlinear Analysis by Means of Finite Element Method

A solid in plane stress condition can be represented with a mesh ofne linear triangular finite elements connected by nl 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).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.
(a)
(b)
(c)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 Figure2(c).The discontinuity path on a contourΓd is not known a priori. Consequently, the contours of triangular elements Γh are assumed as likely discontinuity paths. The linear elastic behavior on the domain outside of possible discontinuity paths Ω⊐Γh is represented by linear triangular elements, while the nonlinear behavior on possible discontinuity paths Γh 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 elemente, the transpose virtual displacement vector δueT=δaeTNeT, where Ne is the shape functions matrix and δae is the nodal virtual displacement vector of the element. Likewise, the transpose virtual strain vector δϵeT=δaeTBeT, where Be is the gradient matrix of the element e [39].The external virtual work of the solid, expressed in (4), can be defined in a mesh of ne triangular finite elements and using matrix notation; thus(9)δWext=∑e=1ne∫ΩeδueTbdΩe+∫ΓteδueTt∗dΓe=∑e=1neδaeTfeext=δaTfext,where the nodal virtual displacement of the mesh corresponds to δa and the external force vector of mesh fext is obtained from assembling of vector feext of ne triangular elements, that is, fext=Ae=1nefeext. The external force vector of a triangular finite element is defined as follows:(10)feext=∫ΩeNeTbdΩe+∫ΓteNeTt∗dΓte.The internal virtual work of the solid, expressed in (5), can be defined in a mesh of ne triangular finite elements; thus(11)δWint=∑e=1ne∫Ωe⊐ΓdδϵeTσedΩe=∑e=1neδaeTfeint=δaTfint,where the internal force vector of the mesh fint is obtained from assembling of the vector feint of the ne triangular elements, that is, fint=Ae=1nefeint. The internal force vector of a triangular element is defined of the following form:(12)feint=∫ΩeBeTσedΩe.The cohesive virtual work in the discontinuity of the solid, expressed in (6), can be defined in a mesh of nl link elements; thus(13)δWcoh=∑l=1nlδalTflcoh=δaTfcoh,where the cohesive force vector of the mesh fcoh is obtained from assembling of the vector flcoh of the nl link elements, that is, fcoh=Al=1nlflcoh.Equations (9), (11), and (13) are substituted into (3), and the transpose nodal virtual displacement vector δaT is factorized and canceled. The obtained equilibrium condition establishes that fint+fcoh-fext=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:(14)Ψ=fint+fcoh-fext.According to Newton-Raphson method, the residual force vectorΨai+Δai evaluated on nodal displacement ai plus the trial incremental displacement vector Δai in the iteration i is equal to(15)Ψai+Δai=Ψai+∂Ψ∂aaiΔai.The derivate of the residual force vector, presented in (14), with respect to nodal displacement is(16)∂Ψ∂a=∂fint∂a+∂fcoh∂a-∂fext∂a=Kint+Kcoh.The external force vector does not change with respect to nodal displacement in the iterative process; then ∂afext=0. However, the derivative of the internal and cohesive force vectors with respect to nodal displacement are equal to the tangent stiffness matrices Kint and Kcoh, respectively.Nonlinear numerical solution method searches the trial incremental nodal displacementΔai for which the residual force vector Ψai+Δai is equal to zero. In each iteration i of this procedure, an equations system is solved; thus(17)Δai=-Kintai+Kcohai-1Ψai.The tangent stiffness matricesKint and Kcoh are obtained from assembling of the tangent stiffness matrices Keint and Klcoh of each finite element, respectively; thus(18a)∂fint∂a=Kint=Ae=1neKeint(18b)∂fcoh∂a=Kcoh=Al=1nlKlcoh,where e identifies the two-dimensional element from 1 to ne and l indicates the link element from 1 to nl.A linear elastic material is assumed in the domain of the triangular elements; thereforeσe=Deϵe, where the constitutive matrix De is constant with respect to the strain state and ϵe=Beae. The internal force vector of a triangular finite element is feint=Keintae, where the nodal displacement vector is ae and the stiffness matrix Keint is defined of the following form [39]:(19)Keint=∫ΩeBeTDeBedΩe.This matrix depends on thickness t of the element and Young’s modulus E and Poisson’s relation ν of material.In contrast, each link elementl exhibits nonlinear behavior, in which the cohesive force vector rate f˙coh corresponds to(20)f˙lcoh=Klcoha˙l,where a˙l is the vector rate of the displacement on the duplicate nodes of link element l. Tangent stiffness matrix Klcoh of this element is defined as follows:(21)Klcoh=klcohnlx2nlxnly-nlx2-nlxnlynlxnlynly2-nlxnly-nly2-nlx2-nlxnlynlx2nlxnly-nlxnly-nly2nlxnlynly2.The longitudinal axis of link element is normal to the side of the neighboring triangular elements and it is defined by directional vector nl=nlxi+nlyj, where nlx=cos⁡θl and nly=sin⁡θl. The angle between the longitudinal axis of the link element and x-global axis is called θl, 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 factorklcoh is equal to zero, when the cohesion is totally lost, and tends to infinity when the displacement compatibility is preserved.


### 3.2. 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. Figure3(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.Figure 3
Finite element meshes: (a) conventional mesh and (b) exploded mesh out of scale.
(a)
(b)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 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.


### 3.3. Evaluation of the Tangent Stiffness of Link Elements

The nodal stress componentsσxx, σyy, and σxy 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 n1 and n2 are computed as indicated in some references for plane stress condition [38, 39]. The positive maximum principal stress σp and its direction np are selected according to (8). The angle between the x-global axis and the directional vector of the positive maximum principal stress np is called θp.If the positive maximum principal stressσp reaches the tensile strength of the material σu, then the tangent stiffness in one of the link elements m,…,nm associated with this set overlapping nodes is modified. If the orientation nl of the link element has the least difference with respect to comparative direction nc, then the tangent stiffness factor klcoh of (21) is equal to zero. Otherwise, this factor preserves the tendency to infinity, that is,(22)klcoh=0ifσp=σu,θlerr=min⁡θmerr,…,θnmerr⟶∞otherwise.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 nc, as shown Figure 4(b).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.
(a)
(b)
(c)Since there is no perfect alignment between the comparative directionnc and the link orientation n• selected because it has the least difference with respect to nc, there is an implicit error θ•err in computing of crack direction. This error is defined in(23)θ•err=θc-θ•,where the angle between x-global axis and the directional vector of a link element n• is called θ• and the angle between x-global axis and the vector of the comparative direction nc is called θc, as shown in Figure 4(c).The angle of the comparative directionθc is equal to the angle of the direction principal θp in the first set of overlapping nodes where σp=σu. 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 vectornl 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 nc, as shown in Figure 4(b).


### 3.4. 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θcr in the set of overlapping nodes r with σpr=σu is defined as follows:(24)θcr=θpr+θlerrr-1,where θlerrr-1 is the error angle computed in (23), which was obtained at the previous set of overlapping nodes r-1 with σpr-1=σu. As was mentioned above, the purpose of this correction is to avoid the implicit error resulting of misalignment of the comparative directions nc 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.Figure 5
Detailed discontinuity path in a mesh of triangular finite elements: (a) with correction and (b) without correction.


## 3.1. Nonlinear Analysis by Means of Finite Element Method

A solid in plane stress condition can be represented with a mesh ofne linear triangular finite elements connected by nl 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).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.
(a)
(b)
(c)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 Figure2(c).The discontinuity path on a contourΓd is not known a priori. Consequently, the contours of triangular elements Γh are assumed as likely discontinuity paths. The linear elastic behavior on the domain outside of possible discontinuity paths Ω⊐Γh is represented by linear triangular elements, while the nonlinear behavior on possible discontinuity paths Γh 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 elemente, the transpose virtual displacement vector δueT=δaeTNeT, where Ne is the shape functions matrix and δae is the nodal virtual displacement vector of the element. Likewise, the transpose virtual strain vector δϵeT=δaeTBeT, where Be is the gradient matrix of the element e [39].The external virtual work of the solid, expressed in (4), can be defined in a mesh of ne triangular finite elements and using matrix notation; thus(9)δWext=∑e=1ne∫ΩeδueTbdΩe+∫ΓteδueTt∗dΓe=∑e=1neδaeTfeext=δaTfext,where the nodal virtual displacement of the mesh corresponds to δa and the external force vector of mesh fext is obtained from assembling of vector feext of ne triangular elements, that is, fext=Ae=1nefeext. The external force vector of a triangular finite element is defined as follows:(10)feext=∫ΩeNeTbdΩe+∫ΓteNeTt∗dΓte.The internal virtual work of the solid, expressed in (5), can be defined in a mesh of ne triangular finite elements; thus(11)δWint=∑e=1ne∫Ωe⊐ΓdδϵeTσedΩe=∑e=1neδaeTfeint=δaTfint,where the internal force vector of the mesh fint is obtained from assembling of the vector feint of the ne triangular elements, that is, fint=Ae=1nefeint. The internal force vector of a triangular element is defined of the following form:(12)feint=∫ΩeBeTσedΩe.The cohesive virtual work in the discontinuity of the solid, expressed in (6), can be defined in a mesh of nl link elements; thus(13)δWcoh=∑l=1nlδalTflcoh=δaTfcoh,where the cohesive force vector of the mesh fcoh is obtained from assembling of the vector flcoh of the nl link elements, that is, fcoh=Al=1nlflcoh.Equations (9), (11), and (13) are substituted into (3), and the transpose nodal virtual displacement vector δaT is factorized and canceled. The obtained equilibrium condition establishes that fint+fcoh-fext=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:(14)Ψ=fint+fcoh-fext.According to Newton-Raphson method, the residual force vectorΨai+Δai evaluated on nodal displacement ai plus the trial incremental displacement vector Δai in the iteration i is equal to(15)Ψai+Δai=Ψai+∂Ψ∂aaiΔai.The derivate of the residual force vector, presented in (14), with respect to nodal displacement is(16)∂Ψ∂a=∂fint∂a+∂fcoh∂a-∂fext∂a=Kint+Kcoh.The external force vector does not change with respect to nodal displacement in the iterative process; then ∂afext=0. However, the derivative of the internal and cohesive force vectors with respect to nodal displacement are equal to the tangent stiffness matrices Kint and Kcoh, respectively.Nonlinear numerical solution method searches the trial incremental nodal displacementΔai for which the residual force vector Ψai+Δai is equal to zero. In each iteration i of this procedure, an equations system is solved; thus(17)Δai=-Kintai+Kcohai-1Ψai.The tangent stiffness matricesKint and Kcoh are obtained from assembling of the tangent stiffness matrices Keint and Klcoh of each finite element, respectively; thus(18a)∂fint∂a=Kint=Ae=1neKeint(18b)∂fcoh∂a=Kcoh=Al=1nlKlcoh,where e identifies the two-dimensional element from 1 to ne and l indicates the link element from 1 to nl.A linear elastic material is assumed in the domain of the triangular elements; thereforeσe=Deϵe, where the constitutive matrix De is constant with respect to the strain state and ϵe=Beae. The internal force vector of a triangular finite element is feint=Keintae, where the nodal displacement vector is ae and the stiffness matrix Keint is defined of the following form [39]:(19)Keint=∫ΩeBeTDeBedΩe.This matrix depends on thickness t of the element and Young’s modulus E and Poisson’s relation ν of material.In contrast, each link elementl exhibits nonlinear behavior, in which the cohesive force vector rate f˙coh corresponds to(20)f˙lcoh=Klcoha˙l,where a˙l is the vector rate of the displacement on the duplicate nodes of link element l. Tangent stiffness matrix Klcoh of this element is defined as follows:(21)Klcoh=klcohnlx2nlxnly-nlx2-nlxnlynlxnlynly2-nlxnly-nly2-nlx2-nlxnlynlx2nlxnly-nlxnly-nly2nlxnlynly2.The longitudinal axis of link element is normal to the side of the neighboring triangular elements and it is defined by directional vector nl=nlxi+nlyj, where nlx=cos⁡θl and nly=sin⁡θl. The angle between the longitudinal axis of the link element and x-global axis is called θl, 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 factorklcoh is equal to zero, when the cohesion is totally lost, and tends to infinity when the displacement compatibility is preserved.


## 3.2. 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. Figure3(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.Figure 3
Finite element meshes: (a) conventional mesh and (b) exploded mesh out of scale.
(a)
(b)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 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.


## 3.3. Evaluation of the Tangent Stiffness of Link Elements

The nodal stress componentsσxx, σyy, and σxy 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 n1 and n2 are computed as indicated in some references for plane stress condition [38, 39]. The positive maximum principal stress σp and its direction np are selected according to (8). The angle between the x-global axis and the directional vector of the positive maximum principal stress np is called θp.If the positive maximum principal stressσp reaches the tensile strength of the material σu, then the tangent stiffness in one of the link elements m,…,nm associated with this set overlapping nodes is modified. If the orientation nl of the link element has the least difference with respect to comparative direction nc, then the tangent stiffness factor klcoh of (21) is equal to zero. Otherwise, this factor preserves the tendency to infinity, that is,(22)klcoh=0ifσp=σu,θlerr=min⁡θmerr,…,θnmerr⟶∞otherwise.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 nc, as shown Figure 4(b).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.
(a)
(b)
(c)Since there is no perfect alignment between the comparative directionnc and the link orientation n• selected because it has the least difference with respect to nc, there is an implicit error θ•err in computing of crack direction. This error is defined in(23)θ•err=θc-θ•,where the angle between x-global axis and the directional vector of a link element n• is called θ• and the angle between x-global axis and the vector of the comparative direction nc is called θc, as shown in Figure 4(c).The angle of the comparative directionθc is equal to the angle of the direction principal θp in the first set of overlapping nodes where σp=σu. 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 vectornl 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 nc, as shown in Figure 4(b).


## 3.4. 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θcr in the set of overlapping nodes r with σpr=σu is defined as follows:(24)θcr=θpr+θlerrr-1,where θlerrr-1 is the error angle computed in (23), which was obtained at the previous set of overlapping nodes r-1 with σpr-1=σu. As was mentioned above, the purpose of this correction is to avoid the implicit error resulting of misalignment of the comparative directions nc 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.Figure 5
Detailed discontinuity path in a mesh of triangular finite elements: (a) with correction and (b) without correction.


## 4. 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.


### 4.1. Example1: Three-Point Beam with Nonconcentric Notch

A simply supported beam of polymethylmethacrylate (PMMA) is subjected to loadP at the upper center part. Figure 6 shows the geometry and load applied on the beam, where the thickness is t=5 mm, and the notch has a of depth and separated b from the loading line on the bottom center face of the beam. The mechanical features of the material are as follows: Young’s modulus E=3100 MPa, Poisson’s relation ν=0.4, and tensile strength σu=76 MPa. Two cases with different values of a and b are studied. In Case I, a=10 mm and b=60 mm are defined and three meshes with different amount of triangular elements are simulated (Table 1).Table 1
Meshes of Case I of Example1.
Number of elements
Number of nodes
Mesh 1
1413
765
Mesh 2
4412
2310
Mesh 3
9432
4869Figure 6
Sketch of the three points beam with nonconcentric notch. The measurements are given in millimeters.Figure7 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.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.
(a)
Mesh 1
(b)
Mesh 2
(c)
Mesh 3The simulation with mesh2 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,a=15 mm and b=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).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.
(a)
(b)


### 4.2. Example2: 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 Figure9. The principal stress and strain are high in the neighborhood of each hole. This produces a curvilinear cracking pattern unique [42]. The depth a and the distance b 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–46].Table 2
Finite element meshes of Example2.
Case I
Case II
Number of elements
Number of nodes
Number of elements
Number of nodes
Mesh 1
2630
1398
2191
1172
Mesh 2
10520
5428
10284
5311
Mesh 3
23200
11846
16296
8360Figure 9
Sketch of the three-point beam with notch and three holes. The measurements are given in millimeters.A deptha=10 mm and the distance b=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–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.Figure 10
Comparison between the discontinuity path in numerical model and the crack path of experimental test in Case I of Example  2.
(a)
Mesh 1
(b)
Mesh 2
(c)
Mesh 3Other 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 saw-tooth form of the proposed model.Figure11 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 a=15 mm and the distance b=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.Figure 11
Evolution of the discontinuity path in the numerical simulation of the finest mesh in Case I of Example  2.
(a)
10th step
(b)
30th step
(c)
Final stepFigure 12
Comparison between the discontinuity path in numerical model and the crack path of experimental test in Case II of Example  2.
(a)
Mesh 1
(b)
Mesh 2
(c)
Mesh 3


### 4.3. Example3: Beam with Notch and Two Loading Points

A concrete beam is subjected to the loads0.13P and P, as shown in Figure 13. The beam has a thickness of t=156 mm and a notch of 82 mm of depth located nearby the support. The material has a Young modulus of E=24800 MPa, Poisson’s relation of ν=0.18, and tensile strength of σu=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.Figure 13
Sketch of the beam with notch and two loading points. The measures are given in millimeters.The notorious relative displacement between nodes of the deformed shape of the numerical model is shown in Figure14(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.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.
(a)
(b)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.


### 4.4. Sensibility of the Discontinuity Path to the Initial Stiffness of Link Elements

The initial stiffness factor of link elementklcoh,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 klcoh,0 is much higher than the maximum coefficient of the stiffness matrix of an equilateral triangular element with linear elastic behavior kemaxint=tE/31-ν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. Table 3 indicates the maximum number of iterations per loading step in column ni and whether (Y) or not (N) the discontinuity path is obtained in column dp, for four values of klcoh,0. The simulations with initial stiffness factor between 1×103 and 1×105 times kemaxint 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×106 times kemaxint, Examples  1 and 2.Table 3
Results of the sensibility analysis of the discontinuity path to the initial stiffness of link elements.
k
l
c
o
h
,
0
/
k
e
(
m
a
x
)
i
n
t
=
1
×
1
0
3
1
×
1
0
4
1
×
1
0
5
1
×
1
0
6
Figuredp
n
i
d
p
n
i
d
p
n
i
d
p
n
i
d
p
Example1, Case I
2
Y
2
Y
2
Y
0
N
Figure7(c)
Example1, Case II
2
Y
2
Y
2
Y
0
N
Figure8(a)
Example2, Case I
2
Y
2
Y
2
Y
0
N
Figure10(c)
Example2, Case II
2
Y
2
Y
2
Y
0
N
Figure12(c)
Example3
2
Y
2
Y
2
Y
2
Y
Figure14(a)Low values of initial stiffness of the link elements were not assigned to the simulations because these do not represent the 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 of1×10-3. The maximum number of iterations per loading step was equal to 2, as shown in Table 3.


## 4.1. Example1: Three-Point Beam with Nonconcentric Notch

A simply supported beam of polymethylmethacrylate (PMMA) is subjected to loadP at the upper center part. Figure 6 shows the geometry and load applied on the beam, where the thickness is t=5 mm, and the notch has a of depth and separated b from the loading line on the bottom center face of the beam. The mechanical features of the material are as follows: Young’s modulus E=3100 MPa, Poisson’s relation ν=0.4, and tensile strength σu=76 MPa. Two cases with different values of a and b are studied. In Case I, a=10 mm and b=60 mm are defined and three meshes with different amount of triangular elements are simulated (Table 1).Table 1
Meshes of Case I of Example1.
Number of elements
Number of nodes
Mesh 1
1413
765
Mesh 2
4412
2310
Mesh 3
9432
4869Figure 6
Sketch of the three points beam with nonconcentric notch. The measurements are given in millimeters.Figure7 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.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.
(a)
Mesh 1
(b)
Mesh 2
(c)
Mesh 3The simulation with mesh2 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,a=15 mm and b=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).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.
(a)
(b)


## 4.2. Example2: 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 Figure9. The principal stress and strain are high in the neighborhood of each hole. This produces a curvilinear cracking pattern unique [42]. The depth a and the distance b 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–46].Table 2
Finite element meshes of Example2.
Case I
Case II
Number of elements
Number of nodes
Number of elements
Number of nodes
Mesh 1
2630
1398
2191
1172
Mesh 2
10520
5428
10284
5311
Mesh 3
23200
11846
16296
8360Figure 9
Sketch of the three-point beam with notch and three holes. The measurements are given in millimeters.A deptha=10 mm and the distance b=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–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.Figure 10
Comparison between the discontinuity path in numerical model and the crack path of experimental test in Case I of Example  2.
(a)
Mesh 1
(b)
Mesh 2
(c)
Mesh 3Other 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 saw-tooth form of the proposed model.Figure11 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 a=15 mm and the distance b=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.Figure 11
Evolution of the discontinuity path in the numerical simulation of the finest mesh in Case I of Example  2.
(a)
10th step
(b)
30th step
(c)
Final stepFigure 12
Comparison between the discontinuity path in numerical model and the crack path of experimental test in Case II of Example  2.
(a)
Mesh 1
(b)
Mesh 2
(c)
Mesh 3


## 4.3. Example3: Beam with Notch and Two Loading Points

A concrete beam is subjected to the loads0.13P and P, as shown in Figure 13. The beam has a thickness of t=156 mm and a notch of 82 mm of depth located nearby the support. The material has a Young modulus of E=24800 MPa, Poisson’s relation of ν=0.18, and tensile strength of σu=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.Figure 13
Sketch of the beam with notch and two loading points. The measures are given in millimeters.The notorious relative displacement between nodes of the deformed shape of the numerical model is shown in Figure14(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.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.
(a)
(b)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.


## 4.4. Sensibility of the Discontinuity Path to the Initial Stiffness of Link Elements

The initial stiffness factor of link elementklcoh,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 klcoh,0 is much higher than the maximum coefficient of the stiffness matrix of an equilateral triangular element with linear elastic behavior kemaxint=tE/31-ν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. Table 3 indicates the maximum number of iterations per loading step in column ni and whether (Y) or not (N) the discontinuity path is obtained in column dp, for four values of klcoh,0. The simulations with initial stiffness factor between 1×103 and 1×105 times kemaxint 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×106 times kemaxint, Examples  1 and 2.Table 3
Results of the sensibility analysis of the discontinuity path to the initial stiffness of link elements.
k
l
c
o
h
,
0
/
k
e
(
m
a
x
)
i
n
t
=
1
×
1
0
3
1
×
1
0
4
1
×
1
0
5
1
×
1
0
6
Figuredp
n
i
d
p
n
i
d
p
n
i
d
p
n
i
d
p
Example1, Case I
2
Y
2
Y
2
Y
0
N
Figure7(c)
Example1, Case II
2
Y
2
Y
2
Y
0
N
Figure8(a)
Example2, Case I
2
Y
2
Y
2
Y
0
N
Figure10(c)
Example2, Case II
2
Y
2
Y
2
Y
0
N
Figure12(c)
Example3
2
Y
2
Y
2
Y
2
Y
Figure14(a)Low values of initial stiffness of the link elements were not assigned to the simulations because these do not represent the 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 of1×10-3. The maximum number of iterations per loading step was equal to 2, as shown in Table 3.


## 5. 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.



---

*Source: 3861526-2017-07-02.xml*