Enriched Element-Free Galerkin Method for Fracture Analysis of Functionally Graded Piezoelectric Materials

A new method using the enriched element-free Galerkin method (EEFGM) to model functionally graded piezoelectric materials (FGPMs) with cracks was presented. To improve the solution accuracy, extended terms were introduced into the approximation function of the conventional element-free Galerkin method (EFGM) to describe the displacement and electric fields near the crack. Compared with the conventional EFGM, the new approach requires smaller domain to describe the crack-tip singular field. Additionally, the domain of the nodes was not affected by the crack. Therefore, the visibility method and the diffraction method were no longer needed. The mechanical response of FGPM was discussed, when its material parameters changed exponentially in a certain direction. The modified J-integrals for FGPM were deduced, whose results were compared with the results of the conventional EFGM and the analytical solution. Numerical example results illustrated that this method is feasible and precise.


Introduction
Functionally graded materials (FGMs) are composite materials formed of two or more constituent phases with a continuously variable composition.During design, the requirements of structural strength, reliability, and lifetime of piezoelectric structures/components call for enhanced mechanical performance, including stress and deformation distribution under multifield loading.In recent years, the emergence of FGMs has demonstrated that they have the potential to reduce stress concentration and to provide improved residual stress distribution, enhanced thermal properties, and higher fracture toughness.Consequently, a new kind of material, functionally graded piezoelectric materials (FGPMs), has been developed to improve the reliability of piezoelectric structures by introducing the concept of the well-known FGM to piezoelectric materials [1].
At present, FGPMs are usually associated with particulate composites where the volume fraction of particles varies in one or several directions.One of the advantages of a monotonous variation of volume fraction of constituent phases is elimination of the stress discontinuities that are often encountered in laminated composites and accordingly avoidance of delamination-related problems.How all these aspects can be improved and what the mechanisms might be are popular topics which have received much attention from researchers.Wang and Noda [2] investigated the thermally induced fracture of a functionally graded piezoelectric layer bonded to a metal.Ueda studied the fracture of an FGPM strip with a normal crack [3,4], a symmetrical FGPM strip with a center crack [5], and an FGPM strip with a two-dimensional crack [6,7].Li and Weng [8] solved the problem of an FGPM strip containing a finite crack normal to boundary surfaces.Hu et al. [9] studied the problem of a crack located in a functionally graded piezoelectric interlayer between two dissimilar homogenous piezoelectric half planes.Rao and Kuna [10] presented an interaction integral method for computing stress intensity factors (SIFs) and the electric displacement intensity factor (EDIF) in FGPM under thermoelectromechanical loading.Borrelli et al. [11] used the energy-decay inequality technique to analyze the decay behavior of end effects in antiplane shear deformation in piezoelectric solids and FGPMs.Zhong and Shang [12] developed an exact solution for a functionally graded piezothermoelectric rectangular plate.Dai et al. [13] conducted a theoretical study of electromagnetoelastic behavior 2 Mathematical Problems in Engineering for an FGPM cylinder and sphere.They then extended their solutions to include thermal effects [14].Zhong and Yu [15] analyzed the FGPM beam with arbitrarily graded material properties along the beam thickness direction.Based on the layerwise finite element model, Shakeri and Mirzaeifar [16] performed a static and dynamic analysis of a thick FGM plate with piezoelectric layers.Wang et al. [17] analytically investigated the axisymmetric bending of circular plates whose material properties varied with the thickness.Using the Fourier transform technique, Chue and Yeh [18] developed a system of singular integral equations for angle cracks in two bonded FGPMs under antiplane shear.Chen and Bian [19] studied the wave propagation characteristics of an axially polarized, functionally graded, piezoceramic cylindrical transducer submerged in an infinite fluid medium.Ueda [20] addressed the problem of two coplanar cracks in an FGPM strip under transient thermal loading.Ben Salah et al. [21] examined the propagation of ultrasonic guided waves in FGPMs.Wang et al. [22] studied the singularity behavior of electroelastic fields in a wedge with angularly graded piezoelectric material (AGPM) under antiplane deformation.Chue and Yeh [23] extended the results of the case of two arbitrarily oriented cracks in two bonded FGM strips.
In the field of engineering technology, accurate solutions computed with analytic methods are only available in problems with relatively simple equations and regular geometry.As for most problems, especially in the case of complex engineering, analytic solutions cannot be computed.Therefore, research has been going on for many years that has led to the development of another approach-the numerical method.With the rapid development and wide application of computers, numerical analysis and theoretical investigation and experimental investigation are considered as the three major research approaches.As one of the most effective tools for the study of mechanics, the finite element and other numerical methods have been widely used in scientific research and engineering practice.FEM analysis of piezoelectric structures with a crack under dynamic electromechanical loading was presented by Enderlein et al. [24].A survey on numerical algorithms for crack analyses in piezoelectric structures to be used with FEM for determining fracture parameters was presented by Kuna [25].Béchet et al. [26] presented an application of XFEM to the analysis of fracture in piezoelectric materials.Nanthakumar et al. [27] analyzed the multiple flaws in piezoelectric structures using XFEM and level sets.Sharma et al. [28] analyzed a subinterface crack in piezoelectric bimaterials with XFEM.Bouvier et al. [29] studied the inverse problems in structural analysis: application to atherosclerotic plaque elasticity reconstruction by using the XFEM.
Compared with the extended finite element method, however, element-free method has a unique feature in solving the problems of crack growth, plastic flow of materials, geometric distortion and phase transition, and singularity.The notable feature of this method is that, in establishing the discrete equation, it does not need mesh but only needs to arrange discrete points in the global domain.Thus, it not only avoids the complicated process of mesh formation, but also greatly reduces the influence of mesh distortion.Various meshless methods have been applied to the analysis of smart materials and structures such as the meshless point collocation method (PCM) [30], the point interpolation meshfree method (PIM) [31], a novel truly hybrid meshlessdifferential order reduction method (hM-DOR) [32], and the local Petrov-Galerkin method (MLPG) [33].Among these meshless methods, the element-free Galerkin method (EFGM) [34], developed by Belytschko et al., has good compatibility and stability and will not have the problem of shear locking of volume even in adopting linear primary function.And fast convergence speed and high accuracy can be achieved.EFGM was widely applied in fracture mechanics.Rabczuk and Belytschko analyzed the problem of a three-dimensional large deformation meshfree method for arbitrary evolving cracks [35].Rabczuk et al. deduced a simplified meshfree method for shear bands with cohesive surfaces [36].
In this study, a type of electromechanical-coupling enriched element-free Galerkin method based on the enriched EFG methods [37,38] is developed.Enriched terms were introduced into the approximation function of the conventional EFGM to describe the displacement and electric fields near the crack.Compared with the conventional EFGM, this method only needs a small domain to describe the crack-tip singular field.Furthermore, the domain of the node is not affected by the crack without using the visibility method and diffraction methods.

