#### Abstract

In the frame of the extended finite element method, the exponent disconnected function is introduced to reflect the discontinuous characteristic of crack and the crack tip enrichment function which is made of triangular basis function, and the linear polar radius function is adopted to describe the displacement field distribution of elastoplastic crack tip. Where, the linear polar radius function form is chosen to decrease the singularity characteristic induced by the plastic yield zone of crack tip, and the triangle basis function form is adopted to describe the displacement distribution character with the polar angle of crack tip. Based on the displacement model containing the above enrichment displacement function, the increment iterative form of elastoplastic extended finite element method is deduced by virtual work principle. For nonuniform hardening material such as concrete, in order to avoid the nonsymmetry characteristic of stiffness matrix induced by the non-associate flowing of plastic strain, the plastic flowing rule containing cross item based on the least energy dissipation principle is adopted. Finally, some numerical examples show that the elastoplastic X-FEM constructed in this paper is of validity.

#### 1. Introduction

The simulation of crack propagation has always been a hot problem of academia. For finite element method, to implement the process simulation of tracking the cracks propagations, the structure containing cracks must be remeshed with the extension of cracks. So at early time, compared with the finite element method, the boundary element method is more broadly used to track the extension of crack. Since the basic ideas and the mathematical foundation of partition of unity finite element method (PUFEM) were discussed by Melenk and Babuska [1] and Duarte and Oden [2], the extended finite element method based on the idea of partition of unity has been gradually constructed and modified. Firstly, Belytschko and Black [3] presented a minimal remeshing finite element method by adding discontinuous enrichment functions to the finite element approximation to account for the presence of a crack. Then, Moës et al. [4] and Dolbow [5] further improved the method and called it the extended finite element method (X-FEM). The core idea of partition of unity in X-FEM mainly embodies its unique displacement model. The displacement model of X-FEM specifically is constructed by introducing the enrichment displacement functions reflecting discontinuous characteristic of crack and singularity at crack tip on the basis of the finite element displacement model. So, the X-FEM not only avoids the complicated mesh regeneration that must adapt the cracks after extension but also inherits the merits of finite element method such as sparse band stiffness matrix. Apparently, the X-FEM is of more broad adaptability compared with the boundary element method. This huge superiority without re-mesh in tracking the crack extension has made the research of X-FEM boom recently, from two-dimensional crack extension problem [6–8] developing to three-dimensional crack extension problem [9], from elastic fracture problem developing to elastoplastic crack extension problem [10], from discontinuous problem developing to weak continuous problem [11] and crack contact problem [12, 13], from crack static extension problem developing to crack dynamics extension problem [14, 15], from crack extension problem developing to interface crack extension problem [16], and from single crack extension problem developing to multicrack problem including branch crack problem [17, 18] and more complicated crack problem [19, 20]. What is more, in order to improve the calculation accuracy further, Qing et al. developed a new X-FEM by incorporating the high precision merit of generalized finite element method—the generalized extended finite element method (GX-FEM) [21].

For elastic-plastic crack extension problem, the crack tip stress field is no longer of singularity because of the production and development of plastic zone. So the enrichment function of crack tip cannot be directly applied to the elastic-plastic crack tip. Elguedj et al. [10] use the triangle basis function based on the Taylor series expansion of power function distribution as the enrichment function of elastoplastic crack tip for Ramberg-Osgood power law hardening material. However, for some other plastic hardening material such as concrete material, the plastic strain distribution of crack tip is difficult to be grasped, especially for complicated stress status. In fact, the plastic zone of crack tip was produced at the crack tip firstly and developed gradually toward outside around induced by the transfer of stress concentration. Therefore, the triangular basis function part with polar angle variable of elastic crack tip can be still used to construct the displacement enrichment function of elastoplastic crack tip, and only the extraction distribution function part with the polar radius reflecting singularity of crack tip needs to be modified to make the singularity disappear. In this paper, the linear distribution function with the polar radius and triangular basis function with polar angle are chosen to construct the displacement enrichment function of elastoplastic crack tip.

The rest of the paper is structured as follows. Section 2 describes the displacement model of X-FEM for elastoplastic crack extension problem. And the strain matrix of elastoplastic X-FEM is presented in Section 3 in detail. In Section 4, the derivation process of the increment finite element iterative form for X-FEM based on the virtual work principle will be presented in detail, where the plastic flowing rule containing cross item based on the least energy dissipation principle and the elastoplastic extension criterion of crack are also introduced. Some numerical examples will be presented in Section 5 to verify the validity of the elastoplastic X-FEM proposed in this paper. Finally, we close with concluding remarks in Section 6.

#### 2. The Displacement Model of Elastoplastic X-FEM

