A Cell-Based Smoothed XFEM for Fracture in Piezoelectric Materials

This paper presents a cell-based smoothed extended finite element method (CS-XFEM) to analyze fractures in piezoelectric materials.Themethod, 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.


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][8][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][22][23][24][25][26][27].Based on generalized gradient smoothing technique a lot of novel and powerful numerical methods have been developed [28][29][30][31][32][33][34][35][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.

Advances in Materials Science and Engineering
In terms of the advantages, smoothing technique was incorporated into the extended FEM [37][38][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 Section 2, 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 -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.

Governing Equations
The electroelastic response of a piezoelectric body of volume Ω and regular boundary surface  is governed by the mechanical and electrostatic equilibrium equations: where   is mechanical body force,  is electric body charge,   is the symmetric Cauchy stress tensor, and   is electric displacement vector components.The constitutive equations for a two-dimensional piezoelectric material in the - plane can be expressed in terms of the stress and the electric field: where   ,   , and   are the strain tensor, the electric displacement vector, and the electric field vector, respectively;   ,   , and   denote elastic stiffness at constant electric field, piezoelectric constants, and dielectric permittivity at constant strain, respectively.
The strains tensor are related to displacements by and the electric field vector is related to electric potential by The piezoelectric body Ω could be subjected to the following essential and natural boundary conditions (Figure 1).] where C  are the elastic compliance constants, e  are piezoelectric constants, and   are the dielectric constants.

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 where    (x),    (x), and    (x) are shape functions of the nodal displacement, while u  are the nodal degrees of freedom associated with node , and a  and b  are additional nodal degrees of freedom corresponding to the Heaviside function (x) and the near-tip functions, {F  } 1≤≤4 , respectively.   (x),    (x), and    (x) are shape functions of the nodal electric potential, while Φ  are the nodal degrees of freedom associated with node , and   and   are additional nodal degrees of freedom corresponding to the Heaviside function (x) and the near-tip functions, {F  } 1≤≤4 , respectively.
Nodes in set  CS- have supports split by crack and nodes in set  CS- 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  CS-FEM is associated with nodes of CS-FEM, shown in Figure 2.
As illustrated in Figure 2, 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 in  CS- whose support domain is split by the crack and tip enriched degrees of freedom are added to nodes in set  CS- whose support domain contains the crack tip.In order to keep the convergence rate as high as possible, a socalled geometric enrichment which is independent from the discretization is used [45][46][47]: where x * is a point on the crack surface; see Figure 3.
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]: where  and  are polar coordinates in the local crack tip coordinate system; see Figure 4.
Employing the strain smoothing operation, the smoothed strain in the domain Ω  from the displacement approximation in (8) can be written in the following matrix form: where B   (x  ) is the smoothed strain gradient matrix for the standard CS-FEM part; B   (x  ) and B   (x  ) 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 Ω  from the displacement approximation in (9) can be written in the following matrix form: where B   (x  ) is the smoothed electric field gradient matrix for the standard CS-FEM part; B   (x  ) and B   (x  ) 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: ] ,  = , , , where where  seg is the number of segments of the boundary Γ   ,  gau is the number of Gauss points used in each segment,  , is the corresponding Gauss weights,   and   are the outward unit normal components to each segment on the smoothing domain boundary, and  , is the th Gaussian point on the th segment of the boundary Γ   .The standard discrete system of equations is obtained: where

