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

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: 𝑢(𝑥)=_𝑖∈𝐼𝑁_𝑖(𝑥)𝑢_𝑖+_{𝑗∈𝐽∩𝐼}𝑁_𝑗(𝑥)𝜙_𝑗(𝑥)𝑎_𝑗+_{𝑚∈𝑀∩𝐼}𝑁_𝑚(𝑥)𝜓_𝑚(𝑥)𝑐_𝑚, (1) 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.

Figure 1 

The nodes containing enrichment degree.

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]: 𝜙_𝑖(𝑥)(𝑑(𝑥=𝑒=ℎ))^{−𝑑(𝑥)}𝑑(𝑥)>0−𝑒^𝑑(𝑥)𝑑(𝑥)<0, (2) where 𝑑(𝑥) is defined as 𝑑(𝑥)=min_{(𝑥)𝑥∈𝑓=0}‖𝑥−𝑥‖sign𝑓(𝑥)⋅𝑓𝑥_𝑖. (3) Here, 𝑓(𝑥)=0 is the curve equation of crack shown in Figure 2.

Figure 2 

The perpendicular distance of point 𝑥 to crack curve 𝑓(𝑥)=0.

Figure 3 

The exponent discontinuous function ℎ(𝑑(𝑥)).

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 as𝜓_𝑚(𝑥)𝑐_𝑚=𝛽_𝑥𝑥′𝛽_𝑥𝑦′𝛽_𝑦𝑥′𝛽_𝑦𝑦′𝑓^𝑚₁₁𝑓^𝑚₁₂𝑓^𝑚₂₁𝑓^𝑚₂₂𝐾^𝑚_I𝐾^𝑚_II𝜓_𝑚,(𝑥)𝑐_𝑚=𝛽_𝑥𝑥′𝛽_𝑥𝑦′𝛽_𝑥𝑧′𝛽_𝑦𝑥′𝛽_𝑦𝑦′𝛽_𝑦𝑧′𝛽_𝑧𝑥′𝛽_𝑧𝑦′𝛽_𝑧𝑧′𝑓^𝑚₁₁𝑓^𝑚₁₂0𝑓^𝑚₂₁𝑓^𝑚₂₂000𝑓^𝑚₃₃𝐾^𝑚_I𝐾^𝑚_II𝐾^𝑚_III, (4a)(4b)where 𝑓^𝑚₁₁=(𝑟/2𝜇√2𝜋)cos(𝜃/2)(𝜅−1+2sin²(𝜃/2)), 𝑓^𝑚₁₂=(𝑟/2𝜇√2𝜋)sin(𝜃/2)(𝜅+1+2cos²(𝜃/2)), 𝑓^𝑚₂₁=(𝑟/2𝜇√2𝜋)sin(𝜃/2)(𝜅+1−2cos²(𝜃/2)), 𝑓^𝑚₂₂=(𝑟/2𝜇√2𝜋)cos(𝜃/2)(𝜅−1−2sin²(𝜃/2)), 𝑓^𝑚₃₃=  (𝑟/𝜇)2/𝜋sin(𝜃/2), 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: 𝜀=𝜀(𝑢)=∇_𝑠𝑢, (5) 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: 𝜀=𝐵_𝑢𝑢^𝑒+𝐵_𝑎𝑎^𝑒+𝐵_𝑐𝑐^𝑒, (6) 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:𝐵_𝑢𝑖=^{𝑁𝑖,𝑥000𝑁𝑖,𝑦000𝑁𝑖,𝑧}_{𝑁𝑖,𝑦𝑁𝑖,𝑥00𝑁𝑖,𝑧𝑁𝑖,𝑦𝑁𝑖,𝑧0𝑁𝑖,𝑥};𝐵_𝑎𝑖=^{𝑁𝑖,𝑥𝜙𝑖+𝑁𝑖𝜙𝑖,𝑥000𝑁𝑖,𝑦𝜙𝑖+𝑁𝑖𝜙𝑖,𝑦000𝑁𝑖,𝑧𝜙𝑖+𝑁𝑖𝜙𝑖,𝑧}_{𝑁𝑖,𝑦𝜙𝑖+𝑁𝑖𝜙𝑖,𝑦𝑁𝑖,𝑥𝜙𝑖+𝑁𝑖𝜙𝑖,𝑥00𝑁𝑖,𝑧𝜙𝑖+𝑁𝑖𝜙𝑖,𝑧𝑁𝑖,𝑦𝜙𝑖+𝑁𝑖𝜙𝑖,𝑦𝑁𝑖,𝑧𝜙𝑖+𝑁𝑖𝜙𝑖,𝑧0𝑁𝑖,𝑥𝜙𝑖+𝑁𝑖𝜙𝑖,𝑥},𝐵^𝑖_𝑐=^{𝑁𝑖,𝑥𝛽𝑥𝑥′𝑓𝑖11+𝛽𝑥𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑥𝑥′𝑓𝑖11,𝑥+𝛽𝑥𝑦′𝑓𝑖21,𝑥𝑁𝑖,𝑥𝛽𝑥𝑥′𝑓𝑖12+𝛽𝑥𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑥𝑥′𝑓𝑖12,𝑥+𝛽𝑥𝑦′𝑓𝑖22,𝑥𝑁𝑖,𝑥𝛽𝑥𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑥𝑧′𝑓𝑖33,𝑥𝑁𝑖,𝑦𝛽𝑦𝑥′𝑓𝑖11+𝛽𝑦𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑦𝑥′𝑓𝑖11,𝑦+𝛽𝑦𝑦′𝑓𝑖21,𝑦𝑁𝑖,𝑦𝛽𝑦𝑥′𝑓𝑖12+𝛽𝑦𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑦𝑥′𝑓𝑖12,𝑦+𝛽𝑦𝑦′𝑓𝑖22,𝑦𝑁𝑖,𝑦𝛽𝑦𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑦𝑧′𝑓𝑖33,𝑦𝑁𝑖,𝑧𝛽𝑧𝑥′𝑓𝑖11+𝛽𝑧𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑧𝑥′𝑓𝑖11,𝑧+𝛽𝑧𝑦′𝑓𝑖21,𝑧𝑁𝑖,𝑧𝛽𝑧𝑥′𝑓𝑖12+𝛽𝑧𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑧𝑥′𝑓𝑖12,𝑧+𝛽𝑧𝑦′𝑓𝑖22,𝑧𝑁𝑖,𝑧𝛽𝑧𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑧𝑧′𝑓𝑖33,𝑧𝑁𝑖,𝑦𝛽𝑥𝑥′𝑓𝑖11+𝛽𝑥𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑥𝑥′𝑓𝑖11,𝑦+𝛽𝑥𝑦′𝑓𝑖21,𝑦𝑁𝑖,𝑦𝛽𝑥𝑥′𝑓𝑖12+𝛽𝑥𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑥𝑥′𝑓𝑖12,𝑦+𝛽𝑥𝑦′𝑓𝑖22,𝑦𝑁𝑖,𝑦𝛽𝑥𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑥𝑧′𝑓𝑖33,𝑦}_{+𝑁𝑖,𝑥𝛽𝑦𝑥′𝑓}^𝑖_{11+𝛽𝑦𝑦′𝑓}^𝑖_{21+𝑁𝑖𝛽𝑦𝑥′𝑓}^𝑖_{11,𝑥+𝛽𝑦𝑦′𝑓}^𝑖_{21,𝑥+𝑁𝑖,𝑥𝛽𝑦𝑥′𝑓}^𝑖_{12+𝛽𝑦𝑦′𝑓}^𝑖_{22+𝑁𝑖𝛽𝑦𝑥′𝑓}^𝑖_{12,𝑥+𝛽𝑦𝑦′𝑓}^𝑖_{22,𝑥+𝑁𝑖,𝑥𝛽𝑦𝑧′𝑓}^𝑖_{33+𝑁𝑖𝛽𝑦𝑧′𝑓}^𝑖_{33,𝑥𝑁𝑖,𝑧𝛽𝑦𝑥′𝑓𝑖11+𝛽𝑦𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑦𝑥′𝑓𝑖11,𝑧+𝛽𝑦𝑦′𝑓𝑖21,𝑧+𝑁𝑖,𝑦𝛽𝑧𝑥′𝑓𝑖11+𝛽𝑧𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑧𝑥′𝑓𝑖11,𝑦+𝛽𝑧𝑦′𝑓𝑖21,𝑦𝑁𝑖,𝑧𝛽𝑥𝑥′𝑓𝑖11+𝛽𝑥𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑥𝑥′𝑓𝑖11,𝑧+𝛽𝑥𝑦′𝑓𝑖21,𝑧+𝑁𝑖,𝑥𝛽𝑧𝑥′𝑓𝑖11+𝛽𝑧𝑦′𝑓𝑖21+𝑁𝑖𝛽𝑧𝑥′𝑓𝑖11,𝑥+𝛽𝑧𝑦′𝑓𝑖21,𝑥𝑁𝑖,𝑧𝛽𝑦𝑥′𝑓𝑖12+𝛽𝑦𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑦𝑥′𝑓𝑖12,𝑧+𝛽𝑦𝑦′𝑓𝑖22,𝑧+𝑁𝑖,𝑦𝛽𝑧𝑥′𝑓𝑖12+𝛽𝑧𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑧𝑥′𝑓𝑖12,𝑦+𝛽𝑧𝑦′𝑓𝑖22,𝑦𝑁𝑖,𝑧𝛽𝑥𝑥′𝑓𝑖12+𝛽𝑥𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑥𝑥′𝑓𝑖12,𝑧+𝛽𝑥𝑦′𝑓𝑖22,𝑧+𝑁𝑖,𝑥𝛽𝑧𝑥′𝑓𝑖12+𝛽𝑧𝑦′𝑓𝑖22+𝑁𝑖𝛽𝑧𝑥′𝑓𝑖12,𝑥+𝛽𝑧𝑦′𝑓𝑖22,𝑥𝑁𝑖,𝑧𝛽𝑦𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑦𝑧′𝑓𝑖33,𝑧+𝑁𝑖,𝑦𝛽𝑧𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑧𝑧′𝑓𝑖33,𝑦𝑁𝑖,𝑧𝛽𝑥𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑥𝑧′𝑓𝑖33,𝑧+𝑁𝑖,𝑥𝛽𝑧𝑧′𝑓𝑖33+𝑁𝑖𝛽𝑧𝑧′𝑓𝑖33,𝑥}. (7)

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: ∇⋅𝜎+𝐹=0,inΩ,𝑛=𝑛,onΓ^𝜏_𝑢,𝜎⋅𝑛=𝑇,onΓ^𝜏_𝑇,𝜎⋅𝑛=0,onΓ^𝜏_𝑐,𝜎open_𝐴⋅𝑛^𝐴+𝜎^𝐵𝑛^𝐵=0,onΓ^𝜏_𝑐(),𝑢close_𝐴⋅𝑛^𝐴+𝑢^𝐵𝑛^𝐵≤0,onΓ^𝜏_𝑐(noembedeachother), (8) 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.

