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

Authors: Xia Xiaozhou; Zhang Qing; Wang Hong; Jiang Qun
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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. Section2 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)u(x)=∑i∈INi(x)ui+∑j∈J∩INj(x)ϕj(x)aj+∑m∈M∩INm(x)ψm(x)cm,
where I is the set of all the nodes of an element, ui the master degree of node i, Ni the interpolation shape function of node i, J the set of the element nodes containing discontinuous enrichment degree, ϕj(x) the discontinuous enrichment displacement function of node j, aj the corresponding enrichment degree of node j presented as the hollow points in Figure 1, M the set of the element nodes containing crack tip enrichment degree, ψm(x) the crack tip enrichment displacement function of node m, and cm the corresponding crack tip enrichment degree of node m presented as the solid points in Figure 1.Figure 1
The nodes containing enrichment degree.For discontinuous enrichment displacement function (see Figure3), considering the fast decline influence with the increment of perpendicular distance to crack, the exponent disconnected function form presented as follows is adopted [22]:
(2)ϕi(x)=h(d(x))={e-d(x)d(x)>0-ed(x)d(x)<0,
where d(x) is defined as
(3)d(x)=minx-∈f(x)=0∥x-x-∥sign(f(x)·f(xi)).
Here, f(x)=0 is the curve equation of crack shown in Figure 2.Figure 2
The perpendicular distance of pointx to crack curve f(x)=0.Figure 3
The exponent discontinuous functionh(d(x)).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(4a)ψm(x)cm=[βxx′βxy′βyx′βyy′][f11mf12mf21mf22m]{KImKIIm},(4b)ψm(x)cm=[βxx′βxy′βxz′βyx′βyy′βyz′βzx′βzy′βzz′][f11mf12m0f21mf22m000f33m]{KImKIImKIIIm},where f11m=(r/2μ2π)cos(θ/2)(κ-1+2sin2(θ/2)), f12m=(r/2μ2π)sin(θ/2)(κ+1+2cos2(θ/2)), f21m=(r/2μ2π)sin(θ/2)(κ+1-2cos2(θ/2)), f22m=(r/2μ2π)cos(θ/2)(κ-1-2sin2(θ/2)), f33m=  (r/μ)2/πsin(θ/2), and [βxx′βxy′βxz′βyx′βyy′βyz′βzx′βzy′βzz′] 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)ε=ε(u)=∇su,
where ∇s 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)ε=Buue+Baae+Bcce,
where the strain matrixes Bu, Ba, and Bc 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:(7)Bui=[Ni,x000Ni,y000Ni,zNi,yNi,x00Ni,zNi,yNi,z0Ni,x];Bai=[Ni,xϕi+Niϕi,x000Ni,yϕi+Niϕi,y000Ni,zϕi+Niϕi,zNi,yϕi+Niϕi,yNi,xϕi+Niϕi,x00Ni,zϕi+Niϕi,zNi,yϕi+Niϕi,yNi,zϕi+Niϕi,z0Ni,xϕi+Niϕi,x],Bci=[Ni,x(βxx′f11i+βxy′f21i)+Ni(βxx′f11,xi+βxy′f21,xi)Ni,x(βxx′f12i+βxy′f22i)+Ni(βxx′f12,xi+βxy′f22,xi)Ni,xβxz′f33i+Niβxz′f33,xiNi,y(βyx′f11i+βyy′f21i)+Ni(βyx′f11,yi+βyy′f21,yi)Ni,y(βyx′f12i+βyy′f22i)+Ni(βyx′f12,yi+βyy′f22,yi)Ni,yβyz′f33i+Niβyz′f33,yiNi,z(βzx′f11i+βzy′f21i)+Ni(βzx′f11,zi+βzy′f21,zi)Ni,z(βzx′f12i+βzy′f22i)+Ni(βzx′f12,zi+βzy′f22,zi)Ni,zβzz′f33i+Niβzz′f33,ziNi,y(βxx′f11i+βxy′f21i)+Ni(βxx′f11,yi+βxy′f21,yi)Ni,y(βxx′f12i+βxy′f22i)+Ni(βxx′f12,yi+βxy′f22,yi)Ni,yβxz′f33i+Niβxz′f33,yi+Ni,x(βyx′f11i+βyy′f21i)+Ni(βyx′f11,xi+βyy′f21,xi)+Ni,x(βyx′f12i+βyy′f22i)+Ni(βyx′f12,xi+βyy′f22,xi)+Ni,xβyz′f33i+Niβyz′f33,xiNi,z(βyx′f11i+βyy′f21i)+Ni(βyx′f11,zi+βyy′f21,zi)+Ni,y(βzx′f11i+βzy′f21i)+Ni(βzx′f11,yi+βzy′f21,yi)Ni,z(βxx′f11i+βxy′f21i)+Ni(βxx′f11,zi+βxy′f21,zi)+Ni,x(βzx′f11i+βzy′f21i)+Ni(βzx′f11,xi+βzy′f21,xi)Ni,z(βyx′f12i+βyy′f22i)+Ni(βyx′f12,zi+βyy′f22,zi)+Ni,y(βzx′f12i+βzy′f22i)+Ni(βzx′f12,yi+βzy′f22,yi)Ni,z(βxx′f12i+βxy′f22i)+Ni(βxx′f12,zi+βxy′f22,zi)+Ni,x(βzx′f12i+βzy′f22i)+Ni(βzx′f12,xi+βzy′f22,xi)Ni,zβyz′f33i+Niβyz′f33,zi+Ni,yβzz′f33i+Niβzz′f33,yiNi,zβxz′f33i+Niβxz′f33,zi+Ni,xβzz′f33i+Niβzz′f33,xi].