For any element, the common form of the displacement model of X-FEM can be presented as follows: where is the set of all the nodes of an element, the master degree of node , the interpolation shape function of node , the set of the element nodes containing discontinuous enrichment degree, the discontinuous enrichment displacement function of node , the corresponding enrichment degree of node presented as the hollow points in Figure 1, the set of the element nodes containing crack tip enrichment degree, the crack tip enrichment displacement function of node , and the corresponding crack tip enrichment degree of node presented as the solid points in Figure 1.

For discontinuous enrichment displacement function (see Figure 3), considering the fast decline influence with the increment of perpendicular distance to crack, the exponent disconnected function form presented as follows is adopted [22]: where is defined as Here, is the curve equation of crack shown in Figure 2.

For the construction of elastoplastic crack tip enrichment displacement function, to decrease the singularity characteristic induced by the plastic yield zone of crack tip, the linear function with polar radius form is chosen, and to describe the displacement distribution character with the polar angle of crack tip, the triangular basis function form derived from the close solution of elastic crack tip for infinite big plate problem is still adopted. So the elastoplastic crack tip enrichment displacement function can be written as follows by combining the linear function polar radius and triangular basis function with polar angle. Formula (4a) is for two-dimensional crack problem and formula (4b) for three-dimensional crack problem aswhere , , , , , and is the transform matrix from local coordinate system established at the crack face fringe to global coordinate.

#### 3. Strain Matrix of X-FEM

Supposing that the problem considered is a small deformation problem, the strain of arbitrary a point can be written as follows in terms of geometry equation: where is the symmetry part of gradient operator. Substituting the displacement model of elastoplastic X-FEM for formula (5), the strain can be gained as follows: where the strain matrixes , , and are derived from common shape function of FEM, discontinuous enrichment function, and crack tip enrichment function, respectively. For three-dimensional problem, the substrain matrixes are written, respectively, as follows:

#### 4. The Increment Iterative Form of Elastic-Plastic X-FEM

##### 4.1. Governing Equation of X-FEM

The equilibrium equation and the boundary condition can be written as follows: where is the body force, the calculation domain, the displacement boundary at time, the force boundary at time, the crack boundary, and and the corresponding outside normal line vector of crack surfaces and , respectively. The sketch of all kinds of boundary conditions of X-FEM at time can be denoted as shown in Figure 4.

##### 4.2. The Plastic Flowing Rule for Kinetic Hardening Material

To avoid the nonsymmetry characteristic of stiffness matrix induced by nonassociated flowing rule presented as formula (9) for kinetic hardening material, the plastic flowing rule [23] containing cross item based on the least energy dissipation principle presented as following formula (10) is adopted. where denotes the th potential function, and denotes the th yield function. For kinetic hardening material, . Here, is the hardening parameters, and is the yield strength parameters as where the first item of equality right means the cross item with tensor product of stress tensor and plastic softness matrix. If the cross item vanishes, the plastic flowing rule is just the associated flowing rule. The denotes the th yield surface. So the constitutive law of kinetic hardening material can be easily deduced as follows [23]:

##### 4.3. The Increment Iterative Form of Elastoplastic X-FEM

The increment displacement of X-FEM can be written easily as follows from formula (1): For crack open problem, the contact condition can be ignored. Therefore, the incremental integral weak form is easily gained as follows from the governing equation (8) and constitutive equation (11): The increment displacement is substituted by its shape function (12) to formula (13), and the integral domains are discretized, and formula (13) becomes as follows: Further combining the same type of increment displacement degree yields Formula (15) can be finally made by simplification as follows: where Here, Because of the randomicity of virtual displacement, the increment X-FEM iterative form based on the least energy dissipation principle can be written as follows:

##### 4.4. The Extension Criterion of Elastoplastic Crack

For the elastoplastic fracture problem, the energy release rate criterion is adopted. Figure 5 shows a planform of a surface shape crack being of virtual extension at the coordinate face. And the definition of the crack local coordinate system is also denoted as in Figure 5. Therefore, the formulas of the strain energy release rates corresponding to the three types of the fracture model are as follows: where , , and are the stress components of virtual crack extension face for actual crack front, and , , and are the open displacement components of virtual crack extension face for the back of the virtual crack front.

For three-dimensional extended finite element with 8 nodes, on the crack front 1 presented in Figure 6, the strain energy release rates of its three types of fracture model are as follows: where , , and are the nodal force of nonnode 1. , , and are the relative displacements between face crack midpoint 3 and midpoint 4 on local coordinate system.

The extension direction of crack for elastoplastic fracture problem can also be determined by the maximum hoop stress criterion [24] based on and . But the values of and should be gained by linear out-interpolating values of and of those points out of the plastic zone presented in Figure 7 as follows:

For the sake of clarity, the algorithm flowchart for elastoplastic X-FEM is presented in Figure 20.

#### 5. Numerical Examples

We analyze several examples to verify the validity of the proposed elastoplastic X-FEM.

##### 5.1. Concrete Plane Thin Plate Containing Crack

