# A Cell-Based Smoothed XFEM for Fracture in Piezoelectric Materials

Authors: Li Ming Zhou; Guang Wei Meng; Feng Li; Shuai Gu
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

This paper presents a cell-based smoothed extended finite element method (CS-XFEM) to analyze fractures in piezoelectric materials. The method, which combines the cell-based smoothed finite element method (CS-FEM) and the extended finite element method (XFEM), shows advantages of both methods. The crack tip enrichment functions are specially derived to represent the characteristics of the displacement field and electric field around the crack tip in piezoelectric materials. With the help of the smoothing technique, integrating the singular derivatives of the crack tip enrichment functions is avoided by transforming interior integration into boundary integration. This is a significant advantage over XFEM. Numerical examples are presented to highlight the accuracy of the proposed CS-XFEM with the analytical solutions and the XFEM results.


---

## Body

## 1. Introduction

Because of their inherent coupling of electric and mechanical behaviors, piezoelectric materials have been widely used in sensors, actuators, signal transmitters and surface acoustic wave devices, aerospace panels, and civil structures. In those applications, piezoelectric materials may experience high mechanical stresses and electric field concentrations. As a result they may fail due to dielectric breakdown or fractures. These materials are usually inhomogeneous and brittle, with low ultimate tensile strength and fracture toughness. Therefore, defects such as cracks and voids should be detected to ensure reliability and durability of the piezoelectric structures. Numerical simulation of fractures in piezoelectric ceramics was conducted in [1], using finite element method (FEM) [2, 3], boundary element method (BEM) [4, 5], meshless method [6], and extended finite element method (XFEM) [7–9].The theoretical fundamentals of piezoelectric fracture mechanics for cracks were presented in [10, 11]. The analytical work to investigate the fracture mechanics of piezoelectric structures was based on Lekhnitskii and Stroh formalism [12]. An elliptic hole with a major axis perpendicular to the polarization direction inside piezoelectric ceramic and the field variables around the cavity were studied in [13]. An overview and critical discussion about the present state in the field of piezoelectric fracture mechanics were given in [14]. A survey on using numerical methods of crack analyses in piezoelectric medium along with FEM to solve fracture parameters is presented in [15]. Recently, XFEM was applied to analyze the 2D crack problems and fully coupled piezoelectric effect of piezoelectric ceramics [16]. The M integral and J integral were used to solve the stress intensity factors and electric displacements intensity factors [17, 18]. The newly developed crack tip enrichment functions of XFEM were found suitable for cracks in piezoelectric materials [19]. An extension of XFEM for dynamic fracture in piezoelectric materials was presented in [20].In recent years, the smoothed FEM has been well adopted to solve fracture mechanics problem [21–27]. Based on generalized gradient smoothing technique a lot of novel and powerful numerical methods have been developed [28–36]. A cell-based smoothed finite element method (CS-FEM) has been developed with the smoothing domains constructed based on cell of the elements. In the method line integration was used along the boundaries of the smoothing cells instead of area integration. Moreover, CS-FEM does not need mapping and derivatives. It also showed that the results are less sensitive to distorted elements.In terms of the advantages, smoothing technique was incorporated into the extended FEM [37–39]. An edge-based smoothed XFEM was developed to combine the advantages of the edge-based smoothed FEM and the XFEM [40]. A node-based smoothed XFEM was applied to linear elastic fracture mechanics [41].Recently, extended finite element method, collocation boundary element, cell-based smoothed finite element method, and so forth were used to solve the problem of fracture in piezoelectric materials. The extended finite element method focuses on the definition of new enrichment functions suitable for cracks in piezoelectric structures and generalized domain integrals are used for the determination of crack tip parameters [42]. The collocation boundary element with subdomain technique is developed, whereby the fundamental solutions are computed by a fast numerical algorithm applying Fourier series [43]. The cell-based smoothed finite element method and VCCT have been used to simulate the fracture mechanics of piezoelectric materials. A piezoelectric element tailored for VCCT was used to study the crack of piezoelectric materials. CS-FEM and VCCT were introduced into fracture mechanics of piezoelectric materials and CSFEM-VCCT for piezoelectric material with cracks was put forward [44]. In this paper cell-based smoothed XFEM (CS-XFEM) is extended to simulate flaws in piezoelectric structures. CS-XFEM combining characteristics of extended finite element method and smoothed finite element method can improve the accuracy of XFEM. Nodal enrichment can model crack propagation without remeshing. CS-XFEM has advantages that the crack tip element does not need fine division; the shape function is simple and no derivatives of shape functions are needed.This paper is outlined as follows. In Section2, the governing equations of piezoelectric materials are introduced. Section 3 focuses on the formulation of cell-based smoothed extended finite element method. Section 4 presents the electromechanical J-integral for 2D crack analysis. In Section 5, numerical examples with the assumption of impermeable crack face boundary conditions are presented to demonstrate the accuracy and efficiency of CS-XFEM. Section 6 is the conclusion.


## 2. Governing Equations