Figure 4 

The sketch of all kinds of boundary of X-FEM at 𝜏 time.

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. 𝑑𝜀^𝑝_𝑖𝑗=^𝑚_𝑝=1𝑑𝜆_𝑝𝜕𝐺_𝑝𝜕𝜎_𝑖𝑗𝐺_𝑃≠𝐹_𝑃, (9) 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 𝑑𝜀^𝑝_𝑖𝑗=−𝜎_𝑘𝑙𝜕𝑑𝜀^𝑝_𝑘𝑙𝜕𝜎_𝑖𝑗+^𝑚_𝑝=1𝑑𝜆_𝑝𝜕𝐹_𝑝𝜕𝜎_𝑖𝑗, (10) 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]: 1𝑑𝜎=𝐷𝑑𝜀−ℎ𝜕𝐹𝐷𝜕𝜎^𝑇𝐷𝜕𝐹𝜕𝜕𝜎𝑑𝜀+2𝐷𝜎𝑑𝜀_𝑝𝜕𝜎=𝐷_𝑒𝑝𝜕𝑑𝜀+2𝐷𝜎𝑑𝜀_𝑝.𝜕𝜎 (11)

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): Δ𝑢(𝑥)=𝑁Δ𝑢^𝑒+𝑁𝜙Δ𝑎^𝑒+𝑁𝜓Δ𝑐^𝑒. (12) 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): _𝑉𝛿∇_𝑠(Δ𝑢)^𝑇𝜏𝐷_𝑒𝑝∇_𝑠(Δ𝑢)𝑑𝑉=_𝑉𝛿(Δ𝑢)^{𝑇𝑡+Δ𝑡}𝐹𝑑𝑉+_𝑆𝜎𝛿(Δ𝑢)^{𝑇𝑡+Δ𝑡}𝑇𝑑𝑆−_𝑉𝛿∇_𝑠(Δ𝑢)^𝑇𝑡𝜎𝑑𝑉−2_𝑉𝛿∇_𝑠(Δ𝑢)^𝑇𝐷𝜕Δ𝜀_𝑝𝜕𝜎^𝑡𝜎𝑑𝑉. (13) 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: _𝑒Ω𝑒𝛿Δ𝑢^𝑒𝑇𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑢Δ𝑢^𝑒+𝛿Δ𝑢^𝑒𝑇𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑎Δ𝑎^𝑒+𝛿Δ𝑢^𝑒𝑇𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑐Δ𝑐^𝑒+𝛿Δ𝑎^𝑒𝑇𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑢Δ𝑢^𝑒+𝛿Δ𝑎^𝑒𝑇𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑎Δ𝑎^𝑒+𝛿Δ𝑎^𝑒𝑇𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑐Δ𝑐^𝑒+𝛿Δ𝑐^𝑒𝑇𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑢Δ𝑢^𝑒+𝛿Δ𝑐^𝑒𝑇𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑎Δ𝑎^𝑒+𝛿Δ𝑐^𝑒𝑇𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑐Δ𝑐^𝑒𝑑𝑣=_𝑒Ω𝑒𝛿𝑁Δ𝑢^𝑒+𝑁𝜙Δ𝑎^𝑒+𝑁𝜓Δ𝑐^{𝑒𝑇𝑡+Δ𝑡}𝑏𝑑𝑣+_𝑒Γ𝑒𝛿𝑁Δ𝑢^𝑒+𝑁𝜙Δ𝑎^𝑒+𝑁𝜓Δ𝑐^{𝑒𝑇𝑡+Δ𝑡}−𝑡𝑑Γ_𝑒Ω𝑒𝛿𝐵_𝑢Δ𝑢^𝑒+𝐵_𝑎Δ𝑎^𝑒+𝐵_𝑐Δ𝑐^𝑒𝑇𝑡𝜎𝑑𝑣−2_𝑒Ω𝑒𝛿𝐵_𝑢Δ𝑢^𝑒+𝐵_𝑎Δ𝑎^𝑒+𝐵_𝑐Δ𝑐^𝑒𝑇𝐷𝜕Δ𝜀_𝑝𝜕𝜎^𝑡𝜎𝑑𝑣. (14) Further combining the same type of increment displacement degree yields _𝑒𝛿Δ𝑢^𝑒𝑇_Ω𝑒𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑢𝑑𝑣⋅Δ𝑢^𝑒+_Ω𝑒𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑎𝑑𝑣⋅Δ𝑎^𝑒+_Ω𝑒𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑐𝑑𝑣⋅Δ𝑐^𝑒+_𝑒𝛿Δ𝑎^𝑒𝑇_Ω𝑒𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑢𝑑𝑣⋅Δ𝑢^𝑒+_Ω𝑒𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑎𝑑𝑣⋅Δ𝑎^𝑒+_Ω𝑒𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑐𝑑𝑣⋅Δ𝑐^𝑒+_𝑒𝛿Δ𝑐^𝑒𝑇_Ω𝑒𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑢𝑑𝑣⋅Δ𝑢^𝑒+_Ω𝑒𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑎𝑑𝑣⋅Δ𝑎^𝑒+_Ω𝑒𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑐𝑑𝑣⋅Δ𝑐^𝑒=_𝑒𝛿Δ𝑢^𝑒𝑇_Ω𝑒𝑁^{𝑇𝑡+Δ𝑡}𝐹𝑑𝑣+_Γ𝑒𝑁^{𝑇𝑡+Δ𝑡}𝑇𝑑Γ+_𝑒𝛿Δ𝑎^𝑒𝑇_Ω𝑒𝑁𝜙^{𝑇𝑡+Δ𝑡}𝐹𝑑𝑣+_Γ𝑒𝑁𝜙^{𝑇𝑡+Δ𝑡}+𝑇𝑑Γ_𝑒𝛿Δ𝑐^𝑒𝑇_Ω𝑒𝑁𝜓^{𝑇𝑡+Δ𝑡}𝐹𝑑𝑣+_Γ𝑒𝑁𝜓^{𝑇𝑡+Δ𝑡}−𝑇𝑑Γ_𝑒𝛿Δ𝑢^𝑒𝑇_Ω𝑒𝐵^𝑇_𝑢𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡−𝜎𝑑𝑣_𝑒𝛿Δ𝑎^𝑒𝑇_Ω𝑒𝐵^𝑇_𝑎𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡𝜎𝑑𝑣−_𝑒𝛿Δ𝑐^𝑒𝑇_Ω𝑒𝐵^𝑇_𝑐𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡𝜎𝑑𝑣. (15) Formula (15) can be finally made by simplification as follows: 𝛿(Δ𝑢)^𝑇𝐾^𝜏_𝑢𝑢Δ𝑢+𝛿(Δ𝑢)^𝑇𝐾^𝜏_𝑢𝑎Δ𝑎+𝛿(Δ𝑢)^𝑇𝐾^𝜏_𝑢𝑐()Δ𝑐+𝛿Δ𝑎^𝑇𝐾^𝜏_𝑎𝑢Δ𝑢+𝛿(Δ𝑎)^𝑇𝐾^𝜏_𝑎𝑎Δ𝑎+𝛿(Δ𝑎)^𝑇𝐾^𝜏_𝑎𝑐()Δ𝑐+𝛿Δ𝑐^𝑇𝐾^𝜏_𝑐𝑢Δ𝑢+𝛿(Δ𝑐)^𝑇𝐾^𝜏_𝑐𝑎Δ𝑎+𝛿(Δ𝑐)^𝑇𝐾^𝜏_𝑐𝑐()Δ𝑐=𝛿Δ𝑢^𝑇𝐹^ext+𝛿(Δ𝑎)^𝑇𝑄^ext+𝛿(Δ𝑐)^𝑇𝐶^ext−𝛿(Δ𝑢)^𝑇_𝑒Ω𝑒𝐵^𝑇_𝑢𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡()𝜎𝑑𝑣−𝛿Δ𝑎^𝑇_𝑒Ω𝑒𝐵^𝑇_𝑎𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡()𝜎𝑑𝑣−𝛿Δ𝑐^𝑇_𝑒Ω𝑒𝐵^𝑇_𝑐𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡𝜎𝑑𝑣, (16) where 𝐾^𝜏_𝑢𝑢=_𝑒𝑔^𝑇_𝑢^𝜏𝐾^𝑒_𝑢𝑢𝑔_𝑢,𝐾^𝜏_𝑢𝑎=_𝑒𝑔^𝑇_𝑢^𝜏𝐾^𝑒_𝑢𝑎𝑔_𝑎,𝐾^𝜏_𝑢𝑐=_𝑒𝑔^𝑇_𝑢^𝜏𝐾^𝑒_𝑢𝑐𝑔_𝑐,𝐾^𝜏_𝑎𝑢=_𝑒𝑔^𝑇_𝑎^𝜏𝐾^𝑒_𝑎𝑢𝑔_𝑢,𝐾^𝜏_𝑎𝑎=_𝑒𝑔^𝑇_𝑎^𝜏𝐾^𝑒_𝑎𝑎𝑔_𝑎,𝐾^𝜏_𝑎𝑐=_𝑒𝑔^𝑇_𝑎^𝜏𝐾^𝑒_𝑎𝑐𝑔_𝑐,𝐾^𝜏_𝑐𝑢=_𝑒𝑔^𝑇_𝑐^𝜏𝐾^𝑒_𝑐𝑢𝑔_𝑢,𝐾^𝜏_𝑐𝑎=_𝑒𝑔^𝑇_𝑐^𝜏𝐾^𝑒_𝑐𝑎𝑔_𝑎,𝐾^𝜏_𝑐𝑐=_𝑒𝑔^𝑇_𝑐^𝜏𝐾^𝑒_𝑐𝑐𝑔_𝑐,𝐹^ext=_𝑒Ω𝑒𝑁^𝑇𝐹𝑑𝑣+_Γ𝑒𝑁^𝑇,𝑄𝑇𝑑Γ_ext=_𝑒Ω𝑒𝑁𝜙^𝑇𝐹𝑑𝑣+_Γ𝑒𝑁𝜙^𝑇,𝐶𝑇𝑑Γ_ext=_𝑒Ω𝑒𝑁𝜓^𝑇𝐹𝑑𝑣+_Γ𝑒𝑁𝜓^𝑇𝑇𝑑Γ. (17) Here, ^𝜏𝐾^𝑒_𝑢𝑢=_Ω𝑒𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑢𝑑𝑣,^𝜏𝐾^𝑒_𝑢𝑎=_Ω𝑒𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑎𝑑𝑣,_𝜏𝐾^𝑒_𝑢𝑐=_Ω𝑒𝐵^𝑇_𝑢^𝜏𝐷_𝑒𝑝𝐵_𝑐𝑑𝑣,^𝜏𝐾^𝑒_𝑎𝑢=_Ω𝑒𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑢𝑑𝑣,_𝜏𝐾^𝑒_𝑎𝑎=_Ω𝑒𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑎𝑑𝑣,^𝜏𝐾^𝑒_𝑎𝑐=_Ω𝑒𝐵^𝑇_𝑎^𝜏𝐷_𝑒𝑝𝐵_𝑐𝑑𝑣,_𝜏𝐾^𝑒_𝑐𝑢=_Ω𝑒𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑢𝑑𝑣,^𝜏𝐾^𝑒_𝑐𝑎=_Ω𝑒𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑎𝑑𝑣,_𝜏𝐾^𝑒_𝑐𝑐=_Ω𝑒𝐵^𝑇_𝑐^𝜏𝐷_𝑒𝑝𝐵_𝑐𝑑𝑣. (18) 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: 𝐾^𝜏_𝑢𝑢𝐾^𝜏_𝑢𝑎𝐾^𝜏_𝑢𝑐𝐾^𝜏_𝑎𝑢𝐾^𝜏_𝑎𝑎𝐾^𝜏_𝑎𝑐𝐾^𝜏_𝑐𝑢𝐾^𝜏_𝑐𝑎𝐾^𝜏_𝑐𝑐=𝐹Δ𝑢Δ𝑎Δ𝑐_ext−∑_𝑒∫_Ω𝑒𝐵^𝑇_𝑢𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡𝑄𝜎𝑑𝑣_ext−∑_𝑒∫_Ω𝑒𝐵^𝑇_𝑎𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡𝐶𝜎𝑑𝑣_ext−∑_𝑒∫_Ω𝑒𝐵^𝑇_𝑐𝜕𝐼+2𝐷Δ𝜀_𝑝𝜕𝜎^𝑡.𝜎𝑑𝑣 (19)

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: 𝐺_I≅lim_Δ𝑎→012Δ𝐴_Δ𝐴𝜎⁽¹⁾_𝑧𝑧(Δ𝑎−𝑟,0)Δ𝑤⁽²⁾()𝐺𝑟,𝜋𝑑𝑥𝑑𝑦,_II≅lim_Δ𝑎→012Δ𝐴_Δ𝐴𝜏⁽¹⁾_𝑥𝑧(Δ𝑎−𝑟,0)Δ𝑢⁽²⁾()𝐺𝑟,𝜋𝑑𝑥𝑑𝑦,_III≅lim_Δ𝑎→012Δ𝐴_Δ𝐴𝜏⁽¹⁾_𝑦𝑧(Δ𝑎−𝑟,0)Δ𝑣⁽²⁾(𝑟,𝜋)𝑑𝑥𝑑𝑦, (20)(21)(22) where 𝜎⁽¹⁾_𝑧𝑧(Δ𝑎−𝑟,0), 𝜏⁽¹⁾_𝑥𝑧(Δ𝑎−𝑟,0), and 𝜏⁽¹⁾_𝑦𝑧(Δ𝑎−𝑟,0) 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.