As shown in Figure 8, a concrete plane thin plate containing a horizontal crack with 0.4445 m length on its left side is subjected to a pair of uniform pressure loads on its top and bottom sides. The size of the plate is 2 m length and 3 m high. The strength parameter cohesive force is given as MPa and internal friction angular is given as . The Drucker-Prager yield criterion is adopted and the uniaxial compression stress-strain curve is given in Figure 9. Suppose that the concrete is a nonassociated flowing material. And the constitutive equation (11) is applied to the crack plate to avoid the production of unsymmetrical stiffness matrix. The elastoplastic X-FEM proposed above is adopted to simulate the crack elastoplastic extension. The results are shown in Figures 10 and 11. It can be seen from Figure 10 that the crack propagates along horizontal direction under the symmetrical loads, though the deformation and plastic zone of crack tip presented in Figure 11 take on some nonsymmetric distribution. The non-symmetrical distribution of deformation and plastic zone is induced by the nonassociated flowing characteristics of concrete material. Since the right outline of the plastic zone perpendicular to the crack, the horizontal propagation direction of crack gained is rational.

**(a)**

**(b)**

**(c)**

##### 5.2. Horizontal Crack on the Right End of Cantilever Beam

A cantilever beam with 10 m length, 1 m width, and 2 m high contains a horizontal crack on the right end of cantilever beam presented in Figure 12. And it is subjected to a pair of uniform tensile stress loads on the top and bottom surfaces of the right end shown in Figure 12. The strength parameters cohesive force MPa and internal friction . The Drucker-Prager yield criterion is adopted and the uniaxial compression stress-strain curve is given in Figure 9. The elastoplastic X-FEM established above is adopted to track the propagation of the horizontal crack face. If the effect of cross item in nonassociated flowing is ignored, the deformation and plastic zone take on the symmetric distribution presented in Figures 13–15 under the symmetric surface loads. Only the plastic zone reaches the saturated status, with the plastic zone moving along extension direction shown in Figures 14 and 15.

**(a) Preextension**

**(b) Postextension**

**(a) Unsaturated plastic zone**

**(b) Saturated plastic zone**

**(a) Preextension**

**(b) Postextension**

##### 5.3. Crack in a Plate with a Hole

The problem shown in Figure 16 is an adaptation of an example presented in [25]. The initial crack length is mm (all dimensions in the sketch of Figure 16 are given in mm), and the parameters are given as elastic modulus GPa, Possion ratio , and plate thickness 16 mm. The load applied to the extended finite element analysis is kN. The extension resistance strength of the aluminum plate is MPa. The Mises yield criterion is adopted and the uniaxial tension stress-strain curve is given in Figure 17. The analysis of quasi-static crack propagation has been carried out with 10 crack increments of mm in each step. The elastoplastic X-FEM established above is adopted to track the propagation of the horizontal crack on the left edge of plate. Since the influence of the hole, the crack propagates toward the hole step by step as presented in Figure 18 which agrees well with the experimental observations shown in Figure 16. What is more, owning to the change of extension direction, the shape of plastic zone of crack tip changes a lot during the crack propagation process presented in Figure 19.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

#### 6. Conclusions

In this paper, we choose the linear polar radius function form to decrease the singularity characteristic induced by the plastic yield zone of crack tip and still adopt the triangle basis function form to describe the displacement distribution character with the polar angle of crack tip. By combining the above two functions, the elastoplastic crack tip enrichment function is constructed. And we choose the exponent disconnected function to reflect the discontinuous characteristic of crack. In the frame of the extended finite element method, the increment iterative form of elastoplastic extended finite element method is deduced by virtual work principle. For nonuniform hardening material such as concrete material, in order to avoid the non-symmetry characteristics of stiffness matrix induced by the nonassociate flowing of plastic strain, the plastic flowing rule containing cross item based on the least energy dissipation principle is adopted. Finally, the simulations of the elastoplastic crack propagation process are implemented through three numerical examples. From the simulation process of the first numerical example, the results show that the elastoplastic X-FEM constructed in this paper can not only track the crack propagation but also model the plastic zone of nonassociated flowing material. It can be seen from the second numerical example results that the elastoplastic X-FEM proposed in this paper can conveniently embody the plastic zone of the associated flowing material when ignoring the influence of the cross item on the flowing rule based on the least energy dissipation principle. From the simulation process of the third numerical example, the calculations simulated by elastoplastic X-FEM proposed in this paper agree well with the experimental observations, which verifies the validity of the construction of displacement enrichment function in X-FEM for elastoplastic crack problem.

#### Acknowledgments

This work was supported by the National Natural ScienceFoundation of China (nos. 51179064, 11372099, and 11132003) and the State Key Laboratory Open Foundation of Hydrology-Water Resources and Hydraulic Engineering at the Hohai University (no. 2011490911).