The electroelastic response of a piezoelectric body of volumeΩ and regular boundary surface S is governed by the mechanical and electrostatic equilibrium equations:(1)σij,j+fi=0inΩ,Di,i-t=0inΩ,where fi is mechanical body force, t is electric body charge, σij is the symmetric Cauchy stress tensor, and Di is electric displacement vector components.The constitutive equations for a two-dimensional piezoelectric material in thex-z plane can be expressed in terms of the stress and the electric field:(2)σij=Cijklεkl-ekijEk,Di=eiklεkl+κikEk,where εkl, Di, and Ek are the strain tensor, the electric displacement vector, and the electric field vector, respectively; Cijkl, ekij, and κik denote elastic stiffness at constant electric field, piezoelectric constants, and dielectric permittivity at constant strain, respectively.The strains tensor are related to displacements by(3)εij=ui,j+uj,i2and the electric field vector is related to electric potential by(4)Ei=-ϕ,i.The piezoelectric body Ω could be subjected to the following essential and natural boundary conditions (Figure 1).Figure 1
Piezoelectric domain with a crack.Essential Boundary Conditions. Consider(5)u=u-orϕ=ϕ-onΓeNatural Boundary Conditions. Consider(6)σijnj=fiorDini=-gonΓn,where u-, ϕ-, fi, g, and ni are mechanical displacement, electric potential, surface force components, surface charge, and outward unit normal vector components, respectively. The crack faces C+ and C- are considered traction-free. The cracks are assumed to be electrically impermeable. Nevertheless, extension to limited permeable cracks is possible.The two-dimensional matrix form of the mechanical and electrical constitutive equations can be given by [15](7)σxσzτxz=C11C130C31C33000C55εxεzγxz-0e310e33e150ExEz,DxDz=00e15e31e330εxεzγxz+κ1100κ33ExEz,where Cij are the elastic compliance constants, eij are piezoelectric constants, and κii are the dielectric constants.


## 3. Cell-Based Smoothed Extended Finite Element Method

