Simulation of cohesive crack growth by a variable-node XFEM

Abstract

A new computational approach that combines the extended finite element method associated with variable-node elements and cohesive zone model is developed. By using a new enriched technique based on sign function, the proposed model using 4-node quadrilateral elements can eliminate the blending element problem. It also allows modeling the equal stresses at both sides of the crack in the crack-tip as assumed in the cohesive model, and is able to simulate the arbitrary crack-tip location. The multiscale mesh technique associated with variable-node elements and the arc-length method further improve the efficiency of the developed approach. The performance and accuracy of the present approach are illustrated through numerical experiments considering both mode-I and mixed-mode fracture in concrete.

Keywords

extended finite element method / cohesive zone model / sign function / crack propagation

Introduction

Quasi-brittle materials such as concrete are of great importance in many engineering applications. Such materials are used for design of key components in civil infrastructures, urban systems, and many others. Cracks are inevitable in concrete structures, and that directly affect the safety and the durability of concrete structures. Hence, reliable predictions of crack growth and crack width in concrete structures are important issues.

The extended finite element method (XFEM) [¹,²] is a powerful numerical method for modeling the evolution of crack as its computational mesh is independent of the crack geometry, so the re-meshing in modeling crack growth is no longer required. In the past decades, there are numerous studies concerning the improvement or application of the XFEM for various discontinuous problems [³–¹⁵]. In particular, the failure in concrete materials is preceded by a gradual development of a nonlinear fracture process zone and a localization of strain. In other words, at the crack-tip region in concrete there exists a fracture process zone (FPZ) [¹⁶], where nonlinear phenomena takes place, while the rest of the body exhibits an elastic behavior. Cohesive crack model is the simplest model to represent the FPZ. The cohesive crack model for analyzing metals was originally proposed by Dugdale [¹⁷] and Barenblatt [¹⁸]. Later, Hillerborg et al. [¹⁹] proposed some concrete cohesive models by introducing fracture energy into the cohesive crack model. Cohesive crack model has been widely used in nonlinear fracture mechanics of quasi-brittle materials [²⁰–²³]. Some studies have been performed for cohesive cracks by using the XFEM. In the cohesive crack model, the crack-opening displacement vanishes at the crack tip without the stress singularity. In general, there are two XFEM enrichment strategies for capturing cohesive cracks, i.e., 1) only the Heaviside step function is used as the enrichment functions, for the first the crack tip is always located at element edge [²⁴], later the position of crack tip may be arbitrary [²⁵,²⁶]; 2) the Heaviside step function and non-singular branch functions are used as the enrichment functions [²⁷], the position of crack tip is arbitrary, but the blending elements exist, which decrease the global convergence rate.

Wells and Sluys [²⁴] developed a cohesive crack model in the framework of the XFEM with the 6-node triangle elements by employing the Heaviside step function as the enrichment function. A straight crack segment is introduced through the entire element when the maximum principal stress at one integration point in the element ahead of the crack exceeds the tensile strength of the material. The limitation is that the crack tip is always located at element edge. Remmers et al. [²⁷] presented a cohesive segment method for crack growth modeling. The crack is represented by a set of overlapping cohesive segments. A new cohesive segment is added when the major principal stress at an integration point within an element reaches the cohesive strength. The limitations are (1) the added cohesive segment passes through the entire element, and (2) the crack path is not continuous and smooth. Moës and Belytschko [²⁸] used the Heaviside step function and non-singular branch functions as the enrichment functions, thus the position of crack tip is arbitrary. The growth of the cohesive zone is governed by requiring the stress intensity factors at the tip of the cohesive zone to be zero. The blending elements exist, and it decreases the overall convergence rate. To avoid the difficulties associated with the branch functions in the tip element, Zi and Belytschko [²⁶] developed a new crack-tip element enriched with the Heaviside step function for the cohesive crack growth modeling, which allows the crack tip to be any location. All cracked elements are enriched by the sign function so that no blending of the local partition of unity is required. However, stresses on both sides of the crack are not equal. Later, Asferg et al. [²⁵] developed a new partly cracked 3-node triangle XFEM element for cohesive crack growth by superposition of the standard nodal shape functions for the element and standard nodal shape functions for a sub-triangle of the cracked element, which can model equal stresses on both sides of the crack. Zhang and Bui [²⁹] developed two new solution algorithms for cohesive crack growth modeling based on Newton-Raphson method, which results in a symmetric tangent matrix and avoid the inconvenience in solving the inversion of an unsymmetrical Jacobian matrix encountered in Ref. [²⁶]. Mougaard et al. [³⁰] developed a new cohesive crack tip element together with a coherent fully cracked element based on the XFEM. The elements are based on a double enriched displacement field of linear strain triangle type, both a discontinuous and a continuous displacement field are used as the enrichment in the two enrichment triangles, this reproduces equal stresses at both sides of the crack at the tip. Cox [³¹] proposed an XFEM with analytical enrichment for cohesive crack problem. These functions can represent displacement gradients in the vicinity of the cohesive crack. Wu and Li [³²] proposed an improved stable XFEM with a novel enrichment function for cohesive crack problem. The Heaviside function stabilized by its linear interpolant is defined as enrichment function, and the improved stable XFEM has sufficient accuracy and well-conditioned system matrix.

