An element-free Galerkin method for electromechanical coupled analysis in piezoelectric materials with cracks

Abstract

We present an element-free Galerkin method for electromechanical coupled fracture analysis in piezoelectric materials. Singularity terms were introduced into the approximation function of the new method to describe the displacement and electric fields near the crack. The new method requires a smaller domain to describe the crack-tip singular field compared with the finite element method. Then, we computed the J-integrals of piezoelectric materials and investigated the effects of crack length on the computational precision. Numerical examples were used to highlight the accuracy of the new method compared with the analytical solutions and finite element method.

Introduction

Piezoelectric materials are a type of functional materials with the piezoelectric effect that converts between electrical and mechanical energies. There are direct and inverse piezoelectric effects. The former, discovered first in quartz crystals by P. Curie and J. Curie in 1880, refers to the appearance of an electric field in a piezoelectric solid under mechanical loading. The second effect refers to the mechanical deformation of piezoelectric solids under the presence of an electric field. The existence of this mechanical behavior was predicted by G. Lippmann in 1881 using energy conservation and electric current conservation and experimentally verified by the Curie brothers a few months later. The discovery of the piezoelectric effect aroused strong interest from the scientific community. So far, abundant achievements have been made in the microscopic and macroscopic theories of piezoelectric devices and their applications in modern science and technology.^1–5

At present, the most widely used piezoelectric materials are piezoelectric ceramics and piezoelectric ceramic–polymer composites. Piezoelectric ceramics are typical brittle materials with low fracture toughness and high defect sensitivity. Thus, the concentration of stress fields and electric fields caused by defects (e.g. cracks, inclusions, and pores) which occurs during manufacture or service often makes these devices dysfunctional and even causes breakdown or fracture damage under electromechanical joint effects. This characteristic restricts further improvement in the performance of related devices and wider application of piezoelectric ceramics. Therefore, in the fields of solid-state mechanics and materials physics, it is urgent to further understand the physiomechanical mechanisms of fractures of piezoelectric materials, to make reliable analysis and prediction, and to uncover the underlying toughening mechanisms.

For >10 years, researchers in the fields of mechanics, physics, and material science had extensively studied the fracture performances of piezoelectric materials by means of theoretical analysis, numerical calculation, and experimental observation and achieved remarkable results. The theoretical framework about the fracture mechanics of piezoelectric materials has been preliminarily established.^6–13

In the field of engineering technology, analytic methods have provided accurate solutions only to problems with relatively simple equations and regular geometry. No analytic solutions can be obtained for most problems, especially in the case of complex engineering. Therefore, years of efforts have been devoted to develop another approach—the numerical method. The rapid development and wide application of computers in recent years have led to the presence of three major research approaches, including numerical analysis, theoretical investigation, and experimental investigation. As one of the most effective tools for mechanics study, the finite element and other numerical methods are widely used in scientific research and engineering practice. For instance, finite element method (FEM) was used to analyze piezoelectric structures with cracks under dynamic electromechanical loading. ¹⁴ Numerical algorithms for fracture analysis in piezoelectric structures were used with FEM to determine fracture parameters. ¹⁵ Extended FEM was applied to analyze fractures in piezoelectric materials. ¹⁶ Partition of unity methods for fractures allows for arbitrary crack growth without remeshing. These methods have been successfully extended to brittle fracture modeling.^17–19

Compared to FEM, however, element-free method is unique in solving the problems of fracture growth, plastic flow, geometric distortion, phase transition, and singularity of materials. Noticeably, in establishing the discrete equation, this method does not need element division but only has to arrange discrete points in the global domain. Thus, it avoids the complicated mesh formation and greatly reduces the impact of mesh distortion.

The moving least-squares (MLS) approximation ²⁰ is an important method to construct the shape function in element-free method. The shape function that is formed with the MLS approximation can obtain a very precise solution. One disadvantage of the MLS approximation is that the final algebra equation system is sometimes ill-conditioned. Thus, sometimes, we cannot obtain a good solution or even correctly obtain a numerical solution with the MLS approximation. Then, improved moving least-squares (IMLS) approximation was discussed by Cheng and colleagues.^21–23 The complex variable moving least-squares (CVMLS) approximation, which is an approximation of a vector function, has been developed and presented. ²⁴ To enhance its computational efficiency, improved complex variable moving least-squares (ICVMLS) approximation was studied. ²⁵ Various meshless methods have been applied to analyze smart materials and structures such as the element-free Galerkin method (EFGM), ²⁶ the meshless point collocation method, ²⁷ point interpolation meshless method ²⁸ and radial point interpolation method,^29,30 a novel truly hybrid meshless-differential-order reduction method, ³¹ and local Petrov–Galerkin method. ³² Among these methods, the EFGM^33–35 is featured with high compatibility, high stability, no shear lock of volume even in use of linear primary function, fast convergence, and high accuracy. And the improved EFGM based on the IMLS approximation was investigated by Zhang et al. ³⁶ Combining CVMLS approximation with the EFGM, the complex variable element-free Galerkin (CVEFG) method for elasto-plasticity problems was presented. ³⁷ Based on the EFGM and the ICVMLS approximation, the improved complex variable element-free Galerkin (ICVEFG) method for two-dimensional (2D) potential problems was presented. ³⁸ EFGM is widely used in fracture mechanics analysis.^39–42 A meshless method was developed for discontinuous mechanical problems.^43–48