In CS-XFEM, the approximation of displacement and electric potential field in a piezoelectric material are given by(8)uhx=∑I∈NCS-FEMNIuxuI+∑J∈NCS-cNJuxHx-HxJaJ+∑K∈NCS-fNKlx∑l=14Flx-FlxKbKl,(9)Φhx=∑I∈NCS-FEMNIϕxΦI+∑J∈NCS-cNJϕxHx-HxJαJ+∑K∈NCS-fNKϕx∑l=14Flx-FlxKβKl,where NIux, NJux, and NKlx are shape functions of the nodal displacement, while uI are the nodal degrees of freedom associated with node I, and aJ and bK are additional nodal degrees of freedom corresponding to the Heaviside function Hx and the near-tip functions, Fl1≤l≤4, respectively. NIϕx, NJϕx, and NKϕx are shape functions of the nodal electric potential, while ΦI are the nodal degrees of freedom associated with node I, and αJ and βK are additional nodal degrees of freedom corresponding to the Heaviside function Hx and the near-tip functions, Fl1≤l≤4, respectively.Nodes in setNCS-c have supports split by crack and nodes in set NCS-f which belong to the smoothing domains contain a crack tip. These nodes are enriched with the Heaviside and asymptotic branch function fields depicted with squares and circles, respectively. The support domain of NCS-FEM is associated with nodes of CS-FEM, shown in Figure 2.Figure 2
Classification of smoothing domains in CS-XFEM.As illustrated in Figure2, four-node quadrilateral cells are used for cell-based strain smoothing operation. The meshing characteristics of CS-XFEM are consistent with XFEM. The complex structure can adopt the fine mesh. For CS-FEM, the number of subcells (SC) would affect the performance of the results. In the case that the solution of SFEM (SC = 1) is overestimated to the exact solution, there exists one optimal value SC > 1 (normally SC = 4) which gives the best results as compared to the exact ones. In the case that the solution of SFEM (SC = 1) is underestimated to the exact solution, it is suggested that we should use SC = 2 to obtain the solution. This solution will be stable and have the smallest displacement and energy norms. In practical calculation, we can use SC = 4 for all problems. The results will always be better than that of standard FEM. And in many cases (not all) this solution (SC = 1) is closest to the exact solution [21]. So we can use four subcells for each quadrilateral.In CS-XFEM, Heaviside enriched degrees of freedom are added to nodes inNCS-c whose support domain is split by the crack and tip enriched degrees of freedom are added to nodes in set NCS-f whose support domain contains the crack tip. In order to keep the convergence rate as high as possible, a so-called geometric enrichment which is independent from the discretization is used [45–47]:(10)Hx=1x-x∗·n≥0-1otherwise,where x∗ is a point on the crack surface; see Figure 3.Figure 3
Normal and tangential coordinates for a crack.For piezoelectric problems, it is advisable to use the regular enrichment functions stemming from the isotropic elasticity. It should be mentioned that similar results have been obtained independently with alternative enrichment functions for cracks in confined plasticity problems [42]. The near-tip enrichment consists of functions which incorporate the radial and angular behaviors of the two-dimensional asymptotic crack tip displacement field [37, 48]:(11)Fl1≤l≤4=rsin⁡θ2cos⁡θ2cos⁡θsin⁡θ2sin⁡θcos⁡θ2,where r and θ are polar coordinates in the local crack tip coordinate system; see Figure 4.Figure 4
Polar coordinate system associated with a crack tip.Employing the strain smoothing operation, the smoothed strain in the domainΩc from the displacement approximation in (8) can be written in the following matrix form:(12)ε-k=∑I∈NCS-FEMB-IuxkuI+∑J∈NCS-cB-JaxkHx-HxJaJ+∑K∈NCS-fB-Kbxk∑l=14Flx-FlxKbKl,where B-Iuxk is the smoothed strain gradient matrix for the standard CS-FEM part; B-Jaxk and B-Kbxk correspond to the enriched parts of the smoothed strain gradient matrix associated with the Heaviside and branch functions, respectively.When employing the electric field smoothing operation, the smoothed electric field in the domainΩc from the displacement approximation in (9) can be written in the following matrix form:(13)E-k=∑I∈NCS-FEMB-IϕxkuI+∑J∈NCS-cB-JαxkHx-HxJαJ+∑K∈NCS-fB-Kβxk∑l=14Flx-FlxKβKl,where B-Iϕxk is the smoothed electric field gradient matrix for the standard CS-FEM part; B-Jαxk and B-Kβxk correspond to the enriched parts of the smoothed electric field gradient matrix associated with the Heaviside and branch functions, respectively. These matrixes can be written as follows:(14)B-Irxk=bIxrxk00bIyrxkbIyrxkbIxrxk,r=u,a,b,B-Itxk=bIxtxkbIytxk,t=ϕ,α,ϕ,where(15)bIhuxk=1Aks∫ΓksnhxkNIudΓs=∑m=1Nseg∑n=1Ngaunhxm,nNIuxm,nwm,n,h=x,z,bIhaxk=1Aks∫ΓksnhxkNIuHx-HxIdΓs=∑m=1Nseg∑n=1Ngaunhxm,nNIuxm,nHxm,n-HxIwm,n,h=x,z,bIhbxk=1Aks∫ΓksnhxkNIu∑l=14Flx-FlxKdΓs=∑m=1Nseg∑n=1Ngaunhxm,nNIuxm,nFlxm,n-FlxIwm,n,h=x,z,l=1,2,3,4,bIhϕxk=1Aks∫ΓksnhxkNIdΓs=∑m=1Nseg∑n=1Ngaunhxm,nNIϕxm,nwm,n,h=x,z,bIhαxk=1Aks∫ΓksnhxkNIϕHx-HxIdΓs=∑m=1Nseg∑n=1Ngaunhxm,nNIϕxm,nHxm,n-HxIwm,n,h=x,z,bIhβxk=1Aks∫ΓksnhxkNIϕ∑l=14Flx-FlxKdΓs=∑m=1Nseg∑n=1Ngaunhxm,nNIϕxm,nFlxm,n-FlxIwm,n,h=x,z,l=1,2,3,4,where Nseg is the number of segments of the boundary Γks, Ngau is the number of Gauss points used in each segment, wm,n is the corresponding Gauss weights, nx and nz are the outward unit normal components to each segment on the smoothing domain boundary, and xm,n is the nth Gaussian point on the mth segment of the boundary Γks.The standard discrete system of equations is obtained:(16)KIJuuKIJuaKIJubKIJauKIJaaKIJabKIJbuKIJbaKIJbbuab+KIJuϕKIJuαKIJuβKIJaϕKIJααKIJαβKIJbϕKIJbαKIJbβΦαβ=ffafb,KIJϕuKIJϕaKIJϕbKIJαuKIJαaKIJαbKIJβuKIJβaKIJβbuab-KIJϕϕKIJϕαKIJϕβKIJαϕKIJααKIJαβKIJβϕKIJβαKIJββΦαβ=ggagb,where (17)KIJuu=∑k=1Ns∫ΩksB-IuTCB-JudΩ=∑k=1NsB-IuTCB-JuAks,KIJua=∑k=1Ns∫ΩksB-IuTCB-JadΩ=∑k=1NsB-IuTCB-JaAks=KIJauT,KIJub=∑k=1Ns∫ΩksB-IuTCB-JbdΩ=∑k=1NsB-IuTCB-JbAks=KIJbuT,KIJaa=∑k=1Ns∫ΩksB-IaTCB-JadΩ=∑k=1NsB-IaTCB-JaAks,KIJab=∑k=1Ns∫ΩksB-IaTCB-JbdΩ=∑k=1NsB-IaTCB-JbAks=KIJbaT,KIJbb=∑k=1Ns∫ΩksB-IbTCB-JbdΩ=∑k=1NsB-IbTCB-JbAks,KIJuϕ=∑k=1Ns∫ΩksB-IuTeB-JϕdΩ=∑k=1NsB-IuTeB-JϕAks=KIJϕuT,KIJuα=∑k=1Ns∫ΩksB-IuTeB-JαdΩ=∑k=1NsB-IuTeB-JαAks=KIJαuT,KIJuβ=∑k=1Ns∫ΩksB-IuTeB-JβdΩ=∑k=1NsB-IuTeB-JβAks=KIJβuT,KIJaϕ=∑k=1Ns∫ΩksB-IaTeB-JϕdΩ=∑k=1NsB-IaTeB-JϕAks=KIJϕaT,KIJaα=∑k=1Ns∫ΩksB-IaTeB-JαdΩ=∑k=1NsB-IaTeB-JαAks=KIJuaT,KIJaβ=∑k=1Ns∫ΩksB-IaTeB-JβdΩ=∑k=1NsB-IaTeB-JβAks=KIJβaT,KIJbϕ=∑k=1Ns∫ΩksB-IbTeB-JϕdΩ=∑k=1NsB-IbTeB-JϕAks=KIJϕbT,KIJbα=∑k=1Ns∫ΩksB-IbTeB-JαdΩ=∑k=1NsB-IbTeB-JαAks=KIJubT,KIJbβ=∑k=1Ns∫ΩksB-IbTeB-JβdΩ=∑k=1NsB-IbTeB-JβAks=KIJβbT,KIJϕϕ=∑k=1Ns∫ΩksB-IϕTκB-JϕdΩ=∑k=1NsB-IϕTκB-JϕAks,KIJϕα=∑k=1Ns∫ΩksB-IϕTκB-JαdΩ=∑k=1NsB-IϕTκB-JαAks=KIJαϕT,KIJϕβ=∑k=1Ns∫ΩksB-IϕTκB-JβdΩ=∑k=1NsB-IϕTκB-JβAks=KIJβϕT,KIJαα=∑k=1Ns∫ΩksB-IαTκB-JαdΩ=∑k=1NsB-IαTκB-JαAks,KIJαβ=∑k=1Ns∫ΩksB-IαTκB-JβdΩ=∑k=1NsB-IαTκB-JβAks=KIJβαT,KIJββ=∑k=1Ns∫ΩksB-IβTκB-JβdΩ=∑k=1NsB-IβTκB-JβAks.