There are some other approaches that have also been introduced in the literature for modeling crack propagation. For instance, a cracking-particle method for modeling discrete cracks was presented in Ref. [³³]. The crack is modeled by a local enrichment of the test and trial functions with a sign function, while the crack growth is described discretely by activation of crack surfaces at individual particles, and the crack's topology is not represented. The cracking-particle method associated with cohesive laws for 3D large deformation problems for arbitrary evolving cracks was reported in Ref. [³⁴]. To save computational cost and enhance the accuracy around the crack tip, h-adaptivity is incorporated in the method. Recently, Rabczuk et al. [³⁵] improved the cracking-particle method by modeling the crack with splitting particles located on opposite sides of the associated crack segments, and without additional unknowns in the variational formulation to capture the displacement discontinuity are required. Peridynamics (PD) is a recently introduced non-local reformulation of classical elasticity theory for modeling materials with discontinuities such as cracks. The partial differential equations of classical solid mechanics are replaced with integro-differential equations, and the internal forces within a body are treated as a network of interactions between material points. Ren et al. [³⁶,³⁷] developed a dual-horizon PD (DH-PD) formulation that naturally includes varying horizon sizes and completely solves the ‘ghost force’ issue. The traditional PD is a special case of the DH-PD, and DH-PD is less sensitive to the spatial than the traditional PD formulation. The cracking-particle method and the DH-PD allow for complex crack patterns without any crack tracking algorithms. The meshfree methods are based on interaction of each node with all its neighbors, and do not require mesh, so they are suitable for cracking simulations. Rabczuk et al. [³⁸] presented a meshfree method for treating fluid-structure interaction of fracturing structures under impulsive loads. The cracks are treated by introducing either discrete or continuous discontinuities into the approximation, and coupling is realized by a master–slave scheme. Rabczuk et al. [³⁹] reviewed different crack tracking techniques in three-dimensions applicable in the context of partition of unity methods, especially meshfree methods. A phantom-node method for three-node shell elements with arbitrary cracks was presented in Ref. [⁴⁰], which are independently of the mesh. In the phantom-node method, the discontinuity is described by superposed elements and phantom nodes by a rearrangement of the XFEM basis and the nodal degrees of freedom. Cracks are treated by adding phantom nodes and superposing elements on the original mesh.

The objective of this study is to develop a computational approach in terms of XFEM with local refined mesh in association with variable-node elements for cohesive crack growth in concrete. The present XFEM is based on the 4-node quadrilateral elements. The multiscale mesh method is adopted in order to save the computational cost. Herein, the fine scale mesh is only used in the concerning area where crack may take place, whereas coarse-scale mesh is for the rest. Then, to connect nonconforming meshes between two distinct scales, the variable-node elements [⁴,⁴¹] is adopted. The number of nodes on the side/surface is arbitrary in the variable-node elements. Hence no special treatment is required for the system matrix under any circumstance. It should be noticed that the present local refined mesh XFEM using multiscale mesh is sometimes termed as multiscale XFEM in the manuscript, and here we neither use error estimations, nor adaptivity algorithm. All cracked elements including the cut elements and tip elements are enriched by the sign function. The proposed enrichment strategy owns the following advantages and characteristics: 1) it is straightforward to handle the blending elements issue; 2) the crack-tip position is arbitrary; 3) the equal stresses on both sides of the crack-tip element hold. In this study, the arc-length method is also used to solve the nonlinear governing equations, which further provides a robust computational framework for nonlinear cohesive crack growth analysis. It is worth mentioning that the implementation of the proposed model is straightforward, thanks to the use of the unified sign function in the enrichment strategies. The efficiency and accuracy of the developed model are demonstrated through three numerical examples including one for mode-I fracture and two for mixed-mode fracture problems. The contributions of this paper are as follows: 1) a new enriched technique based on sign function is developed by which the blending element problem can be eliminated. The equal stresses at both sides of the crack in the crack-tip can be obtained as assumed in the cohesive model, and the crack-tip location is arbitrary; 2) the multiscale mesh technique associated with variable-node elements is developed to improve the efficiency and accuracy.

The outline of this paper is given as follows. In Section 2, a brief description for fundamental equations of cohesive cracks is provided. Section 3 presents a new XFEM formulation for cohesive crack growth. Section 4 analyzes the numerical examples. Finally, the major conclusions are summarized in Section 5.

Fundamental equations

Consider a cracked body \mathrm{\Omega } with its boundary \mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{u}}+{\mathrm{\Gamma }}_{\mathrm{t}}+{\mathrm{\Gamma }}_{\mathrm{c}},{\mathrm{\Gamma }}_{\mathrm{u}}\cap {\mathrm{\Gamma }}_{\mathrm{t}}=\varnothing, where {\mathrm{\Gamma }}_{\mathrm{u}}, {\mathrm{\Gamma }}_{\mathrm{t}}, and {\mathrm{\Gamma }}_{\mathrm{c}} are the displacement boundary, the external traction boundary and the crack face, respectively, as sketched in Fig. 1(a). In this study, the fracture behavior at the region near the crack-tip, called the cohesive zone, is governed by the cohesive tractions law.

Without considering the body force, the boundary value problem is described as follows:

(1)\nabla \cdot \sigma =\mathrm{0},\hspace{0.17em}\mathrm{in}\hspace{0.17em}\mathrm{\Omega },

(2)\sigma \cdot n=\lambda {\tau }_{\mathrm{0}},\hspace{0.17em}\mathrm{on}\hspace{1em}{{\displaystyle \Gamma }}_t,

(3)u=\bar{u},\hspace{0.17em}\mathrm{on}\hspace{1em}{{\displaystyle \Gamma }}_u,

(4){\sigma }^{+}\cdot n^{+}={\sigma }^{-}\cdot n^{-}={\tau }^c,\hspace{0.17em}\hspace{0.17em}\mathrm{on}\hspace{1em}{\mathrm{\Gamma }}_{\mathrm{coh}},

where \sigma is the Cauchy stress, u is the displacement; n is the outward unit normal vector on \mathrm{\Gamma }; \bar{u} is the prescribed displacement; {\tau }_0 is the normalized traction, \lambda is the loading factor; and

