The Numerical Simulation of the Crack Elastoplastic Extension Based on the Extended Finite Element Method

1 State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China 2Department of Engineering Mechanics, Hohai University, Nanjing 210098, China 3 School of Civil and Transportation Engineering, Guangdong University of Technology, Guangzhou 510006, China 4Communication and Electron College, Jiangxi Science & Technology Normal University, Nanchang 330013, China


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 twodimensional crack extension problem [6][7][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 elasticplastic 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.

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, () = 0 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 where ] is the transform matrix from local coordinate system established at the crack face fringe to global coordinate.

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:

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.

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,   (  −    ) −  = 0. Here,    is the hardening parameters, and  is the yield strength parameters as  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]: . (11)

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: .003032 .003638 .004245 .004851 .005457 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:
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  1 ,  1 , and  1 are the nodal force of nonnode 1. Δ 3,4 , ΔV 3,4 , and Δ 3,4 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  I and  II .But the values of  I and  II should be gained by linear out-interpolating values of  I and  II 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.

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

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  = 2 MPa and internal friction angular is given as Φ = 40 ∘ .The Drucker-Prager yield criterion is adopted and the uniaxial compression stressstrain 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.

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  = 2 MPa and internal friction Φ = 40 ∘ .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.

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 0 = 10 mm (all dimensions in the sketch of Figure 16 are given in mm), and the parameters are given as elastic modulus  = 66.7 GPa, Possion ratio  = 0.33, and plate thickness 16 mm.The load applied to the extended finite element analysis is  = 20 kN.The extension resistance strength of the aluminum plate is   = 20 MPa.The Mises yield criterion is adopted and the uniaxial tension stress-strain curve is given in Figure 17.The analysis of quasistatic crack propagation has been carried out with 10 crack increments of Δ = 3 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.

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 of Hydrology-Water Resources and Hydraulic Engineering at the Hohai University (no.2011490911).
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 Mathematical Problems in Engineering

Figure 13 :
Figure 13: The deformation and crack extension of the cantilever beam.

Figure 16 :
Figure 16: Crack in a plate with a hole and the crack growth path of the experimental observations.