## 4. ElectromechanicalJ-Integral

According to research by Rice [49] and Eshelby [50] on pure mechanical applications,(18)Fk=∮sHδkj-σijui,k+DjEknjds,(19)H=12σijεij-DiEi.Equation (18) describes the material force when an electromechanically loaded domain is virtually displaced by the vector δxj. The term in brackets is called the piezoelectric energy momentum tensor:(20)Qkj=Hδkj-σijui,k+DjEk,analogous to the energy momentum tensor, which was introduced by Eshelby [50]. If the integration path S contains neither defects nor source terms bi and ωΩ, the force vanishes and Fk≡0. Now we consider a path Γε enclosing the tip of a crack; see Figure 5. The electromechanical Jkem-integral vector is the material force associated with the crack tip singularity. It is defined as the limit when Γε is shrunk towards r→0(21)Jkem=limr→0∫ΓεQkjnjds=limr→0∫ΓεHδkj-σijui,k+DjEknjds.Figure 5
Combined integration path around a crack tip.For the numerical analysis of theJkem-integral, an equivalent domain integral is easier to handle than a line integral. The condition for the transformation into a domain integral applying the Gaussian integral theorem, a closed integration path S=Γ-Γε+Γ++Γ- with an outward pointing normal vector is given, as can be seen in Figure 5. Then (21) can be written as follows:(22)Jkem=-∫sQkjnjds+∫Γ+Γ++Γ-Qkjnjds.Now a weighting function q is introduced, which should be continuous and satisfy the condition(23)q=0onΓ1onΓε.The calculation of the J-integral is carried out in a ring of elements surrounding the crack tip. The elements within the ring move as a rigid body. q is a constant in these elements; so the derivative of q with respect to xj is zero. For the elements outside the ring, q is zero, and again the derivative of q is zero. For the elements belonging to the ring, the vector q is 1. If the weighting function q is introduced in (22), the contribution of the integral along Γ disappears and (24)Jkem=-∫sQkjnjqds+∫Γ++Γ-Qkjnjqds.The transformation into an equivalent domain integral leads to(25)Jkem=-∫ΩQkjq,jdΩ+∫Γ++Γ-Qkjnjqds.If we use the expression ofQkj for (20), we get the 2-dimensional Jkem-integral as a domain integral: (26)Jkem=-∫ΩHδkj-σijui,k+DjEkq,jdΩ-∫ΩH,kexp-biui,k+ωΩEkqdΩ+∫Γ++Γ-Hnk-Tiui,k-ωsEkqds.If the material is homogenous, no volume forces or charges and no crack face tractions or charges are applied, and the crack face normal points to the z-directions; that is,(27)H,kexp=0,bi=0,ωΩ=0,Ti=0,ωs=0,nk=0,±1T.Equation (26) can be reduced as(28)Jkem=-∫ΩHδkj-σijui,k+DjEkq,jdΩ+∫Γ++Γ-Hnkqds.The x-component of the electromechanical Jkem-integral vector has a physical meaning of the energy release rate G=J1em.The mechanical stresses and the electric fluxes behave singular as1/r, whereas the electric potential and the mechanical displacements show a parabolic shape ~r. The angular functions fij, gij, dj, and vj depend only on material constants. The coefficients KI,KII, and KIII are the well-known mechanical stress intensity factors, which are complemented by the new forth “electric intensity factor” KD, that characterizes the electric field singularity. Mutual interdependence between mechanical and electrical crack tip parameters can be given by(29)σijr,θ=1rKIfijIθ+KIIfijIIθ+KIIIfijIIIθ+KDfijDθ,uijr,θ=1rKIgijIθ+KIIgijIIθ+KIIIgijIIIθ+KDgijDθ,Djr,θ=1rKIdjIθ+KIIdjIIθ+KIIIdjIIIθ+KDdjDθ,φr,θ=1rKIvjIθ+KIIvjIIθ+KIIIvjIIIθ+KDvjDθ.If the limit to an infinitesimal crack growth is considered taking place inside the near-tip solution, its relationship with the intensity factors can be found in [6]. Consider (30)G=12KTYC,e,κK,K=KIIKIKIIIKD.The generalized Irwin-matrix [Y] depends on the elastic, piezoelectric, and dielectric material constants and the relative orientation of the crack with respect to material’s axes and polarization vectors. For a special case when the crack lies perpendicular to the polarization, the crack tip field and the coefficients of [Y] were determined in [10, 11]. Then G is reduced to the following expression, where, CL, CT, CA, e, and κ are material constants:(31)G=GmI+GmII+GmIII+GD=12KI2CT+KII2CL+KIII2CA+KIKDe-KD2κ+KIKDe.The energy release rate consists of the mechanical terms (for each opening mode) and the electric contribution.