In this study, we first introduced the piezoelectric constitutive equations, the MLS approximation, and the formulation of electromechanical coupled EFGM. Then, equivalent-domain J-integral for piezoelectric materials was presented, and numerical examples were used to demonstrate the accuracy and efficiency of the EFGM.

Piezoelectric constitutive equations

The electroelastic response of a piezoelectric body of volume Ω and regular boundary surface S is governed by the mechanical and electrostatic equilibrium equations

$${\mathbf{\sigma }}_{\mathrm{ij},j}+{\mathit{f}}_i=0\left(\mathbf{\Omega }\right)$$

(1)

$${\mathit{D}}_{i,i}-\mathit{q}=0\left(\mathbf{\Omega }\right)$$

(2)

where f _i is the mechanical body force, q is the electric body charge, σ _ij is the symmetric Cauchy stress tensor, and D _i is the electric displacement vector components.

The basic equations for piezoelectric materials include balance laws and constitutive relations. The constitutive relations, usually called piezoelectric equations in the theory of piezoelectricity, describe the interaction effect among stress, electric displacement, and electric field as follows

$${\mathbf{\sigma }}_{\mathrm{ij}}={\mathit{c}}_{\mathrm{ijkl}}^E{\mathbf{\epsilon }}_{\mathrm{kl}}-{\mathit{e}}_{\mathrm{kij}}{\mathit{E}}_k$$

(3)

$${\mathit{D}}_i={\mathit{e}}_{\mathrm{ikl}}{\mathbf{\epsilon }}_{\mathrm{kl}}+{\mathbf{\lambda }}_{\mathrm{ik}}^{\epsilon }{\mathit{E}}_k$$

(4)

where σ _ij and ε _ij are the components of a stress tensor and a strain tensor, respectively; E _k and D _i are the components of an electric field vector and an electric displacement vector, respectively; ${\mathit{C}}_{\mathrm{ijkl}}^E$ , e _kij , and ${\mathbf{\lambda }}_{\mathrm{ik}}^{\epsilon }$ are elastic stiffness, piezoelectric constant, and dielectric constant, respectively; the superscripts E and ε denote constant electric field and strain, respectively.

Strains are related to displacements as follows

$${\epsilon }_{\mathrm{ij}}=\frac{1}{2}(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})$$

(5)

$${\epsilon }_x=\frac{\partial u_x}{\partial x}$$

(6)

$${\epsilon }_z=\frac{\partial u_z}{\partial z}$$

(7)

$${\gamma }_{\mathrm{xz}}=\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}$$

(8)

where u _x and u _z are the displacements in the x- and z-directions, respectively.

The electric field is related to the electric potential as follows

$$E_i=-\frac{\partial \varphi }{\partial x_i}$$

(9)

Essential boundary conditions

$$u=\overline{u}\mathrm{or}\varphi =\overline{\varphi }\mathrm{on}{\Gamma }_e$$

(10)

Natural boundary conditions

$${\mathbf{\sigma }}_{\mathrm{ij}}{n}_j={\mathit{f}}_i\left(\mathrm{or}\right){\mathit{D}}_i{n}_i=-\mathit{g}\mathrm{on}{\Gamma }_n$$

(11)

where $\overline{u}$ , $\overline{\varphi }$ , f _i , g , and n_i are mechanical displacement, electric potential, surface force components, surface charge, and outward unit normal vector components, respectively.

The 2D matrix forms of mechanical and electrical constitutive equations for piezoelectric materials are

$$\left[\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_z}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xz}}}\end{array}\right]=\left[\begin{array}{ccc}{c}_{11}\hfill & c_{13}\hfill & 0\hfill \\ \multicolumn{1}{c}{c_{13}}& c_{33}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & c_{55}\hfill \end{array}\right]\left[\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_z}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xz}}}\end{array}\right]-\left[\begin{array}{cc}0\hfill & e_{31}\hfill \\ \multicolumn{1}{c}{0}& e_{33}\hfill \\ \multicolumn{1}{c}{e_{15}}& 0\hfill \end{array}\right]\left[\begin{array}{c}{E}_x\hfill \\ \multicolumn{1}{c}{E_z}\end{array}\right]$$

(12)