Figure 5 

The crack front and the local coordinate system of the surface shape crack.

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: 𝐺_I≅𝐹_𝑍1Δ𝑤_3,42Δ𝐴,𝐺_II≅𝐹_𝑋1Δ𝑢_3,4,𝐺2Δ𝐴_III≅𝐹_𝑌1Δ𝑣_3,4,2Δ𝐴 (23) where 𝐹_𝑋1, 𝐹_𝑌1, and 𝐹_𝑍1 are the nodal force of nonnode 1. Δ𝑢_3,4, Δ𝑣_3,4, and Δ𝑤_3,4 are the relative displacements between face crack midpoint 3 and midpoint 4 on local coordinate system.

Figure 6 

The sketch of virtual crack close method for extended block element.

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: 𝜃_𝑐=arccos3𝐾²_II+𝐾⁴_I+8𝐾²_I𝐾²_II𝐾²_I+9𝐾²_II. (24)

Figure 7 

The value of 𝐾I by extrapolation of those points out of the plastic zone.

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 𝐶=2MPa and internal friction angular is given as Φ=40^∘. 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.

Figure 8 

The crack contained in the plane thin plate and its mesh.

Figure 9 

The uniaxial compression stress-strain curve of concrete.

(a)

(b)

(c)

Figure 10 

The deformation and crack extension of the thin plate.

Figure 11 

Moving of plastic zone (plastic strain of y-direction).

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 𝐶=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.

Figure 12 

The sketch of crack face in cantilever beam and subjected load.

(a) Preextension

(b) Postextension

Figure 13 

The deformation and crack extension of the cantilever beam.

(a) Unsaturated plastic zone

(b) Saturated plastic zone

Figure 14 

Equivalent plastic strain before expansion.

(a) Preextension

(b) Postextension

Figure 15 

Plastic zone moving (equivalent plastic strain).

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 𝑎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 𝑃=20kN. The extension resistance strength of the aluminum plate is 𝐹_𝑡=20MPa. 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 Δ𝑎=3mm 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)

Figure 16 

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

Figure 17 

The uniaxial extension stress-plastic strain curve of aluminum.

(a)

(b)

(c)

(d)

Figure 18 

The crack growth path predicted using X-FEM.

(a)

(b)

(c)

Figure 19 

The moving of plastic zone (equivalent plastic strain).

Figure 20 

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).