## 5. Numerical Examples

### 5.1. Electromechanical Griffith Crack (Uniaxial Load)

In order to test the accuracy of the CS-XFEM for crack analysis, the methods were applied to a crack in a plane subjected to normal uniaxial tension whenσ=1.0 MPa and electrical flux D=0.001 C/m2. Figure 6 shows the polarization direction as P. The distance of the central crack along the x1-direction is 2a and side length of the plate is 10 m. At the crack faces, impermeable electric boundary conditions are prescribed. Piezoelectric materials PZT-4, P7, and PZT-H5 were adopted for numerical simulation. The material parameters are shown in Table 1.Table 1
Material constant.
Material
Elastic constants/(1010 N/m2)
Piezoelectric constant/(C/m2)
Dielectric constant/(1010 C/V m)
c
11
c
12
c
13
c
44
c
55
e
31
e
33
e
15
d
11
d
33
PZT-4
13.9
7.78
7.43
2.56
11.3
−6.98
13.84
13.44
60.0
54.7
P-7
13.0
8.3
8.3
2.5
11.9
−10.3
14.7
13.5
171.0
186.0
PZT-H5
12.6
7.95
8.41
2.3
11.7
−6.5
23.3
17.44
150.3
130.0Figure 6
Electromechanical Griffith crack (uniaxial load).The exact analytical solution for the crack in the infinite plane under far-field loadsσ33∞, σ13∞, and D3∞ was given by Pak [11]. Consider(32)KI=σ33∞πa,KII=σ13∞πa,KD=D3∞πa.In order to verify the reliability of CS-XFEM, we set crack length as 2a=2 m and use two grid models ((I) 2288 elements, (II) 1521 elements) before the fracture starts, as shown in Figure 7. There are five types of smoothing elements ((i) split smoothing element: there is one Gauss point on each boundary segment for split smoothing element; (ii) split-blending smoothing element: one Gauss point on each boundary segment is sufficient; (iii) tip smoothing element: five Gauss points on a segment of smoothing element are sufficient; (iv) tip-blending smoothing element: five Gauss points on each boundary segment are sufficient; (v) standard smoothing element: one Gauss point on each boundary segment is sufficient) being used for numerical integration as mentioned in [41]. In the simulation every 4-node calculation grid adopts four smoothing elements. The results are compared with those of the XFEM and the theoretical solutions.Figure 7
Mesh of the piezoelectric plate.The results of stress intensity factors and electric displacements factors produced by CS-XFEM and XFEM for three loading cases are listed in Table2. It can be seen that the results of CS-XFEM are closer to the analytical solution than those of XFEM when using the same mesh. This confirms that combining cell-based smoothing technique with the XFEM can improve accuracy.Table 2
Comparison of XFEM results and analytical solutions for electromechanical Griffith crack (uniaxial load).
Loading
Analytical solution
Mesh (I)
Mesh (II)
CS-XFEM
XFEM
CS-XFEM
XFEM
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
Mechanical
1.7724 × 10−1
0
1.7505 × 10−1
0
1.7385 × 10−1
0
1.7494 × 10−1
0
1.7365 × 10−1
0
Electrical
0
1.7724 × 10−4
0
1.7569 × 10−4
0
1.7479 × 10−4
0
1.75481 × 10−4
0
1.7459 × 10−4
Mechanical + electrical
1.7724 × 10−1
1.7724 × 10−4
1.7505 × 10−1
1.7569 × 10−4
1.7385 × 10−1
1.7479 × 10−4
1.7494 × 10−1
1.75481 × 10−4
1.7365 × 10−1
1.7459 × 10−4The normalized mechanical and electrical intensity factors of PZT-4 and P7 with different length of cracks and model (I) grids are listed in Table3. It is obvious that the CS-XFEM can produce more accurate results than XFEM when using the same number of nodes, which indicates that the smoothing technique adopted in this work improve the calculation of normalized mechanical and electrical intensity factors for fracture in piezoelectric materials.Table 3
Normalized intensity factors under different crack lengths for electromechanical Griffith crack (uniaxial load).
a
Normalized intensity factors
PZT-4
P-7
CS-XFEM
XFEM
CS-XFEM
XFEM
0.5 m
K
I/σ33∞πa
0.989
0.965
0.978
0.967
K
D/D33∞πa
0.994
0.974
0.986
0.978
1.0 m
K
I/σ33∞πa
0.941
0.932
0.948
0.925
K
D/D33∞πa
0.963
0.940
0.956
0.948
1.5 m
K
I/σ33∞πa
0.983
0.966
0.988
0.978
K
D/D33∞πa
0.987
0.973
0.979
0.979
2.0 m
K
I/σ33∞πa
0.991
0.969
0.993
0.970
K
D/D33∞πa
0.993
0.987
0.984
0.978


