A study of finite size effects on cracked 2-D piezoelectric media using extended finite element method

Computational Materials Science 50 (2011) 1834–1845

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A study of finite size effects on cracked 2-D piezoelectric media using extended finite element method R.R. Bhargava, Kuldeep Sharma ⇑ Department of Mathematics, Indian Institute of Technology Roorkee, Roorkee 247 667, India

a r t i c l e

i n f o

Article history: Received 29 September 2010 Received in revised form 31 December 2010 Accepted 17 January 2011 Available online 12 February 2011 Keywords: Crack Electric-displacement intensity factor Extended finite element method Interaction integral Piezoelectric ceramics Stroh formalism

a b s t r a c t The extended finite element method is applied on the two-dimensional (2-D) finite piezoelectric media weakened by a crack. The fourfold standard enrichment functions are taken in conjugation with the interaction integral to evaluate the intensity factors. Four sequence of analysis, namely crack–mesh alignment, aspect ratio, mesh with local refinement and domain independency is done on the center and edge crack problems. These four analyses when combined together give an optimum result to study the finite specimen. It is observed that for smaller values of strip width to crack length ratio the finiteness of the specimen size affects the intensity factors. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction Piezoelectric ceramics have been widely used due to the property of mechanical–electrical coupling, as key components in sensors/actuators/transducers. However, piezoelectric ceramics are very brittle in nature, therefore fracture analysis becomes imperative for such materials. The mathematical methods used in linear elastic fracture mechanics (LEFM) such as complex function theory (Lekhnitski, Stroh), integral transformation and singular integral equations have been applied/extended to piezoelectric problems using linear approximation. Many two-dimensional standard crack configurations in infinite domains with different possible electric conditions at the crack surface have been treated successfully. However for practical applications and analysis of fracture test specimens, the bounded geometry, the complex electromechanical boundary conditions and material non-linearity requires the numerical analysis as finite elements method (FEM) or boundary elements method (BEM). Allik and Hughes [1] were the first to apply FEM to investigate vibrations in piezoelectric ceramics. Kumar and Singh [2,3], investigated crack propagation and energy release rate under combined mechanical and electrical loadings. While Kuna et al. [4–6] applied FEM analysis to develop electromechanical J-integral, an equivalent domain integral for a crack piezoelectric 2-D and 3-D problems. They considered the case of ⇑ Corresponding author. Tel.: +91 9456318765. E-mail address: [email protected] (K. Sharma). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.01.026

impermeable, permeable and semi-permeable cracks. Kuna [7,8] gave a detailed survey of different techniques to solve piezoelectric fracture mechanics problems and present status of the work done till date in piezoelectrics. The extended finite element method (X-FEM) proposed by Moës et al. [9], Belytschko and Black [10] have proven to be a very efficient tool for the numerical modeling of cracks in LEFM. In comparison to the standard FEM, the X-FEM provides significant benefits in the numerically modeling of crack propagation. The main advantages are that the finite element mesh needs not to conform to the crack boundaries to account for the geometric discontinuity, and furthermore, mesh regeneration is not needed in crack growth simulations. The X-FEM has also been applied to various fields in engineering [11,12]. More recently X-FEM [13] has been applied to investigate enrichment schemes and their convergence for semi-infinite crack and Griffth–Irwin crack weakening a piezoelectric material. The problem investigated using analytic calculation based on Lekhnitski formalism. It has been concluded that fourfold enrichment functions used in isotropic materials is almost efficient, concerning accuracy both in energy and intensity factors (IFs). Implementation of using these fourfold enrichment functions becomes simpler and also involves less computation. To the best of the author’s knowledge, there is paucity in literature on X-FEM solutions for finite piezoelectric domains, particularly those encountered in practical situations. We attempt this paucity, to investigate the effect of finite dimension of specimen on a cracked piezoelectric material.

1835

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

X-FEM has been applied using fourfold enrichment functions under in-plane tensile and electric loading. Two benchmark problems: center and edge crack have been discussed in this paper. The opening mode stress intensity factor (KI) and electric-displacement intensity factor (KIV) have been evaluated using an interaction integral technique [7,8] and the crack tip behavior has been given by Stroh formalism [14,15]. To reduce the computational cost a mesh with local refinement has been taken and also the results of KI and KIV are compared with those of the structured mesh.

2. Basic equations for piezoelectric media As is well-known the governing equations and the boundary conditions which form the foundation of piezoelectric media are given below.

3. Crack tip fields in homogeneous piezoelectric media For cracks in homogeneous piezoelectric media the asymptotic behavior of the field quantities has been given by Sosa [14] and Pak [15]. The electromechanical stress and electrical displacement fields can be written in polar coordinates (r, h) (with the origin at the crack tip) as

K N fijN ðhÞ;

ð9Þ

1 X K N g Ni ðhÞ; Di ðr; hÞ ¼ pffiffiffiffiffiffiffiffiffi 2pr N

ð10Þ

2.1. Field equations

 Constitutive equations

rij ¼ C ijks eks  esij Es ; Di ¼ eiks eks þ jis Es ;

ð1Þ

 Kinematic equations

1 2

eij ¼ ðui;j þ uj;i Þ; Ei ¼ /;i ;

ð2Þ

 Equilibrium equations

rij;j ¼ 0; Di;i ¼ 0;

ð3Þ

where rij, eij, Di and Ei denote the components of the stress, strain, electric displacement, and electric field; Cijks and eiks denote the elastic and piezoelectric constants, respectively; jis denotes the dielectric permitivities. In Eqs. (2) and (3), comma denotes partial differentiation with respect to argument following it; ui is the component of the elastic displacement vector u; / is the electric potential; where i, j, k and s = 1, 2, 3. 2.2. Boundary conditions Consider a piezoelectric medium occupying the space X enclosed by surface S. On the boundaries Sr and SD, the resultant of stresses and electric displacements are respectively