## 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:(8)∇·σ+F=0,inΩ,n=n-,onΓuτ,σ·n=T-,onΓTτ,σ·n=0,onΓcτ(open),σA·nA+σBnB=0,onΓcτ(close),uA·nA+uBnB≤0,onΓcτ(noembedeachother),
where F is the body force, Ω the calculation domain, Γuτ the displacement boundary at τ time, ΓTτ the force boundary at τ time, Γcτ the crack boundary, and nA and nB the corresponding outside normal line vector of crack surfaces A and B, 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.
(9)dεijp=∑p=1mdλp∂Gp∂σij(GP≠FP),
where GP denotes the pth potential function, and FP denotes the pth yield function. For kinetic hardening material, FP(σij-αijp)-K=0. Here, αijp is the hardening parameters, and K is the yield strength parameters as
(10)dεijp=-σkl∂(dεklp)∂σij+∑p=1mdλp∂Fp∂σij,
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 Fp denotes the pth yield surface. So the constitutive law of kinetic hardening material can be easily deduced as follows [23]:
(11)dσ=D(dε-1h(∂F∂σD)TD∂F∂σdε)+2Dσ∂(dεp)∂σ=Depdε+2Dσ∂(dεp)∂σ.


### 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)Δu(x)=NΔue+(Nϕ)Δae+(Nψ)Δce.
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):
(13)∫Vδ(∇s(Δu))TDepτ∇s(Δu)dV=∫Vδ(Δu)TFt+ΔtdV+∫Sσδ(Δu)TT-t+ΔtdS-∫Vδ(∇s(Δu))TtσdV-2∫Vδ(∇s(Δu))TD∂(Δεp)∂σtσdV.
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:
(14)∑e∫Ωe[δ(Δce)TBcTDepτBcΔceδ(Δue)TBuTDepτBuΔue+δ(Δue)TBuTDepτBaΔae+δ(Δue)TBuTDepτBcΔce+δ(Δae)TBaTDepτBuΔue+δ(Δae)TBaTDepτBaΔae+δ(Δae)TBaTDepτBcΔce+δ(Δce)TBcTDepτBuΔue+δ(Δce)TBcTDepτBaΔae+δ(Δce)TBcTDepτBcΔce]dv=∑e∫Ωeδ(NΔue+(Nϕ)Δae+(Nψ)Δce)Tt+Δtbdv+∑e∫Γeδ(NΔue+(Nϕ)Δae+(Nψ)Δce)Tt+Δtt-dΓ-∑e∫Ωeδ(BuΔue+BaΔae+BcΔce)Tσtdv-2∑e∫Ωeδ(BuΔue+BaΔae+BcΔce)TD∂(Δεp)∂σσtdv.
Further combining the same type of increment displacement degree yields
(15)∑eδ(Δue)T[∫ΩeBuTDepτBudv·Δue+∫ΩeBuTDepτBadv·Δae+∫ΩeBuTDepτBcdv·Δce]+∑eδ(Δae)T[∫ΩeBaTDepτBudv·Δue+∫ΩeBaTDepτBadv·Δae+∫ΩeBaTDepτBcdv·Δce]+∑eδ(Δce)T[∫ΩeBcTDepτBudv·Δue+∫ΩeBcTDepτBadv·Δae+∫ΩeBcTDepτBcdv·Δce]=∑eδ(Δue)T[∫ΩeNTFt+Δtdv+∫ΓeNTT-t+ΔtdΓ]+∑eδ(Δae)T[∫Ωe(Nϕ)TFt+Δtdv+∫Γe(Nϕ)TT-t+ΔtdΓ]+∑eδ(Δce)T[∫Ωe(Nψ)TFt+Δtdv+∫Γe(Nψ)TT-t+ΔtdΓ]-∑eδ(Δue)T∫ΩeBuT(I+2D∂(Δεp)∂σ)σtdv-∑eδ(Δae)T∫ΩeBaT(I+2D∂(Δεp)∂σ)σtdv-∑eδ(Δce)T∫ΩeBcT(I+2D∂(Δεp)∂σ)σtdv.
Formula (15) can be finally made by simplification as follows:
(16)δ(Δu)TKuuτΔu+δ(Δu)TKuaτΔa+δ(Δu)TKucτΔc+δ(Δa)TKauτΔu+δ(Δa)TKaaτΔa+δ(Δa)TKacτΔc+δ(Δc)TKcuτΔu+δ(Δc)TKcaτΔa+δ(Δc)TKccτΔc=δ(Δu)TFext+δ(Δa)TQext+δ(Δc)TCext-δ(Δu)T∑e∫ΩeBuT(I+2D∂(Δεp)∂σ)σtdv-δ(Δa)T∑e∫ΩeBaT(I+2D∂(Δεp)∂σ)σtdv-δ(Δc)T∑e∫ΩeBcT(I+2D∂(Δεp)∂σ)σtdv,
where
(17)Kuuτ=∑eguTKτuuegu,Kuaτ=∑eguTKτuaega,Kucτ=∑eguTKτucegc,Kauτ=∑egaTKτauegu,Kaaτ=∑egaTKτaaega,Kacτ=∑egaTKτacegc,Kcuτ=∑egcTKτcuegu,Kcaτ=∑egcTKτcaega,Kccτ=∑egcTKτccegc,Fext=∑e(∫ΩeNTFdv+∫ΓeNTT-dΓ),Qext=∑e(∫Ωe(Nϕ)TFdv+∫Γe(Nϕ)TT-dΓ),Cext=∑e(∫Ωe(Nψ)TFdv+∫Γe(Nψ)TT-dΓ).
Here,
(18)Kτuue=∫ΩeBuTDτepBudv,Kτuae=∫ΩeBuTDτepBadv,Kτuce=∫ΩeBuTDτepBcdv,Kτaue=∫ΩeBaTDτepBudv,Kτaae=∫ΩeBaTDτepBadv,Kτace=∫ΩeBaTDτepBcdv,Kτcue=∫ΩeBcTDτepBudv,Kτcae=∫ΩeBcTDτepBadv,Kτcce=∫ΩeBcTDτepBcdv.
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:
(19)[KuuτKuaτKucτKauτKaaτKacτKcuτKcaτKccτ]{ΔuΔaΔc}={Fext-∑e∫ΩeBuT(I+2D∂(Δεp)∂σ)σtdvQext-∑e∫ΩeBaT(I+2D∂(Δεp)∂σ)σtdvCext-∑e∫ΩeBcT(I+2D∂(Δεp)∂σ)σtdv}.