### 5.2. Electromechanical Griffith Crack (Shear Load)

In the second example of straight electromechanical crack in the plane, the normal mechanical load is replaced by a shear load. The far-field loadings areσ33∞=0, σ13∞=1 MPa, and D3∞=0.001 C/m2; see Figure 8.Figure 8
Electromechanical Griffith crack (uniaxial load).From Table4, we can observe that the result of the mechanical and electrical intensity factors of PZT-H5 has a high precision in two models. The accuracy of CS-XFEM is higher than that of XFEM. The results prove that the CS-XFEM can decrease stiffness of the system and improve solution accuracy. Also model (I) of PZT-H5 was used under the CS-XFEM and XFEM using gauss integral calculation efficiency. CS-XFEM takes 46.247 seconds, while XFEM takes 48.252 second with the following CPU setup: Intel Core i5-3470 3.20 GHz, RAM: 8 G. The efficiency of CS-XFEM has been improved but is not obvious.Table 4
Comparison of XFEM results and analytical solutions for electromechanical Griffith crack (shear load).
Loading
Analytical solution
Mesh (I)
Mesh (II)
CS-XFEM
XFEM
CS-XFEM
XFEM
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
Mechanical
1.7724 × 10−1
0
1.7452 × 10−1
0
1.7315 × 10−1
0
1.7441 × 10−1
0
1.7295 × 10−1
0
Electrical
0
1.7724 × 10−4
0
1.7534 × 10−4
0
1.7392 × 10−4
0
1.7513 × 10−4
0
1.7372 × 10−4
Mechanical + electrical
1.7724 × 10−1
1.7724 × 10−4
1.7452 × 10−1
1.7534 × 10−4
1.7315 × 10−1
1.7392 × 10−4
1.7441 × 10−1
1.7513 × 10−4
1.7295 × 10−1
1.7372 × 10−4The normalized mechanical and electrical intensity factors of PZT-4 and P7 with different length of cracks when using model (I) grids are listed in Table5. It is obvious that the CS-XFEM can produce more accurate results than XFEM using the same number of nodes.Table 5
Normalized intensity factors under different crack lengths for electromechanical Griffith crack (shear load).
a
Normalized intensity factors
PZT-4
P-7
CS-XFEM
XFEM
CS-XFEM
XFEM
0.5 m
K
II/τ13∞πa
0.988
0.979
0.977
0.964
K
D/D33∞πa
0.994
0.972
0.985
0.971
1.0 m
K
II/τ13∞πa
0.968
0.948
0.957
0.937
K
D/D33∞πa
0.986
0.961
0.975
0.958
1.5 m
K
II/τ13∞πa
0.958
0.940
0.948
0.935
K
D/D33∞πa
0.973
0.955
0.963
0.940
2.0 m
K
II/τ13∞πa
0.994
0.987
0.981
0.972
K
D/D33∞πa
0.998
0.991
0.996
0.989


### 5.3. Piezoelectric Model with a Hole with Cracks

A piezoelectric model, with a center circular holeR=10 mm and horizontal cracks a on the left and right, is subjected to unidirectional uniform tensile σ∞=10 MPa and electric displacement D∞=10-3 C/m2 at infinity. The geometric model is expressed as in Figure 9, the length of side l=200 mm, P is the direction of polarization, and material is PZT-H5. The situation of meshing when a=4 mm was given by Figure 10.Figure 9
Piezoelectric model with a hole with cracks.Figure 10
Discrete model.The stress intensity factor and electric displacement intensity factor on the left side of the crack tip obtained by CS-XFEM and XFEM in different crack lengths are shown in Table6. The maximum error of calculation in CS-XFEM is 2.9%, and the maximum error of calculation in XFEM is 4.1%. The results verified the accuracy of CS-XFEM.Table 6
Intensity factor under different crack lengths.
a
Intensity factor
CS-XFEM
XFEM
Analytical solution
1 mm
K
I/107  × Pa⋅mm0.5
0.2995
0.2950
0.3076
K
D/10−9  × C⋅mm−1.5
0.2986
0.2983
0.3076
2 mm
K
I/107  × Pa⋅mm0.5
0.5322
0.5263
0.5416
K
D/10−9  × C⋅mm−1.5
0.5300
0.5230
0.5416
4 mm
K
I/107  × Pa⋅mm0.5
0.8513
0.8368
0.8681
K
D/10−9  × C⋅mm−1.5
0.8454
0.8351
0.8681
6 mm
K
I/107  × Pa⋅mm0.5
1.0517
1.0560
1.0800
K
D/10−9  × C⋅mm−1.5
1.0603
1.0432
1.0800


## 5.1. Electromechanical Griffith Crack (Uniaxial Load)