/ðr; hÞ ¼

p

N

p

N

rffiffiffiffiffi 2r X

N

K N di ðhÞ;

ð11Þ

K N mN ðhÞ;

ð12Þ

( ) M ia NaN pa Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; cos h þ pa sin h a¼1 ( ) 4 X M ia N aN fi2N ¼ Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; cos h þ pa sin h a¼1 fi1N ¼ 

4 X

ð13Þ

) M 4a N aN pa ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ¼ Re ; cos h þ pa sin h a¼1 ( ) 4 X M 4a N a N g N2 ¼ Re pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; cos h þ pa sin h a¼1 4 X

g N1

N

4 X

a¼1

ð4Þ ð5Þ

rffiffiffiffiffi 2r X

where i, j = 1,2, and the summation over N = {II, I, III, IV} comprises the fracture opening mode intensity factors KII and KIII denote mode-II, mode-III stress intensity factor, respectively. In case of two-dimensional piezoelectric structure, KIII = 0. The functions N fijN ðhÞ; g Ni ðhÞ; di ðhÞ, and mN(h), are the standard angular functions for a crack in a homogeneous piezoelectric–elastic medium, which depend only on the material properties and can be determined by means of the extended Stroh formalism and semi-analytical calculations. They can be expressed in terms of complex material eigenvalues pa, eigenvectors AMa and matrices MMa and NaN (Kuna [7,8]; Rao and Kuna [16]) as

di ¼

rij nj ¼ t0j ; on Sr ; Dj nj ¼ x0 ; on SD ;

N

and the near tip displacement field and electrical potential can be obtained from

ui ðr; hÞ ¼

In a fixed rectangular coordinate system xj (j = 1, 2, 3), the field equations for a linear piezoelectric medium subjected to electromechanical loads in the absence of body forces and charges are

X

1

rij ðr; hÞ ¼ pffiffiffiffiffiffiffiffiffi 2pr

mN ¼

4 X

