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.

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 [58]. Other authors propose to enhance the trial function of the finite element in order to represent a jump of the displacement field [911]. 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 [1214].

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 [2430]. 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 [2527], 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 [3437].

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 vector at the domain and surface forces at the contour with normal vector , as shown in Figure 1(a). The prescribed displacements are indicated by the vector at the contour . The state stress in each material point is expressed by the tensor . Equations (1a) to (1c) present the equilibrium and boundary conditions in the solid:

In Figure 1(b), a discontinuity surface with normal vector preserves continuity of the traction vector , 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 is produced by the body and surface forces and and is expressed as follows:where is the virtual displacement vector. The internal virtual work is generated by stress tensor and the virtual strain tensor 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 , where is the linear elastic constitutive tensor of the material. The cohesive virtual work is produced by traction vector between the faces of the discontinuity and the virtual displacement jump vector ; thusThe 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 ; thuswhere 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 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 , where , , 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.

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 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 , where is the shape functions matrix and is the nodal virtual displacement vector of the element. Likewise, the transpose virtual strain vector , where 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; thuswhere the nodal virtual displacement of the mesh corresponds to and the external force vector of mesh is obtained from assembling of vector 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; thuswhere the internal force vector of the mesh is obtained from assembling of the vector 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; thuswhere the cohesive force vector of the mesh is obtained from assembling of the vector of the link elements, that is, .

Equations (9), (11), and (13) are substituted into (3), and the transpose nodal virtual displacement vector is factorized and canceled. The obtained equilibrium condition establishes that . 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, . The residual force vector is defined as follows:

According to Newton-Raphson method, the residual force vector evaluated on nodal displacement plus the trial incremental displacement vector in the iteration is equal toThe derivate of the residual force vector, presented in (14), with respect to nodal displacement isThe external force vector does not change with respect to nodal displacement in the iterative process; then . However, the derivative of the internal and cohesive force vectors with respect to nodal displacement are equal to the tangent stiffness matrices and , respectively.

Nonlinear numerical solution method searches the trial incremental nodal displacement for which the residual force vector is equal to zero. In each iteration of this procedure, an equations system is solved; thus

The tangent stiffness matrices and are obtained from assembling of the tangent stiffness matrices and of each finite element, respectively; thuswhere identifies the two-dimensional element from 1 to and indicates the link element from 1 to .

A linear elastic material is assumed in the domain of the triangular elements; therefore , where the constitutive matrix is constant with respect to the strain state and . The internal force vector of a triangular finite element is , where the nodal displacement vector is and the stiffness matrix 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 corresponds towhere is the vector rate of the displacement on the duplicate nodes of link element . Tangent stiffness matrix 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 , where and . 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 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. 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 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 , , and on each set of overlapping nodes are obtained from the average among the stress components of surrounding linear triangular elements. The principal stresses and and their directions and are computed as indicated in some references for plane stress condition [38, 39]. The positive maximum principal stress and its direction are selected according to (8). The angle between the -global axis and the directional vector of the positive maximum principal stress 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 of the link element has the least difference with respect to comparative direction , then the tangent stiffness factor 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 , as shown Figure 4(b).

Since there is no perfect alignment between the comparative direction and the link orientation selected because it has the least difference with respect to , there is an implicit error in computing of crack direction. This error is defined inwhere the angle between -global axis and the directional vector of a link element is called and the angle between -global axis and the vector of the comparative direction 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 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 , 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 in the set of overlapping nodes with is defined as follows:where is the error angle computed in (23), which was obtained at the previous set of overlapping nodes with . As was mentioned above, the purpose of this correction is to avoid the implicit error resulting of misalignment of the comparative directions 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.

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. 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  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  MPa, Poisson’s relation , and tensile strength  MPa. Two cases with different values of and are studied. In Case I,  mm and  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 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 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 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,  mm and  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).

4.2. 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 [4246].

A depth  mm and the distance  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 [4347]. 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 saw-tooth 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  mm and the distance  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.

4.3. Example  3: Beam with Notch and Two Loading Points

A concrete beam is subjected to the loads and , as shown in Figure 13. The beam has a thickness of  mm and a notch of 82 mm of depth located nearby the support. The material has a Young modulus of  MPa, Poisson’s relation of , and tensile strength of  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.

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

The initial stiffness factor of link element 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 is much higher than the maximum coefficient of the stiffness matrix of an equilateral triangular element with linear elastic behavior . 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 and whether (Y) or not (N) the discontinuity path is obtained in column , for four values of . The simulations with initial stiffness factor between and times 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 times , 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 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 . 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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.