### 4.4. The Extension Criterion of Elastoplastic Crack

For the elastoplastic fracture problem, the energy release rate criterion is adopted. Figure5 shows a planform of a surface shape crack being of virtual extension Δa at the X-Y coordinate face. And the definition of the crack local coordinate system X-Y-Z- 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:
(20)GI≅limΔa→012ΔA∫ΔAσ-zz(1)(Δa-r,0)Δw-(2)(r,π)dxdy,(21)GII≅limΔa→012ΔA∫ΔAτ-xz(1)(Δa-r,0)Δu-(2)(r,π)dxdy,(22)GIII≅limΔa→012ΔA∫ΔAτ-yz(1)(Δa-r,0)Δv-(2)(r,π)dxdy,
where σ-zz(1)(Δa-r,0), τ-xz(1)(Δa-r,0), and τ-yz(1)(Δa-r,0) are the stress components of virtual crack extension face for actual crack front, and Δu-(2)(r,π), Δv-(2)(r,π), and Δw-(2)(r,π) 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 Figure6, the strain energy release rates of its three types of fracture model are as follows:
(23)GI≅F-Z1Δw-3,42ΔA,GII≅F-X1Δu-3,42ΔA,GIII≅F-Y1Δv-3,42ΔA,
where F-X1, F-Y1, and F-Z1 are the nodal force of nonnode 1. Δu-3,4, Δv-3,4, and Δw-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 KI and KII. But the values of KI and KII should be gained by linear out-interpolating values of KI and KII of those points out of the plastic zone presented in Figure 7 as follows:
(24)θc=arccos3KII2+KI4+8KI2KII2KI2+9KII2.Figure 7
The value ofKI 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 Figure20.