n^{+}=-n^{-}=n^c

are the outward unit normal on the crack faces, {\sigma }^{\pm } and {\tau }^c are the stress and normal traction at the crack faces, as depicted in Fig. 1(b).

As small strain is assumed, the following relation holds

(5)\epsilon (u)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[\nabla u+{(\nabla u)}^T].

Along with the constitutive equations

(6a)\sigma =C\epsilon ,\hspace{0.17em}\mathrm{in}\hspace{0.17em}\Omega ,

(6b){\tau }^c={\tau }^c(w),\hspace{0.17em}\hspace{0.17em}\mathrm{on}\hspace{1em}{{\displaystyle \Gamma }}_{\mathrm{coh}},

where C is the matrix form of the fourth rank elasticity tensor; and w is the crack opening, which is defined as

(6c)w=u^{-}-u^{+},\hspace{0.17em}\hspace{0.17em}\mathrm{on}\hspace{1em}{{\displaystyle \mathrm{\Gamma }}}_{\mathrm{coh}},

with u^{-} and u^{+} are the displacements on the crack faces, respectively.

Applying the principle of virtual work to Eq. (6), the weak-form of the equilibrium equations can then be obtained by

(7){\displaystyle {\int }_{\Omega }\sigma }:\epsilon (v)\mathrm{d\Omega }={\displaystyle {\int }_{{\Gamma }_t}\lambda {\tau }_{\mathrm{0}}}\cdot v\mathrm{ds}-{\displaystyle {\int }_{{\Gamma }_{\mathrm{coh}}^{+}}{\sigma }^{+}{n}^c}\cdot v^{+}\mathrm{ds\; +}{\displaystyle {\int }_{{\Gamma }_{\mathrm{coh}}^{-}}{\sigma }^{-}{n}^c}\cdot v^{-}\mathrm{ds,}

where v is the admissible virtual displacement.

Considering the continuity condition in Eq. (4) and w(v)=v^{-}-v^{+}, Eq. (7) can be rewritten as

(8){\displaystyle {\int }_{\mathrm{\Omega }}\sigma }:\epsilon (v)\mathrm{d\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{coh}}^{}}{\tau }^c\cdot w(v)}\mathrm{ds}={\displaystyle {\int }_{{\mathrm{\Gamma }}_t}\lambda }{\tau }_0\cdot v\mathrm{ds}.

In general, the cohesive traction force {\tau }^{\mathrm{c}} on the crack face is a function of the crack-opening displacement in both normal and tangential directions. For the sake of simplicity, only the normal traction mode is considered in this study, and the shear traction is neglected as addressed in Refs. [²⁸,⁴²]. It is worth stressing out again here that this choice is not an intrinsic limitation of the present method. The mixed-mode cohesive law [⁴³,⁴⁴] can be easily implemented in such a framework.

The linear cohesive law [⁴⁵], which is schematically sketched in Fig. 2, is adopted here for this analysis. Herein, the normal cohesive traction {\tau }^{\mathrm{c}} is thus expressed by