$$\left[\begin{array}{c}{D}_x\hfill \\ \multicolumn{1}{c}{D_z}\end{array}\right]=\left[\begin{array}{ccc}0\hfill & 0\hfill & e_{15}\hfill \\ \multicolumn{1}{c}{e_{31}}& e_{33}\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_z}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xz}}}\end{array}\right]-\left[\begin{array}{cc}{\lambda }_{11}^{\epsilon }\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\lambda }_{33}^{\epsilon }\hfill \end{array}\right]\left[\begin{array}{c}{E}_x\hfill \\ \multicolumn{1}{c}{E_z}\end{array}\right]$$

(13)

By substituting equations (5)–(9) into equations (12) and (13), we can recast equations (3) and (4) using mechanical displacement and electric potential. For a 2D piezoelectric problem, the equation governing electrical equilibrium is

$$\begin{array}{l}\left(e_{31}+e_{15}\right)u_{x,xz}+e_{15}{u}_{z,xx}+e_{33}{u}_{z,zz}\\ -{\lambda }_{11}^{\epsilon }{\varphi }_{,xx}-{\lambda }_{33}^{\epsilon }{\varphi }_{,zz}=0\end{array}$$

(14)

and the equations governing mechanical equilibrium are

$$\begin{array}{l}{c}_{11}{u}_{x,xx}+c_{55}{u}_{x,zz}+\left(c_{13}+c_{55}\right)u_{z,xz}\\ +\left(e_{31}+e_{15}\right){\varphi }_{,xz}=0\end{array}$$

(15)

$$\begin{array}{l}\left(c_{13}+c_{55}\right)u_{x,xz}+c_{33}{u}_{z,zz}+c_{55}{u}_{z,xx}\\ +e_{33}{\varphi }_{,zz}+e_{15}{\varphi }_{,xx}=0\end{array}$$

(16)

MLS approximation

We consider function $u(\mathbf{x})$ over the domain $\Omega \subseteq {\Re }^K$ , where K = 1, 2, or 3. Let ${\Omega }_{\mathbf{x}}\subseteq \Omega$ be a subdomain describing the neighborhood of a point, $\mathbf{x}\subseteq {\Re }^K$ , located in $\Omega$ . According to MLS, the approximation $u^h(\mathbf{x})$ of $u(\mathbf{x})$ is

$$u\left(\mathbf{x}\right)\cong u^h\left(\mathbf{x}\right)=\sum _{i=1}^m{p}_i\left(\mathbf{x}\right)a_i\left(\mathbf{x}\right)={\mathit{p}}^{\mathrm{T}}\left(\mathbf{x}\right)\mathit{a}\left(\mathbf{x}\right)$$

(17)

where $p_i(\mathbf{x})$ is a given monomial in the polynomial basis function in the space coordinates ${\mathbf{x}}^{\mathrm{T}}=[x,y]$ , m is the number of terms of monomials in the basis function, and $a_i(\mathbf{x})$ is a coefficient to be determined for $p_i(\mathbf{x})$ . For one-dimensional (1D) and 2D spaces, the linear basis functions are

$${\mathit{p}}^{\mathrm{T}}\left(\mathbf{x}\right)=\left[\begin{array}{cc}1\hfill & x\hfill \end{array}\right]m=2,\left(1\mathrm{D}\right)$$

(18)

$${\mathit{p}}^{\mathrm{T}}\left(\mathbf{x}\right)=\left[\begin{array}{ccc}1\hfill & x\hfill & y\hfill \end{array}\right]m=3,\left(2\mathrm{D}\right)$$

(19)

and the quadratic basis functions are

$${\mathit{p}}^{\mathrm{T}}\left(\mathbf{x}\right)=\left[\begin{array}{ccc}1\hfill & x\hfill & x^2\hfill \end{array}\right]m=3,\left(1\mathrm{D}\right)$$

(20)

$${\mathit{p}}^{\mathrm{T}}\left(\mathbf{x}\right)=\left[\begin{array}{cccccc}1\hfill & x\hfill & y\hfill & x^2\hfill & \mathrm{xy}\hfill & y^2\hfill \end{array}\right]m=6,\left(2\mathrm{D}\right)$$

(21)

For piezoelectric problems, it is advisable to use the regular enriched basis functions stemming from the isotropic elasticity. It should be mentioned that similar results have been obtained independently using alternative regular enriched basis functions for cracks in confined plasticity problems. ⁴⁹ The regular enriched basis functions incorporate the radial and angular behaviors of the 2D asymptotic crack-tip displacement field

$$\begin{array}{l}{\mathit{p}}^{\mathrm{T}}\left(\mathit{x}\right)=[1xy\sqrt{r}\cos \frac{\theta }{2}\sqrt{r}\sin \frac{\theta }{2}\sqrt{r}\sin \frac{\theta }{2}\sin \theta \sqrt{r}\cos \frac{\theta }{2}\sin \theta ]\\ m=7\end{array}$$

(22)

where r and θ are polar coordinates with the crack tip as the origin.

$a_i(\mathbf{x})$ in equation (17) is a function of $\mathbf{x}$ , which can be determined by minimizing a function of weighted residual