## 4.1. Governing Equation of X-FEM

The equilibrium equation and the boundary condition can be written as follows:(8)∇·σ+F=0,inΩ,n=n-,onΓuτ,σ·n=T-,onΓTτ,σ·n=0,onΓcτ(open),σA·nA+σBnB=0,onΓcτ(close),uA·nA+uBnB≤0,onΓcτ(noembedeachother),
where F is the body force, Ω the calculation domain, Γuτ the displacement boundary at τ time, ΓTτ the force boundary at τ time, Γcτ the crack boundary, and nA and nB the corresponding outside normal line vector of crack surfaces A and B, 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.
(9)dεijp=∑p=1mdλp∂Gp∂σij(GP≠FP),
where GP denotes the pth potential function, and FP denotes the pth yield function. For kinetic hardening material, FP(σij-αijp)-K=0. Here, αijp is the hardening parameters, and K is the yield strength parameters as
(10)dεijp=-σkl∂(dεklp)∂σij+∑p=1mdλp∂Fp∂σij,
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 Fp denotes the pth yield surface. So the constitutive law of kinetic hardening material can be easily deduced as follows [23]:
(11)dσ=D(dε-1h(∂F∂σD)TD∂F∂σdε)+2Dσ∂(dεp)∂σ=Depdε+2Dσ∂(dεp)∂σ.