(9)\{\begin{array}{ll}{\tau }^{\mathrm{c}}=f_t\left(1-{\displaystyle \frac{w}{w_0}}\right),& 0\le w\le w_0,\\ {\tau }^{\mathrm{c}}\mathrm{=0,}& w>w_0,\end{array}

where f_{\mathrm{t}} is the tensile strength; w and

w_0

are the crack-opening displacement and the critical crack-opening displacement in the normal direction, respectively. The cohesive fracture energy is then defined by G={\displaystyle \frac{1}{2}}{f}_t{w}_0.

A new XFEM scheme for cohesive crack growth

Displacement approximation

A 4-node quadrilateral element shown in Fig. 3 containing a crack-tip cutting edge 1-4 is considered. The main objective of the proposed approximation is to construct an enriched field that disappears at the edges 1-2, 2-3, and 3-4, and that continuous at the edge 1-4.

An auxiliary line is plotted through crack-tip point p, the intersection points between the line and the element boundaries are points 5 and 6. Thus, the crack-tip enriched element is element 1-5-6-4, the nodes 1 and 4 are enriched. The discontinuous displacement field of the original crack-tip element is thus expressed as

(10)u_{\mathrm{disc}}=b_1{N}_{\mathrm{1}}(x^{\ast })(H(x^{\ast })-H(x_1))+b_4{N}_{\mathrm{4}}(x^{\ast })(H(x^{\ast })-H(x_4)),

where N_i(x^{*}) are the standard nodal shape function value at x^{\ast } which is located in the element 1-5-6-4, as shown in Fig. 3(a); the local coordinates of nodes 6 and 5 in the isoparametric element 1-2-3-4 are shown in Fig. 3(b); H(x)is a modified Heaviside step function which takes on the value+ 1 above the crack and -1 below the crack.

Similarly, if the crack cuts edge 1-2, the discontinuous displacement field of the original crack tip element is described as

(11)u_{\mathrm{disc}}=b_1{N}_{\mathrm{1}}(x^{\ast })(H(x^{\ast })-H(x_1))+b_2{N}_{\mathrm{2}}(x^{\ast })(H(x^{\ast })-H(x_2)).

The enriched displacement for the crack-tip is expressed as

(12)u_{\mathrm{disc}}={\displaystyle \underset{i\in N^k}{\mathrm{\Sigma }}{b}_i{N}_i(x^{\ast })(H(x^{\ast })-H(x_i))},

where N^k is the node set on the edge of crack tip element cut by the crack.

In XFEM, the displacement field for a cracked element is the sum of the continuous and the discontinuous displacement fields. Hence, the displacement approximation of XFEM with 4-node quadrilateral element for cohesive crack model is written as

(13)u^h(x)={\displaystyle \underset{i\in N^s}{\mathrm{\Sigma }}{u}_i}{N}_i(x)+{\displaystyle \underset{i\in N^{\mathrm{cut}}\cup N^k}{\mathrm{\Sigma }}{a}_i}{N}_i(x)(H(x)-H(x_i)),

where N^s is the set of all nodes in the discretization, N^{\mathrm{cut}} is the set of nodes whose basis function support is entirely split by the crack; u_i and a_iare the vector of nodal degrees of freedom defined in standard finite elements and the vector of nodal enrichment variable, respectively.

From Eq. (13), it can be found that a unified displacement enrichment function approximation is used in the elements entirely split by the crack and the crack-tip element, it provides convenience for programming. The crack-opening displacement w can be obtained with the following equation.

(14)w={\displaystyle \underset{i\in N^{\mathrm{cut}}\cup N^k}{\mathrm{\Sigma }}\mathrm{2}{a}_i}{N}_i(x).

Discretized governing equations

Applying the XFEM discretization described in subsection 3.1 to Eq. (8), implying the following system of equations

(15)K\delta =R,

where K and R are the global stiffness matrix and external nodal force vector, respectively; \delta is the vector of nodal unknowns.

The local stiffness matrix Kfor each element is given by

(16k_{ij}=\left[\begin{array}{ll}{k}_{ij}^{uu}& k_{ij}^{ua}\\ k_{ij}^{au}& k_{ij}^{aa}\end{array}\right],

and

(17)k_{ij}^{rs}={\displaystyle {\int }_{{\mathrm{\Omega }}^e}{(B_i^r)}^{\mathrm{T}}C{B}_j^s\mathrm{d}\mathrm{\Omega }}\begin{array}{ll}& \end{array}(r,s=u,a),

in which the derivative of shape function is defined as

(18a)B_i^u=\left[\begin{array}{ll}{N}_{i,x}& 0\\ 0& N_{i,y}\\ N_{i,y}& N_{i,x}\end{array}\right],

(18b)B_i^a=(H(x)-H(x_i))\left[\begin{array}{ll}{N}_{i,x}& 0\\ 0& N_{i,y}\\ N_{i,y}& N_{i,x}\end{array}\right].

The external nodal force R for each element is formulated by

(19)r_i^{}={\left[\begin{array}{ll}{r}_i^u& r_i^a\end{array}\right]}^{\mathrm{T}},

and

(20a)r_i^u={\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{N}_i\lambda {\tau }_{\mathrm{0}}\mathrm{d}\mathrm{\Gamma }},

(20b)r_i^a={\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{t}}}{\int }{N}_i(H(x)-H(x_i))\lambda {\tau }_{\mathrm{0}}\mathrm{d}\mathrm{\Gamma }}+{\displaystyle \underset{{\mathrm{\Gamma }}_{\mathrm{coh}}}{\int }\mathrm{2}{N}_i\tau \mathrm{d}\mathrm{\Gamma }}.

In practice, the crack-opening displacement w is determined through Eq. (14), then the crack opening displacement (COD) in the normal direction \mathit{w} is computed through coordinate transformation. In the case that \mathrm{0}\le w\le w_c at all integral points on the crack surface, the crack segment is here considered to be located in the fracture process zone or the cohesive zone. Herein the cohesive law (depending onw) is used to describe the decohesion process. On the contrary, if the normal crack opening at any integral point verifiesw\ge w_c, the crack segment is considered to be located in the real crack section or fully cracked area.

Multiscale mesh and variable-node elements

It has been well-known that the accuracy of XFEM analysis can be significantly improved by using a fine-scale mesh, but it obviously requires more computational efforts. A possible technique can be adopted to avoid this problem is to use a multiscale mesh, by which only the critical region (expected cracking) is discretized by the fine-scale mesh, whereas coarse-scale mesh is for the rest of the body. A major difficulty induced by this technique is how to effectively/accurately connect two different mesh-scales. As represented by the crosshatched lines in Fig. 4, there exists one layer of variable-node elements between two different scale elements.

In this study, the variable-node elements [⁴¹] are adopted, which enable us to connect the nonconforming meshes between two zones. The shape functions of variable-node elements are developed based on the concept of generic point interpolation with special bases that have slope discontinuities in 2D domains [⁴]. Within this framework, the approximated displacements u^h(\xi ) are expressed as

(21)u^h(\xi )={\displaystyle \underset{i=1}{\overset{N_p}{\mathrm{\Sigma }}}{N}_i(\xi )u_i}=a^{\mathrm{T}}p(\xi ),

where N_p is the number of sample points in the point interpolation; N_i is the shape function matrix of the ith node; u_i is the nodal variable vector; a^{\mathrm{T}} is the matrix of the unknown coefficients; and p(\xi ) is the column vector of the polynomial basis.

Figure 5 shows a so-called (4+ k + m)-node element, where k is the number of extra nodes on the top and bottom edges of element, while m is the number of extra nodes on the left and right edges. The polynomial basis following [⁴] is given as

(22)p(\xi )={[\mathrm{1},\xi ,\eta ,\xi \eta ,\left|\xi -{{\displaystyle \xi }}_{\mathrm{5}}\right|(\eta +\sin ({{\displaystyle \xi }}_{\mathrm{5}})),\mathrm{...},\left|\xi -{{\displaystyle \xi }}_{\mathrm{4}+k}\right|(\eta +\mathrm{sign}({{\displaystyle \xi }}_{\mathrm{4}+k})),\left|\eta -{{\displaystyle \eta }}_{\mathrm{4}+k+\mathrm{1}}\right|(\xi +\mathrm{sign}({{\displaystyle \xi }}_{\mathrm{4}+k})),\mathrm{...},\left|\eta -{{\displaystyle \eta }}_{\mathrm{4}+k+m}\right|(\xi +\mathrm{sign}({{\displaystyle \xi }}_{\mathrm{4}+k+m}))]}^T,