In order to test the accuracy of the CS-XFEM for crack analysis, the methods were applied to a crack in a plane subjected to normal uniaxial tension whenσ=1.0 MPa and electrical flux D=0.001 C/m2. Figure 6 shows the polarization direction as P. The distance of the central crack along the x1-direction is 2a and side length of the plate is 10 m. At the crack faces, impermeable electric boundary conditions are prescribed. Piezoelectric materials PZT-4, P7, and PZT-H5 were adopted for numerical simulation. The material parameters are shown in Table 1.Table 1
Material constant.
Material
Elastic constants/(1010 N/m2)
Piezoelectric constant/(C/m2)
Dielectric constant/(1010 C/V m)
c
11
c
12
c
13
c
44
c
55
e
31
e
33
e
15
d
11
d
33
PZT-4
13.9
7.78
7.43
2.56
11.3
−6.98
13.84
13.44
60.0
54.7
P-7
13.0
8.3
8.3
2.5
11.9
−10.3
14.7
13.5
171.0
186.0
PZT-H5
12.6
7.95
8.41
2.3
11.7
−6.5
23.3
17.44
150.3
130.0Figure 6
Electromechanical Griffith crack (uniaxial load).The exact analytical solution for the crack in the infinite plane under far-field loadsσ33∞, σ13∞, and D3∞ was given by Pak [11]. Consider(32)KI=σ33∞πa,KII=σ13∞πa,KD=D3∞πa.In order to verify the reliability of CS-XFEM, we set crack length as 2a=2 m and use two grid models ((I) 2288 elements, (II) 1521 elements) before the fracture starts, as shown in Figure 7. There are five types of smoothing elements ((i) split smoothing element: there is one Gauss point on each boundary segment for split smoothing element; (ii) split-blending smoothing element: one Gauss point on each boundary segment is sufficient; (iii) tip smoothing element: five Gauss points on a segment of smoothing element are sufficient; (iv) tip-blending smoothing element: five Gauss points on each boundary segment are sufficient; (v) standard smoothing element: one Gauss point on each boundary segment is sufficient) being used for numerical integration as mentioned in [41]. In the simulation every 4-node calculation grid adopts four smoothing elements. The results are compared with those of the XFEM and the theoretical solutions.Figure 7
Mesh of the piezoelectric plate.The results of stress intensity factors and electric displacements factors produced by CS-XFEM and XFEM for three loading cases are listed in Table2. It can be seen that the results of CS-XFEM are closer to the analytical solution than those of XFEM when using the same mesh. This confirms that combining cell-based smoothing technique with the XFEM can improve accuracy.Table 2
Comparison of XFEM results and analytical solutions for electromechanical Griffith crack (uniaxial load).
Loading
Analytical solution
Mesh (I)
Mesh (II)
CS-XFEM
XFEM
CS-XFEM
XFEM
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
K
I/(MNm−3/2)
K
D/(CM−3/2)
Mechanical
1.7724 × 10−1
0
1.7505 × 10−1
0
1.7385 × 10−1
0
1.7494 × 10−1
0
1.7365 × 10−1
0
Electrical
0
1.7724 × 10−4
0
1.7569 × 10−4
0
1.7479 × 10−4
0
1.75481 × 10−4
0
1.7459 × 10−4
Mechanical + electrical
1.7724 × 10−1
1.7724 × 10−4
1.7505 × 10−1
1.7569 × 10−4
1.7385 × 10−1
1.7479 × 10−4
1.7494 × 10−1
1.75481 × 10−4
1.7365 × 10−1
1.7459 × 10−4The normalized mechanical and electrical intensity factors of PZT-4 and P7 with different length of cracks and model (I) grids are listed in Table3. It is obvious that the CS-XFEM can produce more accurate results than XFEM when using the same number of nodes, which indicates that the smoothing technique adopted in this work improve the calculation of normalized mechanical and electrical intensity factors for fracture in piezoelectric materials.Table 3
Normalized intensity factors under different crack lengths for electromechanical Griffith crack (uniaxial load).
a
Normalized intensity factors
PZT-4
P-7
CS-XFEM
XFEM
CS-XFEM
XFEM
0.5 m
K
I/σ33∞πa
0.989
0.965
0.978
0.967
K
D/D33∞πa
0.994
0.974
0.986
0.978
1.0 m
K
I/σ33∞πa
0.941
0.932
0.948
0.925
K
D/D33∞πa
0.963
0.940
0.956
0.948
1.5 m
K
I/σ33∞πa
0.983
0.966
0.988
0.978
K
D/D33∞πa
0.987
0.973
0.979
0.979
2.0 m
K
I/σ33∞πa
0.991
0.969
0.993
0.970
K
D/D33∞πa
0.993
0.987
0.984
0.978


## 5.2. Electromechanical Griffith Crack (Shear Load)