## 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)Δu(x)=NΔue+(Nϕ)Δae+(Nψ)Δce.
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):
(13)∫Vδ(∇s(Δu))TDepτ∇s(Δu)dV=∫Vδ(Δu)TFt+ΔtdV+∫Sσδ(Δu)TT-t+ΔtdS-∫Vδ(∇s(Δu))TtσdV-2∫Vδ(∇s(Δu))TD∂(Δεp)∂σtσdV.
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:
(14)∑e∫Ωe[δ(Δce)TBcTDepτBcΔceδ(Δue)TBuTDepτBuΔue+δ(Δue)TBuTDepτBaΔae+δ(Δue)TBuTDepτBcΔce+δ(Δae)TBaTDepτBuΔue+δ(Δae)TBaTDepτBaΔae+δ(Δae)TBaTDepτBcΔce+δ(Δce)TBcTDepτBuΔue+δ(Δce)TBcTDepτBaΔae+δ(Δce)TBcTDepτBcΔce]dv=∑e∫Ωeδ(NΔue+(Nϕ)Δae+(Nψ)Δce)Tt+Δtbdv+∑e∫Γeδ(NΔue+(Nϕ)Δae+(Nψ)Δce)Tt+Δtt-dΓ-∑e∫Ωeδ(BuΔue+BaΔae+BcΔce)Tσtdv-2∑e∫Ωeδ(BuΔue+BaΔae+BcΔce)TD∂(Δεp)∂σσtdv.
Further combining the same type of increment displacement degree yields
(15)∑eδ(Δue)T[∫ΩeBuTDepτBudv·Δue+∫ΩeBuTDepτBadv·Δae+∫ΩeBuTDepτBcdv·Δce]+∑eδ(Δae)T[∫ΩeBaTDepτBudv·Δue+∫ΩeBaTDepτBadv·Δae+∫ΩeBaTDepτBcdv·Δce]+∑eδ(Δce)T[∫ΩeBcTDepτBudv·Δue+∫ΩeBcTDepτBadv·Δae+∫ΩeBcTDepτBcdv·Δce]=∑eδ(Δue)T[∫ΩeNTFt+Δtdv+∫ΓeNTT-t+ΔtdΓ]+∑eδ(Δae)T[∫Ωe(Nϕ)TFt+Δtdv+∫Γe(Nϕ)TT-t+ΔtdΓ]+∑eδ(Δce)T[∫Ωe(Nψ)TFt+Δtdv+∫Γe(Nψ)TT-t+ΔtdΓ]-∑eδ(Δue)T∫ΩeBuT(I+2D∂(Δεp)∂σ)σtdv-∑eδ(Δae)T∫ΩeBaT(I+2D∂(Δεp)∂σ)σtdv-∑eδ(Δce)T∫ΩeBcT(I+2D∂(Δεp)∂σ)σtdv.
Formula (15) can be finally made by simplification as follows:
(16)δ(Δu)TKuuτΔu+δ(Δu)TKuaτΔa+δ(Δu)TKucτΔc+δ(Δa)TKauτΔu+δ(Δa)TKaaτΔa+δ(Δa)TKacτΔc+δ(Δc)TKcuτΔu+δ(Δc)TKcaτΔa+δ(Δc)TKccτΔc=δ(Δu)TFext+δ(Δa)TQext+δ(Δc)TCext-δ(Δu)T∑e∫ΩeBuT(I+2D∂(Δεp)∂σ)σtdv-δ(Δa)T∑e∫ΩeBaT(I+2D∂(Δεp)∂σ)σtdv-δ(Δc)T∑e∫ΩeBcT(I+2D∂(Δεp)∂σ)σtdv,
where
(17)Kuuτ=∑eguTKτuuegu,Kuaτ=∑eguTKτuaega,Kucτ=∑eguTKτucegc,Kauτ=∑egaTKτauegu,Kaaτ=∑egaTKτaaega,Kacτ=∑egaTKτacegc,Kcuτ=∑egcTKτcuegu,Kcaτ=∑egcTKτcaega,Kccτ=∑egcTKτccegc,Fext=∑e(∫ΩeNTFdv+∫ΓeNTT-dΓ),Qext=∑e(∫Ωe(Nϕ)TFdv+∫Γe(Nϕ)TT-dΓ),Cext=∑e(∫Ωe(Nψ)TFdv+∫Γe(Nψ)TT-dΓ).
Here,
(18)Kτuue=∫ΩeBuTDτepBudv,Kτuae=∫ΩeBuTDτepBadv,Kτuce=∫ΩeBuTDτepBcdv,Kτaue=∫ΩeBaTDτepBudv,Kτaae=∫ΩeBaTDτepBadv,Kτace=∫ΩeBaTDτepBcdv,Kτcue=∫ΩeBcTDτepBudv,Kτcae=∫ΩeBcTDτepBadv,Kτcce=∫ΩeBcTDτepBcdv.
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:
(19)[KuuτKuaτKucτKauτKaaτKacτKcuτKcaτKccτ]{ΔuΔaΔc}={Fext-∑e∫ΩeBuT(I+2D∂(Δεp)∂σ)σtdvQext-∑e∫ΩeBaT(I+2D∂(Δεp)∂σ)σtdvCext-∑e∫ΩeBcT(I+2D∂(Δεp)∂σ)σtdv}.


## 4.4. The Extension Criterion of Elastoplastic Crack