where \xi and \eta describe the local coordinates in the isoparametric element.

The point interpolation is given by

(23)u^{\mathrm{h}}(\xi )=a^{\mathrm{T}}p(\mathrm{\xi })=U^{\mathrm{T}}{q}^{\mathrm{-1}}p(\mathrm{\xi }),

which defines the shape functions of the (4+ k + m)-node element as

(24){[N_1,\cdots ,N_{4+k+m}]}^T=q^{\mathrm{-1}}p(\xi ),

with

(25a)q=[p({\xi }_1),\mathrm{...},p({\xi }_{4+k+m})],

(25b)\begin{array}{l}p({\xi }_1)=[\mathrm{1},{\xi }_{\mathrm{1}},{\eta }_{\mathrm{1}},{\xi }_{\mathrm{1}}{\eta }_{\mathrm{1}},\left|{\xi }_{\mathrm{1}}-{\xi }_5\right|({\eta }_1+\mathrm{sign}({\xi }_5)),\mathrm{...},\left|{\xi }_1-{\xi }_{\mathrm{4+}k}\right|({\eta }_1+\mathrm{sign}({\xi }_{\mathrm{4+}k})),\\ {\left|{\eta }_1-{\eta }_{\mathrm{4+}k+1}\right|({\xi }_1+\mathrm{sign}({\eta }_{\mathrm{4+}k+1})),\mathrm{...},\left|{\eta }_1-{\eta }_{\mathrm{4+}k+m}\right|({\xi }_1+\mathrm{sign}({\eta }_{\mathrm{4+}k+m}))]}^T,\end{array}

(25c)\begin{array}{l}p({\xi }_{4+k+m})=[\mathrm{1},{\xi }_{\mathrm{4}+k+m},{\eta }_{\mathrm{4}+k+m},{\xi }_{\mathrm{4}+k+m}{\eta }_{\mathrm{4}+k+m},\left|{\xi }_{4+k+m}-{\xi }_5\right|({\eta }_{4+k+m}+\mathrm{sign}({\xi }_5)),\mathrm{...},\left|{\xi }_{4+k+m}-{\xi }_{4+k}\right|({\eta }_{4+k+m}+\mathrm{sign}({\xi }_{4+k})),\\ {\left|{\eta }_{4+k+m}-{\eta }_{\mathrm{4+}k+1}\right|({\xi }_{4+k+m}+\mathrm{sign}({\eta }_{\mathrm{4+}k+1})),\mathrm{...},\left|{\eta }_{4+k+m}-{\eta }_{\mathrm{4+}k+m}\right|({\xi }_{4+k+m}+\mathrm{sign}({\eta }_{\mathrm{4+}k+m}))]}^T,\end{array}

(26)U^T=\left[\begin{array}{lllll}{u}_1& \mathrm{...}& u_i& \mathrm{...}& u_{4+k+m}\\ v_1& \mathrm{...}& v_i& \mathrm{...}& v_{4+k+m}\end{array}\right],

The variable-node elements retain the linear interpolation between any two neighboring nodes, so they can naturally connect different-scale elements.

Crack growth modeling

The crack growth modeling is governed by the maximum normal tensile stress criterion. This theory states that the crack propagation will take place when the projection of the fictitious crack-tip stress tensor in the normal direction reaches the tensile strength of the material [²⁹], as follows

(27)n_2\cdot C\cdot B(x_{\mathrm{t}})\delta =S\cdot \delta =f_{\mathrm{t}},

with

(28)n_2=[l^{\mathrm{2}} m^{\mathrm{2}} \mathrm{2}lm],

(29)S=n_2\cdot C\cdot B(x_{\mathrm{t}}),

where l and m are components of the normal vector n at the fictitious crack-tip in the local coordinate system; x_t is the location of the fictitious crack-tip.

Furthermore, the crack growth direction is determined by using the maximum hoop stress criterion [²] by

