- Email: [email protected]
Boundary element modeling of cracks in piezoelectric solids
Engineering Analysis with Boundary Elements 25 (2001) 771±781
www.elsevier.com/locate/enganabound
Boundary element modeling of cracks in piezoelectric solids R.K.N.D. Rajapakse*, X.-L. Xu Department of Civil and Geological Engineering, University of Manitoba, Winnipeg, Canada R3T 5V6 Received 19 June 2000; accepted 1 March 2001
Abstract This study focuses on the application of boundary element methods for linear fracture mechanics of two-dimensional piezoelectric solids. A complete set of piezoelectric Green's functions, based on the extended Lekhnitskii's formalism and distributed dislocation modeling, are presented. Special Green's functions are obtained for an in®nite medium containing a conducting crack or an impermeable crack. The numerical solution of the boundary integral equation and the computation of fracture parameters are discussed. The concept of crack closure integral is utilized to calculate energy release rates. Accuracy of the boundary element solutions is con®rmed by comparing with analytical solutions reported in the literature. The present scheme can be applied to study complex cracks such as branched cracks, forked cracks and microcrack clusters. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Piezoelectric; Impermeable crack; Boundary element; Conducting crack; Fracture mechanics
1. Introduction Piezoelectric actuators are widely used in the development of adaptive (smart) structures. The main drawback of piezoelectrics is their brittleness. Defects due to manufacturing process and complex electromechanical loading encountered in modern applications increase the likelihood of fracture and failure of piezoceramic actuators. A thorough understanding of fracture behavior of piezoelectric materials is therefore vital to the advancement of adaptive structures technology. A number of researchers have examined crack problems in piezoelectric materials using analytical methods. Different types of electrical boundary conditions on crack faces have been proposed. Mikhailov and Parton [12] considered the crack faces as electrically permeable, which led to a controversial conclusion that an electric ®eld has no contribution to crack propagation. Suo et al. [26] argued that this condition is not realistic, as there is an electric potential drop across the crack. Using an impermeable crack model, Sosa and Suo et al. [24,26] added an electric intensity factor to the well known stress intensity factors. [13,15] used the energy release rates to study fracture. A conducting crack model has also been considered [25]. Cracks in piezoelectric * Corresponding author. Address: Department of Mechanical Engineering, University of British Columbia, 2324 Main Mall, Vancouver, Canada V6T1Z4. Tel.: 11-604-822-0497; fax: 11-604-822-0944. E-mail address: [email protected] (R.K.N.D. Rajapakse).
actuators may emanate from an electrode edge. A conducting crack model is suitable when a conducting species migrates on the crack surfaces. Qin and Mai and Xu and Rajapakse [19,29], among others, considered crack branching problems. The above studies provided the theoretical basis for fracture mechanics of piezoelectrics. However, the study of fracture problems related to practical situations involving ®nite geometries and complex loading requires the development of numerical techniques. Some recent studies have been concerned with the development of boundary element methods for piezoelectrics. The ef®ciency of boundary element method is mainly attributed to the reduction in the dimensionality of the problem. Moreover, compared to the ®nite element method, the boundary element method is highly accurate and computationally ef®cient when dealing with linear fracture mechanics [2,3]. Indeed, fracture mechanics problems are one of the principal areas of application of boundary element method [4]. Green's functions (fundamental solutions) are essential for all boundary element methods. Superior numerical ef®ciency can be achieved when closed form Green's functions are available. Deeg [5] presented a comprehensive theoretical study of Green's functions, dislocations and cracks in piezoelectric solids under general anisotropy. However, his Green's functions are not in closed form due to the consideration of general anisotropy. Lee and Jiang and Rajapakse [9,20] derived two-dimensional piezoelectric Green's functions for materials with hexagonal symmetry by using double Fourier transformation and
0955-7997/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0955-799 7(01)00058-3
772
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Fourier transformation, respectively. Pan [14] used Deeg's formulation to rederive two-dimensional Green's functions for in®nite planes, half planes and bi-material planes. Dunn and Wienecke [6] derived three-dimensional Green's functions. Parton and Kudryavtsev [16] presented a boundary integral equation for linear piezoelectric problems. Lee and Jiang, Lu and Mahrenholtz, Qin and Qin and Lu [9,10,17,18], among others, also considered boundary element formulations for piezoelectric solids. Rajapakse [21] discussed three boundary integral formulations (the direct and indirect methods, and the ®ctitious stress-charge method) for linear piezoelectricity. However, only a few studies presented numerical implementations. Hill and Farris [8] conducted a three-dimensional boundary element analysis of piezoelectric solids by a numerically evaluating Green's functions presented by Deeg [5]. Xu and Rajapakse [28] developed a two-dimensional boundary element code based on the closed form Green's functions [20] to examine stress and electric ®eld concentration around voids. Pan [14] proposed a single-domain boundary element method for two-dimensional crack problems. Numerical results for ®eld intensity factors, crack opening displacement and electric potential jump across crack faces were presented for an inclined crack. The objective of this study is to extend the application of the conventional boundary element method to fracture mechanics of piezoelectric solids. A complete set of piezoelectric Green's functions are derived. Conventional Green's functions, which are electroelastic solutions due to concentrated line forces and electric charges in an in®nite homogeneous plane, are rederived by using the extended Lekhnitskii's formalism. Crack Green's functions are then derived by considering an in®nite piezoelectric plane containing a center crack. The crack faces can be either electrically impermeable or conducting. The application of these Green's functions in boundary element analysis and the computation of fracture parameters are then discussed. The accuracy of the boundary element solutions is con®rmed and selected numerical solutions for crack problems are presented.
2. Basic equations A plane piezoelectric medium is considered in Fig. 1. The medium is polarized along the z 0 -direction, which makes angle b with the z-direction. In the absence of body forces and body electric charges, the governing equations for the medium can be expressed as
s ij;j 0;
Di;i 0
1
where s ij i; j x; z and Di denote the stress tensor and electric displacement in the i-direction, respectively. The constitutive relations in the xz system can be
z
z’ (pol. direc.) β
Γ x
y (y’) Fig. 1. A plane piezoelectric medium.
expressed as 8 9 0 exx > a11 > > > < = B ezz B @ a12 > > > > : ; 2exz a13 (
Ex
)
Ez
2
a12 a22 a23
b11
b12
b21
b22
9 0 18 a13 > s xx > b11 > > = B C< B a23 C A> s zz > 1 @ b12 > > : ; a33 s xz b13 8 9 s xx > !> > > < = b13 d11 s zz 1 > > > b23 > d12 : ; s xz
1 b21 ( ) C Dx b22 C A Dz b23
d12 d22
!(
Dx
)
Dz
2 where e ij and Ei denote the strain tensor, and electric ®eld in the i-direction, respectively; aij, bij and dij denote twodimensional elastic, piezoelectric and dielectric constants, respectively. These material constants are functions of polarization angle b and material constants in Eq. (A1), and are different for plane stress and plane strain cases. The following boundary conditions are admissible on the boundary G of the medium u i u~ i
or
s ij nj t~ i
3
f f~
or
Di ni 2q~
4
where ui and f denote the displacement in the i-direction and electric potential, respectively; u~ i ; t~ i ; ni, f~ and q~ denote speci®ed displacement and traction components in the idirection, components of outward unit normal vector, speci®ed electric potential and speci®ed surface charge respectively. By utilizing the extended Lekhnitskii's formalism, the general solutions for plane piezo electrics can be expressed as [29] {ux ; uz ; f}T 2Re
3 X n1
{pn ; qn ; sn }T wn zn
5
where Re denotes real part of a complex valued quantity; pn, qn and sn are complex constants de®ned in Appendix A that depend on material properties and polarization angle b ; a
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
superscript T denotes transpose of a vector; zn x 1 mn z where mn n 1; 2; 3 are the roots of Eq. (A2); wn zn are complex potential functions to be determined from the loading and boundary conditions.
Fig. 2b involves a homogeneous medium subjected to a line force/charge, and has been solved previously. In Fig. 2c, a continuous generalized dislocation ®eld (distributed discontinuities in displacements and electric potential) with densities bj ( j 1,2,3) is utilized to model the impermeable crack. The resultant electroelastic ®elds on the crack line due to line force/charge and distributed dislocation ®eld should satisfy the boundary conditions given by Eq. (7). The solution for a single generalized edge dislocation is sought ®rst. The dislocation ®eld bi is then obtained by enforcing boundary conditions of Eq. (7) and making use of Eq. (5). For the sake of brevity, the details of derivation are omitted. The ®nal results for the potential functions wn zn corresponding to the subproblem in Fig. 2c are obtained as
3. Green's functions 3.1. Green's functions for line force/charge Several previous studies have considered the derivation of conventional Green's functions for a piezoelectric in®nite plane. The present study utilizes the extended Lekhnitskii's formalism. Omitting the details, the ®nal results are given below. Consider a concentrated line force/charge P(P1, P2, P3), i.e. horizontal force P1, vertical force P2 and electric charge P3, applied at point z 0(x0, z0). The complete Green's functions are obtained by substituting the following potentials into general solutions given by Eq. (5)
wn zn
1 p ln zn 2 z0n kni Pi 2pi
w 0n zn
1 1 kp P 2pi zn 2 z0n ni i
kni Im 2pi " # pp 3 X z0n 2 2 a2 z2n 2 a2 1 z0n zn 2 a2 p £ vni ln knj Pj zn 2 z0n n1
wn zn
6
kni Im 2p i " # p ! 3 X z0n 2 2 a2 p 1 £ 1 2 p knj Pj vni zn 2 z0n z2n 2 a2 n1
w 0n zn
p where complex constants kni n; i 1; 2; 3 de®ned in Appendix A are functions of material properties and polarization angle b .
8
3.2. Crack Green's functions
where coef®cients kni and vni n; i 1; 2; 3 given in the Appendix A are functions of material properties and polarization angle b . Superposition of Eqs. (6) and (8) yields the complex potential functions corresponding to impermeable crack Green's functions. Now assume that the crack faces in Fig. 2a are conducting. Therefore
Availability of crack Green's functions allows the development of computationally ef®cient boundary element methods for crack problems. For the case of purely elastic solids, Ting [27] derived crack Green's functions based on Stroh's formalism. In this subsection, Green's functions for impermeable cracks and conducting cracks are derived based on the extended Lekhnitskii's formalism. The case of an impermeable crack is considered ®rst (Fig. 2a). The following boundary conditions hold on the crack faces
s zz 0;
s xz 0;
Dz 0; 2a # x # a
s zz 0;
7
z
Ex 0; 2a # x # a
z
z P
(a)
s xz 0;
P x
=
x (b)
9
The decomposition scheme and dislocation modeling are utilized again. Similar to the case of an impermeable crack, the generalized dislocation ®eld bi which satis®es the conducting boundary conditions can be obtained. The ®nal results for the potential functions wn zn corresponding to
The problem in Fig. 2a can be decomposed into two subproblems shown in Fig. 2b and c. The subproblem in
2a
773
+
bj (c)
Fig. 2. Decomposition used to derive crack Green's.
x
774
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Fig. 2c are k^ wn zn ni Im "2pi pp # 3 X z0n 2 2 a2 z2n 2 a2 1 z0n zn 2 a2 p £ v^ni ln knj Pj zn 2 z0n n1 k^ni Im "2pi # p ! 3 X z0n 2 2 a2 p 1 £ 1 2 p knj Pj v^ni zn 2 z0n z2n 2 a2 n1
w 0n zn
10 where coef®cients k^ni and v^ni n; i 1; 2; 3 given in the Appendix A are functions of the material properties and polarization angle b . Superposition of Eqs. (6) and (10) yields the complex potential functions corresponding to conducting crack Green's functions. 4. Boundary element method Based on the reciprocal relations for a piezoelectric solid, the following integral equation can be derived in the absence of body forces and charges [16] Z c x 0 ui x 0 Gji x; x 0 tj x 2 Hji x; x 0 uj x dG 11
where N is the total number of boundary elements, n is the total number of nodes on the l-th boundary element; tkj ; ukj ; qk ; fk are j-th component of nodal traction and displacement, electric charge and electric potential respectively of the k-th local node of the l-th boundary element; the p Jacobian J h is equal to 2x=2h2 1 2z=2h2 : Quadratic shape functions M k h are used in the present study. When the ®eld point x coincides with the node k, attention should be paid to evaluate singular integrals arising from Gij and Hij in Eq. (12). For conventional Green's functions in homogeneous piezoelectric media, based on the solutions in Section 3.2, Gij have ln r type singularity and Hij have 1=r type singularity. The inspection of crack Green's functions indicates that no additional singularity exists. Therefore, piezoelectric Green's functions have the same type of singularity as the elastic counterpart, and the integration scheme for elastic problems can be directly applied to calculate the singular integrals in Eq. (12). The numerical implementation and computation of boundary and internal ®eld variables are similar to the case of elastic materials. Details of numerical scheme for elastic problems are given by Banerjee [2]. While the present study focuses on the application of traditional boundary element method to piezoelectric solids, it is worth noting that Pan [14] proposed a single-domain boundary element method for piezoelectric crack problems.
G
where i x; z; q; u q 2f; tq 2q; and summation is implied on the index j j x; z; q: Gji and Hji are the Green's functions representing the electroelastic ®elds due to concentrated line force/charge in an in®nite medium: Gji x; x 0 and Hji x; x 0 denote the displacement and traction component at ®eld point x [ G in the j-direction j x; z due to a unit concentrated line load applied in the i-direction i x; z at point x 0 respectively; G qi x; x 0 and Hqi x; x 0 denote the electric potential and electric charge at ®eld point x due to unit concentrated line load applied in the i-direction at point x 0 respectively; Gjq x; x 0 and Hjq x; x 0 denote the displacement and traction components in the j-direction at ®eld point x due to unit electric line charge at point x 0 respectively; Gqq x; x 0 and Hqq x; x 0 denote the electric potential and electric charge at ®eld point x due to the unit electric charge at load point x 0 respectively. The coef®cient c x 0 can be calculated by the concept of rigid body movement. Upon discretizing the boundary G , Eq. (11) becomes " # ( N n Z X X 0 0 0 k Gji x; x M hJ h c x u i x DG l
l1
)
(
£ dG h t jk 2 " £
0
Hji x; x
n X k1
k1
Z
12
DG l
# k
M hJ h dG h
!
) u jk
5. Computation of fracture parameters Field intensity factors and energy release rates are commonly used fracture parameters. It has been shown that, for both impermeable and conducting cracks, p classical inverse square root type singularity 1= r exists for stresses andpelectric displacements (®elds) at the crack tip, and r type singularity exists for displacements and electric potential [30]. Therefore, methods of calculating stress intensity factors and energy release rates for elastic problems can be directly extended to piezoelectric problems. In the application of conventional Green's functions, the popular technique of quarter-point elements can avoid mesh re®nement at the crack tip [2]. By moving the mid-side node of a quadratic elementpto a quarter point position from the crack tip, the desired r type singularity for displacements and electric potential can be achieved. Stress and electrical displacement (®eld) intensity factors can then be determined from nodal tractions and electrical charge or from the crack opening displacements and electrical potential jump. In the application of crack Green's functions, intensity factors are computed from the solution of internal points near the crack tip [23]. Mechanical (strain) energy and electrical energy coexist in piezoelectric solids, and their sum is the total energy. The direct computation of energy release rates involves two boundary element analyses: one for the original crack
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
775
Fig. 3. Crack front boundary elements.
shape and the other for the crack shape after a small extension [1]. The energy release rates can then be evaluated by a ®nite difference scheme of computed strain or electrical energy. In the case of conventional Green's functions, the crack closure integral provides an alterative way to compute energy release rates on the crack line [7], which can effectively make up for the disadvantage of crack
face discretization. The crack closure calculation of energy release rates in terms of quadratic boundary element solutions is discussed below. It is noted that the corresponding formulae are not reported in the literature for the elastic case. Suppose a crack extends by a small amount da, the total energy release rate can be expressed in the following form z (polar.)
z (polar.) a 2a
x
x
w b). An edge-cracked half plane
a). A center cracked plane
z (polar.)
z (polar.)
b a
c). A branched crack
b1
ω
x
a
b2
ω1 ω2 x
d). A double-branched (forked) crack Fig. 4. Crack con®gurations.
776
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
quadratic interpolation functions Mj h are
Table 1 Strain energy release rate GM I (N/m) in PZT-4 for various l=a l/a [15] Present study Error (%)
0.10 2.749 2.681 2.47
0.12 2.710 1.42
0.15 2.741 0.29
0.16 2.748 0.04
0.18 2.758 2 0.32
0.20 2.765 2 0.58
by using a polar coordinate system with the origin at the crack tip.
M1 h 2 h 2 0:5 h 2 1 M3 h 2h h 2 0:5
GM I
where i x; z; ui da 2 x; ^p; f da 2 x; ^p and Dz da 2 x; ^p denote displacement jump, electric potential jump and normal electric displacement jump across the crack faces, respectively; the ®rst term of the integral contributes to the mechanical energy release rate (G M), the second and third terms to the electric energy release rate (G E). Similar to the elastic case, i z leads to Mode I cracks and i x for Mode II. It is noted that the third term in Eq. (13) vanishes for an impermeable crack, and the second term vanishes for a conducting crack. Fig. 3 shows crack front elements connected to the crack fronts. The initial crack front is at point B (Fig. 3a). After extension (Fig. 3b), the crack front is at point C. The combination of these two boundary element solutions yields the energy release rate by using Eq. (13). When the crack extension d l is very small compared to the crack length, boundary element solutions on element BC after crack extension are approximately equal to those on AB before extension. Therefore, only one boundary element analysis of the case shown in Fig. 3a is required. In case of E an impermeable crack, GM I and GI can be expressed in terms of boundary element nodal data in the following manner: GM I
2 3 3 3 X 1 Z1 4X Mj h uzj1 2 uzj2 5 Mj hs z z j dh 2 0 j1 j1
16 1 4Dz 1Df 1 1 2Dz 2Df 1 1 16Dz 2Df 2 60 1 2Dz 3Df 2 2 Dz 3Df 1 1 2Dz 1Df 2 where Du u1 2 u2 ; Df f1 2 f2 and the numbers in the parentheses are the appropriate node numbers. GM II can be obtained in the same manner. Similar expressions can be obtained for conducting cracks. 6. Numerical results and discussion Numerical examples are presented in this section for crack con®gurations shown in Fig. 4 to show the effectiveness of the boundary element method for piezoelectric crack problems. The accuracy of boundary element results is veri®ed by comparison with available analytical solutions in the literature. Plane strain conditions are assumed. PZT-4 is used in the computation, and its properties based on Eq. (A1) are given below [15] c11 13:9 £ 1010 N=m2 ;
c12 7:78 £ 1010 N=m2 ;
c13 7:43 £ 1010 N=m2 c33 11:3 £ 1010 N=m2 ; e31 26:98 C=m2 ;
c44 2:56 £ 1010 N=m2 e33 13:84 C=m2 ;
e15 13:44 C=m2 (14)
GEI
1 4s zz 1Du 1 1 2s zz 2Du 1 1 16s zz 2Du 2 60 1 2s zz 3Du 2 2 s zz 3Du 1 1 2s zz 1Du 2
GEI 13
15
The substitution of Eq. (15) into Eq. (14) and subsequent manipulations result in
1 Zda G lim {s iz x; 0ui da 2 x; ^p 1 Dz x; 0f da da!0 2da 0 2 x; ^p 1 f x; 0Dz da 2 x; ^p} dx
M2 h 24h h 2 1
2 3 3 3 X 1 Z1 4X Mj h fj1 2 fj2 5 Mj hDz j dh 2 0 j1 j1
where the stresses and electric displacements are taken from element BC and the displacements are taken from element AB. The subscripts j 1 and j 2 denote the two nodes belonging to the upper and lower surfaces of the crack which have identical coordinates before deformation. The
e11 6:00 £ 1029 CV=m;
e33 5:47 £ 1029 CV=m
where subscript 1, 2, 3 correspond to x 0 , y 0 , z 0 in Eq. (A1), respectively. A center crack in an in®nite PZT-4 medium is considered ®rst, as shown in Fig. 4a W ! 1: Assuming the crack surfaces are impermeable, Park and Sun [15] presented an analytical solution for strain energy release rate GM I and total energy release rate GI for this problem. [13] presented the analytical solution for GI of the same problem when the medium is PZT-5. Conventional Green's functions are used to solve this problem, and energy release rates are obtained from crack closure calculation. Due to symmetry,
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
half of the medium is discretized. Five elements are used to model the crack surface, and altogether 19 elements are employed. Let s zz1 10 MPa; a 1 mm: Table 1 shows the comparison of GM I corresponding to different ratios of the length of the crack front element to the half crack length l=a: The present results approach the results given by Park and Sun [15]. Note that l=a 0:16 is optimum for twodimensional crack problems in piezoelectricity. To ®nd a suitable l=a ratio for general 2-D crack problems, the twodimensional elastic case is also investigated. The results show that the optimal l=a ratio is 0.16 for two-dimensional boundary element crack closure calculations. Note Farris and Liu [7] found that 0.10 is optimal for three-dimensional elastic crack problems. Now consider combined mechanical and electric loading at far-®eld. A comparison of the total energy release rate and the strain energy release rate with the analytical results [15] is shown in Fig. 5a and b, respectively. Fig. 5a presents the
777
variation of total energy release rate with mechanical loading for three different electric loading. Very good agreement is observed. A negative electric loading tends to decrease the total energy release rate whereas a positive electric loading can either increase or decrease energy release rate depending on the magnitude of mechanical loading. Fig. 5b shows the variation of strain energy release rate with mechanical loading for three different magnitudes of applied electric ®eld. A signi®cant dependence of the strain energy release on electric ®eld is noted. A positive electric ®eld tends to increase the strain energy release rate while a negative ®eld tends to decrease it. Note the realistic energy release rate in Fig. 5 should be set to zero when it becomes negative, due to the closure of crack surfaces. Additional numerical results show that present results agree with the analytical solution by Pak [13] for PZT-5. It can be concluded that the boundary element crack closure scheme is numerically ef®cient and accurate. Note that Sack [22]
Fig. 5. Energy release rates for a centre-cracked PZT-4 plane under combined mechanical and electrical loading.
778
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Fig. 6. Strain energy release rates for an edge-cracked PZT-4 half plane under combined mechanical and electrical loading..
used 5000 elements for ®nite element crack closure calculation of the above problem. The case of a crack perpendicular to the edge of a semiin®nite PZT-4 medium is shown in Fig. 4b. The study of this problem has attained importance since cracks extending inward from a free edge are frequently found in practice. Impermeable crack faces are assumed. Again ®ve elements are employed to model the crack surface a 1 mm: Mode I s zz1 and Ez1 and Mode II s xz1 and D1 z loading are considered in Fig. 6a and b, respectively. A signi®cant change in the strain energy release rate is observed under different electric loading. For a given mechanical load, the strain energy release rate is higher for a positive electric displacement (®eld) and lower for a negative electric displacement (®eld). Consider a ®nite value of W in Fig. 4a. The boundary effect of a center-cracked specimen on fracture parameters is portrayed in Figs. 7 and 8. Assume uniform tensile stress loading and impermeable crack faces. The variation of
strain energy release rate with relative crack length 2a=W is shown in Fig. 7. The energy release rates are calculated by crack closure scheme. The variation of stress intensity factor with 2a=W is shown in Fig. 8. The impermeable crack Green's functions are used and eight elements are employed to model the outside boundary. The stress intensity factor is computed from internal point solutions. The accuracy of the numerical scheme is justi®ed by a comparison with the analytical solution for an in®nite medium (error , 4%). Strong boundary effects are noted from the solutions presented in Figs. 7 and 8. These boundary element results can be useful when dealing with ®nitewidth specimens of piezoceramics (e.g. practical situations and experimental studies). The case of conducting crack faces is also examined by using conducting crack Green's functions. The corresponding stress intensity factors are found to be virtually identical to those shown in Fig. 8. An impermeable branched crack in an in®nite medium is
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781 25
779
∞
σzz=8 MPa
σ∞ =10 MPa 20
zz
σ∞ =12 MPa
15
I
GM (N/m)
zz
10
5
0 0.1
0.2
0.3 2a/W
0.4
0.5
Fig. 7. Variation of strain energy release rate with cracked length for a centre-cracked PZT-4 specimen under uniform tension.
shown in Fig. 4c. Branched cracks are commonly encountered in the fracture of brittle materials. The technique of multi-domain boundary element modeling is used. The entire medium is divided into two domains so that each domain contains one crack face. Forty-eight elements are used to model each domain. Quarter-point elements are employed at the branch tip. Consider the case where the ratio of branch-main crack length, b=a 0:25; and the remote electrical displacement loading. The corresponding ®eld intensity factors are shown in Table 2, where KD denotes the electric displacement intensity factor for the main crack (without branch). The corresponding analytical results given by [29] are also shown. Fairly good agreement is observed. The case of a double-branched (forked) crack is considered in Fig. 4d. [11] observed that an impermeable crack bifurcates and has a feathered appearance in a piezoceramic sample. A theoretical study of forked cracks in piezoelectrics has not appeared in the literature. The entire medium is divided into three domains that join along crack lines. The
®eld intensity factors at the branch tip are shown in Table 3 for the case of symmetrical branches and remote tension loading. KI denotes the stress intensity factor for the main crack. It is reasonable to observe that KIb decreases when the branch angle increases. The stress intensity factors for asymmetrical branches are shown in Table 4. It is noted that an increase in branch length results in a stronger KIb at the corresponding branch tip. 7. Conclusions A solution scheme based on the boundary integral equation method is developed to analyze plane cracks in piezoelectric solids. A complete set of piezoelectric Green's functions are presented in closed form, based on the extended Lekhnitskii's formalism. Special Green's functions are obtained for a medium containing an impermeable or conducting crack. The accurate calculation of ®eld intensity factors and energy release rates is discussed and
2
KI (MPam
1/2
)
1.5
1
σ∞ =8 MPa
0.5
zz
σ∞ =10 MPa zz
σ∞ =12 MPa zz
0 0.1
0.2
0.3 2a/W
0.4
0.5
Fig. 8. Variation of stress intensity factor with crack length for a centre-cracked PZT-4 specimen under unifor tension.
780
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Table 2 Field intensity factors at branch tip in PZT-4 b 0; b=a 0:25 Loading
D1 z 1
v (deg)
30 40 50 60
Present study
Table 4 Field intensity factors at tips of an asymmetrically branched crack in PZT-4 (b 0; v1 30; v2 45 in Fig. 4d)
Ref. [29]
KIb =KD 107
KDb =KD
7 b I =KD 10
KDb =KD
0.09 0.21 0.40 0.71
1.00 0.94 0.88 0.82
0.10 0.22 0.44 0.75
1.05 1.00 0.94 0.87
illustrated through numerical examples. Comparison with analytical solutions reported in the literature con®rms the accuracy and effectiveness of the boundary element schemes. The present results con®rm that the boundary element methods are computationally ef®cient alternatives for the ®nite element method when applied to piezoelectric crack problems. A logical extension to the present study is the derivation of multi-crack Green's functions and its boundary element implementation, which can be used to ef®ciently analyze micro-crack clusters.
Branch length
KIb1 =KI
KIb2 =KI
s zz1 1
b1 =a 0:25; b2 =a 0:25 b1 =a 0:25; b2 =a 0:30 b1 =a 0:25; b2 =a 0:40
0.78 0.73 0.67
0.50 0.58 0.65
mn n 1; 2; 3 in Eq. (5) are chosen such that their imaginary parts are positive. The constants pn, qn and sn in Eq. (5) are pn a11 m2n 1 a12 2 a13 mn 1 dn b11 mn 2 b21 ; sn b11 m2n 1 b12 2 b13 mn 2 dn d11 mn 2 d12
The work presented in this paper was supported by the Natural Sciences and Engineering Research Council of Canada grant A-6507.
dn l2 m=l1 m p The constants kni n; i 1; 2; 3 appearing in Eq. (6) are
gp1n mn 1 m i finp ;
A4 gp2n 1 1 finp ;
gp3n dn 1 d i finp
p f1n qn p2 s3 2 p 3 s2 2 pn q2 s3 2 q 3 s2 2 sn p 2 q3
2 p3 q 2 =Dp1
Appendix A The constitutive equations for z 0 -polarized ceramics in the x 0 y 0 z 0 system can be expressed as [16] 0
T
0
s ce 2 e E ;
0
0
0
D ee 1 eE A1
where a prime denotes variables with respect to x 0 y 0 z 0 system; cij, eij and e ij denote elastic, piezoelectric and dielectric constants, respectively. mn are the roots of the characteristic equation P m l1 ml3 m 1 l22 m 0
A2
where l1(m) d11m2 2 2d12m 1 d22, l2(m) b11m3 2 (b21 1 b13)m 2 1 (b12 1 b23)m 2 b22, l3(m) a11m4 2 2a13m 3 1 (2a121 a33)m 2 2 2a23m 1 a22. The roots (m) of Eq. (A2) are complex with three pair wise conjugates, and they are generally distinct. Note Table 3 Field intensity factors at tips of a symmetrically branched crack in PZT-4 (b 0; v1 v2 v; b1 =a b2 =a 0:25 in Fig. 4d) Loading
v (deg)
KIb =KI
b II =KI
s zz1 1
30 45 60
0.68 0.59 0.31
0.23 0.43 0.55
KDb =KI 2 0.25 £ 10 210 2 0.23 £ 10 210 2 0.22 £ 10 210
(A3)
qn a12 m2n 1 a22 2 a23 mn 1 dn b12 mn 2 dn b22 =mn ;
kp gp 21
Acknowledgements
0
Loading
(A5)
p f2n pn q1 s3 2 q 3 s1 2 qn p1 s3 2 p 3 s1 1 sn p 1 q3
2 p3 q 1 =Dp1 p f3n qn p1 s2 2 p 2 s1 2 pn q1 s2 2 q 2 s1 2 sn p 1 q2
2 p2 q 1 =Dp1
Dp1 p 1 q 2 s3 2 q3 s2 2 p 2 q 1 s3 2 q3 s1 1 p3 q 1 s2 2 q2 s1 and an overbar denotes the complex conjugate of a complex-valued quantity. The constants kni n; i 1; 2; 3 appearing in Eq. (8) are k g21 g1n pn 1 pi fin ;
A6 g2n qn 1 q i fin ; g3n sn 1 si fin
f1n m 3 d2 2 m 2 d3 1 mn d 3 2 mn d 2 2 dn m 3 1 dn m 2 =D1 f2n m 1 d 3 2 m3 d 1 2 mn d 3 1 mn d 1 2 dn m 1 1 dn m 3 =D1 f3n m 2 d 1 2 m 1 d 2 2 mn d 1 1 mn d 2 2 dn m 2 1 dn m 1 =D1
D1 m 1 d 2 2 d 3 1 m 2 d 3 2 d 1 1 m 3 d 1 2 d 2 vni n; i 1; 2; 3
(A7)
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
appearing in Eq. (8) are vni ui1 2 mn ui2 2 dn ui3
A8
where u s21 ; and {s1j ; s2j ; s3j }T Im
3 X n1
{1; 2mn ; 2dn }T knj :
The constants k^ni n; i 1; 2; 3 appearing in Eq. (10) are ^ g ^ 21 k
A9 g2n qn 1 qi f^in ; g3n dn 1 d i f^in
g^ 1n pn 1 p n f^in ;
f^1n m 3 s2 2 m 2 s3 1 mn s3 2 mn s2 2 sn m 3 1 sn m 2 =D^ 1 f^2n m 1 s3 2 m 3 s1 2 mn s3 1 mn s1 2 sn m 1 1 sn m 3 =D^ 1 f^3n m 2 s1 2 m 1 s2 2 mn s1 1 mn s2 2 sn m 2 1 sn m 1 =D^ 1
D^ 1 m 1 s2 2 s3 1 m 2 s3 2 s1 1 m 3 s1 2 s2 A10 v^ni n; i 1; 2; 3 appearing in Eq: 10 are v^ni u^i1 2 mn u^ i2 2 sn u^ i3
A11
where ^ ^s21 ; u and {^s1j ; s^2j ; s^3j }T Im
3 X n1
{1; 2mn ; 2sn }T k^nj :
References [1] Aliabadi MH, Rooke DP. Numerical fracture mechanics. Southampton: Kluwer Academic Publishers, 1991. [2] Banerjee PK. The boundary element methods in engineering. New York: McGraw-Hill, 1992. [3] Cruse TA. Boundary element analysis in computational fracture mechanics. Amsterdam: Kluwer Academic Publishers, 1988. [4] Cruse TA. BIE fracture mechanics analysis: 25 years of development. Comput Mech 1996;18:1±11.
781
[5] Deeg WF. The analysis of dislocation, crack and inclusion problems in piezoelectric solids. PhD, Thesis, Stanford University, USA 1980. [6] Dunn ML, Wienecke HA. Green's functions for transversely isotropic piezoelectric solids. Int J Solids Struct 1996;30:4571±81. [7] Farris TN, Liu M. Boundary element crack closure calculation of three-dimensional stress intensity factors. Int J Fracture 1993;60:33±47. [8] Hill LR, Farris TN. Three-dimensional piezoelectric boundary element method. Proc SPIE Conf Smart Struct Mater 1997;3039:406±17. [9] Lee KJ, Jiang LZ. A boundary integral formulation and 2D fundamental solutions for piezoelectric media. Mech Res Comm 1994;21:47±54. [10] Lu P, Mahrenholtz O. A variational boundary element formulation for piezoelectricity. Mech Res Comm 1994;21:605±11. [11] Lynch CS, Yang W, Collier L, Suo Z, McMeeking RM. Electric ®eld induced cracking in ferroelectric ceramics. Ferroelectrics 1995;166:11±30. [12] Mikhailov GK, Parton VZ. Electromagnetoelasticity. New York: Hemisphere, 1990. [13] Pak YE. Linear electro-elastic fracture mechanics of piezoelectric materials. Int J Fracture 1992;54:79±100. [14] Pan E. A BEM analysis of fracture mechanics in 2D anisotropic piezoelectric solids. Engng Analysis Boundary Elem 1999;23:67±76. [15] Park SB, Sun CT. Fracture criteria for piezoelectric ceramics. J Am Ceramic Soc 1995;78:1475±80 (and correction (1996)). [16] Parton VZ, Kudryavtsev BA. Electromagnetoelasticity. New York: Gordon and Breach Science Publishers, 1988. [17] Qin QH. Thermoelectroelastic analysis of cracks in piezoelectric halfplane by BEM. Comput Mech 1999;23:353±60. [18] Qin QH, Lu M. BEM for crack-inclusion problems of plane thermopiezolectric solids. Int J Numer Meth Engng 2000;48:1071±88. [19] Qin QH, Mai YW. Crack branch in piezoelectric bi-material system. Int J Engng Sci 2000;38:667±87. [20] Rajapakse RKND. Plane strain/stress solutions for piezoelectric solids. Composites Engng: part B 1997;28:385±96. [21] Rajapakse RKND. Boundary element methods for piezoelectric materials. Proc SPIE Conf Smart Struct Mater 1997;3039:418±28. [22] Sack EE. Evaluation of ®nite analysis of piezoelectric materials with holes and cracks. M.Sc Thesis, Purdue University, USA 1994. [23] Snyder MD, Cruse TA. Boundary-integral equation analysis of cracked anisotropic plates. Int J Fracture 1975;11:315±28. [24] Sosa HA. On the fracture mechanics of piezoelectric solids. Int J Solids Struct 1992;29:2613±22. [25] Suo Z. Models for breakdown-resistant dielectric and ferroelectric ceramics. J Mech Phys Solids 1993;41:1155±76. [26] Suo Z, Kuo CM, Barnett DM, Willis JR. Fracture mechanics for piezoelectric ceramics. J Mech Phys Solids 1992;40:739±65. [27] Ting TCT. Anisotropy elasticity: theory and applications. Oxford: Oxford University Press, 1996. [28] Xu X-L, Rajapakse RKND. Boundary element analysis of piezoelectric solids with defects. Composites Engng: part B 1998;29:655±69. [29] Xu X-L, Rajapakse RKND. A theoretical study of branched cracks in piezoelectrics. Acta Materialia 2000;48:1865±82. [30] Xu X-L, Rajapakse RKND. On singularities in composite piezoelectric wedges and junctions. Int J Solids Struct 2000;37:3253±75.
www.elsevier.com/locate/enganabound
Boundary element modeling of cracks in piezoelectric solids R.K.N.D. Rajapakse*, X.-L. Xu Department of Civil and Geological Engineering, University of Manitoba, Winnipeg, Canada R3T 5V6 Received 19 June 2000; accepted 1 March 2001
Abstract This study focuses on the application of boundary element methods for linear fracture mechanics of two-dimensional piezoelectric solids. A complete set of piezoelectric Green's functions, based on the extended Lekhnitskii's formalism and distributed dislocation modeling, are presented. Special Green's functions are obtained for an in®nite medium containing a conducting crack or an impermeable crack. The numerical solution of the boundary integral equation and the computation of fracture parameters are discussed. The concept of crack closure integral is utilized to calculate energy release rates. Accuracy of the boundary element solutions is con®rmed by comparing with analytical solutions reported in the literature. The present scheme can be applied to study complex cracks such as branched cracks, forked cracks and microcrack clusters. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Piezoelectric; Impermeable crack; Boundary element; Conducting crack; Fracture mechanics
1. Introduction Piezoelectric actuators are widely used in the development of adaptive (smart) structures. The main drawback of piezoelectrics is their brittleness. Defects due to manufacturing process and complex electromechanical loading encountered in modern applications increase the likelihood of fracture and failure of piezoceramic actuators. A thorough understanding of fracture behavior of piezoelectric materials is therefore vital to the advancement of adaptive structures technology. A number of researchers have examined crack problems in piezoelectric materials using analytical methods. Different types of electrical boundary conditions on crack faces have been proposed. Mikhailov and Parton [12] considered the crack faces as electrically permeable, which led to a controversial conclusion that an electric ®eld has no contribution to crack propagation. Suo et al. [26] argued that this condition is not realistic, as there is an electric potential drop across the crack. Using an impermeable crack model, Sosa and Suo et al. [24,26] added an electric intensity factor to the well known stress intensity factors. [13,15] used the energy release rates to study fracture. A conducting crack model has also been considered [25]. Cracks in piezoelectric * Corresponding author. Address: Department of Mechanical Engineering, University of British Columbia, 2324 Main Mall, Vancouver, Canada V6T1Z4. Tel.: 11-604-822-0497; fax: 11-604-822-0944. E-mail address: [email protected] (R.K.N.D. Rajapakse).
actuators may emanate from an electrode edge. A conducting crack model is suitable when a conducting species migrates on the crack surfaces. Qin and Mai and Xu and Rajapakse [19,29], among others, considered crack branching problems. The above studies provided the theoretical basis for fracture mechanics of piezoelectrics. However, the study of fracture problems related to practical situations involving ®nite geometries and complex loading requires the development of numerical techniques. Some recent studies have been concerned with the development of boundary element methods for piezoelectrics. The ef®ciency of boundary element method is mainly attributed to the reduction in the dimensionality of the problem. Moreover, compared to the ®nite element method, the boundary element method is highly accurate and computationally ef®cient when dealing with linear fracture mechanics [2,3]. Indeed, fracture mechanics problems are one of the principal areas of application of boundary element method [4]. Green's functions (fundamental solutions) are essential for all boundary element methods. Superior numerical ef®ciency can be achieved when closed form Green's functions are available. Deeg [5] presented a comprehensive theoretical study of Green's functions, dislocations and cracks in piezoelectric solids under general anisotropy. However, his Green's functions are not in closed form due to the consideration of general anisotropy. Lee and Jiang and Rajapakse [9,20] derived two-dimensional piezoelectric Green's functions for materials with hexagonal symmetry by using double Fourier transformation and
0955-7997/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0955-799 7(01)00058-3
772
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Fourier transformation, respectively. Pan [14] used Deeg's formulation to rederive two-dimensional Green's functions for in®nite planes, half planes and bi-material planes. Dunn and Wienecke [6] derived three-dimensional Green's functions. Parton and Kudryavtsev [16] presented a boundary integral equation for linear piezoelectric problems. Lee and Jiang, Lu and Mahrenholtz, Qin and Qin and Lu [9,10,17,18], among others, also considered boundary element formulations for piezoelectric solids. Rajapakse [21] discussed three boundary integral formulations (the direct and indirect methods, and the ®ctitious stress-charge method) for linear piezoelectricity. However, only a few studies presented numerical implementations. Hill and Farris [8] conducted a three-dimensional boundary element analysis of piezoelectric solids by a numerically evaluating Green's functions presented by Deeg [5]. Xu and Rajapakse [28] developed a two-dimensional boundary element code based on the closed form Green's functions [20] to examine stress and electric ®eld concentration around voids. Pan [14] proposed a single-domain boundary element method for two-dimensional crack problems. Numerical results for ®eld intensity factors, crack opening displacement and electric potential jump across crack faces were presented for an inclined crack. The objective of this study is to extend the application of the conventional boundary element method to fracture mechanics of piezoelectric solids. A complete set of piezoelectric Green's functions are derived. Conventional Green's functions, which are electroelastic solutions due to concentrated line forces and electric charges in an in®nite homogeneous plane, are rederived by using the extended Lekhnitskii's formalism. Crack Green's functions are then derived by considering an in®nite piezoelectric plane containing a center crack. The crack faces can be either electrically impermeable or conducting. The application of these Green's functions in boundary element analysis and the computation of fracture parameters are then discussed. The accuracy of the boundary element solutions is con®rmed and selected numerical solutions for crack problems are presented.
2. Basic equations A plane piezoelectric medium is considered in Fig. 1. The medium is polarized along the z 0 -direction, which makes angle b with the z-direction. In the absence of body forces and body electric charges, the governing equations for the medium can be expressed as
s ij;j 0;
Di;i 0
1
where s ij i; j x; z and Di denote the stress tensor and electric displacement in the i-direction, respectively. The constitutive relations in the xz system can be
z
z’ (pol. direc.) β
Γ x
y (y’) Fig. 1. A plane piezoelectric medium.
expressed as 8 9 0 exx > a11 > > > < = B ezz B @ a12 > > > > : ; 2exz a13 (
Ex
)
Ez
2
a12 a22 a23
b11
b12
b21
b22
9 0 18 a13 > s xx > b11 > > = B C< B a23 C A> s zz > 1 @ b12 > > : ; a33 s xz b13 8 9 s xx > !> > > < = b13 d11 s zz 1 > > > b23 > d12 : ; s xz
1 b21 ( ) C Dx b22 C A Dz b23
d12 d22
!(
Dx
)
Dz
2 where e ij and Ei denote the strain tensor, and electric ®eld in the i-direction, respectively; aij, bij and dij denote twodimensional elastic, piezoelectric and dielectric constants, respectively. These material constants are functions of polarization angle b and material constants in Eq. (A1), and are different for plane stress and plane strain cases. The following boundary conditions are admissible on the boundary G of the medium u i u~ i
or
s ij nj t~ i
3
f f~
or
Di ni 2q~
4
where ui and f denote the displacement in the i-direction and electric potential, respectively; u~ i ; t~ i ; ni, f~ and q~ denote speci®ed displacement and traction components in the idirection, components of outward unit normal vector, speci®ed electric potential and speci®ed surface charge respectively. By utilizing the extended Lekhnitskii's formalism, the general solutions for plane piezo electrics can be expressed as [29] {ux ; uz ; f}T 2Re
3 X n1
{pn ; qn ; sn }T wn zn
5
where Re denotes real part of a complex valued quantity; pn, qn and sn are complex constants de®ned in Appendix A that depend on material properties and polarization angle b ; a
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
superscript T denotes transpose of a vector; zn x 1 mn z where mn n 1; 2; 3 are the roots of Eq. (A2); wn zn are complex potential functions to be determined from the loading and boundary conditions.
Fig. 2b involves a homogeneous medium subjected to a line force/charge, and has been solved previously. In Fig. 2c, a continuous generalized dislocation ®eld (distributed discontinuities in displacements and electric potential) with densities bj ( j 1,2,3) is utilized to model the impermeable crack. The resultant electroelastic ®elds on the crack line due to line force/charge and distributed dislocation ®eld should satisfy the boundary conditions given by Eq. (7). The solution for a single generalized edge dislocation is sought ®rst. The dislocation ®eld bi is then obtained by enforcing boundary conditions of Eq. (7) and making use of Eq. (5). For the sake of brevity, the details of derivation are omitted. The ®nal results for the potential functions wn zn corresponding to the subproblem in Fig. 2c are obtained as
3. Green's functions 3.1. Green's functions for line force/charge Several previous studies have considered the derivation of conventional Green's functions for a piezoelectric in®nite plane. The present study utilizes the extended Lekhnitskii's formalism. Omitting the details, the ®nal results are given below. Consider a concentrated line force/charge P(P1, P2, P3), i.e. horizontal force P1, vertical force P2 and electric charge P3, applied at point z 0(x0, z0). The complete Green's functions are obtained by substituting the following potentials into general solutions given by Eq. (5)
wn zn
1 p ln zn 2 z0n kni Pi 2pi
w 0n zn
1 1 kp P 2pi zn 2 z0n ni i
kni Im 2pi " # pp 3 X z0n 2 2 a2 z2n 2 a2 1 z0n zn 2 a2 p £ vni ln knj Pj zn 2 z0n n1
wn zn
6
kni Im 2p i " # p ! 3 X z0n 2 2 a2 p 1 £ 1 2 p knj Pj vni zn 2 z0n z2n 2 a2 n1
w 0n zn
p where complex constants kni n; i 1; 2; 3 de®ned in Appendix A are functions of material properties and polarization angle b .
8
3.2. Crack Green's functions
where coef®cients kni and vni n; i 1; 2; 3 given in the Appendix A are functions of material properties and polarization angle b . Superposition of Eqs. (6) and (8) yields the complex potential functions corresponding to impermeable crack Green's functions. Now assume that the crack faces in Fig. 2a are conducting. Therefore
Availability of crack Green's functions allows the development of computationally ef®cient boundary element methods for crack problems. For the case of purely elastic solids, Ting [27] derived crack Green's functions based on Stroh's formalism. In this subsection, Green's functions for impermeable cracks and conducting cracks are derived based on the extended Lekhnitskii's formalism. The case of an impermeable crack is considered ®rst (Fig. 2a). The following boundary conditions hold on the crack faces
s zz 0;
s xz 0;
Dz 0; 2a # x # a
s zz 0;
7
z
Ex 0; 2a # x # a
z
z P
(a)
s xz 0;
P x
=
x (b)
9
The decomposition scheme and dislocation modeling are utilized again. Similar to the case of an impermeable crack, the generalized dislocation ®eld bi which satis®es the conducting boundary conditions can be obtained. The ®nal results for the potential functions wn zn corresponding to
The problem in Fig. 2a can be decomposed into two subproblems shown in Fig. 2b and c. The subproblem in
2a
773
+
bj (c)
Fig. 2. Decomposition used to derive crack Green's.
x
774
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Fig. 2c are k^ wn zn ni Im "2pi pp # 3 X z0n 2 2 a2 z2n 2 a2 1 z0n zn 2 a2 p £ v^ni ln knj Pj zn 2 z0n n1 k^ni Im "2pi # p ! 3 X z0n 2 2 a2 p 1 £ 1 2 p knj Pj v^ni zn 2 z0n z2n 2 a2 n1
w 0n zn
10 where coef®cients k^ni and v^ni n; i 1; 2; 3 given in the Appendix A are functions of the material properties and polarization angle b . Superposition of Eqs. (6) and (10) yields the complex potential functions corresponding to conducting crack Green's functions. 4. Boundary element method Based on the reciprocal relations for a piezoelectric solid, the following integral equation can be derived in the absence of body forces and charges [16] Z c x 0 ui x 0 Gji x; x 0 tj x 2 Hji x; x 0 uj x dG 11
where N is the total number of boundary elements, n is the total number of nodes on the l-th boundary element; tkj ; ukj ; qk ; fk are j-th component of nodal traction and displacement, electric charge and electric potential respectively of the k-th local node of the l-th boundary element; the p Jacobian J h is equal to 2x=2h2 1 2z=2h2 : Quadratic shape functions M k h are used in the present study. When the ®eld point x coincides with the node k, attention should be paid to evaluate singular integrals arising from Gij and Hij in Eq. (12). For conventional Green's functions in homogeneous piezoelectric media, based on the solutions in Section 3.2, Gij have ln r type singularity and Hij have 1=r type singularity. The inspection of crack Green's functions indicates that no additional singularity exists. Therefore, piezoelectric Green's functions have the same type of singularity as the elastic counterpart, and the integration scheme for elastic problems can be directly applied to calculate the singular integrals in Eq. (12). The numerical implementation and computation of boundary and internal ®eld variables are similar to the case of elastic materials. Details of numerical scheme for elastic problems are given by Banerjee [2]. While the present study focuses on the application of traditional boundary element method to piezoelectric solids, it is worth noting that Pan [14] proposed a single-domain boundary element method for piezoelectric crack problems.
G
where i x; z; q; u q 2f; tq 2q; and summation is implied on the index j j x; z; q: Gji and Hji are the Green's functions representing the electroelastic ®elds due to concentrated line force/charge in an in®nite medium: Gji x; x 0 and Hji x; x 0 denote the displacement and traction component at ®eld point x [ G in the j-direction j x; z due to a unit concentrated line load applied in the i-direction i x; z at point x 0 respectively; G qi x; x 0 and Hqi x; x 0 denote the electric potential and electric charge at ®eld point x due to unit concentrated line load applied in the i-direction at point x 0 respectively; Gjq x; x 0 and Hjq x; x 0 denote the displacement and traction components in the j-direction at ®eld point x due to unit electric line charge at point x 0 respectively; Gqq x; x 0 and Hqq x; x 0 denote the electric potential and electric charge at ®eld point x due to the unit electric charge at load point x 0 respectively. The coef®cient c x 0 can be calculated by the concept of rigid body movement. Upon discretizing the boundary G , Eq. (11) becomes " # ( N n Z X X 0 0 0 k Gji x; x M hJ h c x u i x DG l
l1
)
(
£ dG h t jk 2 " £
0
Hji x; x
n X k1
k1
Z
12
DG l
# k
M hJ h dG h
!
) u jk
5. Computation of fracture parameters Field intensity factors and energy release rates are commonly used fracture parameters. It has been shown that, for both impermeable and conducting cracks, p classical inverse square root type singularity 1= r exists for stresses andpelectric displacements (®elds) at the crack tip, and r type singularity exists for displacements and electric potential [30]. Therefore, methods of calculating stress intensity factors and energy release rates for elastic problems can be directly extended to piezoelectric problems. In the application of conventional Green's functions, the popular technique of quarter-point elements can avoid mesh re®nement at the crack tip [2]. By moving the mid-side node of a quadratic elementpto a quarter point position from the crack tip, the desired r type singularity for displacements and electric potential can be achieved. Stress and electrical displacement (®eld) intensity factors can then be determined from nodal tractions and electrical charge or from the crack opening displacements and electrical potential jump. In the application of crack Green's functions, intensity factors are computed from the solution of internal points near the crack tip [23]. Mechanical (strain) energy and electrical energy coexist in piezoelectric solids, and their sum is the total energy. The direct computation of energy release rates involves two boundary element analyses: one for the original crack
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
775
Fig. 3. Crack front boundary elements.
shape and the other for the crack shape after a small extension [1]. The energy release rates can then be evaluated by a ®nite difference scheme of computed strain or electrical energy. In the case of conventional Green's functions, the crack closure integral provides an alterative way to compute energy release rates on the crack line [7], which can effectively make up for the disadvantage of crack
face discretization. The crack closure calculation of energy release rates in terms of quadratic boundary element solutions is discussed below. It is noted that the corresponding formulae are not reported in the literature for the elastic case. Suppose a crack extends by a small amount da, the total energy release rate can be expressed in the following form z (polar.)
z (polar.) a 2a
x
x
w b). An edge-cracked half plane
a). A center cracked plane
z (polar.)
z (polar.)
b a
c). A branched crack
b1
ω
x
a
b2
ω1 ω2 x
d). A double-branched (forked) crack Fig. 4. Crack con®gurations.
776
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
quadratic interpolation functions Mj h are
Table 1 Strain energy release rate GM I (N/m) in PZT-4 for various l=a l/a [15] Present study Error (%)
0.10 2.749 2.681 2.47
0.12 2.710 1.42
0.15 2.741 0.29
0.16 2.748 0.04
0.18 2.758 2 0.32
0.20 2.765 2 0.58
by using a polar coordinate system with the origin at the crack tip.
M1 h 2 h 2 0:5 h 2 1 M3 h 2h h 2 0:5
GM I
where i x; z; ui da 2 x; ^p; f da 2 x; ^p and Dz da 2 x; ^p denote displacement jump, electric potential jump and normal electric displacement jump across the crack faces, respectively; the ®rst term of the integral contributes to the mechanical energy release rate (G M), the second and third terms to the electric energy release rate (G E). Similar to the elastic case, i z leads to Mode I cracks and i x for Mode II. It is noted that the third term in Eq. (13) vanishes for an impermeable crack, and the second term vanishes for a conducting crack. Fig. 3 shows crack front elements connected to the crack fronts. The initial crack front is at point B (Fig. 3a). After extension (Fig. 3b), the crack front is at point C. The combination of these two boundary element solutions yields the energy release rate by using Eq. (13). When the crack extension d l is very small compared to the crack length, boundary element solutions on element BC after crack extension are approximately equal to those on AB before extension. Therefore, only one boundary element analysis of the case shown in Fig. 3a is required. In case of E an impermeable crack, GM I and GI can be expressed in terms of boundary element nodal data in the following manner: GM I
2 3 3 3 X 1 Z1 4X Mj h uzj1 2 uzj2 5 Mj hs z z j dh 2 0 j1 j1
16 1 4Dz 1Df 1 1 2Dz 2Df 1 1 16Dz 2Df 2 60 1 2Dz 3Df 2 2 Dz 3Df 1 1 2Dz 1Df 2 where Du u1 2 u2 ; Df f1 2 f2 and the numbers in the parentheses are the appropriate node numbers. GM II can be obtained in the same manner. Similar expressions can be obtained for conducting cracks. 6. Numerical results and discussion Numerical examples are presented in this section for crack con®gurations shown in Fig. 4 to show the effectiveness of the boundary element method for piezoelectric crack problems. The accuracy of boundary element results is veri®ed by comparison with available analytical solutions in the literature. Plane strain conditions are assumed. PZT-4 is used in the computation, and its properties based on Eq. (A1) are given below [15] c11 13:9 £ 1010 N=m2 ;
c12 7:78 £ 1010 N=m2 ;
c13 7:43 £ 1010 N=m2 c33 11:3 £ 1010 N=m2 ; e31 26:98 C=m2 ;
c44 2:56 £ 1010 N=m2 e33 13:84 C=m2 ;
e15 13:44 C=m2 (14)
GEI
1 4s zz 1Du 1 1 2s zz 2Du 1 1 16s zz 2Du 2 60 1 2s zz 3Du 2 2 s zz 3Du 1 1 2s zz 1Du 2
GEI 13
15
The substitution of Eq. (15) into Eq. (14) and subsequent manipulations result in
1 Zda G lim {s iz x; 0ui da 2 x; ^p 1 Dz x; 0f da da!0 2da 0 2 x; ^p 1 f x; 0Dz da 2 x; ^p} dx
M2 h 24h h 2 1
2 3 3 3 X 1 Z1 4X Mj h fj1 2 fj2 5 Mj hDz j dh 2 0 j1 j1
where the stresses and electric displacements are taken from element BC and the displacements are taken from element AB. The subscripts j 1 and j 2 denote the two nodes belonging to the upper and lower surfaces of the crack which have identical coordinates before deformation. The
e11 6:00 £ 1029 CV=m;
e33 5:47 £ 1029 CV=m
where subscript 1, 2, 3 correspond to x 0 , y 0 , z 0 in Eq. (A1), respectively. A center crack in an in®nite PZT-4 medium is considered ®rst, as shown in Fig. 4a W ! 1: Assuming the crack surfaces are impermeable, Park and Sun [15] presented an analytical solution for strain energy release rate GM I and total energy release rate GI for this problem. [13] presented the analytical solution for GI of the same problem when the medium is PZT-5. Conventional Green's functions are used to solve this problem, and energy release rates are obtained from crack closure calculation. Due to symmetry,
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
half of the medium is discretized. Five elements are used to model the crack surface, and altogether 19 elements are employed. Let s zz1 10 MPa; a 1 mm: Table 1 shows the comparison of GM I corresponding to different ratios of the length of the crack front element to the half crack length l=a: The present results approach the results given by Park and Sun [15]. Note that l=a 0:16 is optimum for twodimensional crack problems in piezoelectricity. To ®nd a suitable l=a ratio for general 2-D crack problems, the twodimensional elastic case is also investigated. The results show that the optimal l=a ratio is 0.16 for two-dimensional boundary element crack closure calculations. Note Farris and Liu [7] found that 0.10 is optimal for three-dimensional elastic crack problems. Now consider combined mechanical and electric loading at far-®eld. A comparison of the total energy release rate and the strain energy release rate with the analytical results [15] is shown in Fig. 5a and b, respectively. Fig. 5a presents the
777
variation of total energy release rate with mechanical loading for three different electric loading. Very good agreement is observed. A negative electric loading tends to decrease the total energy release rate whereas a positive electric loading can either increase or decrease energy release rate depending on the magnitude of mechanical loading. Fig. 5b shows the variation of strain energy release rate with mechanical loading for three different magnitudes of applied electric ®eld. A signi®cant dependence of the strain energy release on electric ®eld is noted. A positive electric ®eld tends to increase the strain energy release rate while a negative ®eld tends to decrease it. Note the realistic energy release rate in Fig. 5 should be set to zero when it becomes negative, due to the closure of crack surfaces. Additional numerical results show that present results agree with the analytical solution by Pak [13] for PZT-5. It can be concluded that the boundary element crack closure scheme is numerically ef®cient and accurate. Note that Sack [22]
Fig. 5. Energy release rates for a centre-cracked PZT-4 plane under combined mechanical and electrical loading.
778
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Fig. 6. Strain energy release rates for an edge-cracked PZT-4 half plane under combined mechanical and electrical loading..
used 5000 elements for ®nite element crack closure calculation of the above problem. The case of a crack perpendicular to the edge of a semiin®nite PZT-4 medium is shown in Fig. 4b. The study of this problem has attained importance since cracks extending inward from a free edge are frequently found in practice. Impermeable crack faces are assumed. Again ®ve elements are employed to model the crack surface a 1 mm: Mode I s zz1 and Ez1 and Mode II s xz1 and D1 z loading are considered in Fig. 6a and b, respectively. A signi®cant change in the strain energy release rate is observed under different electric loading. For a given mechanical load, the strain energy release rate is higher for a positive electric displacement (®eld) and lower for a negative electric displacement (®eld). Consider a ®nite value of W in Fig. 4a. The boundary effect of a center-cracked specimen on fracture parameters is portrayed in Figs. 7 and 8. Assume uniform tensile stress loading and impermeable crack faces. The variation of
strain energy release rate with relative crack length 2a=W is shown in Fig. 7. The energy release rates are calculated by crack closure scheme. The variation of stress intensity factor with 2a=W is shown in Fig. 8. The impermeable crack Green's functions are used and eight elements are employed to model the outside boundary. The stress intensity factor is computed from internal point solutions. The accuracy of the numerical scheme is justi®ed by a comparison with the analytical solution for an in®nite medium (error , 4%). Strong boundary effects are noted from the solutions presented in Figs. 7 and 8. These boundary element results can be useful when dealing with ®nitewidth specimens of piezoceramics (e.g. practical situations and experimental studies). The case of conducting crack faces is also examined by using conducting crack Green's functions. The corresponding stress intensity factors are found to be virtually identical to those shown in Fig. 8. An impermeable branched crack in an in®nite medium is
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781 25
779
∞
σzz=8 MPa
σ∞ =10 MPa 20
zz
σ∞ =12 MPa
15
I
GM (N/m)
zz
10
5
0 0.1
0.2
0.3 2a/W
0.4
0.5
Fig. 7. Variation of strain energy release rate with cracked length for a centre-cracked PZT-4 specimen under uniform tension.
shown in Fig. 4c. Branched cracks are commonly encountered in the fracture of brittle materials. The technique of multi-domain boundary element modeling is used. The entire medium is divided into two domains so that each domain contains one crack face. Forty-eight elements are used to model each domain. Quarter-point elements are employed at the branch tip. Consider the case where the ratio of branch-main crack length, b=a 0:25; and the remote electrical displacement loading. The corresponding ®eld intensity factors are shown in Table 2, where KD denotes the electric displacement intensity factor for the main crack (without branch). The corresponding analytical results given by [29] are also shown. Fairly good agreement is observed. The case of a double-branched (forked) crack is considered in Fig. 4d. [11] observed that an impermeable crack bifurcates and has a feathered appearance in a piezoceramic sample. A theoretical study of forked cracks in piezoelectrics has not appeared in the literature. The entire medium is divided into three domains that join along crack lines. The
®eld intensity factors at the branch tip are shown in Table 3 for the case of symmetrical branches and remote tension loading. KI denotes the stress intensity factor for the main crack. It is reasonable to observe that KIb decreases when the branch angle increases. The stress intensity factors for asymmetrical branches are shown in Table 4. It is noted that an increase in branch length results in a stronger KIb at the corresponding branch tip. 7. Conclusions A solution scheme based on the boundary integral equation method is developed to analyze plane cracks in piezoelectric solids. A complete set of piezoelectric Green's functions are presented in closed form, based on the extended Lekhnitskii's formalism. Special Green's functions are obtained for a medium containing an impermeable or conducting crack. The accurate calculation of ®eld intensity factors and energy release rates is discussed and
2
KI (MPam
1/2
)
1.5
1
σ∞ =8 MPa
0.5
zz
σ∞ =10 MPa zz
σ∞ =12 MPa zz
0 0.1
0.2
0.3 2a/W
0.4
0.5
Fig. 8. Variation of stress intensity factor with crack length for a centre-cracked PZT-4 specimen under unifor tension.
780
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
Table 2 Field intensity factors at branch tip in PZT-4 b 0; b=a 0:25 Loading
D1 z 1
v (deg)
30 40 50 60
Present study
Table 4 Field intensity factors at tips of an asymmetrically branched crack in PZT-4 (b 0; v1 30; v2 45 in Fig. 4d)
Ref. [29]
KIb =KD 107
KDb =KD
7 b I =KD 10
KDb =KD
0.09 0.21 0.40 0.71
1.00 0.94 0.88 0.82
0.10 0.22 0.44 0.75
1.05 1.00 0.94 0.87
illustrated through numerical examples. Comparison with analytical solutions reported in the literature con®rms the accuracy and effectiveness of the boundary element schemes. The present results con®rm that the boundary element methods are computationally ef®cient alternatives for the ®nite element method when applied to piezoelectric crack problems. A logical extension to the present study is the derivation of multi-crack Green's functions and its boundary element implementation, which can be used to ef®ciently analyze micro-crack clusters.
Branch length
KIb1 =KI
KIb2 =KI
s zz1 1
b1 =a 0:25; b2 =a 0:25 b1 =a 0:25; b2 =a 0:30 b1 =a 0:25; b2 =a 0:40
0.78 0.73 0.67
0.50 0.58 0.65
mn n 1; 2; 3 in Eq. (5) are chosen such that their imaginary parts are positive. The constants pn, qn and sn in Eq. (5) are pn a11 m2n 1 a12 2 a13 mn 1 dn b11 mn 2 b21 ; sn b11 m2n 1 b12 2 b13 mn 2 dn d11 mn 2 d12
The work presented in this paper was supported by the Natural Sciences and Engineering Research Council of Canada grant A-6507.
dn l2 m=l1 m p The constants kni n; i 1; 2; 3 appearing in Eq. (6) are
gp1n mn 1 m i finp ;
A4 gp2n 1 1 finp ;
gp3n dn 1 d i finp
p f1n qn p2 s3 2 p 3 s2 2 pn q2 s3 2 q 3 s2 2 sn p 2 q3
2 p3 q 2 =Dp1
Appendix A The constitutive equations for z 0 -polarized ceramics in the x 0 y 0 z 0 system can be expressed as [16] 0
T
0
s ce 2 e E ;
0
0
0
D ee 1 eE A1
where a prime denotes variables with respect to x 0 y 0 z 0 system; cij, eij and e ij denote elastic, piezoelectric and dielectric constants, respectively. mn are the roots of the characteristic equation P m l1 ml3 m 1 l22 m 0
A2
where l1(m) d11m2 2 2d12m 1 d22, l2(m) b11m3 2 (b21 1 b13)m 2 1 (b12 1 b23)m 2 b22, l3(m) a11m4 2 2a13m 3 1 (2a121 a33)m 2 2 2a23m 1 a22. The roots (m) of Eq. (A2) are complex with three pair wise conjugates, and they are generally distinct. Note Table 3 Field intensity factors at tips of a symmetrically branched crack in PZT-4 (b 0; v1 v2 v; b1 =a b2 =a 0:25 in Fig. 4d) Loading
v (deg)
KIb =KI
b II =KI
s zz1 1
30 45 60
0.68 0.59 0.31
0.23 0.43 0.55
KDb =KI 2 0.25 £ 10 210 2 0.23 £ 10 210 2 0.22 £ 10 210
(A3)
qn a12 m2n 1 a22 2 a23 mn 1 dn b12 mn 2 dn b22 =mn ;
kp gp 21
Acknowledgements
0
Loading
(A5)
p f2n pn q1 s3 2 q 3 s1 2 qn p1 s3 2 p 3 s1 1 sn p 1 q3
2 p3 q 1 =Dp1 p f3n qn p1 s2 2 p 2 s1 2 pn q1 s2 2 q 2 s1 2 sn p 1 q2
2 p2 q 1 =Dp1
Dp1 p 1 q 2 s3 2 q3 s2 2 p 2 q 1 s3 2 q3 s1 1 p3 q 1 s2 2 q2 s1 and an overbar denotes the complex conjugate of a complex-valued quantity. The constants kni n; i 1; 2; 3 appearing in Eq. (8) are k g21 g1n pn 1 pi fin ;
A6 g2n qn 1 q i fin ; g3n sn 1 si fin
f1n m 3 d2 2 m 2 d3 1 mn d 3 2 mn d 2 2 dn m 3 1 dn m 2 =D1 f2n m 1 d 3 2 m3 d 1 2 mn d 3 1 mn d 1 2 dn m 1 1 dn m 3 =D1 f3n m 2 d 1 2 m 1 d 2 2 mn d 1 1 mn d 2 2 dn m 2 1 dn m 1 =D1
D1 m 1 d 2 2 d 3 1 m 2 d 3 2 d 1 1 m 3 d 1 2 d 2 vni n; i 1; 2; 3
(A7)
R.K.N.D. Rajapakse, X.-L. Xu / Engineering Analysis with Boundary Elements 25 (2001) 771±781
appearing in Eq. (8) are vni ui1 2 mn ui2 2 dn ui3
A8
where u s21 ; and {s1j ; s2j ; s3j }T Im
3 X n1
{1; 2mn ; 2dn }T knj :
The constants k^ni n; i 1; 2; 3 appearing in Eq. (10) are ^ g ^ 21 k
A9 g2n qn 1 qi f^in ; g3n dn 1 d i f^in
g^ 1n pn 1 p n f^in ;
f^1n m 3 s2 2 m 2 s3 1 mn s3 2 mn s2 2 sn m 3 1 sn m 2 =D^ 1 f^2n m 1 s3 2 m 3 s1 2 mn s3 1 mn s1 2 sn m 1 1 sn m 3 =D^ 1 f^3n m 2 s1 2 m 1 s2 2 mn s1 1 mn s2 2 sn m 2 1 sn m 1 =D^ 1
D^ 1 m 1 s2 2 s3 1 m 2 s3 2 s1 1 m 3 s1 2 s2 A10 v^ni n; i 1; 2; 3 appearing in Eq: 10 are v^ni u^i1 2 mn u^ i2 2 sn u^ i3
A11
where ^ ^s21 ; u and {^s1j ; s^2j ; s^3j }T Im
3 X n1
{1; 2mn ; 2sn }T k^nj :
References [1] Aliabadi MH, Rooke DP. Numerical fracture mechanics. Southampton: Kluwer Academic Publishers, 1991. [2] Banerjee PK. The boundary element methods in engineering. New York: McGraw-Hill, 1992. [3] Cruse TA. Boundary element analysis in computational fracture mechanics. Amsterdam: Kluwer Academic Publishers, 1988. [4] Cruse TA. BIE fracture mechanics analysis: 25 years of development. Comput Mech 1996;18:1±11.
781
[5] Deeg WF. The analysis of dislocation, crack and inclusion problems in piezoelectric solids. PhD, Thesis, Stanford University, USA 1980. [6] Dunn ML, Wienecke HA. Green's functions for transversely isotropic piezoelectric solids. Int J Solids Struct 1996;30:4571±81. [7] Farris TN, Liu M. Boundary element crack closure calculation of three-dimensional stress intensity factors. Int J Fracture 1993;60:33±47. [8] Hill LR, Farris TN. Three-dimensional piezoelectric boundary element method. Proc SPIE Conf Smart Struct Mater 1997;3039:406±17. [9] Lee KJ, Jiang LZ. A boundary integral formulation and 2D fundamental solutions for piezoelectric media. Mech Res Comm 1994;21:47±54. [10] Lu P, Mahrenholtz O. A variational boundary element formulation for piezoelectricity. Mech Res Comm 1994;21:605±11. [11] Lynch CS, Yang W, Collier L, Suo Z, McMeeking RM. Electric ®eld induced cracking in ferroelectric ceramics. Ferroelectrics 1995;166:11±30. [12] Mikhailov GK, Parton VZ. Electromagnetoelasticity. New York: Hemisphere, 1990. [13] Pak YE. Linear electro-elastic fracture mechanics of piezoelectric materials. Int J Fracture 1992;54:79±100. [14] Pan E. A BEM analysis of fracture mechanics in 2D anisotropic piezoelectric solids. Engng Analysis Boundary Elem 1999;23:67±76. [15] Park SB, Sun CT. Fracture criteria for piezoelectric ceramics. J Am Ceramic Soc 1995;78:1475±80 (and correction (1996)). [16] Parton VZ, Kudryavtsev BA. Electromagnetoelasticity. New York: Gordon and Breach Science Publishers, 1988. [17] Qin QH. Thermoelectroelastic analysis of cracks in piezoelectric halfplane by BEM. Comput Mech 1999;23:353±60. [18] Qin QH, Lu M. BEM for crack-inclusion problems of plane thermopiezolectric solids. Int J Numer Meth Engng 2000;48:1071±88. [19] Qin QH, Mai YW. Crack branch in piezoelectric bi-material system. Int J Engng Sci 2000;38:667±87. [20] Rajapakse RKND. Plane strain/stress solutions for piezoelectric solids. Composites Engng: part B 1997;28:385±96. [21] Rajapakse RKND. Boundary element methods for piezoelectric materials. Proc SPIE Conf Smart Struct Mater 1997;3039:418±28. [22] Sack EE. Evaluation of ®nite analysis of piezoelectric materials with holes and cracks. M.Sc Thesis, Purdue University, USA 1994. [23] Snyder MD, Cruse TA. Boundary-integral equation analysis of cracked anisotropic plates. Int J Fracture 1975;11:315±28. [24] Sosa HA. On the fracture mechanics of piezoelectric solids. Int J Solids Struct 1992;29:2613±22. [25] Suo Z. Models for breakdown-resistant dielectric and ferroelectric ceramics. J Mech Phys Solids 1993;41:1155±76. [26] Suo Z, Kuo CM, Barnett DM, Willis JR. Fracture mechanics for piezoelectric ceramics. J Mech Phys Solids 1992;40:739±65. [27] Ting TCT. Anisotropy elasticity: theory and applications. Oxford: Oxford University Press, 1996. [28] Xu X-L, Rajapakse RKND. Boundary element analysis of piezoelectric solids with defects. Composites Engng: part B 1998;29:655±69. [29] Xu X-L, Rajapakse RKND. A theoretical study of branched cracks in piezoelectrics. Acta Materialia 2000;48:1865±82. [30] Xu X-L, Rajapakse RKND. On singularities in composite piezoelectric wedges and junctions. Int J Solids Struct 2000;37:3253±75.