For the elastoplastic fracture problem, the energy release rate criterion is adopted. Figure5 shows a planform of a surface shape crack being of virtual extension Δa at the X-Y coordinate face. And the definition of the crack local coordinate system X-Y-Z- 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:
(20)GI≅limΔa→012ΔA∫ΔAσ-zz(1)(Δa-r,0)Δw-(2)(r,π)dxdy,(21)GII≅limΔa→012ΔA∫ΔAτ-xz(1)(Δa-r,0)Δu-(2)(r,π)dxdy,(22)GIII≅limΔa→012ΔA∫ΔAτ-yz(1)(Δa-r,0)Δv-(2)(r,π)dxdy,
where σ-zz(1)(Δa-r,0), τ-xz(1)(Δa-r,0), and τ-yz(1)(Δa-r,0) are the stress components of virtual crack extension face for actual crack front, and Δu-(2)(r,π), Δv-(2)(r,π), and Δw-(2)(r,π) 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 Figure6, the strain energy release rates of its three types of fracture model are as follows:
(23)GI≅F-Z1Δw-3,42ΔA,GII≅F-X1Δu-3,42ΔA,GIII≅F-Y1Δv-3,42ΔA,
where F-X1, F-Y1, and F-Z1 are the nodal force of nonnode 1. Δu-3,4, Δv-3,4, and Δw-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 KI and KII. But the values of KI and KII should be gained by linear out-interpolating values of KI and KII of those points out of the plastic zone presented in Figure 7 as follows:
(24)θc=arccos3KII2+KI4+8KI2KII2KI2+9KII2.Figure 7
The value ofKI 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 Figure20.


## 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 Figure8, 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 C=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.The deformation and crack extension of the thin plate.
(a)
(b)
(c)Figure 11
Moving of plastic zone (plastic strain ofy-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 Figure12. 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 C=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.The deformation and crack extension of the cantilever beam.
(a)
Preextension
(b)
PostextensionEquivalent plastic strain before expansion.
(a)
Unsaturated plastic zone
(b)
Saturated plastic zonePlastic zone moving (equivalent plastic strain).
(a)
Preextension
(b)
Postextension


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

The problem shown in Figure16 is an adaptation of an example presented in [25]. The initial crack length is a0=10 mm (all dimensions in the sketch of Figure 16 are given in mm), and the parameters are given as elastic modulus E=66.7 GPa, Possion ratio μ=0.33, and plate thickness 16 mm. The load applied to the extended finite element analysis is P=20kN. The extension resistance strength of the aluminum plate is Ft=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 Δa=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.Crack in a plate with a hole and the crack growth path of the experimental observations.
(a)
(b)Figure 17
The uniaxial extension stress-plastic strain curve of aluminum.The crack growth path predicted using X-FEM.
(a)
(b)
(c)
(d)The moving of plastic zone (equivalent plastic strain).
(a)
(b)
(c)Figure 20


## 5.1. Concrete Plane Thin Plate Containing Crack

As shown in Figure8, 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 C=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.The deformation and crack extension of the thin plate.
(a)
(b)
(c)Figure 11
Moving of plastic zone (plastic strain ofy-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 Figure12. 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 C=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.The deformation and crack extension of the cantilever beam.
(a)
Preextension
(b)
PostextensionEquivalent plastic strain before expansion.
(a)
Unsaturated plastic zone
(b)
Saturated plastic zonePlastic zone moving (equivalent plastic strain).
(a)
Preextension
(b)
Postextension


## 5.3. Crack in a Plate with a Hole

The problem shown in Figure16 is an adaptation of an example presented in [25]. The initial crack length is a0=10 mm (all dimensions in the sketch of Figure 16 are given in mm), and the parameters are given as elastic modulus E=66.7 GPa, Possion ratio μ=0.33, and plate thickness 16 mm. The load applied to the extended finite element analysis is P=20kN. The extension resistance strength of the aluminum plate is Ft=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 Δa=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.Crack in a plate with a hole and the crack growth path of the experimental observations.
(a)
(b)Figure 17
The uniaxial extension stress-plastic strain curve of aluminum.The crack growth path predicted using X-FEM.
(a)
(b)
(c)
(d)The moving of plastic zone (equivalent plastic strain).
(a)
(b)
(c)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.



---

*Source: 157130-2013-12-07.xml*