$$J\left(\mathbf{x}\right)=\sum _{I=1}^n{w}_I\left(\mathbf{x}\right){[{\mathit{p}}^{\mathrm{T}}\left({\mathbf{x}}_I\right)\mathit{a}\left(\mathbf{x}\right)-u_I]}^2={\mathit{C}}^{\mathrm{T}}\mathit{WC}$$

(23)

where $\mathit{C}=\mathit{Pa}-\mathit{u}$ ; $w_I(\mathbf{x})$ is the weight function associated with node I, such that $w_I(\mathbf{x})>0$ for all $\mathbf{x}$ in the support ${\Omega }_{\mathbf{x}}$ of $w_I(\mathbf{x})$ and otherwise, $w_I(\mathbf{x})=0$ ; n is the number of nodes in ${\Omega }_{\mathbf{x}}$ for which $w_I(\mathbf{x})>0$ .

Moreover

$$\mathit{P}=\left[\begin{array}{cccc}{p}_1\left({\mathbf{x}}_1\right)\hfill & p_2\left({\mathbf{x}}_1\right)\hfill & \cdots \hfill & p_m\left({\mathbf{x}}_1\right)\hfill \\ \multicolumn{1}{c}{p_1\left({\mathbf{x}}_2\right)}& p_2\left({\mathbf{x}}_2\right)\hfill & \cdots \hfill & p_m\left({\mathbf{x}}_2\right)\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{p_1\left({\mathbf{x}}_n\right)}& p_2\left({\mathbf{x}}_n\right)\hfill & \cdots \hfill & p_m\left({\mathbf{x}}_n\right)\hfill \end{array}\right]=\left[\begin{array}{c}{\mathit{p}}^{\mathrm{T}}\left({\mathbf{x}}_1\right)\hfill \\ \multicolumn{1}{c}{{\mathit{p}}^{\mathrm{T}}\left({\mathbf{x}}_2\right)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{p}}^{\mathrm{T}}\left({\mathbf{x}}_n\right)}\end{array}\right]$$

(24)

$$\mathit{W}=\left[\begin{array}{cccc}{w}_1\left(\mathbf{x}\right)\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& w_2\left(\mathbf{x}\right)\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & w_n\left(\mathbf{x}\right)\hfill \end{array}\right]$$

(25)

$$\mathit{a}\left(\mathbf{x}\right)={\left[\begin{array}{cccc}{a}_1\left(\mathbf{x}\right)\hfill & a_2\left(\mathbf{x}\right)\hfill & \cdots \hfill & a_m\left(\mathbf{x}\right)\hfill \end{array}\right]}^{\mathrm{T}}$$

(26)

$$\mathit{u}={\left[{\begin{array}{cccc}{u}_1\hfill & u_2\hfill & \cdots \hfill & u\hfill \end{array}}_n\right]}^{\mathrm{T}}$$

(27)

The stationarity of $J(\mathbf{x})$ with respect to $a_j(\mathbf{x})$ can be written as

$$\frac{\partial J\left(\mathbf{x}\right)}{\partial a_j\left(\mathbf{x}\right)}=2\frac{\partial {\mathit{C}}^{\mathrm{T}}}{\partial a_j\left(\mathbf{x}\right)}\mathit{WC}=2{\mathit{P}}^{\mathrm{T}}\mathit{WC}=0$$

(28)

which results in the following equation system

$${\mathit{P}}^{\mathrm{T}}\mathit{WPa}\left(\mathbf{x}\right)={\mathit{P}}^{\mathrm{T}}\mathit{Wu}$$

(29)

where

$$\mathit{A}\left(\mathbf{x}\right)={\mathit{P}}^{\mathrm{T}}\mathit{WP}=\sum _{I=1}^n{w}_I\left(\mathbf{x}\right){\mathit{p}}^{\mathrm{T}}\left({\mathbf{x}}_I\right)\mathit{p}\left({\mathbf{x}}_I\right)$$

(30)

$$\mathit{B}\left(\mathbf{x}\right)={\mathit{P}}^{\mathrm{T}}\mathit{W}=\left[\begin{array}{ccccc}{\mathit{B}}_1\left(\mathbf{x}\right)\hfill & \cdots \hfill & {\mathit{B}}_I\left(\mathbf{x}\right)\hfill & \cdots \hfill & {\mathit{B}}_n\left(\mathbf{x}\right)\hfill \end{array}\right]$$

(31)

$${\mathit{B}}_I\left(\mathbf{x}\right)=\mathit{p}\left({\mathbf{x}}_I\right)w_I\left(\mathbf{x}\right)$$

(32)

Then, $\mathit{a}(\mathbf{x})$ from equation (29) is solved as follows

$$\mathit{a}\left(\mathbf{x}\right)={\mathit{A}}^{-1}\left(\mathbf{x}\right)\mathit{B}\left(\mathbf{x}\right)\mathit{u}$$