Electromechanical 𝐽-Integral
According to research by Rice [49] and Eshelby [50] on pure mechanical applications, Equation ( 18) describes the material force when an electromechanically loaded domain is virtually displaced by the vector   .The term in brackets is called the piezoelectric energy momentum tensor: analogous to the energy momentum tensor, which was introduced by Eshelby [50].If the integration path  contains neither defects nor source terms   and  Ω , the force vanishes and F  ≡ 0. Now we consider a path Γ  enclosing the tip of a crack; see Figure 5.The electromechanical  em  -integral vector is the material force associated with the crack tip singularity.It is defined as the limit when Γ  is shrunk towards  → 0 For the numerical analysis of the  em  -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  = Γ − Γ  + Γ + + Γ − with an outward pointing normal vector is given, as can be seen in Figure 5. Then ( 21) can be written as follows: Now a weighting function  is introduced, which should be continuous and satisfy the condition The calculation of the -integral is carried out in a ring of elements surrounding the crack tip.The elements within the ring move as a rigid body. is a constant in these elements; so the derivative of  with respect to   is zero.For the elements outside the ring,  is zero, and again the derivative of  is zero.
For the elements belonging to the ring, the vector  is 1.If the weighting function  is introduced in ( 22), the contribution of the integral along Γ disappears and The transformation into an equivalent domain integral leads to If we use the expression of   for (20), we get the 2dimensional  em  -integral as a domain integral: 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 -directions; that is, Equation ( 26) can be reduced as The -component of the electromechanical  em  -integral vector has a physical meaning of the energy release rate  =  em 1 .The mechanical stresses and the electric fluxes behave singular as 1/√, whereas the electric potential and the mechanical displacements show a parabolic shape ∼√.The angular functions   ,   ,   , and V  depend only on material constants.The coefficients  I , II , and  III are the well-known mechanical stress intensity factors, which are complemented by the new forth "electric intensity factor"  D , that characterizes the electric field singularity.Mutual interdependence between mechanical and electrical crack tip parameters can be given by 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 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  is reduced to the following expression, where,   ,   ,   , , and  are material constants: The energy release rate consists of the mechanical terms (for each opening mode) and the electric contribution.

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  = 0.001 C/m 2 .Figure 6 shows the polarization direction as .The distance of the central crack along the  1 -direction is 2 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.The exact analytical solution for the crack in the infinite plane under far-field loads  ∞ 33 ,  ∞ 13 , and  ∞ 3 was given by Pak [11].Consider In order to verify the reliability of CS-XFEM, we set crack length as 2 = 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.The results of stress intensity factors and electric displacements factors produced by CS-XFEM and XFEM for three loading cases are listed in Table 2.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.
The normalized mechanical and electrical intensity factors of PZT-4 and P7 with different length of cracks and model (I) grids are listed in Table 3.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.

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  ∞ 3 = 0.001 C/m 2 ; see Figure 8.
From Table 4, 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.
The normalized mechanical and electrical intensity factors of PZT-4 and P7 with different length of cracks when using model (I) grids are listed in Table 5.It is obvious that    the CS-XFEM can produce more accurate results than XFEM using the same number of nodes.

Piezoelectric Model with a Hole with Cracks.
A piezoelectric model, with a center circular hole  = 10 mm and horizontal cracks  on the left and right, is subjected to unidirectional uniform tensile  ∞ = 10 MPa and electric displacement  ∞ = 10 −3 C/m 2 at infinity.The geometric model is expressed as in Figure 9, the length of side  = 200 mm,  is the direction of polarization, and material is PZT-H5.The situation of meshing when  = 4 mm was given by Figure 10.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 Table 6.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.

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

Figure 1 :
Figure 1: Piezoelectric domain with a crack.

Figure 3 :Figure 4 :
Figure 3: Normal and tangential coordinates for a crack.

Figure 5 :
Figure 5: Combined integration path around a crack tip.

Figure 7 :
Figure 7: Mesh of the piezoelectric plate.

Figure 9 :
Figure 9: Piezoelectric model with a hole with cracks.

Table 2 :
Comparison of XFEM results and analytical solutions for electromechanical Griffith crack (uniaxial load).

Table 3 :
Normalized intensity factors under different crack lengths for electromechanical Griffith crack (uniaxial load).

Table 4 :
Comparison of XFEM results and analytical solutions for electromechanical Griffith crack (shear load).

Table 5 :
Normalized intensity factors under different crack lengths for electromechanical Griffith crack (shear load).

Table 6 :
Intensity factor under different crack lengths.