Basic Equations for Two-Dimensional FGPM
The governing equations and the boundary conditions of FGPM are briefly given.

Constitutive Equations.
Consider in which C  , e  , and   are the elastic, piezoelectric, and dielectric constants, respectively.  ,   , E  , and D  are the stress tensor, strain tensor, electric field, and electrical displacements.
Strains are related to displacements by the expression where u  and u  are, respectively, the displacements in and -directions.
The electric field is related to the electric potential by the expression Types of Gradation.The properties of FGM are usually assumed to have the same functions of certain space coordinates.Exponential material gradation is commonly used graded forms.
All the material constants including elastic constants, piezoelectric parameters, and dielectric constants follow the exponential law: where M  represents material constants such as C  , e  , or   , M 0  is the corresponding value at the plane  3 = 0, and  denotes a material graded parameter.
The field equations of electroelasticity are reduced to two-dimensional form in the special cases: plane strain.Considering a transversely isotropic FGPM, according to (1), the - plane is the isotropic plane, and one can employ either the - or - plane for the study of plane electromechanical phenomena.The plane strain conditions require that By substitution of ( 10) into (1), we have or inversely where   ,   , and   are the reduced material constants.They are related to the elastic compliance tensor   , the piezoelectric tensor   , and the dielectric impermeability tensor   by the following relations [39]: ( Boundary Conditions.In electroelasticity theory, mechanical boundary conditions are formulated just as in classical elasticity theory.The electric boundary conditions are, however, still controversial.The first attempt to define the electric boundary conditions over crack faces was done by Parton [40].He assumed that although the magnitude of the normal electrical displacement component at the crack face was very small, the electrical displacement was continuous across the crack faces.He used the following electric boundary conditions: Later, Hao and Shen [41] improved the above assumption by taking the electric permeability of air in the crack gap into consideration.In addition to (14), they presented an equation for the boundary condition at crack faces: where   is the permittivity of air.However, (15) has remained disregarded for a long time due to its complex mathematical treatment.
As pointed out by Suo et al. [42], the above assumption is not physically practical as there will clearly be a potential drop across the lower capacitance crack.This is particularly true for those piezoelectric ceramics with permittivity 10 3 times higher than that in the air.For this reason, Deeg [43] Crack line proposed another set of electric boundary conditions over crack faces: Equation ( 16) is derived from the constitutive equation   =     .This is equivalent to having the crack surfaces free of surface charge which is the electrical boundary condition.Thus the electric displacement vanishes in the environment.

Electromechanical-Coupling Enriched Element-Free Galerkin Method
Displacement u and electric potential  were adopted as basic field quantities for the solution of the enriched element-free Galerkin method, and the displacement and electric potential interpolation of a typical point x were conducted as follows: where   (x) is the MLS shape function at the point x.Ω is the support domain.Ω Γ is a set of all nodes whose support is cut by the crack.The set Ω  is a set of all nodes that lie within a fixed region around the crack tip.  and   are the additional degrees of freedom of the displacement.  and   are the additional degrees of freedom of the electric potential.
The first terms of ( 15) and ( 16) are the conventional EFGM approximation which simulates the displacement and electric potential fields.The second term is the displacement and electric potential approximation function of the nodes in Ω Γ .(x) is the Heaviside function: where x * is the projection of point x on the crack.
The last term in (15) and ( 16) is the displacement and electric potential approximation function of the nodes in Ω  .() is the branch functions given by where  and  are polar coordinates in the local coordinate system.
The set Ω Γ includes the nodes whose support contains point x or is cut by the crack.The set Ω  is nodes whose support contains point x and crack tip x tip (see Figure 1).
The expression for the potential energy in the domain Ω is given by where f  = { f q }, T  = {T/q}, u  = {u/}, and   = {     } are the generalized volume force vector, the given generalized boundary force vector, the generalized boundary displacement vector, and the generalized penalty function coefficient vector, respectively.
In the expressions above, are the generalized strain matrix and generalized elastic matrices.
When the system reaches equilibrium, the potential energy of the system is minimized.The use of the variation principle on ( 9) and ( 10) leads to the equilibrium equation: The standard discrete equations is obtained, For the enriched terms, ( , ] , and K   = (K   )  , ,  = , , , , , .

A Modified 𝐽-Integral for Functionally Graded Piezoelectric Materials
Owing to the influence of material inhomogeneity in functionally graded piezoelectric materials, the standard integral loses the path-independence.This section presents a modified -integral, which retains the features of the pathindependence.
The -integral for piezoelectric material was deduced by Pak [44] as where ℎ is the electric enthalpy, Γ is an arbitrary contour around the crack, and   is the unit outward normal component along path Γ as shown in Figure 2: Crack face (−) Note that the crack is considered to be impermeable and traction-free.
Numerical determination of the -integral is very difficult because the contour, surface, and volume integrals should be defined precisely in the grid.The closed integration path is transformed into an equivalent domain integral (EDI) by applying Gauss's integral theorem, as can be seen in Figure 3.According to the path independence, (24) could be written as where  = Γ + Γ + + Γ − − Γ  ,   is an outward normal vector, and  is a weighting function which must be continuous and meet the following conditions: For the functionally graded piezoelectric material, the influence of material inhomogeneity causes the standard integral to lose the characteristic of the path-independence.A modified -integral for calculating the energy release rates of functionally graded piezoelectric material was presented.Substitute ( 7) into (17), Let   u  = u  0 , and ( 30) is found that where (33) represents the equilibrium equation of the piezoelectric plate whose material constants are the corresponding values at the plane  = 0 in FGPMs.From ( 32) and ( 33), it can be observed that u  0 /u  =   .Equations ( 3)-( 6) could be rewritten into Substituting ( 32)-( 35) into ( 28), it is found that where  0  denotes the -integral of the piezoelectric plate whose material constants are the corresponding values at the plane  3 = 0 in FGPMs.The standard -integral has the characteristic of the path-independence.So does the  FGPM .

Numerical Example
Case 1.The centrally cracked infinite FGPM plate is shown in Figure 4.The height (2ℎ) and width (2) are both 0.4 m and the crack length (2) is 0.02 m.The width of the plate is 20 times the crack length which can be considered to be large enough to simulate an infinite plate with a crack.The positive  3 direction is defined as the polarization direction.The remote uniform distributed stress load is  ∞ = 10 6 Pa.The remote uniform distributed electric field intensity is E ∞ = −10 5 V/m.The crack is considered to be permeable and traction-free.The positive  3 direction is defined as the gradient direction of FGPM.The electron micrograph picture of the crack is shown in Figure 5.The material gradation function:   =  0    ,  = −1, −0.5, 0, 0.5, 1.
The material properties at  3 = 0 were as follows: As the consequence of the symmetry of the problem, symmetric boundary conditions are imposed; only half of the plate needs to be modeled.In this process, 41 × 81 nodes are uniformly distributed in the computational domain, using 40×80 background grids.In each grid, 4×4 Gaussian integral and cubic spline weight function are adopted.If a background mesh is present, nodes and the vertices of the integration usually coincide (as in conventional FEM meshes).When cell structures are utilized, a regular array of domains is created, independently of the particle position.It is worth noting that the  matrix at a Gauss point is composed of two parts: the standard and the enriched part, where the standard part is always computed and the enriched part is only computed if, in the nodes whose supports cover Gauss point, there exist enriched nodes.
Owing to the existence of additional degrees of freedom, the assembly procedure of stiffness matrix should be revised.Virtual nodes are used to deal with these additional degrees of freedom.One virtual node is added at a Heaviside function enriched node and four virtual nodes are added at a tip enriched node, as depicted in Figure 6.The number of these virtual nodes starts from the total number of true nodes.For example, if there are six nodes numbered from one to six, the third node and fifth node are near tip enriched nodes, and the second node, fourth node, and sixth node are enriched with the Heaviside function.Then, we have 6 + 3 × 1 + 2 × 4 = 17 nodes.And then we add a virtual node numbered 7 at the second node, add four virtual nodes numbered 8, 9, 10, and 11 at the third node, add a virtual node numbered 12 at the fourth node, add four virtual nodes numbered 13, 14, 15, and 16 at the fifth node, and add a virtual node numbered 17 at the sixth node.
The variation of displacements,  1 and  3 , electric potential, , stresses,  1 and  3 , electric fields,  1 and  3 , and electric displacements,  1 and  3 , are shown in Figure 7.The material property gradient index  takes five values: −1, −0.5, 0, 0.5, and 1.As shown in Figure 7, Figure 7(a) presents the displacement  3 of  line.For different Gauss integral, the modified -integral for FGPM has been computed.The results are listed in Table 1.The largest difference is only 2.08% when compared with the analytical solution which had been obtained from (36).Table 3 illustrates the accuracy of EFGM and EEFGM for calculating the crack of different lengths instead of the different gauss integral.
For a different number of nodes, the modified J-integrals for FGPM have been computed.21 × 41, 41 × 81, and 81 × 161 nodes are uniformly distributed in the computational domain, using 20×40, 40×80, and 80×160 background grids.In each grid, 4 × 4 Gaussian integrals are adopted.The results are listed in Table 2. From Table 2, it can be observed that the accuracy of the results by EEFGM has been greatly improved compared with the results by the conventional EFGM.In this process, 41 × 81 nodes are uniformly distributed in the computational domain, using 40 × 80 background grids.In each grid, 4 × 4 Gauss integral and cubic spline weight function are adopted.
The modified -integral for FGPM has been computed as the crack is under different length.The results are listed in Table 2.The largest difference is only 3% when compared with the analytical solution.From the results listed in Table 2,  it can be observed that the J-integral results obtained by EEFGM have higher accuracy than those calculated by the conventional EFGM.
It is important to note that it can approximately simulate quasi-static crack propagation process with the increase of the crack length.

Conclusions
An alternative electromechanical-coupling enriched element-free Galerkin method is proposed by introducing extended terms into the approximation function of conventional element-free Galerkin method to describe the displacement and electric fields near the crack.The major advantage of the present method is that it only needs a small domain to describe the crack-tip singular field, and the domain of the node is not affected by the crack.It can improve the computational efficiency without the introduction of the visibility and the diffraction methods, compared with the conventional EFGM.A series of numerical examples for infinite square plate with central cracks of FGPM is solved, and the J-integrals are calculated.EEFGM has high precision compared with the conventional EFGM with the material property gradient index  taking five values.The values of -integrals agree well with the analytical solutions.EEFGM can be considered as an alternative numerical method for 0.1 0.3 0.5 0.7 0.9 0.0   electromechanical-coupling problems.It may be potentially attractive for extensions to dynamic fracture analysis of piezoelectric materials.

Figure 7 (
b) presents the displacement  1 of  line.

Figure 7 (
c) presents electric potential  of  line.

Figure 7 (
f) presents electric field  3 of  line.Figure7(g) presents electric field  3 of  line.Figure7(h) presents electric displacement  3 of  line.Figure7(i) presents electric displacement  1 of  line.

Figure 5 :Figure 6 :
Figure 5: The electron micrograph of the  3 direction and the crack.

Figure 7 :
Figure 7: Variation of physical quantities of  and  line.

Table 1 :
The modified -integral under different Gauss integral (N/m).

Table 2 :
The modified -integral under different number of nodes (N/m).

Table 3 :
The modified -integral under different Gauss integral (N/m).