(33)

The EFGM for electromechanical coupled problems

Displacement u and electric potential $\mathbf{\varphi }$ were adopted as basic field quantities for the solution of the EFGM, and the interpolation of u and $\mathbf{\varphi }$ at a typical point x were conducted as

$$u_{\mathrm{iI}}\left(\mathbf{x}\right)=\sum _{I\in \Omega }{\Phi }_I\left(\mathbf{x}\right)u_{\mathrm{iI}}=\mathbf{\Phi }\left(\mathbf{x}\right)\mathit{u}$$

(34)

$${\varphi }_I\left(\mathbf{x}\right)=\sum _{I\in \Omega }{\Phi }_I\left(\mathbf{x}\right){\varphi }_I=\mathbf{\Phi }\left(\mathbf{x}\right)\mathbf{\varphi }$$

(35)

where Ω is the support domain and Φ(x) is the MLS shape function at point x

$$\Phi \left(\mathrm{x}\right)={\left[\begin{array}{cccc}{\phi }_1\left(\mathrm{x}\right)& {\begin{array}{cc}\cdots & \phi \end{array}}_I\left(\mathrm{x}\right)& \cdots & {\phi }_n\left(\mathrm{x}\right)\end{array}\right]}^{\mathrm{T}}=p^{\mathrm{T}}\left(\mathrm{x}\right)A^{-1}\left(\mathrm{x}\right)B\left(\mathrm{x}\right)$$

(36)

$${\phi }_I\left(\mathbf{x}\right)={\mathit{p}}^{\mathrm{T}}\left(\mathbf{x}\right){\mathit{A}}^{-1}\left(\mathbf{x}\right){\mathit{B}}_I\left(\mathbf{x}\right)$$

(37)

The potential energy in the domain Ω is expressed as follows

$$\prod ={\displaystyle \underset{\Omega }{\int }\frac{1}{2}{{\mathrm{\epsilon }}^{\prime }}^{\mathrm{T}}\mathbf{G}{\mathrm{\epsilon }}^{\prime }\mathrm{d}\Omega }-{\displaystyle \underset{\Omega }{\int }{\mathbf{f}}^{\prime }{\mathbf{u}}^{\prime }\mathrm{d}\Omega }-{\displaystyle {\int }_{{\Omega }^{\sigma }}{\mathbf{T}}^{\prime }{\mathbf{u}}^{\prime }\mathrm{d}\Gamma }+{\displaystyle \underset{{\Omega }^u}{\int }\frac{1}{2}{\alpha }^{\prime }\left({\mathbf{u}}^{\prime }-\overline{{\mathbf{u}}^{\prime }}\right)\mathrm{d}\Gamma }$$

(38)

where

$$\mathit{f}\prime =\left\{\begin{array}{c}\mathit{f}\hfill \\ \multicolumn{1}{c}{\mathit{q}}\end{array}\right\}$$

$$\mathit{T}\prime =\left\{\begin{array}{c}\overline{\mathit{T}}\hfill \\ \multicolumn{1}{c}{\overline{\mathit{q}}}\end{array}\right\}$$

$$\overline{\mathit{u}\prime }=\left\{\begin{array}{c}\overline{\mathit{u}}\hfill \\ \multicolumn{1}{c}{\overline{\varphi }}\end{array}\right\}$$

$$\mathbf{\alpha }\prime =\left\{\begin{array}{c}{\alpha }^{\mathrm{u}}\hfill \\ \multicolumn{1}{c}{{\alpha }^{\varphi }}\end{array}\right\}$$

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. Due to the great difference between mechanical field and electric field coefficient matrices, the two penalty functions must be selected separately.

In equation (38)

$$\mathbf{\epsilon }\prime =\left\{\begin{array}{c}{\epsilon }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{E_x}\\ \multicolumn{1}{c}{E_z}\end{array}\right\}$$

and

$$\mathit{G}=\left[\begin{array}{ccccc}{c}_{11}\hfill & c_{13}\hfill & 0\hfill & 0\hfill & e_{31}\hfill \\ \multicolumn{1}{c}{c_{13}}& c_{33}\hfill & 0\hfill & 0\hfill & e_{33}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & c_{44}\hfill & e_{15}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & e_{15}\hfill & {\lambda }_{11}\hfill & 0\hfill \\ \multicolumn{1}{c}{e_{31}}& e_{33}\hfill & 0\hfill & 0\hfill & {\lambda }_{33}\hfill \end{array}\right]$$

are the generalized strain matrix and generalized elastic matrix, respectively.

When the system reaches equilibrium, its potential energy is minimized. Using the variation principle in equation (38) leads to the following equilibrium equation

$$\mathit{Ku}\prime =\mathit{f}\prime$$

(39)

where $\mathit{K}$ is a global stiffness matrix, $\mathit{u}\prime$ is a generalized displacement vector, and $\mathit{f}\prime$ is a generalized load vector.