(30)\theta =\mathrm{2}\arctan {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(K_{{\rm I}}^{\mathrm{ext}}/K_{\Pi }^{\mathrm{ext}}\pm \sqrt{{(K_{{\rm I}}^{\mathrm{ext}}/K_{\Pi }^{\mathrm{ext}})}^{\mathrm{2}}+\mathrm{8}}),

in which K_{{\rm I}}^{\mathrm{ext}} and K_{\Pi }^{\mathrm{ext}} are the stress intensity factors (SIFs) in mode I and mode II under the external force, respectively.

Solving the nonlinear governing equations

Due to the nonlinearity of the cohesive traction law, the discrete governing equation described in Eq. (15) is nonlinear. An iterative method is therefore required to solve the aforementioned nonlinear-problem. In this study, the arc length method, similar to one proposed in Ref. [²⁹], is adopted to that aim. The load factor \lambda and the nodal variables \delta are considered as unknowns. The unbalanced force \tilde{R}(\delta ,\lambda ) is formulated by

(31)\tilde{R}(\delta ,\lambda )=K\delta -\lambda f_0^{\mathrm{ext}}-f^{\mathrm{coh}}(\delta ),

where f_0^{\mathrm{ext}} and f^{\mathrm{coh}}(\delta ) are the equivalent nodal load vectors for the external tractions {\tau }_{\mathrm{0}} and for the cohesive tractions {\tau }^c, respectively.

Using a Taylor expansion at first order of the unbalanced force \tilde{R}(\delta ,\lambda ) at the iteration ith we can express

(32)\begin{array}{l}\tilde{R}={\tilde{R}}^i+{\left({\displaystyle \frac{\partial \tilde{R}}{\partial \delta }}\right)}^i\Delta {\delta }^i+{\displaystyle \frac{\partial \tilde{R}}{\partial \lambda }}\Delta {\lambda }^i={\tilde{R}}^i+{\left(K-{\displaystyle \frac{\partial f^{\mathrm{coh}}(\delta )}{\partial \delta }}\right)}^i\Delta {\delta }^i-f_0^{\mathrm{ext}}\Delta {\lambda }^i\\ = {\tilde{R}}^i+K_T^i\Delta {\delta }^i-f_0^{\mathrm{ext}}\Delta {\lambda }^i.\end{array}

Setting \tilde{R}=\mathrm{0}, it implies

(33)\mathrm{\Delta }{\delta }^i=-{(K_T^i)}_{}^{-\mathrm{1}}{\tilde{R}}^i+{(K_T^i)}_{}^{-\mathrm{1}}{f}_0^{\mathrm{ext}}\mathrm{\Delta }{\lambda }^i=\mathrm{\Delta }{\overline{\delta }}^i+{\delta }_e^i\mathrm{\Delta }{\lambda }^i.

Substituting Eq. (33) into {\delta }^{i+\mathrm{1}}={\delta }^i+\Delta {\delta }^i and along with Eq. (28), yield

(34)S\cdot \left({\delta }^i+\mathrm{\Delta }{\overline{\delta }}^i+{\delta }_e^i\mathrm{\Delta }{\lambda }^i\right)=f_{\mathrm{t}}.

The incremental loading factor is here defined

(35)\mathrm{\Delta }{\lambda }^i={\displaystyle \frac{f_t-S\cdot \left({\delta }^i+\mathrm{\Delta }{\overline{\delta }}^i\right)}{S\cdot {\delta }_e^i}}.

In summary, the incremental unknowns are deduced by Eqs. (33), (35). Then the loading and displacement at the (i + 1)th iteration are given by {\lambda }^{i+\mathrm{1}}={\lambda }^i+\Delta {\lambda }^iand {\delta }^{i+\mathrm{1}}={\delta }^i+\Delta {\delta }^i.

When the ratio between the norm of current unbalanced force and that of the previous step is less than a certain value, it is considered that the equilibrium condition is satisfied, i.e., the solution is convergent. In the numerical examples, the value is adopted as 1%.

Major steps of the numerical implementation

The main solution procedure of modeling cohesive crack growth using the variable-node XFEM with local refined mesh is summarized as follows:

1) Give the information on the mesh, crack, external load, and boundary condition, etc;

2) Solve the displacement field under the external load;

3) Compute the SIFs and the crack growth direction;

4) The fictitious crack grows a given length in the crack growth direction;

5) Solve the nonlinear governing equations and simultaneously obtain the FPZ size;

6) Compute the SIFs for the current fictitious crack-tip under the current external load and the new crack growth direction;

7) Go to step 4 and repeat the computation for the next step.

Numerical examples

In this section, three numerical examples, including one mode-I fracture and two mixed mode fracture problems, are presented to demonstrate the accuracy and performance of the proposed approach. Note that, in all numerical experiments, the plane strain condition is assumed.

Three-point bending test of a simply supported beam with mode-I fracture

A simply supported beam containing an initial crack under three-point bending as schematically sketched in Fig. 6 is investigated. The dimensions of the considered beam are L× H× B = 2000 mm × 500 mm × 200 mm, and the initial crack length a₀ = 270 mm. The material properties of concrete are chosen as follows [⁴⁶], with Young’s modulus E = 24.94 GPa, Poisson’s ratio \nu = 0.2, tensile strength f_t = 1.43 N/mm², and the cohesive fracture energy G = 102.22 N/m. The mesh sensitivity is analyzed by performing the numerical simulations on three different meshes as depicted in Fig. 7, showing a coarse scale mesh, a fine scale mesh, and a multiscale mesh.

The load-CMOD (crack mouth opening displacement) curves for different meshes are represented in Fig. 8. Herein, we also provide the reference result obtained from the experiment [⁴⁶] for our validation purpose. A quantitative comparison on the computational time and the correlation with reference data for different cases is reported in Table 1. In more detail, the correlation between the model prediction using the coarse scale mesh and the reference data are about 5.44% with the computational time 181.8 s. A much better agreement with literature observation is demonstrated for the fine scale mesh, where the difference is 1.50%. However, as said, the fine scale mesh task obviously requires much more computational cost, and with that the calculation time is estimated by 343.6 s. More interestingly, the multiscale mesh with the computational time is only 260.6 s but it still offers good accuracy, with the difference in comparison with the reference data are 2.26%. The obtained results demonstrate one important issue that a balance between the accuracy and efficiency could be gained on the sense that the multiscale mesh holds. In particular, we here can save 24.2% of the computational time, while the accuracy is still ensured.

The crack propagation of the model using multiscale mesh for four different loading forces is thus represented in Fig. 9. A post-processing is adopted to highlight the crack opening by multiply the displacement by a scale factor. The mode-I fracture is clearly captured, demonstrating the performance of the proposed technique.

Three-point bending test with mixed-mode fracture

Next, a pre-cracked concrete beam with dimension D = 80 mm as depicted in Fig. 10 subjected to a mixed mode fracture test is considered. The initial crack length is a = 20 mm, and its distance with the prescribed force (at the middle) is e = 25 mm. The material properties are chosen as follows: Young’s modulus E = 38 GPa, Poisson’s ratio\nu = 0.2, tensile strength f_t = 3.5 N/mm², and cohesive fracture energy G= 80 N/m. As usual, the simulation is also performed on three different meshes, i.e., coarse scale mesh, fine scale mesh, and multiscale mesh, as shown in Fig. 11.