(

ð14Þ

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Re Aia NaN cos h þ pa sin h ;  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Re A4a NaN cos h þ pa sin h ;

ð15Þ

a¼1

where t0j is prescribed on surface Sr and x0 is given on surface SD. n is the outward drawn unit normal vector on S. The displacement vector u is prescribed on surface Su and the electric potential / is prescribed on surface S/ as

uj ¼ u0j ;

on Su ;

ð6Þ

/ ¼ /0 ;

on S/

ð7Þ

2.3. Impermeable crack conditions In case of impermeable crack conditions on the upper and lower crack faces Sc,

Dj nj ¼ 0;

on Sc :

ð8Þ

where Re{. . .} denote the real part, of the quantity in brackets. The four conjugate pairs of eigenvalues pa can be obtained by solving the following



  C i2k1 þ C i1k2 ei21 þ ei12 þ p e1k1 j11 e2k1 þ e1k2 j12  j21     Ai C i2k2 ei22 p2 ¼ 0: þ e2k2 j22 A4 C i1k1

ei11



ð16Þ

Only the four eigenvalues having positive imaginary part and the corresponding eigenvectors are used as a column vectors in the matrix A. The matrices MMa and NaN are calculated by

N1 aN ¼ M M a ¼



ðC i2k1 þ C i2k2 pa ÞAka

ðe1i2 þ e2i2 pa ÞA4a

ðe2k1 þ e2k2 pa ÞAka

ðj21  j22 pa ÞA4a

 :

ð17Þ

1836

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

4. Extended finite element approximation Modeling of cracks using FEM is cumbersome in 2-D for complex structures or crack geometries. But X-FEM models crack(s) geometry independent of the mesh leading to a simplification in mesh generation and avoids remeshing as the crack grows. The X-FEM exploits the partition of unity property of FEM first identified in [9], which allows local enrichment functions to be easily incorporated into a finite element approximation. A standard approximation is thus enriched in a region of interest by the local functions in conjunction with additional degrees of freedom. For crack problems the enrichment functions are the near tip asymptotic fields and a discontinuous function to represent the jump in the displacement across the crack boundary. Fig. 1. The set of nodes elected for enrichment.

4.1. Description of crack geometry To represent a crack, Level Set Method (LSM) has been applied which is a numerical method proposed by Osher and Sethian [17]. The principle of the method is to represent an interface by the zero of a function, called level set function, and to update this function with Hamilton–Jacobi equations knowing the speed of the interface in the direction normal to this interface. In two-dimensions, the crack is an open curve bounded by the crack tip(s) and in three-dimensions, it is an open surface bounded by the crack front. An extension of the LSM is proposed by Stolarska et al. [18] to represent open boundaries using two level set functions  One the normal level set function, w1, that is the signed distance to the union of the crack and the tangent extension from its front,  And other the tangent level set function, w2, that is the signed distance function to a surface that passes by the crack boundary and that is normal to the crack. The crack surface is defined as the subset of the zero level set of

parameters. Fk(r, h) is the basis for the Westergaard field for the crack tip, which are defined as

F k ðr; hÞ ¼

where r is the kx  xTIPk; N the set of all nodes in the discretization; NTIP the set of all nodes that are connected to elements containing crack tip(s); and Ncr is the set of nodes that are connected to elements containing the crack but not in NTIP. The set of nodes elected for enrichment are shown in Fig. 1. Substituting the approximate displacement from Eq. (18) and the electric potential from Eq. (19) into the weak form illustrated in [19], the standard discrete system of equations is obtained

K s d ¼ f ext

2

bU kij

Considering the impermeable boundary conditions defined in Section 2.3 at the crack surface, the extended finite element approximation for the displacement and electric potential can be written as follows:

X

NI ðxÞuI þ

X I2N

/h ðxÞ ¼

X

TIP

4 X

X I2NTIP

k

ðF k ðr; hÞ  F k ðxI ÞÞbI ;

ð18Þ

k¼1

NI ðxÞ/I þ

I2N

þ

NI ðxÞ

X I2N

NI ðxÞðHðf h ðxÞÞ  HðfI ÞÞcI



ð19Þ

k¼1

1 if y < 0 þ1 if y > 0

ab 7 kij 7 5

ð22Þ

bb kij

e

UU

ð23Þ

While the element contribution to the global element force vector fext is

ð24Þ

where k k

U ¼ fu/gT ; a ¼ faI cI gT ; b ¼ fbI dI g and Z  r T s  rs Bi C Bj dX ðr; s ¼ U; a; bÞ; i; j ¼ 1; 2; 3; 4; kij ¼ fiU ¼

ð25Þ

ð20Þ

And the shape functions, NI(x), are isoparametric linear quadrilateral element shape functions that construct the partition of unity. The column matrices uI and /I are the nodal displacements and k k electric potential respectively, and aI, bI and cI, dI are the additional

Z Z @ Xe

NitdC þ

Z

Ni fdC;

Xe

Z Ni ðHðf h ðxÞÞ  Hðfi ÞÞtdC þ Ni ðHðf h ðxÞÞ  Hðfi ÞÞfdC; @ Xe Xe Z Z N i ðF k ðxÞ  F k ðxi ÞÞtdC þ Ni ðF k ðxÞ  F k ðxi ÞÞfdC; fib ¼

fia ¼

 :

ba kij

3

Xe

where H(f(x)) is a modified Heaviside step function

HðyÞ ¼

aa

kij

Ub

kij

and in case of non-enriched element kij ¼ kij :

cr

4 X k NI ðxÞ ðF k ðr; hÞ  F k ðxI ÞÞdI ;

Ua

kij



T fie ¼ fiU fia fib and fie ¼ fiU for the non-enriched element:

NI ðxÞðHðf h ðxÞÞ  HðfI ÞÞaI

I2Ncr

I2N

þ

X

UU

kij

6 aU e kij ¼ 6 4 kij

4.2. Enriched approximation

uh ðxÞ ¼

ð21Þ

where fext is the vector of external nodal forces and Ks the stiffness matrix. The elementary stiffness matrix for an enriched element is defined as

w1 where w2 is negative. The crack front is defined as the intersection of the two zero level sets.

  pffiffiffi pffiffiffi h pffiffiffi h pffiffiffi h h r sin ; r cos ; r sin sin h; r cos sin h ; 2 2 2 2

@ Xe

Xe

ð26Þ Ni is the standard finite element shape function defined at node i, and B0i s are the nodal matrices of the shape derivatives. The detailed of this method is given in the following flow chart

1837

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

Percentage error in IFs ¼

ðIFsX-FEM  IFsAnalytically Þ  100; IFsAnalytically

ð30Þ

and

Normalized IFs ¼

IFsX-FEM ; IFsAnalytically

ð31Þ

respectively. In this section we present numerical results for the computation of the intensity factors in piezoelectric material using interaction integral. The interaction integral method is an effective tool for calculating the SIFs and EDIF in the homogeneous piezoelectric materials [16,22]. The path independent electromechanical J-integral for a homogeneous piezoelectric cracked body is given by Cherepanov [23].

J¼

 Z  @ui @/ nj dC; Wd1j  rij  Dj @x1 @x1 C

ð32Þ

R R where W = rijdeij  DidEi is the electric enthalpy density, nj is the jth component of the outward unit vector normal to an arbitrary contour C enclosing the crack tip and dij is the Kronecker delta. For linear piezoelectric material, an equivalent domain form is given as

J¼

Z 

rij

A

 @ui @/ @q þ Dj  Wd1j dA; @x1 @x1 @xj

ð33Þ

where A is the area inside the contour C and q is a smooth weight function chosen such that it has a value of unity at the crack tip, zero along the boundary of the domain C, and smoothly interpolated in between. Consider two independent equilibrium states of the cracked body: let state 1 correspond to the actual state for the given boundary conditions, and let state 2 correspond to an auxiliary state, which can be near tip electromechanical fields of any of the fracture opening modes I, II, III and IV. Superposition of these two states leads to another equilibrium state called (state s) for which the domain form of the J-integral is

J

ðsÞ

¼

Z A

! ð1Þ

 þ /ð2Þ Þ @q ð1Þ ð2Þ @ð/ ðsÞ þ Dj þ Dj  W d1j dA; @x1 @xj

5. Interaction integral For an infinite domain problem in piezoelectric fracture mechanics under tensile stress r1 and electric-displacement D1 on the remote boundary, the mode-I stress intensity factor KI and electric-displacement intensity factor (EDIF) KIV at the crack tip pffiffiffiffiffiffi pffiffiffiffiffiffi are r1 pa and D1 pa respectively. In case of a finite dimensional specimen, a dimensionless correction factor [20,21] that depends on the ratio wa is introduced for the IFs as

K I ¼ C r1

pffiffiffiffiffiffi pffiffiffiffiffiffi pa and K IV ¼ CD1 pa:

ð27Þ

where

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

pa ffi

a 2

a 4  ; 1  0:025 C ¼ sec þ 0:06 2w w w

ð28Þ

a is the semi-crack length and w is the semi-width of the specimen for the center crack problem and

a

a 2

a 3

a 4 C ¼ 1:12  0:231  21:72 þ 30:39 ; þ 10:55 w w w w ð29Þ a is the edge crack length and w is the width of the specimen in case of edge crack problem. The percentage error and normalized value in IFs are calculated by

 0

 @ uið1Þ þ uið2Þ ð1Þ ð2Þ @ r þr ij ij @x1 ð34Þ

where superscript i = 1, 2, and s indicate fields and quantities associated with state i and

W ðsÞ ¼

n







ð2Þ ð2Þ ð1Þ ð2Þ ð1Þ ð2Þ rð1Þ eð1Þ  Dj þ Dj Ej þ Ej ij þ rij ij þ eij

o. 2:

By expanding Eq. (34), we may write

J ðsÞ ¼ J ð1Þ þ J ð2Þ þ Ið1;2Þ ; where J(1) and J(2) are the electromechanical J-integrals for states 1 and 2 respectively, and interaction integral I(1,2) is ! Z ð2Þ ð1Þ ð2Þ ð1Þ @ui @q ð1Þ @/ ð2Þ @ui ð2Þ @/ ð1;2Þ Ið1;2Þ ¼ rð1Þ þ D r þ D  W d dA; 1j ij j j @x1 @x1 ij @x1 @x1 @xj A

In Eq. (35), W

ð1;2Þ

¼

ð1Þ ð2Þ ij ij

ð2Þ ð1Þ ij ij

r e þr e 

ð1Þ ð2Þ Dj Ej



ð2Þ ð1Þ Dj Ej

 =2:

ð35Þ

As is well-known for piezoelectric solids under mixed-mode loading conditions; the electromechanical J-integral can be written (using Kuna [7]) as

J¼

1 T K YK; 2

ð36Þ

1838

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

where K ¼ fK II K I K III K IV gT is the vector of the intensity factors, and Y is the 4  4 generalized Irwin matrix, which depends on the elastic, piezoelectric, and dielectric material constants and is given by

Y MN ¼ ImfAMa NaN g:

ð37Þ

Therefore, the interaction integral for two-dimensional case is given by ð1Þ

ð2Þ

ð1Þ

ð2Þ

ð1Þ

ð2Þ

Ið1;2Þ ¼ K II K II Y 11 þ K I K I Y 22 þ K IV K IV Y 44



 ð1Þ ð2Þ ð1Þ ð2Þ ð1Þ ð2Þ ð1Þ ð2Þ þ K I K II þ K II K I Y 12 þ K II K IV þ K IV K II Y 14

 ð1Þ ð2Þ ð1Þ ð2Þ þ K I K IV þ K IV K I Y 24 :

ryy = 10 MPa, electric displacement Dy = 0.001 C/m2. The material considered for the analysis in this paper is PZT-5H and for which material constants are given in Table 1 taken from [24]. The IFs are calculated on a domain radius equal to 0.8a. 6.1.1. Crack–mesh alignment vs. Non-alignment The analysis of crack–mesh alignment vs. non-alignment has been done on the structured mesh Figs. 4 and 5. The effect of crack–mesh alignment is considered for two possible cases: full alignment and non-alignment. The elements used in the mesh con-

ð38Þ

The individual IFs for the actual state can be obtained by judiciously choosing the auxiliary state (state 2). If state 2 is taken to be (i) For fracture opening mode-I ð2Þ

ð2Þ

ð2Þ

then, K I ¼ 0; K II ¼ 1; and K IV ¼ 0: Hence, Eq. (38) can be reduced to ð1Þ

ð1Þ

ð1Þ

Ið1;IÞ ¼ K I Y 22 þ K II Y 12 þ K IV Y 24 :

ð39Þ

(ii) For fracture opening mode-II then Eq. (38) reduces to ð1Þ

ð1Þ

ð1Þ

Ið1;IIÞ ¼ K I Y 12 þ K II Y 11 þ K IV Y 14 ;

ð40Þ

and (iii) For fracture opening mode-IV ð1Þ

ð1Þ

ð1Þ

Ið1;IVÞ ¼ K I Y 24 þ K II Y 14 þ K IV Y 44 ;

ð41Þ

respectively. ð1Þ

ð1Þ

Solving the simultaneous Eqs. (39)–(41), the IFs K I ; K II and are obtained. The interaction integrals I(1,I), I(1,II) and I(1,IV) are evaluated from Eq. (35). ð1Þ K IV

6. Numerical applications

Fig. 2. Model geometry of the center crack.

Two benchmarks problems in fracture mechanics of piezoelectric materials are considered with the following objectives  Comparison /trends of IFs in case of crack line and tip alignment (full alignment) to the mesh with neither the crack line nor the crack tip align to the mesh (non-alignment).  Influence of aspect ratio on IFs accuracy.  To show the accuracy obtainable on mesh with local refinement having full alignment case.  To study the domain independence in the IFs computations.  Explore finite size specimen effects and trends of the IFs. 6.1. Center-crack in-plane tension and electric loading First case considered is of a center-cracked test specimen under in-plane tension and electrical loading. The model geometry shown in Fig. 2 consists of a homogeneous poled piezoelectric plate of length L = 2h and width 2w. The plate is weakened by a center crack occupies the interval [a, a]. The crack is perpendicular to the polarization direction. The plate is subjected to tensile stress r1 or ryy and electric-displacement D1 or Dy on the remote boundary. Due to symmetry, only half the specimen is considered as shown in Fig. 3 with appropriate displacement boundary conditions to remove rigid-body motion and the edge singularity. The plain strain case has been discussed with tensile loading

Fig. 3. Half plate symmetry for center crack problem.

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

1839

Table 1 Material constants of PZT-5H. Properties

PZT-5H

Elastic constants

c11 = 126 GPa c13 = 53 GPa c44 = 35.3 GPa e15 = 17 C/m2 e33 = 23.3 C/m2 j11 = 15.1 nC/V m

Piezoelectric constants Permittivity

c12 = 55 GPa c33 = 117 GPa e31 = 6.5 C/m2

j33 = 13.1 nC/V m

Fig. 6. Enrichment scheme in full alignment case.

Fig. 4. Structured mesh with full alignment.

Fig. 7. Enrichment scheme in non-alignment case.

Fig. 5. Structured mesh with non-alignment.

sist of uniform four noded quadrilateral isoparametric element. The enrichment scheme used in both the cases is shown in Figs. 6 and 7. Plate dimensions of semi-width w = 2 and length L = 4, respectively, and a half-crack length a = 0.4 are used. Figs. 8 and

9 show the convergence of IFs with respect to the number of elements per unit length. The results of IFs are presented for both the cases, i.e. full alignment and non-alignment to the mesh using X-FEM, and further these results are compared with the FEM or without enrichment. The results of FEM are calculated on the same structured mesh and number of elements. It clearly shows that the results of X-FEM are better than the conventional FEM results. The convergence in case of full alignment and non-alignment has been obtained for 1/he > 50. It is observed that behavior of convergence in KI and KIV is the same vis-a-vis for all the three cases. It is also observed that the results with the non-alignment case are the best one and as the mesh density increases the results of full alignment case converges to non-alignment case in X-FEM. The trends of KI in both the cases are in support of the work [25].

1840

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

for the same number of elements. It is observed that the percentage of increment in Dofs is very small in comparison to the total number of elements employed in the mesh. Fig. 10 shows the computational time taken by the CPU to run the Matlab code for X-FEM and FEM with respect to the mesh density. It is found that full alignment case takes more computational time than the FEM and this difference increases further with the mesh density. Even using the less number of elements in non-alignment case takes almost the same computational time as in the case of FEM. But still the difference is not more than 10% for higher mesh density. The accuracy

Fig. 8. Convergence study of KI w.r.t 1/he for center crack problem.

Fig. 10. Computational comparison of X-FEM and FEM w.r.t 1/he for center crack problem.

Fig. 9. Convergence study of KIV w.r.t 1/he for center crack problem.

The results of IFs are tabulated in Table 2 for different specimen sizes and for a constant value of a/w = 0.2. It is found that the trends of IFs are independent of the specimen size. The results are only deviating with respect to the values of he/a. This may allow the authors to work on a higher dimensions with variation in the values of he/a. As the application of X-FEM involves additional degree of freedoms due to enrichment, so it takes more computational cost than the FEM on the same structured mesh. The comparison of degree of freedoms (Dofs) used in X-FEM and FEM are tabulated in Table 3

Fig. 11. Influence of aspect ratio in the center crack problem.

Table 2 Comparison of IFs w.r.t no. of elements for different specimen sizes in the center crack problem. No. of elements

20  40 40  80 60  120 80  160 100  200

he/a

0.25 0.125 0.0833 0.0625 0.0250

w = 2, L = 4, a = 0.4

w = 5, L = 10, a = 1.0

w = 8, L = 16, a = 1.6

he

% Error in KI

% Error in KIV

he

% Error in KI

% Error in KIV

he

% Error in KI

% Error in KIV

0.1 0.05 0.033 0.025 0.01

0.32 0.43 0.46 0.47 0.48

2.08 2.24 2.34 2.36 2.41

0.25 0.125 0.0833 0.0625 0.05

0.41 0.43 0.46 0.47 0.48

2.13 2.19 2.33 2.36 2.39

0.4 0.2 0.1333 0.1 0.08

0.36 0.43 0.46 0.47 0.48

2.11 2.25 2.34 2.36 2.41

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

1841

Fig. 12. Mesh with local refinement.

Fig. 14. Structured mesh and mesh with local refinement vs. w for the center crack problem.

Fig. 13. Comparison between structured mesh and MESH with local refinement for the center crack problem.

of the results with X-FEM on the structured mesh in comparison to the FEM is always preferable. This will also avoid the computational cost involved in using higher order elements and the proper meshing techniques in conventional FEM.

Fig. 15. Computational comparison between the structured mesh and mesh with local refinement.

6.1.2. Aspect ratio analysis Influence of the plate aspect ratio (i.e. h/w), on the IFs calculations is investigated. The semi-width of the plate is taken as w = 5, the half-crack length is a = 1 and an element length near the tip (he) equal to 0.25. The semi-length of the specimen is then varied from 8 to 30. Fig. 11 shows the percentage error in KI and KIV with respect to aspect ratio and is calculated on the structured mesh for full alignment case. It is observed in Fig. 11, a dramatic reduction in the measured error can be seen around h/w = 1.0, beyond it the percentage error stabilizes to 0.8% in case of KI and 1.0% in case of KIV. The deviation of the X-FEM results from the theory at low aspect ratio is most likely due to the fact that the applied loads can no longer be considered as far-field loadings. The behavior of KI is validated by the work of McNary [25].

The size of the specimen and crack length is taken as w = 5, L = 10 and a = 1.0. The results depicted in Fig. 13 show good accuracy and agreement of IFs in both type meshes. Fig. 14 depicts the variation of KI and KIV with respect to w for structured mesh and mesh with local refinement considering he = 0.1. It is observed in this case also that the results for each of KI and KIV in both the cases are almost identical. Graphs are plotted in Fig. 15 for the total number of elements required in structured mesh and mesh with local refinement with respect to 1/he. It is amply clear from the graphs drawn that mesh with local refinement needs almost one third of the elements as compared to that in case of structured mesh. Consequently the computation cost is drastically reduced with use of mesh with local refinement. And this is one of the objectives of present work. In further computational analysis for center crack the mesh with local refinement has been used.

6.1.3. Structured mesh and Mesh with local refinement Results for KI and KIV are compared for structured mesh shown in Fig. 4 with the mesh having local refinement shown in Fig. 12. The case considered for the crack and its tip align with the mesh.

6.1.4. Domain independency The accuracy of interaction integral on numerical solution for different domain sizes has been tabulated in Table 4. The sample

1842

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

Table 3 Comparison of degree of freedoms used in X-FEM and FEM. No. of elements

20  40 40  80 60  120 80  160 100  200

Dofs

Increment in Dofs due to enrichment (in %)

Fullalignment

FEM

2700 10,092 22,284 39,276 61,068

2595 9987 22,179 39,171 60,963

No. of elements

Dofs

Non-alignment 4.04 1.05 0.47 0.27 0.17

19  39 39  79 59  119 79  159 99  199

2466 9690 21,714 38,538 60,162

Table 4 Domain independence study for center crack. he/a

Radius (rd)

Normalized KI

% Error in KI

Normalized KIV

% Error in KIV

0.25

0.5 0.6 0.7 0.8 0.9

1.0028 1.0030 1.0031 1.0032 1.0032

0.28 0.30 0.31 0.32 0.32

1.0226 1.0229 1.0238 1.0235 1.0235

2.26 2.29 2.38 2.35 2.35

0.4 0.5 0.6 0.7 0.8

1.0042 1.0045 1.0046 1.0044 1.0044

0.42 0.45 0.46 0.44 0.44

1.0232 1.0234 1.0237 1.0235 1.0235

2.32 2.34 2.37 2.35 2.35

0.1

 For the study of finite specimen size, the width of the specimen is the only varying parameter. For a fixed length of an element near the crack tip, the higher values of w require more number of elements in case of structured mesh. Therefore mesh with local refinement is preferred.  Domain independency analysis suggested us to take the domain radius equal to 0.8a for better results of far off crack tip elements. The IFs KI and KIV have been tabulated in Table 5. The normalized values of KI and KIV are calculated with the help of Eqs. (27), (28), and (31). As expected, with an increase in w/a, the numerically computed results obtained for both KI and KIV are collectively best and approximate to the exact solution with less than 0.23% error. For w/a P 20, finite specimen effects are negligible, and hence specimen dimensions of that order or larger can be used with confidence to model the infinite domain problem. 6.2. Edge crack in-plane tension and electric loading

size chosen is the same as in Section 6.1.3 with full alignment case. Results are listed for two different values of he/a = 0.25 and 0.1. The numerically calculated IFs KI and KIV are normalized by using Eqs. (27), (28), and (31). The domain for computation of interaction integral is selected as the domain that falls within a ball of radius rd. In Table 4, the results of the domain independence study are listed. It has been observed that the normalized IFs are approximately the same, consequently domain independency for rd > 0.6 (in case he = 0.25) and rd > 0.4 (in case of he = 0.1) is established.

In case of edge crack problem, a specimen of width w, length L and a crack of length a is taken. A tensile loading r1 = ryy = 10 Mpa and an in-plane electric loading D1 = Dy = 0.001 C/m2 is prescribed. The model geometry and the load specified are shown in Fig. 16. The plain strain case has been considered. The domain, mesh and the data for the analysis of objectives mentioned in Section 6 are the same as in the center crack problem, respectively. The IFs are calculated on a domain radius equal to 0.8a. The only change is in the symmetric boundary conditions and the geometric correction factor [21]. The analysis in this case is almost parallel to that for the center crack when considered half-plate symmetry. Here too study for crack–mesh alignment, aspect ratio, mesh with local refinement and domain independency has been carried out. 6.2.1. Crack–mesh alignment vs. Non-alignment To study the convergence of IFs a specimen having width w = 2, length L = 4 and an edge crack length a = 0.4 has been considered.

6.1.5. Finite size specimen effects Effect of finite size of specimen is carried out including.  The length of the element near the tip is kept fixed which requires the full alignment in case of structured mesh.  As is observed in the aspect ratio analysis and for optimal computational cost, a square specimen gives the best result.

Table 5 Finite specimen effects for center crack. he/a

w/a

Normalized KI

% Error in KI

Normalized KIV

% Error in KIV

0.25 rd = 0.8

5 7 10 20 30 40

1.0032 1.0004 0.9991 0.9980 0.9977 0.9977

0.32 0.04 0.09 0.20 0.23 0.23

1.0235 1.0125 1.0061 1.0011 1.0001 1.0001

2.35 1.25 0.61 0.11 0.01 0.01

0.1 rd = 0.8

5 10 20 30

1.0044 1.0009 0.9997 0.9997

0.44 0.09 0.03 0.03

1.0235 1.0057 1.0005 1.0005

2.35 0.57 0.05 0.05

0.05 rd = 0.8

5 10

1.0045 1.0011

0.45 0.11

1.0238 1.0061

2.38 0.61

Fig. 16. Model geometry of edge crack.

1843

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

Figs. 17 and 18 depict the variation of percentage error of KI and KIV with respect to the number of elements per unit length in full alignment and non-alignment case. The results are also compared with the FEM or without enrichment on the same structured mesh and number of elements. From the graphs as could be observed that the values of KI show a good convergence rate whereas KIV has almost the constant values and have a very low convergence rate in respect to all the cases for analysis. This is somewhat different from the center crack problem.

It is also found that the values of IFs in full alignment case and non-alignment case are not differing significantly vis-a-vis each other. The values of KI are better in full alignment case than the non-alignment case. This is not observed in the center crack problem. The trends of KI for edge crack problem are in agreement with [25]. Table 6 shows the results of IFs with respect to the number of elements per unit length for different specimen sizes. It is observed that as in the case of center crack here too the trends of IFs are same and they are independent of the specimen size. The only variation in IFs are found with respect to he/a. So, in further sections a structured mesh with full alignment case having specimen size of w = 5, L = 10 and an edge crack length a = 1.0 has been taken for the analysis. Fig. 19 shows the computational time taken by the CPU to run the Matlab code for X-FEM and FEM with respect to the number of elements per unit length. It is found that in case of edge crack problem the computational time taken is less than the center crack. But the behavior of computational time taken in all the cases is similar to the center crack problem. Again the key issue is the accuracy of the results obtained by X-FEM on a structured mesh. 6.2.2. Aspect ratio analysis Fig. 20 illustrates the influence of aspect ratio on the IFs KI and KIV. It is observed that for aspect ratio, h/w P 1, the variation in the results for KI and KIV is almost constant. It is noted that this behavior matches with that of the center crack case. Again, the more percentage error at low aspect ratios is most likely due to the decreased validity of far-field loadings. It has been also supported by the analysis done in [25].

Fig. 17. Convergence study of KI w.r.t 1/he for edge crack problem.

6.2.3. Structured mesh and mesh with local refinement Figs. 21 and 22 show the behavior of structured mesh and mesh with local refinement with respect to mesh density and width of the specimen, respectively. The results for each of KI and KIV are in good agreement for structured and mesh with local refinement. And as the mesh density is increased the results for two cases match with each other. It has been concluded that mesh with local refinement is independent of the nature of the problems considered. 6.2.4. Domain independency Normalized values of IFs KI and KIV are shown in Table 7. Two cases of he/a are considered. Results show that for rd > 0.5 (in case of he = 0.25) and rd > 0.4(in case of he = 0.1) are very near to each other and therefore may be considered to be domain independent. The domain independence was also observed in the center crack case. 6.2.5. Finite size specimen effects As is done in the case of center crack case, here too the first four objectives implemented for edge crack problem have been applied on the finite specimen effect analysis. The results are given in Table

Fig. 18. Convergence study of KIv w.r.t 1/he for edge crack problem.

Table 6 Comparison of IFs w.r.t no. of elements for different specimen sizes in the edge crack problem. No. of elements

20  40 40  80 60  120 80  160 100  200

he/a

0.25 0.125 0.0833 0.0625 0.0250

w = 2, L = 4, a = 0.4

w = 5, L = 10, a = 1.0

w = 8, L = 16, a = 1.6

he

% Error in KI

% Error in KIV

he

% Error in KI

% Error in KIV

he

% Error in KI

% Error in KIV

0.1 0.05 0.033 0.025 0.01

0.43 0.21 0.18 0.17 0.14

45.15 45.87 45.99 46.00 46.12

0.25 0.125 0.0833 0.0625 0.05

0.51 0.26 0.18 0.17 0.14

44.77 45.7 45.99 46.00 46.12

0.4 0.2 0.1333 0.1 0.08

0.43 0.21 0.18 0.17 0.14

45.15 45.87 45.99 46.01 46.13

1844

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845

Fig. 19. Computational comparison of X-FEM and FEM w.r.t 1/he for the edge crack problem.

Fig. 20. Influence of aspect ratio in the edge crack.

Fig. 22. Structured mesh and mesh with local refinement vs. w for the edge crack problem.

8. In case of edge crack the results are quite interesting. For he = 0.25, with the increase in the values of w/a, results of IFs KI and KIV approximate to the exact solution for infinite domain. As the length of the element near the tip is decreased, the percentage error in KI increases and in KIV decreases with the increasing values w/a. It is to be noted that the amount of increment in KI is very less than the decrement in KIV. It may be remembered that for piezoelectric ceramics the collective behaviors of KI and KIV is very important. It is also to be noted that for w/a > 20, the percentage error decreases with increasing values of w/a, for both IF’s KI and KIV. In nutshell, it may be emphasized that the results are collectively better as w/a ratio increased. For still higher values of w/a P 20, the domain of the specimen may be considered as an infinite one. As the value of percentage error in KIV is higher than the center crack case, so the results of X-FEM are also compared with the Wang and Mai [26] for different electrical loadings and a/w ratio. The specimen size and crack length has been taken as in Section 6.2.1. A structured mesh with full alignment case and an element length he = 0.02 is also considered for the X-FEM analysis. Figs. 23 and 24 ensure that the results of X-FEM are similar to the results of [26] for all the applied loadings. Fig. 23 and Table 8 show that the normalized values of KIV increases with the increasing values of a/w and the applied electrical loadings. It is amply clear from Fig. 24 that KI is independent of the applied electric loadings. It may also conclude that in edge crack problem KIV has more affect of finite specimen and applied electrical loadings than the center crack.

Table 7 Domain independence study for edge crack.

Fig. 21. Comparison between structured mesh and mesh with local refinement for the edge crack problem.

he/a

Radius (rd)

Normalized KI

% Error in KI

Normalized KIV

% Error in KIV

0.25

0.5 0.6 0.7 0.8 0.9

0.9871 0.9869 0.9872 0.9875 0.9875

1.29 1.31 1.28 1.25 1.25

1.4293 1.4265 1.4275 1.4288 1.4288

42.93 42.65 42.75 42.88 42.88

0.1

0.4 0.5 0.6 0.7 0.8

0.9964 0.9965 0.9965 0.9964 0.9964

0.36 0.35 0.35 0.36 0.36

1.4547 1.4550 1.4549 1.4548 1.4548

45.47 45.50 45.49 45.48 45.48

R.R. Bhargava, K. Sharma / Computational Materials Science 50 (2011) 1834–1845 Table 8 Finite specimen effects for edge crack.

1845

7. Conclusions

he/a

w/a

Normalized KI

% Error in KI

Normalized KIV

% Error in KIV

0.25 rd = 0.8

5 10 20 30 40 50

0.9875 0.9982 1.0014 1.0009 1.0003 0.9998

1.25 0.18 0.14 0.09 0.03 0.02

1.4288 1.2700 1.2155 1.2025 1.1972 1.1944

42.88 27.00 21.55 20.25 19.72 19.44

0.1 rd = 0.8

5 10 20 30

0.9964 1.0047 1.0072 1.0066

0.36 0.47 0.72 0.66

1.4548 1.2897 1.2333 1.2202

45.48 28.97 23.33 22.02

0.05 rd = 0.8

5 10

0.9980 1.0059

0.20 0.59

1.4596 1.2936

45.96 29.36

The conclusions are  It is noted that the results obtained with mesh having local refinement can predict the behavior extremely well and having a good computational advantage.  A finite specimen effect is observed in both the center and edge crack problems. Even the latter case has a stronger effect than the center crack. From our computation we observe that (i) For center cracked specimen a ratio w/a P 10 gives a good accuracy. (ii) For edge cracked specimen a ratio w/a P 20 gives a good accuracy.  The percentage error in KIV is higher than the percentage error in KI in case of edge crack problem. This behavior is also validated with Kuna [13]. The results presented for an edge crack using X-FEM are also similar to the Wang and Mai [26].  Finally, our analysis confirms the statement of Kuna [13] that the fourfold enrichment functions used in the X-FEM are efficient and accurate to provide the IFs and energy release rate for piezoelectric materials.

Acknowledgements The authors are grateful to Prof. R.D. Bhargava (Senior Professor and Head, Retd., Indian Institute of Technology Bombay, Mumbai, India) for the encouragement throughout the course of this work. References

Fig. 23. Variation of normalized value of KIV w.r.t a/w for different electric loading conditions in case of edge crack problem.

Fig. 24. Variation of normalized value of KI w.r.t a/w for different electric loading conditions in case of edge crack problem.

[1] H. Allik, T.J.R. Hughes, International Journal for Numerical Methods in Engineering 2 (1970) 151–157. [2] S. Kumar, R.N. Singh, Acta Materialia 44 (1996) 173–200. [3] S. Kumar, R.N. Singh, Acta Materialia 45 (1997) 849–868. [4] M. Abendroth, U. Groh, M. Kuna, A. Ricoeur, International Journal of Fracture 114 (2002) 359–378. [5] K. Wippler, A. Ricoeur, M. Kuna, Engineering Fracture Mechanics 71 (2004) 2567–2587. [6] A. Ricoeur, M. Enderlein, M. Kuna, Archive of Applied Mechanics 74 (2005) 536–549. [7] M. Kuna, Archives in Applied Mechanics 76 (2006) 725–745. [8] M. Kuna, Engineering Fracture Mechanics 77 (2010) 309–326. [9] N. Moës, J. Dolbow, T. Belytschko, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. [10] T. Belytschko, T. Black, International Journal for Numerical Methods in Engineering 45 (1999) 601–620. [11] A. Yazid, N. Abdelkader, H. Abdelmadjid, Applied Mathematical Modelling 33 (2009) 4269–4282. [12] T. Belytschko, R. Gracie, G. Ventura, Modelling and Simulation in Materials Science and Engineering 17 (2009) 043001. [13] E. Bechet, M. Scherzer, M. Kuna, International Journal for Numerical Methods in Engineering 77 (2009) 1535–1565. [14] H. Sosa, International Journal of Solids and Structures 29 (1992) 2613– 2622. [15] Y.E. Pak, International Journal of Fracture 54 (1992) 79–100. [16] B.N. Rao, M. Kuna, International Journal of Solids and Structures 45 (2008) 5237–5257. [17] S. Osher, J.A. Sethian, Journal of Computational Physics 79 (1988) 12–49. [18] M. Stolarska, D.L. Chopp, N. Moës, T. Belytschko, International Journal for Numerical Methods in Engineering 51 (2001) 943–960. [19] V. Piefort, A. Preumont, Finite Element Modeling of Piezoelectric Structures, Active Structures Laboratory, ULB-CP 165/42. [20] Y. Mukamai, Stress Intensity Factors Handbook, Pergamon Press, 1987. [21] H. Ewalds, R. Wanhill, Fracture Mechanics, Edward Arnold, 1984. [22] M. Enderlein, A. Ricoeur, M. Kuna, International Journal of Fracture 134 (2005) 191–208. [23] G.P. Cherepanov, Journal of Applied Mathematics and Mechanics 41 (1977) 397–410. [24] T.J.C. Liu, Theoretical and Applied Fracture Mechanics 51 (2009) 102–110. [25] M.J. McNary, Implementation of the XFEM in the ABAQUS Software Package, M.Sc. Thesis, Georgia Institute of Technology, 2009. [26] B.L. Wang, Y.-W. Mai, International Journal of Solids and Structures 39 (2002) 4501–4524.