In this work, a weight function based on Gaussian distribution is used

$$w_I\left(\mathbf{x}\right)=\{\begin{array}{cc}\frac{{(1+{\beta }^2{r}^2)}^{-\left(\frac{1+\beta }{2}\right)}-{(1+{\beta }^2)}^{-\left(\frac{1+\beta }{2}\right)}}{1-{(1+{\beta }^2)}^{-\left(\frac{1+\beta }{2}\right)}},\hfill & r\le 1\hfill \\ \multicolumn{1}{c}{0,}& r>1\hfill \end{array}$$

(40)

where $\beta$ is a shape-controlling parameter and $r=d_I/d_{\mathrm{mI}}$ . Moreover, $d_I=\Vert \mathbf{x}-{\mathbf{x}}_I\Vert$ is the distance from a sampling point $\mathbf{x}$ to a node ${\mathbf{x}}_I$ ; $d_{\mathrm{mI}}$ is the influence domain of node I such that

$$d_{\mathrm{mI}}=\mathrm{scale}\times c_I$$

(41)

where scale is a scaling parameter and $c_I$ is a characteristic nodal spacing distance such that node I has sufficient neighbors for regularity of $\mathit{A}(\mathbf{x})$ in equation (30).

To avoid any fracture-caused discontinuity in the shape functions, we used a diffraction method to modify $d_I$ in the weight function (Figure 1) as follows

$$d_I={\left(\frac{s_1+s_2\left(\mathbf{x}\right)}{s_0\left(\mathbf{x}\right)}\right)}^{\lambda }{s}_0\left(\mathbf{x}\right)$$

(42)

where $s_1=\Vert {\mathbf{x}}_I-{\mathbf{x}}_c\Vert$ , $s_2(\mathbf{x})=\Vert \mathbf{x}-{\mathbf{x}}_c\Vert$ , and $s_0(\mathbf{x})=\Vert \mathbf{x}-{\mathbf{x}}_I\Vert$ ; ${\mathbf{x}}_c={\left[\begin{array}{cc}{x}_{1c}\hfill & y_{1c}\hfill \end{array}\right]}^{\mathrm{T}}$ is a vector of coordinates at the crack tip; $1\le \lambda \le 2$ is a diffraction parameter.

Figure 1.

Diffraction method.

Equivalent-domain J-integral for piezoelectric materials

The J-integral defined in Rice ⁵⁰ is one parameter of the nonlinear fracture problem and also a path-independent integral (Figure 2)

$$J=\underset{\Gamma }{\int }(W{n}_1-{\sigma }_{\mathrm{ij}}{n}_j\frac{\partial u_i}{\partial x})\mathrm{d}\Gamma$$

(43)

where $W=(1/2){\sigma }_{\mathrm{ij}}{\epsilon }_{\mathrm{ij}}$ is the energy density, Γ is an arbitrary contour around the fracture, and n_i is the unit outward normal component along the path Γ.

Figure 2.

J-integral path around a crack tip.

The J-integral for piezoelectric materials was deduced as follows ⁵¹

$$J=\underset{\Gamma }{\int }(H{n}_1-{\sigma }_{\mathrm{ij}}{n}_j{u}_{i,1}+D_i{n}_i{E}_1)\mathrm{d}s$$

(44)

where H is the electric enthalpy

$$H=\frac{1}{2}{C}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{ij}}{\epsilon }_{\mathrm{kl}}-\frac{1}{2}{\lambda }_{\mathrm{ij}}{E}_i{E}_j-e_{\mathrm{ikl}}{\epsilon }_{\mathrm{kl}}{E}_i$$

(45)

$${\sigma }_{\mathrm{ij}}=\frac{\partial H}{\partial {\epsilon }_{\mathrm{ij}}}=C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}-e_{\mathrm{kij}}{E}_k$$

(46)

$$D_i=-\frac{\partial H}{\partial E_i}=e_{\mathrm{ikl}}{\epsilon }_{\mathrm{kl}}+{\lambda }_{\mathrm{ik}}{E}_k$$

(47)

Note that the crack is considered as impermeable and traction free. Numerical determination of the J-integral is very difficult because the contour, surface, and volume integrals in the grids should be defined precisely. The closed integration path is transformed to an equivalent-domain integral (EDI) by applying Gauss’ integral theorem (Figure 3). According to the rule of path independence, equation (44) could be rewritten as

$$\begin{array}{l}J=-{\displaystyle \underset{S}{\int }\left(h{n}_1-{\sigma }_{ij}{n}_j{u}_{i,1}+D_i{n}_i{E}_1\right)q\mathrm{d}\Gamma }\\ +{\displaystyle \underset{\Gamma +{\Gamma }^{+}+{\Gamma }^{-}}{\int }\left(h{n}_1-{\sigma }_{ij}{n}_j{u}_{i,1}+D_i{n}_i{E}_1\right)q\mathrm{d}\Gamma }\end{array}$$