In the second example of straight electromechanical crack in the plane, the normal mechanical load is replaced by a shear load. The far-field loadings areσ33∞=0, σ13∞=1 MPa, and D3∞=0.001 C/m2; see Figure 8.Figure 8
Electromechanical Griffith crack (uniaxial load).From Table4, we can observe that the result of the mechanical and electrical intensity factors of PZT-H5 has a high precision in two models. The accuracy of CS-XFEM is higher than that of XFEM. The results prove that the CS-XFEM can decrease stiffness of the system and improve solution accuracy. Also model (I) of PZT-H5 was used under the CS-XFEM and XFEM using gauss integral calculation efficiency. CS-XFEM takes 46.247 seconds, while XFEM takes 48.252 second with the following CPU setup: Intel Core i5-3470 3.20 GHz, RAM: 8 G. The efficiency of CS-XFEM has been improved but is not obvious.Table 4
Comparison of XFEM results and analytical solutions for electromechanical Griffith crack (shear load).
Loading
Analytical solution
Mesh (I)
Mesh (II)
CS-XFEM
XFEM
CS-XFEM
XFEM
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
K
II/(MNm−3/2)
K
D/(Cm−3/2)
Mechanical
1.7724 × 10−1
0
1.7452 × 10−1
0
1.7315 × 10−1
0
1.7441 × 10−1
0
1.7295 × 10−1
0
Electrical
0
1.7724 × 10−4
0
1.7534 × 10−4
0
1.7392 × 10−4
0
1.7513 × 10−4
0
1.7372 × 10−4
Mechanical + electrical
1.7724 × 10−1
1.7724 × 10−4
1.7452 × 10−1
1.7534 × 10−4
1.7315 × 10−1
1.7392 × 10−4
1.7441 × 10−1
1.7513 × 10−4
1.7295 × 10−1
1.7372 × 10−4The normalized mechanical and electrical intensity factors of PZT-4 and P7 with different length of cracks when using model (I) grids are listed in Table5. It is obvious that the CS-XFEM can produce more accurate results than XFEM using the same number of nodes.Table 5
Normalized intensity factors under different crack lengths for electromechanical Griffith crack (shear load).
a
Normalized intensity factors
PZT-4
P-7
CS-XFEM
XFEM
CS-XFEM
XFEM
0.5 m
K
II/τ13∞πa
0.988
0.979
0.977
0.964
K
D/D33∞πa
0.994
0.972
0.985
0.971
1.0 m
K
II/τ13∞πa
0.968
0.948
0.957
0.937
K
D/D33∞πa
0.986
0.961
0.975
0.958
1.5 m
K
II/τ13∞πa
0.958
0.940
0.948
0.935
K
D/D33∞πa
0.973
0.955
0.963
0.940
2.0 m
K
II/τ13∞πa
0.994
0.987
0.981
0.972
K
D/D33∞πa
0.998
0.991
0.996
0.989


## 5.3. Piezoelectric Model with a Hole with Cracks

A piezoelectric model, with a center circular holeR=10 mm and horizontal cracks a on the left and right, is subjected to unidirectional uniform tensile σ∞=10 MPa and electric displacement D∞=10-3 C/m2 at infinity. The geometric model is expressed as in Figure 9, the length of side l=200 mm, P is the direction of polarization, and material is PZT-H5. The situation of meshing when a=4 mm was given by Figure 10.Figure 9
Piezoelectric model with a hole with cracks.Figure 10
Discrete model.The stress intensity factor and electric displacement intensity factor on the left side of the crack tip obtained by CS-XFEM and XFEM in different crack lengths are shown in Table6. The maximum error of calculation in CS-XFEM is 2.9%, and the maximum error of calculation in XFEM is 4.1%. The results verified the accuracy of CS-XFEM.Table 6
Intensity factor under different crack lengths.
a
Intensity factor
CS-XFEM
XFEM
Analytical solution
1 mm
K
I/107  × Pa⋅mm0.5
0.2995
0.2950
0.3076
K
D/10−9  × C⋅mm−1.5
0.2986
0.2983
0.3076
2 mm
K
I/107  × Pa⋅mm0.5
0.5322
0.5263
0.5416
K
D/10−9  × C⋅mm−1.5
0.5300
0.5230
0.5416
4 mm
K
I/107  × Pa⋅mm0.5
0.8513
0.8368
0.8681
K
D/10−9  × C⋅mm−1.5
0.8454
0.8351
0.8681
6 mm
K
I/107  × Pa⋅mm0.5
1.0517
1.0560
1.0800
K
D/10−9  × C⋅mm−1.5
1.0603
1.0432
1.0800


## 6. Conclusions

With the growing applications of piezoelectric structures in innovative technical areas, problems of strength and reliability become important and have to be carefully investigated. In order to quantitatively assess fracture and fatigue, sophisticated analysis of cracks’ electromechanical properties is needed. The fracture mechanics approach for crack-like defects in piezoelectric materials reveals coupled electrical and mechanical field singularities. Effective numerical methods are needed to evaluate fracture behavior of cracks in arbitrary piezoelectric structures subjected to combined electromechanical loading.In this paper CS-XFEM with crack tip enrichment functions was presented for electromechanical crack analyses. Two examples are used to verify the accuracy of CS-XFEM. Through the numerical simulation conclusions can be drawn as follows:(1)
The cell-based smoothing technique is extended into XFEM to combine the advantages of XFEM and CS-FEM. Thus, CS-XFEM has high accuracy and convergence rate. The present method also simplifies the integration of discontinuous approximation by transforming interior integration into boundary integration. More importantly, no derivatives of shape functions are needed to compute the stiffness matrix.(2)
The examples show that CS-XFEM can perform better in accuracy and convergence rate compared with XFEM in terms of stress intensity factors and electric displacements factors.(3)
CS-XFEM is superior to XFEM in regard to the calculation of the stiffness matrix and singularity in the integrand. It also avoids the mapping process, which will increase the complexity of the calculation.



---

*Source: 4125307-2016-01-13.xml*