The result of crack propagation for the multiscale mesh at different loading is represented in Fig. 12, which shows a typical mixed-mode fracture behavior. The final crack path obtained from the simulation for all three meshes, along with the reference data [⁴⁷], and experimental observation in Ref. [⁴⁷] is visualized in Fig. 13. A good agreement of the numerical prediction derived from the present model and literature reference results is obtained. In particular, the crack trajectory obtained by fine mesh and multiscale mesh shows a high correlation with the reference result reported in Ref. [⁴⁷]. Consequently, that results exhibit the accuracy of the proposed model.

The results of load-CMOD curves obtained from the numerical simulation on three different meshes are depicted in Fig. 14. The comparison of the numerical performance (computational time, accuracy), among different cases is provided in the Table 2, in which the accuracy of the present framework is estimated by comparing with the reference data in Ref. [⁴⁷].

Once again, the robustness of the multiscale mesh is proved. An acceptable error 1.78% is obtained, with the saving of computational time is about 26.78% compared to the fine scale mesh.

To show the performance of the developed method using irregular mesh, two irregular meshes as shown in Fig. 15 are additionally considered. The computed results of the load-CMOD curves obtained by the numerical simulation for two different meshes are thus depicted in Fig. 16. It can be observed from Fig. 16 that the present approach using irregular meshes also offers acceptable solutions as compared with the regular ones. In other words, the obtained results by the present approach are insensitive to the element shape.

A four-point shear specimen with mixed mode fracture

The last numerical example deals with a more complex fracture behavior. More specifically, a pre-cracked concrete beam under the mixed mode fracture is analyzed using the developed method. The details of geometry and boundary condition are depicted in Fig. 17 (length unit is in m, and force unit is in kN). An initial crack length is typically chosen a = 82.4 mm. The material properties are taken as follows: Young’s modulus E = 24.8 GPa, Poisson’s ratio\nu = 0.18, tensile strength f_t = 4 MPa, and cohesive fracture energy G = 125 N/m. Three kinds of mesh mentioned above, as shown in Fig. 18, are also analyzed.

The validation of the final crack paths obtained from different meshes with the reference result reported in Ref. [⁴⁷] is shown in Fig. 19. The complex crack trajectories are reproduced by the present model. The convergence to the reference solution is obtained when refining mesh. Moreover, to quantitatively assess the accuracy of the proposed model, we present the comparison of the mechanical response for different cases in Fig. 20, showing a high correlation. All of that have demonstrated the performance of the developed method, which can ensure the accuracy with the significant improvement in the computational cost.

Conclusions and outlook

In this paper, we have presented a new computational approach based on the combination of the XFEM with local refined mesh and the cohesive zone model. The cracked elements are enriched by the sign function. The blending element problem is treated by the shifted sign enrichment function. Hence, the equal stresses at both sides of the crack in the crack-tip are effectively modeled. The multiscale mesh using variable-node elements along with the arc-length technique is adopted to further enhance the efficiency of the proposed model.

The new model is applied to study both mode-I and mixed-mode fractures in concrete material. The very promising result is obtained. As numerically demonstration, it can be found that by using a multiscale mesh, a balance between accuracy and efficiency can be obtained. In other words, the computational time is significantly reduced, while the accuracy is still acceptable. In addition, the proposed model is extremely robust and easy to implement, which provides a new performance computational tool to investigate the fracture behavior of materials.

A potential extension of the proposed framework in the future would be devoted to adopt the mixed-mode cohesive law and then apply to study the crack propagation in highly heterogeneous material, like cement-based materials at the microscopic scale. The proposed methodology has the potential to solve crack nucleation and 3D structures, which have been the direction for our future research.

In the present study, the size of the refined zone (i.e., the zone of interest) is artificially defined. A more scientific and reasonable method is to determine the fine-scale meshes by a posteriori error estimation [⁴]. In addition, the robustness of the developed methodology need to be tested by applying the present method to study complex crack growth like out of plane crack propagation, crack growth in composite structures.

References

[1] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[2] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[3] Yi G, Yu T T, Bui T Q, Tanaka S. Bi-material V-notched SIFs analysis by XFEM and conservative integral approach. Computers & Structures, 2018, 196: 217–232

[4] Yu T T, Bui T Q. Numerical simulation of 2-D weak and strong discontinuities by a novel approach based on XFEM with local mesh refinement. Computers & Structures, 2018, 196: 112–133

[5] Ma C P, Yu T T, Lich L V, Bui T Q. An effective computational approach based on XFEM and a novel three-step detection algorithm for multiple complex flaw clusters. Computers & Structures, 2017, 193: 207–225

[6] Yi G, Yu T T, Bui T Q, Ma C P, Hirose S. SIFs evaluation of sharp V-notched fracture by XFEM and strain energy approach. Theoretical and Applied Fracture Mechanics, 2017, 89: 35–44

[7] Wang Z, Yu T T, Bui T Q, Tanaka S, Zhang C, Hirose S, Curiel-Sosa J L. 3-D local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 375–405

[8] Wang Z, Yu T T, Bui T Q, Trinh N A, Luong N T H, Duc N D, Doan D H. Numerical modeling of 3-D inclusions and voids by a novel adaptive XFEM. Advances in Engineering Software, 2016, 102: 105–122

[9] Yu T T, Bui T Q, Yin S H, Doan D H, Wu C T, Do T V, Tanaka S. On the thermal buckling analysis of functionally graded plates with internal defects using extended isogeometric analysis. Composite Structures, 2016, 136: 684–695