(48)

where S = Γ + Γ⁺ + Γ⁻ − Γ _ε , n_i is an outward normal vector, and q is a weighting function that must be continuous and meet the following condition

$$q=\{\begin{array}{cc}0\hfill & \mathrm{on}\Gamma \hfill \\ \multicolumn{1}{c}{1}& \mathrm{on}{\Gamma }_{\epsilon }\hfill \end{array}$$

(49)

Figure 3.

Equivalent-domain J-integral.

Numerical examples

Infinite PZT-4 square plate with a central crack

The centrally cracked infinite PZT-4 plate (height 2h = 40 m, width 2w = 40 m, and crack length 2a = 2 m) is shown in Figure 4. The plate width which is 20 times of crack length can be considered as large enough for simulation of an infinite cracked plate. As shown in Figure 4, the positive x₃ direction is defined as the polarization direction, σ ^∞ is the remote uniformly distributed stress load, and E ^∞ is the remote uniformly distributed electric field intensity. The crack is considered as impermeable and traction free.

Figure 4.

Piezoelectric plate with central crack.

The material properties of PZT-4 are as follows: c₁₁ = 13.9 × 10¹⁰ N/m, c₁₂ = 7.78 × 10¹⁰ N/m, c₁₃ = 7.43 × 10¹⁰ N/m, c₃₃ = 11.3 × 10¹⁰ N/m, c₄₄ = 2.56 ×10¹⁰ N/m, e₃₁ = −6.98 C/m², e₃₃ = 13.84 C/m², e₁₅ = 13.44 C/m², ε₁₁ = 6.00 × 10⁻⁹ C/V m, and ε₃₃ = 5.47 × 10⁻⁹ C/V m.

Due to the symmetry of the plate, only half of the plate needs to be modeled. In this process, 41 × 81 nodes were uniformly distributed in the computational domain, and here, we used 40 × 80 background grids. In each grid, an 8 × 8 Gaussian integral was adopted.

To improve the precision, we added four nodes at the crack tip and added the number of nodes with the background grid coordinates at the crack (Figure 5).

Figure 5.

Additional nodes around the crack.

For different combinations of stress and electric field intensity, we computed the J-integral (EDI) with the integral area as shown in Figure 6. The results are listed in Table 1. The largest difference is only 0.12% compared with the real energy release rate that had been determined for a centrally cracked infinite PZT-4 plate. When the numbers of nodes were approximately equal and the crack-tip elements of FEM were locally refined, EFGM was still more precise than FEM (Table 1). Figure 7 shows the finite element mesh by FEM, which includes 1424 elements

$$J={10}^{-10}a(0.363{\sigma }_{\infty }^2+0.373{\sigma }_{\infty }{E}_{\infty }-138.3E_{\infty }^2)$$

(50)

Figure 6.

Equivalent-domain J-integral.

Table 1.

J-integrals for centrally cracked infinite PZT-4 plate (N/m).

σ ∞ (MPa)	E ∞ (100 kV/m)	Analytical	FEM 52	EFGM
1.0	1.0	−98.27	−98.10	−98.23
1.0	−1.0	−105.73	−105.60	−105.76
1.0	0.1	35.29	34.93	35.32
0.1	1.0	−137.60	−137.0	−137.71
0.1	−1.0	−138.31	−137.7	−138.48

FEM: finite element method; EFGM: element-free Galerkin method.

Figure 7.

Finite element mesh for the quadrant of the piezoelectric plate.

Infinite PZT-5H square plate with a central crack

The material of the plate in this section is PZT-5H, where the other conditions are the same as in section “Infinite PZT-4 square plate with a central crack.” The plate size, mechanics, electric load, and boundary conditions are shown in Figure 4.

The material properties of PZT-5H are as follows: c₁₁ = 12.6 × 10¹⁰ N/m, c₁₂ = 5 × 10¹⁰ N/m, c₁₃ = 8.41 × 10¹⁰ N/m, c₃₃ = 11.7 × 10¹⁰ N/m, c₄₄ = 2.30 ×10¹⁰ N/m, e₃₁ = −6.50 C/m², e₃₃ = 23.3 C/m², e₁₅ = 17.44 C/m², ε₁₁ = 150.3 × 10⁻¹⁰ C/V m, and ε₃₃ = 130 × 10⁻¹⁰ C/V m.

The analytical solution of J-integral in this plate was calculated as follows

$$J={10}^{-10}a(0.4248{\sigma }_{\infty }^2+0.6952{\sigma }_{\infty }{E}_{\infty }-389.44E_{\infty }^2)$$

(51)

The J-integral results obtained by EFGM are more accurate compared with those determined by FEM (Table 2).

Table 2.

J-integrals for centrally cracked infinite PZT-5H plate (N/m).

σ ∞ (MPa)	E ∞ (100 kV/m)	Analytical	FEM	EFGM
1.0	1.0	−340.01	−338.67	−341.31
1.0	−1.0	−353.91	−351.26	−352.82
1.0	0.1	39.28	38.34	40.21
0.1	1.0	−388.32	−387.45	−387.75
0.1	−1.0	−389.71	−388.25	−390.47

FEM: finite element method; EFGM: element-free Galerkin method.

Infinite PZT-5H square plate with central cracks of different lengths

The PZT-5H square plate is 20 m high (2h) and 20 m wide (2w) with the crack length (2a) of 0.5, 0.75, 1, 1.25, and 1.5 m (Figure 4). The J-integral was calculated for each crack. The computational domain involved 81 × 161 uniformly distributed nodes, of which 80 × 160 background grids were used. In each grid, a 4 × 4 gauss Gaussian integral was adopted.

The J-integral calculations for different crack lengths and different scales are shown in Figure 8. The results show that the errors at scale = 1.3 are large due to the small radius of the supporting zone. At scale = 1.7, 2.3, and 3.0, the EFGM calculated results agree well with the analytical solutions no matter how the crack length changes. Moreover, as the size of the influence domain increases, the accuracy of the results is improved but changes very little in general, which indicates that the method has high stability.

Figure 8.

J-integral of different crack lengths at (a) scale = 1.3, (b) scale = 1.7, (c) scale = 2.3, and (d) scale = 3.0.

A rectangular plate with central inclined crack

Here, a finite rectangular PZT-4 plate with a central crack inclined by 45° to the horizontal direction is considered (Figure 9). Other parameters are as follows: uniform tension σ = 1 GPa and electric displacement D = 1 C/m² in the y-direction; the ratio of crack length to width a/h = 0.1 (h = l m, a = 0.1 m).

Figure 9.

A rectangular plate with a central and inclined crack.

The high accuracy of EFGM is observed in Table 3. The comparison of relative errors shows that the EFGM results are more accurate than the FEM results.

Table 3.

Error comparison between two methods of the mechanical and electrical intensity factors for a rectangular plate with an inclined crack under coupled σ = 1 GPa and D = 1 C/m².

Method	Analytical	FEM	EFGM
K I (GPa/m1/2)	0.28025	0.2780	0.2787
Error (%)	–	0.81	0.55
K II (GPa/m1/2)	0.28025	0.2781	0.2789
Error (%)	–	0.78	0.48
KD (C/m3/2)	0.39633	0.3931	0.3941
Error (%)	–	0.82	0.56

FEM: finite element method; EFGM: element-free Galerkin method.

Electromechanical kinked crack

In this example, a kinked crack problem is analyzed. The intensity factors of a kinked crack in the infinite plate under electrical far-field load D = 1 C/m² are calculated and compared with the analytical solution and FEM results. The piezoelectric material is PZT-4. The kinked crack model is shown in Figure 10. A straight main crack of the length a continues into a straight crack branch of the length c = 0.25a, deviating from the main crack by the angle α = 30°. The poling axis of the piezoelectric material is oriented perpendicular to the main crack. The K_II corresponding analytical results are not given in Xu and Rajapakse. ⁵³

Figure 10.

Kinked crack in PZT-4.

Table 4 lists the corresponding field intensity factors, analytical results, and FEM results. The comparison of relative errors demonstrates that the EFGM results are more accurate than the FEM results.

Table 4.

The mechanical and electrical intensity factors at branch tip in PZT-4.

Method	Analytical	FEM	EFGM
K I (MPa/m1/2)	1.2533	1.2270	1.2345
Error (%)	–	2.1	1.5
KD (C/m3/2)	1.3160	1.2609	1.2712
Error (%)	–	4.2	3.4

FEM: finite element method; EFGM: element-free Galerkin method.

Conclusion

The wide application of piezoelectric structures into innovative technical areas has highlighted the problems of strength and reliability, which have to be carefully investigated. Sophisticated analysis of electromechanical properties of cracks is needed to quantitatively assess fracture and fatigue. The fracture mechanics approach for crack-like defects in piezoelectric materials reveals the electromechanical coupled field singularities. Effective numerical methods are needed to evaluate fracture behavior of cracks in arbitrary piezoelectric structures under combined electromechanical loading.

In this article, the EFGM was proposed for electromechanical crack analyses, and its accuracy was verified by four examples. The conclusions from numerical simulation are listed in the following:

1. Singularity terms which could describe the displacement and electric fields near the crack were introduced to the approximation function of the EFGM.

2. The new method is unique that it only uses a small domain to describe the crack-tip singular field.

3. A series of numerical examples for infinite square plate with cracks of piezoelectric materials were solved, and the J-integrals were calculated. EFGM is more precise than FEM. With the increase in crack length, the calculations of J-integrals agree well with the analytical solutions. EFGM can be an alternative numerical method for analysis of electromechanical coupling problems. It may be potentially extended to dynamic fracture analysis of piezoelectric materials.