[10] Yu T T, Bui T Q, Liu P, Zhang C, Hirose S. Interfacial dynamic impermeable cracks analysis in dissimilar piezoelectric materials under coupled electromechanical loading with the extended finite element method. International Journal of Solids and Structures, 2015, 67–68: 205–218

[11] Shi L Y, Yu T T, Bui T Q. Numerical modeling of hydraulic fracturing in rock mass by XFEM. Soil Mechanics and Foundation Engineering, 2015, 52(2): 74–83

[12] Liu P, Bui T Q, Zhu D, Yu T T, Wang J W, Yin S H, Hirose S. Buckling failure analysis of cracked functionally graded plates by a stabilized discrete shear gap extended 3-node triangular plate element. Composites. Part B, Engineering, 2015, 77: 179–193

[13] Liu P, Yu T T, Bui T Q, Zhang C, Xu Y P, Lim C W. Transient thermal shock fracture analysis of functionally graded piezoelectric materials by the extended finite element method. International Journal of Solids and Structures, 2014, 51(11–12): 2167–2182

[14] Yu T T. The extended finite element method (XFEM) for discontinuous rock masses. Engineering Computations, 2011, 28(3): 340–369

[15] Liu S, Yu T, Van Lich L, Yin S, Bui T Q. Size effect on cracked functional composite micro-plates by an XIGA-based effective approach. Meccanica, 2018, 53(10): 2637–2658

[16] Bažant Z P, Planas J. Fracture and Size Effect in Concrete and Other Quasibrittle Materials. Boca Raton: CRC Press, 1998

[17] Dugdale D S. Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids, 1960, 8(2): 100–104

[18] Barenblatt G I. The mathematical theory of equilibrium of cracks in brittle fracture. Advances in Applied Mechanics, 1962, 7: 55–129

[19] Hillerborg A, Modéer M, Petersson P E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6(6): 773–781

[20] Barpi F, Valente S. The cohesive frictional crack model applied to the analysis of the dam-foundation joint. Engineering Fracture Mechanics, 2010, 77(11): 2182–2191

[21] Yang B, Ravi-Chandar K. A single-domain dual-boundary-element formulation incorporating a cohesive zone model for elastostatic cracks. International Journal of Fracture, 1998, 93(1/4): 115–144

[22] Belytschko T, Organ D, Gerlach C. Element-free galerkin methods for dynamic fracture in concrete. Computer Methods in Applied Mechanics and Engineering, 2000, 187(3–4): 385–399

[23] Sancho J M, Planas J, Fathy A M, Gálvez J C, Cendón D A. Three-dimensional simulation of concrete fracture using embedded crack elements without enforcing crack path continuity. International Journal for Numerical and Analitical Methods in Geomechanics, 2007, 31(2): 173–187

[24] Wells G N, Sluys L J. A new method for modeling cohesive cracks using finite elements. International Journal for Numerical Methods in Engineering, 2001, 50(12): 2667–2682

[25] Asferg J L, Poulsen P N, Nielsen L O. A consistent partly cracked XFEM element for cohesive crack growth. International Journal for Numerical Methods in Engineering, 2007, 72(4): 464–485

[26] Zi G, Belytschko T. New crack-tip elements for XFEM and applications to cohesive cracks. International Journal for Numerical Methods in Engineering, 2003, 57(15): 2221–2240

[27] Remmers J J C, Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. Computational Mechanics, 2003, 31(1–2): 69–77

[28] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics, 2002, 69(7): 813–833

[29] Zhang X D, Bui T Q. A fictitious crack XFEM with two new solution algorithms for cohesive crack growth modeling in concrete structures. Engineering Computations, 2015, 32(2): 473–497

[30] Mougaard J F, Poulsen P N, Nielsen L O. A partly and fully cracked triangular XFEM element for modeling cohesive fracture. International Journal for Numerical Methods in Engineering, 2011, 85(13): 1667–1686

[31] Cox J V. An extended finite element method with analytical enrichment for cohesive crack modeling. International Journal for Numerical Methods in Engineering, 2009, 78(1): 48–83

[32] Wu J Y, Li F B. An improved stable XFEM (Is-XFEM) with a novel enrichment function for the computational modeling of cohesive cracks. Computer Methods in Applied Mechanics and Engineering, 2015, 295: 77–107

[33] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[34] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[35] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[36] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[37] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[38] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[39] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23–24): 1391–1411

[40] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92–93: 242–256

[41] Lim J H, Sohn D, Im S. Variable-node element families for mesh connection and adaptive mesh computation. Structural Engineering and Mechanics, 2012, 43(3): 349–370

[42] Roth S N, Léger P, Soulaïmani A. A combined XFEM-damage mechanics approach for concrete crack propagation. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 923–955

[43] Xu Y, Yuan H. Computational modeling of mixed-mode fatigue crack growth using extended finite element methods. International Journal of Fracture, 2009, 159(2): 151–165

[44] Liao F, Huang Z. An extended finite element model for modelling localised fracture of reinforced concrete beams in fire. Computers & Structures, 2015, 152: 11–26

[45] Yu T T. The Extended Finite Element Method-Theory, Application and Program. Beijing: Science Press, 2014

[46] Zhang X F, Xu S L. Determination of fracture energy of three-point bending concrete beam using relationship between load and crack-mouth opening displacement. Journal of Hydraulic Engineering, 2008, 39(6): 714–719

[47] Cendón D A, Gálvez J C, Elices M, Planas J. Modelling the fracture of concrete under mixed loading. International Journal of Fracture, 2000, 103(3): 293–310

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature