Application of The X-FEM To The Fracture of Piezoelectric Materials

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2009; 77:1535–1565

Published online 21 September 2008 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme.2455
Application of the X-FEM to the fracture of

piezoelectric materials
E. Béchet1, ∗, † , M. Scherzer2, ‡ and M. Kuna2

1 Laboratoire de Physique et Mécanique des Matériaux, Université de Metz, UMR CNRS 7554, Ile du Saulcy,
57045 Metz Cedex 1, France
2 Institut für Mechanik und Fluiddynamik, TU Bergakademie Freiberg, Lampadiusstraße 4,
09596 Freiberg, Germany
SUMMARY
This paper presents an application of the extended finite element method (X-FEM) to the analysis of
fracture in piezoelectric materials. These materials are increasingly used in actuators and sensors. New
applications can be found as constituents of smart composites for adaptive electromechanical structures.
Under in service loading, phenomena of crack initiation and propagation may occur due to high electrome-
chanical field concentrations. In the past few years, the X-FEM has been applied mostly to model cracks
in structural materials. The present paper focuses at first on the definition of new enrichment functions
suitable for cracks in piezoelectric structures. At second, generalized domain integrals are used for the
determination of crack tip parameters. The approach is based on specific asymptotic crack tip solutions,
derived for piezoelectric materials. We present convergence results in the energy norm and for the stress
intensity factors, in various settings. Copyright q 2008 John Wiley & Sons, Ltd.
Received 30 August 2007; Revised 8 July 2008; Accepted 22 July 2008
KEY WORDS: finite elements; extended finite element method; piezoelectric materials; convergence;
crack
1. INTRODUCTION
Materials showing a strong piezoelectric effect are used in many applications, where they serve
as sensors, actuators or transducers. These applications range from sub-millimeter length scales
in micro-electro-mechanical systems up to large scales in the design of smart electromechanical
structures as, for example, wings in the aerospace industry. As for regular materials subjected to
high mechanical stresses, the knowledge of fracture behavior for these smart materials is often
∗ Correspondence to: E. Béchet, Laboratoire de Physique et Mécanique des Matériaux, Université de Metz, UMR
CNRS 7554, Ile du Saulcy, 57045 Metz Cedex 1, France.
†
E-mail: eric.bechet@univ-metz.fr
‡
M. Scherzer has passed away while we were working on this paper after the review process. We wish pay tribute
and dedicate this article to his living memory.
Copyright q 2008 John Wiley & Sons, Ltd.

1536 E. BÉCHET, M. SCHERZER AND M. KUNA
crucial within the design of parts under high electrical and mechanical loading. The extended finite
element method (X-FEM) [1] has been originally designed for crack growth analysis in isotropic
elastic materials. In association with level sets [2] as a mean for representing the crack geometry,
it is a powerful way to get rid of the costly constrained remeshing needed with conventional
finite element techniques. Under certain circumstances, this method is also able to achieve regular
convergence rates in the energy error even for a cracked domain [3]. The X-FEM has also been
applied to various fields in engineering, one can cite phase transformations [4], modeling of
dislocations [5], fracture in shells [6], boundary layers [7] and many others.
Piezoelectric materials exhibit a coupling between mechanical and electrical variables, as well
as transversely isotropic material characteristics. This leads to completely different near-crack
tip fields of the mechanical and electrical quantities. Moreover, the electrical variables show a
singularity as well. The fundamentals of piezoelectric fracture mechanics can be found in [8–11].
For the treatment of cracks in piezoelectric structures with the finite element method (FEM), many
of the ‘classical’ linear elastic techniques have been extended to coupled electromechanical crack
analysis, see [12]. To treat this new problem in the X-FEM, specific crack analysis tools and
changes in the enrichment functions are needed, which allow to represent the crack—both in the
mechanical and electrical functional spaces. It is worth noting that the methodology presented here
may be applied to anisotropic elastic materials such as composites as well.
2. PIEZOELECTRIC MODEL
We recall the linear material equations of piezoelectricity under small strain assumptions. The
material (PZT4) considered here for the numerical experiments may be either part of an actuator
or a sensor (eg. see [13]); however, devices that usually withstand higher loads are actuators. The
cracked domain is denoted by and its boundary *∪C + ∪C − by B. The boundary B is subjected
to boundary conditions: At the part Bn tractions or the electric field is imposed (natural boundary
conditions), whereas at part Be the displacement u i or the electrical potential is prescribed
(essential boundary conditions). The crack faces C + and C − are assumed to be both traction-free
and electrically impermeable. In the following, units are those from the international system of
units (SI)

i j n j = 0
on C − ∪C + (1)
Di n j = 0

i j n j = t j
on Bn (2)
Di n j = q

u i = u i0
in Be (3)
= 0
The Cauchy stress tensor i j and the electric displacement D j must satisfy the equilibrium
equations:

i j, j = 0
on (4)
D j, j = 0
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2009; 77:1535–1565
DOI: 10.1002/nme
APPLICATION OF THE X-FEM 1537
The strain εi j and electric field E j are derived from the displacement u i and the electric potential :
εi j = 12 (u i, j +u j,i ) (5)
E j = −, j (6)
The constitutive law is defined by
i j = ci jkl εkl −eki j E k (7)
Di = eikl εkl +ik E k (8)
where ci jkl is Hooke’s elasticity tensor, i j is the dielectric tensor and ei jk is the piezoelectric
tensor. These tensors are known experimentally for various kinds of piezoelectric materials. They
are usually not isotropic.
3. ANALYTICAL SOLUTIONS
In the following, analytical solutions will be recalled for a semi-infinite crack in an infinite
piezoelectric domain. The asymptotic expansion will be used to generate enrichment functions
that span the exact solution at crack tips for every loading case. In addition, we will present the
analytical solution for a Griffith–Irwin crack in an infinite domain. The latter will be used in
the sequel to compare enrichment strategies, by evaluating the convergence rate for an increasing
density of the finite element mesh.
3.1. Asymptotic expansion for a semi-infinite crack

In order to make the paper self-contained, we will briefly look back on the extraction of the
asymptotic eigenfunction expansion at crack tips in piezoelectric materials. The details can be
found, for instance, in [10, 14]. Thus, the concrete forms of the general equations (1)–(8) for
their application in the X-FEM will be derived. We will restrict our analysis to plane strain
conditions (εx z = ε yz = εzz = E z = 0) for transversely isotropic piezoelectric materials with respect
to the material axes x and y (poling direction) as considered in [15]. Thus, the constitutive equations
(7) and (8) are specified to
⎛ ⎞ ⎛ ⎞
⎧ ⎫ a11 a12 0 ⎧ ⎫ 0 b21
⎪ ⎬ε x x ⎪ ⎜ ⎪ x x ⎪ ⎜
⎨ ⎟⎨ ⎬ ⎟
⎜ a12 a22 0 ⎟ ⎜ 0 b22 ⎟ Dx
ε yy = ⎜ ⎟ yy + ⎜ ⎟ (9)
⎪
⎩ ⎪ ⎭ ⎝ a33 ⎠ ⎪
⎩ ⎪
⎭ ⎝ b13 ⎠ Dy
εx y 0 0 x y 0
2 2
⎧ ⎫
⎪ x x ⎪
Ex 0 0 b13 ⎨ ⎬ 11 0 Dx
=− yy + (10)
Ey b21 b22 0 ⎪ ⎩ ⎪
⎭ 0 22 Dy
x y
The equilibrium equations are
*x x *x y *x y * yy *Dx *D y
+ = 0, + = 0, + =0 (11)
*x *y *x *y *x *y
DOI: 10.1002/nme
and the compatibility equations for the deformations and electric field components read as
2 2 2
* εx x * ε yy * ε yx *E x *E y
+ =2 , − =0 (12)
*y 2 *x 2 *x*y *y *x
The relations (9)–(12) represent the basic equations for our next analytical developments. Equations
(12) are the integrability of (5) and (6) regarding the displacements and the electric potential. The
material coefficients a11 , . . . , a33 , b21 , . . . , 11 and 22 follow from the general three-dimensional
properties by means of the formulae given in Appendix A of this paper. Let us now consider the
general solution of the piezoelectric boundary value problem within an infinite domain based on
the material coordinates x and y (poling axis). Introducing the potentials U (x, y) and (x, y) [15]
with the definitions
x x = U (x, y),yy , yy =U (x, y),x x , x y = −U (x, y),x y (13)
Dx = (x, y),y , D y = −(x, y),x (14)
the equilibrium equations (11) are fulfilled automatically. Using the material equations (9) and
(10), in which the stresses and the dielectric displacements are expressed through U (x, y) and
(x, y), the compatibility equations (12) can be reduced to a partial differential equation of sixth
order for U (x, y) (see [15]):
L4 [L2 (U )]+L3 [L3 (U )] = 0 (15)
Equation (15) is the extension of the classical two-dimensional Goursat’s equation of elasticity
(see e.g. [16]) to piezoelectric analysis. The operators L2 , L3 and L4 are defined by means of
2 2 3 3
* * * *
L2 = 22 +11 2 , L3 = b22 +(b12 +b13 )
*x 2 *y *x 3 *x*y 2
(16)
4 4 4
* * *
L4 = a22 4 +a11 4 +(2a12 +a33 ) 2 2
*x *y *x *y
in [15]. The solution of (15) can be given in the form
U (x, y) =U (x +ay) with a = a (re) +ia (im) (17)
√
where a (re) and a (im) represent the real and imaginary part of a. i (= −1) is the imaginary unit.
After inserting U (x +ay) into (15), the characteristic equation of the differential equation (15) is
obtained in terms of a:
a11 11 a 6 +[a11 22 +(2a12 +a33 )11 +b12 (b12 +2b13 )+b13
2
]a 4
+[a22 11 +(2a12 +a33 )22 +2b22 (b12 +b13 ]a 2 +(a22 22 +b22
2
)=0 (18)
Each root of (18) ai (i = 1, 2, . . . , 6) for a realizes the fulfillment of (15) based on (17). Thus, the
general solution for U (x, y) is built-up by means of the six roots a1 , a2 , . . . , a6 , following from
the characteristic equation for a through the linear combination:

6
U (x, y) = Ui (x +ai y) (19)
i=1
DOI: 10.1002/nme
The Ui (z i ) in (19) are functions of z i = x +ai y. It is clear that (x, y) has the same structure
as U (x, y) in (19) with the corresponding functions i (z i ). The potential (x, y) follows from
U (x, y) by usual integration, analogous to the deductions given in [15], that is, using the potential
representations (13) and (14) together with (10), the second equation of (12) results for each Ui (z i )
and i (z i ) in
*i (z i ) (b21 +b13 )ai2 +b22 *2 Ui (z i )

=− (20)
*z i 11 ai2 +22 *z i2
Based on (19) and (20), all mechanical and electrical variables can be calculated only by means
of the potential U (x, y) (Figure 1).
In the next step we search for the solution around the tip of a straight crack in an infinite plane,
see Figure 2. To allow for arbitrary poling directions, the material coordinate system (x, y) is
rotated by an angle with respect to the crack axis (, ). In order to fulfill the boundary conditions
at crack faces, it is necessary and convenient to expand the general solution in Laurent-like series
Ω n
+
C
C n Be
u, ϕ
0 0
Bn
t,q
Figure 1. Cracked body with boundary conditions.
α
poling
y
r x
ω θ
α
Figure 2. Definition of the material axes around the crack tip.
DOI: 10.1002/nme
(see [16, 17]) using general power functions for Ui (z i ). The origin of the coordinates x and y is
taken at the crack tip:

6
U (x, y) = di (
k )(x +ai y)
k +2 (21)
k i=1
The complex variables di (

k ) are free coefficients of the series expansion at the crack tip. They can
be determined only from the overall solution of the considered boundary value problem. The values
k are generally complex. They represent the infinite discrete number of roots of the solvability
equation for the given conditions at the crack tip. Within this paper, we will restrict ourselves
to classical homogeneous conditions (vanishing normal and shear stresses and vanishing normal
dielectric component at the crack surfaces). We know that these boundary conditions do not reflect
the reality of piezoelectric crack surfaces. In general, these conditions are non-linear with respect
to dielectric displacements and the electric potential [18], but within this paper these features are
of minor importance. Our intentions are directed towards the application of the X-FEM. First
numerical experiences based on these complicated boundary conditions can be found, for instance,
in [12, 19]. This subject should be analyzed within another publication.
In order to fulfill the classical homogeneous boundary conditions, one has to represent the stress
function, the stresses and the dielectric displacements in a polar coordinate system (r, ) with = 0
at the crack prolongation ligament as it is usually assumed. Another point consists in the fact that
the material properties considered give six characteristic roots a1 , a2 , a3 , a1 , a2 , a3 representing
three conjugate complex pairs. Owing to the real property of U (x, y) =U (r, ), the corresponding
terms of ai and ai in (21) have to be combined mutually for real
k . Thus, we make an ansatz for
U (r, ) with respect to (21) in the form:

3
U (r, ) = di (
)(r cos(−)+ai r sin(−))
+2
i=1

3
+ di (
)(r cos(−)+ai r sin(−))
+2 (22)
i=1
In (22)
is assumed to be real and it represents one of
k in (21). Equation (22) contains six
real constants based on the three complex constants di (
). The real representation of each term in
(22) for each pair of ai and ai (i = 1, 2, 3) gives expressions of the form ( = −):

ei p
+2 cos (
+2) + + f i p
+2 sin (
+2) +
2 2
with

(re) (im) 2 (re)
p = r (cos )2 +ai sin(2)+(sin )2 [(ai ) +(ai )2 ]
(re) (23)
cos()+ai sin() (re) (im)
= arctan (im)
, di (
) = ei (
)+i f i (
), ai = ai +iai
ai sin()
Based on the ansatz (22), which fulfills the equilibrium equations, the homogeneous boundary
conditions at the crack flanks ( (r, = ± ) = r (r, = ± ) = D (r, = ± ) = 0) define a linear
homogeneous equation system for the six unknown coefficients ei , f i . For the existence of a
DOI: 10.1002/nme
non-vanishing solution of this system, one gets the corresponding solvability equation leading to
the unknown
.§ Generally, an infinite discrete number of
, i.e.
k (k = 1, 2, . . .), can be calculated.
In addition, the appropriate eigenvectors must be determined for the solution representation (21).
Thus, the general series solution for the crack tip except the coefficients of the eigenfunctions, i.e.
the K-factors, etc., is constructed. For the considered homogeneous crack face boundary conditions,
one gets for
k the values:

1 = − 12 ,
2 = 0,
3 = 12 ,
4 = 1,
5 = 32 , . . . (24)
The root
1 = − 12 generates three independent eigenvectors based on the coefficients e1 (
1 ), f 1 (
1 ),
e2 (
1 ), f 2 (
1 ), e3 (
1 ) and f 3 (
1 ) as a consequence of the boundary conditions. In this manner,
one can construct three independent
√ singular eigenfunctions based on (13), (14), (20), (22) and
(23) having the classical 1/ r -crack tip singularity. It is interesting to note that the eigenfunctions
depend on the angle = − and not on and explicitly, and so do six linear independent
functions as ingredients of these eigenfunctions, which can be identified and used as enrichment
functions for the X-FEM. This fact is important, because the enrichment functions can be calculated
one and for all for every poling directions with respect to the crack path.
3.2. Analytical solutions for a Griffith–Irwin crack in an infinite medium

In order to assess the influence of the piezoelectric enrichment functions and to study their
convergence, we need a complete analytical solution for a crack within a certain boundary value
problem. The use of the eigensolutions presented is not advantageous in this context. This is
because these solutions, or a part of them, will be applied as enrichment functions for the X-FEM.
Thus, for proving the developed numerical techniques it is necessary to use other independent
analytical solutions. Here, solutions for an Griffith–Irwin crack in an infinite piezoelectric medium
will be used. The loads are given at infinity as shown in Figure 3. We apply transversely isotropic
piezoelectric material (9) with arbitrary poling direction with regard to the crack, i.e. the material
axes 0x y are in the same manner inclined to the crack tip as in the section above.
In Appendix A the general solution procedure is presented for the symmetric piezoelectric
Griffith problem without shear at infinity and poling direction perpendicular to the crack. The
general cases with arbitrary poling direction and arbitrary stresses and dielectric displacements at
infinity can be solved in the same manner as shown in Appendix A, but the algebraic manipulations
are more complicated. Regarding these general cases one can find in the literature [15, 20, 21]
the application of the classical approach of Lekhnitskii [22], whereby the solution is given for
an elliptical void with corresponding conform mapping and by letting the small axis of the void
approach zero. In these publications it is stated that the rigid body displacements can be assigned
to any of the coefficients of the linear analytic function part used. In order to prove the numerical
techniques developed in this paper we found that this statement cannot be accepted and will be
shown evidently in Appendix A for the considered example there. This fact is important for the
general correct solution representation even if the mentioned coefficients do not have essential
influence on the K-factors. They are connected with the values at infinity and involved in the specific
kinematic conditions at the crack tips and thus have to be approached in a correct physical sense.
In Section 3.1 the material axes (x, y) were used for the general solution representation. In order
to get an asymptotic expansion at crack tips, this solution was transformed into the usual polar
§ If
is a complex number one has to combine the terms with the corresponding conjugate complex terms.
DOI: 10.1002/nme
Figure 3. Griffith–Irwin crack inside an infinite piezoelectric solid.
co-coordinate system at the crack tip. For the crack problem considered here with the boundary
conditions shown in Figure 3 an optimal solution deduction can be given in the coordinate system
0 closely connected with the crack (Figure 3). In this manner, one has to rewrite the material
equations (9) for this coordinate system (0 ). Independent from the fact that within this coordinate
system the material behavior is anisotropic, i.e. more complicated as the transversely isotropic
case, the corresponding potentials for the stresses (13) and for electric displacements (fields) (14),
i.e. U (, ) and (, ) with the definitions,
= U (, ), , =U (, ), , = −U (, ),
(25)
D = (, ), , D = −(, ),
can be used in order to fulfill the equilibrium equations. In the same manner as considered in
Section 3.1, we get one partial differential equation of sixth order for the stress function U (, ) as
a consequence of the compatibility equations. Further, one applies U (, ) in the form U (, ) =
U (+a ) and the resulting partial differential equation gives a more complicated characteristic
equation for a having six roots: a1 , a 1 , a2 , a 2 , a3 and a 3 . Exploiting the real value of U (+a ),
the stresses and the dielectric displacements result in

3 3 3
= 2RE ak2 k (k ) , = 2RE k (k ) , = −2RE ak k (k )
k=1 k=1 k=1
(26)

3
3
D = 2RE ak k k (k ) , D = −2RE k k (k )
k=1 k=1
DOI: 10.1002/nme
RE denotes the real part of a complex number. In (26), the dependencies k (k ) stand for holomor-
phic functions [16] of the corresponding complex arguments k = +ak with the characteristic
roots ak . It is clear that the analytic functions k (k ) are closely related to the stress function
U , that is, they represent the second derivatives of the corresponding functions of Ui in (19).
For the deduction of (26), one starts with the sum of three analytic functions k (k ) for and
accordingly, , , D and D follow from (25) and from the extended relations of the form
(20). k are constants depending on the material parameters, the angle and the characteristic
roots. Because of the complexity, the precise expressions of these constants will not be given here.
Knowing the functions k (k ) the boundary value problem is solved. The deformations ε , ε ,
ε and the electric field E and E are defined from the stresses and the dielectric displacements
by means of the material equations. The displacements and the electric potential follow after
integration from ε , ε , ε , E and E .
To guarantee the fulfillment of the boundary conditions at the crack flanks and the conditions at
infinity (constant stresses and constant dielectric displacements), it is obvious that the holomorphic
functions k (z k ) must have the form:
Ak (re) (im) (re) (im)

k (k ) = k + Bk , Ak = Ak +iAk , Bk = Bk +iBk (27)
2k −l 2
l means the half crack length in (27). The representations (27) are in accordance with the results
achieved by means of the conform mapping technique used in [15, 20, 21]. It is worthwhile to note
that the solutions of (27) are guessed in the same manner as it was done for the Griffith crack
in classical elasticity. The radical analytic functions are able to fulfill the boundary condition on
the crack faces and the other part of the solution give constant values of stresses and dielectric
displacements at infinity. Clearly, the equilibrium equations are fulfilled. Thus, the 12 arbitrary real
(re) (im) (re) (im)
constants (Ak , Ak , Bk , Bk , k = 1, 2, 3) in (27) result from the boundary conditions given
in Figure 3. The determination of these constants requires some tireless algebraic manipulations,
within which one has to take care for the exclusion of the rigid body motions and for the fixing of the
electric potential with respect to a constant value. Because of the complexity, the concrete relations
will not be given here, but they will be used in comparing numerical solutions obtained with the
X-FEM. As already mentioned, a special case of these computations is given in Appendix A.
Regarding the solution of the boundary value problem just explained, another interesting point
consists in the fact that the values at infinity cannot be given arbitrarily, otherwise it is possible
to get physically not acceptable overlapping of the crack faces, which was already mentioned
in [20], and piercing shear deformations. These facts are connected with the already mentioned
determination of the coefficients representing the rigid displacements and will be discussed in a
forthcoming publication.
3.3. Results of K-factor (SIF) calculations for inclined piezoelectric Griffith crack configurations
In order to prove the efficiency of the following numerical technique in calculating stress intensity
factors (SIFS) for cracks within piezoelectric materials, K-factors are calculated for different
inclined cracks (angle ) with respect to the ∞ ∞ ∞
yy − D y -loading (x y = 0) along the poling direction
as shown in Figure 4. The results of the calculations are based on the solutions of Section 3.2
and presented in Table I. Thus, the solution is given in the co-ordinate system 0 for which the
DOI: 10.1002/nme
Figure 4. Griffith crack inclined to the load and poling directions.
Table I. K-factors of the inclined crack.

0 0.1 0.2 /4 0.3 0.4 /2
√ KI ∞ 1.772454 1.603199 1.160086 0.8862269 0.6123677 0.1692543 0
m yy
√ K II∞ 0 −0.5209111 −0.8428519 −0.8862269 −0.8428519 −0.5209111 0
m yy
√ K IV∞ 1.772454 1.685704 1.433945 1.253314 1.041822 0.5477184 0
m Dy
The values coincide for both tips. l = 1 m.
values at infinity can be calculated by means of
∞ ∞
= yy ·(cos ) ,
2
∞ ∞
= − yy ·sin cos , ∞ ∞
= yy ·(sin )
2
(28)
D ∞ = D ∞
y ·cos , D∞ = −D ∞
y ·sin
The K-factors presented in Table I are defined in the classical form

√
K I = lim 2 r · ( =l +r, = 0)
r →0
√
K II = lim 2 r · ( =l +r, = 0) (29)
r →0
√
K IV = lim 2 r · D ( =l +r, = 0)
r →0
DOI: 10.1002/nme
They can be calculated using the solution presented above by means of the explicit expressions
(29). To our knowledge, the K-factors for arbitrary values of are not found in the literature up to
now. We provide those for = 0, 0.1 , 0.2 , /4, 0.3 , 0.4 and /2. The cases = 0 and − /4
are given, for instance, in [23]. It is interesting to note that, on the one hand, ∞
yy influences K I
∞
and K II only. On the other hand, K IV depends only on D y for all values of . Another point
consists in the fact that the values presented in Table I do not depend on the piezoelectric material
parameters. In isotropic elasticity the full solution does not depend on the material properties
(Young’s modulus, Poisson ratio) because the characteristic roots (±i) are not depending on these
properties. In piezoelectric problems, the solution depends on the material properties, except at the
ligament lines of the crack tips. Hence, we come to the interesting conclusion that the K-factors
(29) do not depend on the material properties as well.
4. THE CLASSICAL X-FEM FORMULATION
The idea of partition of unity methods (PUM) [24–26] is to enrich an approximation field coming
from the standard finite elements with functions that are able to describe some specific behavior
of the exact solution field. The X-FEM is based on the local PUM [1, 27, 28]. In the sequel, we
explain first the application to cracks in linear isotropic materials. In this case, one uses elementary
displacement fields obtained from an asymptotic solution at the crack tip.
4.1. Crack description

As in the original article from Moës et al. [1], the crack is described by a pair of level set functions,
as shown in Figure 5. In the following, the local polar coordinates (r, ) at the crack tip, introduced
in Section 3.1, are used. They are computed from the values of the level set functions lsn and
lst in the vicinity of the crack. These level sets are interpolated on the same mesh used for the
finite element computation. This is of course not mandatory but may ease the fast computation of
r and .
lsn
lst
Figure 5. Cracked body with a level set representation of the crack.
DOI: 10.1002/nme
The dependence of r and is the following:

r = lsn2 +lst2
(30)
= atan2 (lst, lsn)
where atan2 is the four-quadrant inverse tangent function.
4.2. Enrichment functions

The classical enrichment function used for isotropic elastic materials are as follows. They are
derived from the asymptotic expansion for a semi-infinite crack in an infinite domain.
There are four basic enrichment functions (C1 = 1, 2, . . . , 4) needed to span the Westergaard–
Williams–Sneddon solution [27, 29, 30] at the crack tip:

√ √ √ √
gi (r, ) =
e
r sin , r cos , r sin sin , r cos sin (31)
2 2 2 2
Only elements immediately around the crack tip are usually enriched with these functions, but
in the sequel we will deal with geometrical and topological enrichment schemes as in [3]. For
the elements that are cut by the crack, the enrichment is chosen to be the modified Heaviside
discontinuous function so that the resulting displacement field contains a discontinuity at the
location of the crack [1]:

+1 if lst0
h(lst) = (32)
−1 if lst<0
In this function, lst is the level set representing the signed distance to the crack plane. In the
remaining portion of the domain, there is no enrichment, hence, the conventional finite element
basis is the only one existing. In the sequel, the conventional finite element used is a linear
three-node triangle. The resulting displacement field is

uh = Ni i + Ni g ej i j + Ni hi (33)
i∈R i∈T j=1,...,4 i∈H
Here R, T and H are, respectively, the regular set of nodes having the regular finite element
basis, the set of nodes, that need to be enriched for the tip functions, and the set of nodes for the
Heaviside enrichment.
5. NEW ENRICHMENT FUNCTIONS FOR CRACKS IN PIEZOELECTRIC MATERIALS
The asymptotic expansion for a crack in an infinite piezoelectric medium has been derived in
Section 3.1. From this expansion, we can determine a set of elementary functions that span the same
functional space, for any supposed orientation of the crack and any loading. For plane problems,
there are three independent eigenmodes (only two in the isotropic elasticity). These modes are
combined linearly for each spatial variable x or y. It should be noted that these functions depend
on the orientation of the crack with respect to the material axes , contrary to the isotropic case.
Altogether, one needs six functions to describe all the possible states at the crack tip, as opposed to
DOI: 10.1002/nme
the only four functions used to span the Westergaard–Williams–Sneddon solution in the isotropic
elastic case. These independent functions are extracted from the asymptotic singular solution (21),
(22) and (23), derived in Section 3.1:
p √ √ √ √ √ √
gi (r, ) = { r f 1 (), r f 2 (), r f 3 (), r f 4 (), r f 5 (), r f 6 ()} (34)
where
(re) (im)
f i () = ((, ), ai , ai )
⎧ (re) (im)
⎪
⎪ (re) (im) (, ai , ai ) (im)
⎪
⎨ (, a , a ) cos if ai >0
i i
2
= (35)
⎪
⎪ (re) (im)
(, ai , ai )
⎪
⎩ (re) (im) (im)
(, ai , ai ) sin if ai 0
2
(re) (im)
The complex numbers ai = ai +iai are the six roots of the characteristic equation following
from (15), i.e. (18), and is the orientation of the material axes with respect to the crack. The
angles and as well as the modified radius have the form
= − (36)

(re) (im)
(, ai , ai )= + int
2
⎛ ⎞
(re)
cos − int +ai sin − int
⎜ ⎟
− arctan ⎝ ⎠ (37)
(im)
|ai | sin − int

(re) (im)
(, ai ai )

1 4 (re) (im) (re) (re) (im)
= √ (ai )2 +(ai )2 +ai sin 2−[(ai )2 +(ai )2 −1] cos 2 (38)
2
In addition to this, the same enrichment function as in the case of isotropic elasticity is used to
represent the jump in the displacement and in the electric potential. This whole set is then combined
with regular FEM shape functions to build the approximation space for both the displacements u
and the potential :
p
uh = Ni i + Ni g j i j + Ni hi (39)
i∈R i∈T j=1,...,6 i∈H
p
h = Ni i + Ni g j i j + Ni hi (40)
i∈R i∈T j=1,...,6 i∈H
It is worth noting that the same enrichment functions are used for the displacements and the electric
potential. In Figure 6, we show the transformation between the pseudo-polar coordinates (, )
used in the six-fold enrichment, and the raw geometrical polar coordinates (r, ) at the crack tip
as a function of .
DOI: 10.1002/nme
Figure 6. Comparison between () (dashed line) and r () = 1 (dotted–dashed line) appearing,

respectively, in the six-fold piezoelectric enrichment and the usual four-fold enrichment, and
between () (dotted line) and (solid line).
6. CONVERGENCE STUDY IN THE ENERGY NORM
The convergence is studied by comparing exact solutions for a reference problem to the actual
results of a numerical model. The total (internal) energy of the system is given by

1
W= (ci jkl εi j εkl +i j E i E j ) dV (41)
2
The error in the (total) energy norm is then given by

1
errW = ci jkl (i j −iexj )(kl −ex
kl )+i j (E i − E i )(E j − E j ) dV
ex ex (42)
2
where εiexj and E ex
j are obtained from the exact displacement and electric fields u
ex and ex that
follow from the solution in Section 3.2. The error energy norm errW does not contain terms
involving the piezoelectric coupling ei jk , because of symmetries. This would not have been the
case if we had used the enthalpy H as reference, but convergence results shown in the sequel are
also valid for H .
The example problem is that of a Griffith–Irwin crack in an infinite medium, discussed in
Section 3.2. Boundary conditions are applied at infinity. Those are D ∞ = 10−3 and ∞ = 10
7
(case 1) and D ∞ = 10−3 and ∞ = 0 (case 2). The computations were made with an FEM model,
whose mesh is structured and has been gradually refined. Of course the computational domain
is finite, hence, we apply the exact solution for an infinite domain as boundary condition for the
finite computational domain. The length of the crack is 2 and the size of the domain is 8-by-8 (see
Figure 7). As stated before, all units are SI. The material axes have been rotated by an angle of 30◦
with respect to the plane of the crack to avoid trivial symmetry of the solution. Except for enriched
DOI: 10.1002/nme
Figure 7. Example of mesh used for the convergence study. Round dots are crack tip enrichment, square
dots are Heaviside enrichment. These are examples in the case of a geometrical enrichment, for h = 16 .
The length of the crack is 2 and the domain is 8-by-8. There are 625 nodes and 1152 elements.
degrees of freedom (dofs), the finite element discretization is linear for both displacements and
potential and we use three nodes linear triangles as the finite element basis.
We performed five sets of simulations with various enrichment strategies, for every case:
1. a topological enrichment done with six functions;
2. a geometrical enrichment done with six functions (radius of the enriched domain is 0.24);
3. a geometrical enrichment done with six functions and a bigger enriched domain (radius of
the enriched domain is 0.48);
4. a geometrical enrichment done with the four functions used in linear elasticity (radius of the
enriched domain is 0.24);
5. no enrichment at all.
Figures 8 (for case 1) and 9 (for case 2) show the relationship between the mesh density and the
error in the energy norm in a log–log scale. It should be stressed that any enrichment strategy allows
better results than conventional finite elements without enrichment. This fact has been explained
in earlier papers extensively, hence, we will not discuss it further. Thus, these results without
enrichment are for reference only. What we would like to show here is the relative performance
between (a) a topological enrichment and a geometrical enrichment and (b) an enrichment made
with the six functions that are specifically designed for piezoelectric materials and the simple use
of the four standard functions that are theoretically valid only for linear isotropic elasticity. For
case 2, with only electrical boundary conditions and a vanishing mechanical load, the idea is to
see whether the relative effectiveness of each enrichment depends on the nature of the applied
loads. Figure 9 clearly shows us that the answer is ‘no’. This is coherent with the fact that in the
current setting (an electrically impermeable crack) the electrical load induces a singularity both in
the stresses and in the electric displacements. As can be seen from Section 3, the mathematical
DOI: 10.1002/nme
100
topological enrichment (slope 0.5)
geometrical enrichment (slope 1.0)
big geometrical enrichment (slope 1.0)
geo std 4
no enrichment
error
10
1
1 10 100 1000
1/h
Figure 8. Convergence for various enrichment strategies, with a mixed

loading (D ∞ = 10−3 , ∞ = 10 ) (case 1).
7
topological enrichment (slope ~ 0.5)

geometrical enrichment (slope ~ 1.0)
big geometrical enrichment (slope ~ 1.0)
10 geo std 4
no enrichment
error
1 10 100 1000
1/h
Figure 9. Same as Figure 8, with a vanishing mechanical loading (D ∞ = 10−3 , ∞

= 0) (case 2).
form of the asymptotic solution at the crack tip represents an inherent coupling of electrical and
mechanical fields. It is characteristic for that class of material and does not depend on type or
magnitudes of external loadings. Therefore, no difference is obtained between both loading cases
with respect to convergence and enrichments. We would like to put emphasis on the fact that the
material axes are rotated arbitrarily. Thus, the exact solution cannot be represented by any of the
six enrichment functions taken alone. Only a linear combination of those may represent exactly
the solution, and the parameters of the linear combination are unknown a priori.
6.1. Relative performance: geometrical enrichment vs topological enrichment

Figures 8 and 9 show that the convergence rate for geometrical enrichment is roughly double
than that of the topological enrichment. This is what was expected (as shown in [3, 31]) when the
DOI: 10.1002/nme
enrichment function well represents the physics of the phenomenon. Here, the six enriched function
that were specifically developed for the piezoelectric crack are ‘doing the job’ and hence allow
to recover the classical convergence rate for finite elements with smooth problems (1.0 for linear
elements). The drawback is an increased computational effort due to the bad conditioning of the
resulting matrices. This has been discussed in [3] and solutions have been proposed to address this
issue. For very coarse meshes, the geometrical enrichment lead to less enriched elements than the
topological enrichment, therefore it is outperformed. In practice, when the geometrical enrichment
domain is smaller than the topological one, one should retain the latter as the minimal set of
enriched elements.
6.2. Relative performance: six-fold enrichment vs four-fold enrichment

Comparing results obtained using a specifically designed six-fold enrichment with the standard
four-fold enrichment of the isotropic elasticity in similar settings shows almost no difference. It
means that the classical four-fold enrichment spans a functional space that is large enough to
represent well the more complex crack tip behavior found in piezoelectric materials. Although
this may be surprising, it must be tempered by the fact that piezoelectric materials are not very
far from having a isotropic behavior in both the mechanical and electrical variables. Furthermore,
the coupling, which is clearly anisotropic, has a very weak magnitude. Thus, the additional dof
and shape of the six-fold enrichment may be of little advantage for piezoelectric materials. This is
not a bad news, since it means that reliable piezoelectric simulations can be made using the much
simpler four-fold enrichment.
7. COMPUTATION OF THE GENERALIZED SIFS
In order to calculate the SIFs for the stresses and the dielectric displacements, the common technique
via the J -integral can be generalized to cracks in piezoelectric materials. In the context of the
FEM, J is generally computed via an equivalent domain integral, see [32] for isotropic elasticity
in three dimensions. This procedure has been successfully adapted to the case of a piezoelectric
media [33].
Let qi = n i be a virtual crack extension field, normal to the crack front and parallel to the
crack surface. The scalar is equal to one on the crack front and vanishes at the boundary S of
the integration domain (noted V in the sequel).
Then, the J -integral represents the energy release rate for the virtual crack extension qi :

J =− (qi Q i j ), j dV − qi Q i j n j dS (43)
V SC
Now, this definition is extended to a piezoelectric material (small strains and linear material
law). Instead of the classical Eshelby tensor for isotropic elastic materials, Q i j is now replaced
by the piezoelectric energy momentum tensor [8]:
Q i j = 12 (kl εkl − Dk E k )i j −k j u k,i − D j ,i (44)
DOI: 10.1002/nme
We choose to cancel the integral over the crack surface SC (term in brackets), because a load-free
crack is considered with no surface charges or tractions.
One disadvantage of the J -integral technique consists in the problem that it gives solely the
energy release rate G, which cannot be trivially discriminated into the K -factors for the individual
crack opening modes. To overcome this drawback in the case of mixed mode crack problems, the
interaction integral approach has been suggested by Yau et al. [34] for two dimension and lately
by Gosz and Moran [35] for three-dimensional elastic crack problems. The interaction integral
was applied in the context of X-FEM for elastic crack analyzes as well [28].
Obviously, any good piezoelectric crack analysis needs SIFs. These are comprised from dielectric
K IV and mechanical field intensity factors K I , K II , K III . To extend the interaction integral approach
to piezoelectric materials, two loading states of the crack configuration will be considered:
• The actual loading (marked by ), for which the K -factors are wanted, which is analyzed by
finite elements.
• An auxiliary loading acting as reference state, denoted by aux .
Both loading states will now be superimposed to give a fictitious state, labeled by ∗ . As a
consequence, all field quantities are added, i.e. u ∗ = u +u aux , ε ∗ = ε +ε aux , ∗ = +aux , ∗ = +
aux , D ∗ = D + D aux and E ∗ = E + E aux . Also, the K -factors are superimposed: K N∗ = K N + K Naux ,
where N ∈ {I, II, III, IV}. Putting the composed field variables into the energy momentum tensor
Q i j of (44), the following relation comes out for the J -integral (43):
J ∗ = J + J aux + I (45)
Here, J and J aux contain all terms stemming from field variables of either the actual state or
the auxiliary state aux , respectively. The last term I groups all mixed expressions, due to the fact
that the J -integral is not a linear form with respect to the field variables. Therefore, it is called
the interaction integral.
Developing the above relation, one obtains
I = J ∗ − J − J aux

= (qi (k j u aux
k,i +k j u k,i + D j ,i + D j ,i )), j dV
aux aux aux
V

1
− qi (kl εkl
aux
+aux
kl εkl − Dk E k − Dk E k )i j
aux aux
dV (46)
V 2 ,j
The main symmetries of the material tensors allow us to write that kl εkl
aux = aux ε and D E aux =
kl kl k k
aux
Dk E k , thus, we have

I= (qi (k j u aux
k,i +k j u k,i + D j ,i + D j ,i )), j dV
aux aux aux
V

− (qi (kl εkl
aux
− Dk E kaux )i j ), j dV (47)
V
DOI: 10.1002/nme
Expanding this derivation yields

I = qi, j (k j u aux
k,i +k j u k,i + D j ,i + D j ,i ) dV
aux aux aux
V

+ qi (k j u aux
k,i +k j u k,i + D j ,i + D j ,i ), j dV
aux aux aux
V

− qi, j ((kl εkl
aux
− Dk E kaux )i j ) dV
V

− qi ((kl εkl
aux
− Dk E kaux )i j ), j dV (48)
V
In a two-dimensional setting, one can cancel the terms containing additional derivatives of the
auxiliary fields because they fulfill the equilibrium equations (4). This is not the case in a three-
dimensional setting, because complete solutions are usually not known for the auxiliary fields.
The constitutive relations are satisfied and can be used to show that the following relation holds:
kl εkl,i
aux
− Dk E k,i
aux
= εkl aux
kl,i − E k Dk,i
aux
(49)
Thus, we obtain an expression of the interaction integral, which depends on fewer auxiliary
variables. In the next expression, only u auxk,i , k j , ,i , D j , u k,i j , kl,i , ,i j , Dk,i , εkl , E k
aux aux aux aux aux aux aux aux aux
appear. The variables εkl,i

aux and E aux are eliminated.
k,i

I = qi, j (k j u aux
k,i +k j u k,i + D j ,i + D j ,i ) dV
aux aux aux
V

− qi, j (kl εkl
aux
− Dk E kaux )i j dV (50)
V
This final form of the interaction integral is valid for a load-free crack surface and in the absence
of volumetric forces and charges. It is identical with the derivations, made for two-dimensional
dynamic piezoelectric crack analyses in [19].
The relationship between the interaction integral and the energy release rate is the following in
a two-dimensional setting:
I = G iaux qi (51)
Since qi is normal to the crack front, we do not need to know other components of G iaux than the
one normal to the crack front, denoted by G aux . Thus, we have the following relation:
I = G aux = Y M N K M K Naux (52)
The matrix Y M N is the generalized Irwin matrix, depending on the piezoelectric material
properties.¶ All field intensity factors K M and K Naux are arranged in a vector with N ∈ {I, II, III, IV}.
¶ The energy release rate is still G = 12 Y M N K M K N as usual. The lack of factor 12 in Equation (52) comes from the
fact that G aux come from the superposition of two states: G(A, B) = G(A)+ G(B)+ G(A, B). Since G(X ) ∼ X 2 , I
or G(A, B) incorporates a factor 2 because (A + B)2 = A2 + B 2 +2AB.
DOI: 10.1002/nme
Now, a set of auxiliary fields is chosen, for which the K -factors are already known, from which
the vector K Naux is preset. For isotropic elastic materials, the following relations exists, which is
valid for any crack orientation:
⎡ ⎤
2(1−2 )
⎢ 0 0 ⎥ ⎧ aux ⎫
⎢ E ⎥ ⎪ KI ⎪
⎢ ⎥⎨ ⎬
⎢ 2(1− 2) ⎥ aux
Y M N K M K N = {K I K II K III } ⎢
aux
0 ⎥
0 ⎥ ⎪ II ⎪K (53)
⎢
⎢ E ⎥ ⎩ K aux ⎭
⎣ 2(1+) ⎦ III
0 0
E
For a piezoelectric material, however, the Irwin matrix is not so easy to determine. For a crack
perpendicular to the poling direction, this relation holds
⎡ 1 1 ⎤
0 0
⎢ cT e ⎥
⎢ ⎥ ⎧ K aux ⎫
⎢ ⎥⎪⎪ I ⎪
⎢ 0 1
0 ⎥ ⎪ aux ⎪ ⎪
⎢ 0 ⎥⎪⎨K ⎪ ⎬
⎢ cL ⎥ II
Y M N K M K Naux = {K I K II K III K IV } ⎢ ⎥ (54)
⎢ 1 ⎥⎪⎪
aux ⎪
K III ⎪
⎢ 0 0 ⎥ ⎪ ⎪
⎢ 0 ⎪
⎥ ⎩ aux ⎪ ⎭
⎢ cA ⎥ K IV
⎣ ⎦
1 1
0 0 −
e
The parameters e, cT , c L , c A and depend on material characteristics and may be found in the
literature [8] for common piezoelectric materials. For more general situations, the Irwin matrix
is not constant and depends on the orientation of the crack relative to the material axes and the
poling direction. In the following, we propose a way to compute the coefficients of the Irwin
matrix. For this purpose, we consider two auxiliary fields for which the K -factors are known.
These two auxiliary fields belong to any of the four eigenfunctions that respect the crack geometry
and material laws. Using these auxiliary fields, one can compute their interaction integral with the
procedure detailed above. The following relation holds:
G aux1,aux2 = Y M N K M K N = I aux1,aux2
aux1 aux2
(55)
By combining two auxiliary functions (each taken from the set of four available eigenfunctions)
as test functions in the interaction integral, one can determine numerically every component of the
Irwin matrix. For instance, by choosing
{K Iaux1 = 1 K IIaux1 = 0 K III

aux1
= 0 K IV
aux1
= 0} (56)
and
{K Iaux2 = 0 K IIaux2 = 1 K III

aux2
= 0 K IV
aux2
= 0} (57)
we get the following relation, which allows to determine Y12 and Y21 :
I aux1,aux2 = Y M N K M K N = Y12 = Y21

aux1 aux2
(58)
DOI: 10.1002/nme
The Irwin matrix is symmetric, hence, only 10 such integrals have to be computed in the three-
dimensional case to determine it completely for a given angle between the normal of the crack
and the poling direction. In a two-dimensional in-plane problem, K III has no meaning, hence, only
six integrals are needed. In order to identify completely the material law with regard to fracture
parameters, this has to be done for every . One can advantageously interpolate the -dependency
as a Fourier series, because it allows an efficient use in crack propagation algorithms without
requiring to compute the integrals each time the Irwin matrix is needed. It is expected that the
-dependency is well approximated even by few terms in the Fourier series.
Once this matrix is known, we can compute the interaction integrals I Naux between the actual
FEM-solution and the auxiliary fields selected for each mode N . Then, the following expression
gives the K -factors:
⎧ ⎫ ⎧ aux ⎫
⎪ KI ⎪ ⎪ I
⎨ ⎬ ⎨ I ⎪ ⎬
−1
K II = Y M N IIIaux (59)
⎪
⎩ ⎪
⎭ ⎪
⎩ aux ⎪ ⎭
K IV IIV
In the sequel, the Irwin matrix is computed numerically with the same quadrature rule as the
interaction integrals. Afterwards, it was used to extract the K-factors for a given orientation of the
material axes.
8. CONVERGENCE OF THE SIFS
We consider two cases: (i) a semi-infinite crack in an infinite piezoelectric medium and (ii) a
Griffith–Irwin crack in an infinite piezoelectric medium. The computational domain is actually
finite (a 4-by-4 square) but loaded along its boundaries with those fields, corresponding to the
exact solution that is sought numerically. In both cases, the numerical values of the K-factors are
determined using the equivalent integrals for a given contour (the radius of this domain is set
to 0.8). The radius of the enriched domain in the case of a geometrical enrichment is 0.24. In
Figures 10–14, a comparison is made between topological and geometrical piezoelectric enrichment
(six-fold) as well as with the classical four-fold enrichment. The results without enrichment (regular
finite elements) are also shown. We should stress that, when looking at Figures 10–14, results for
1/ h<50 are not very interesting because the mesh is very coarse with regard to the characteristic
length of the problem and the area where the domain integrals are computed. However, we kept
these results for the sake of completeness.
8.1. Semi-infinite crack in an infinite domain

Here, we consider three independent problems, for which the exact values of the K-factors are,
respectively, (K I = 1, K II = 0, K IV = 0), (K I = 0, K II = 1, K IV = 0) and (K I = 0, K II = 0, K IV = 1).
These cases are obtained by loading the boundaries of the domain with the relevant eigenfunctions
(auxiliary fields used to determine the K-factors) that are detailed in Section 3.1. As for the energy
error, the performance of the four-fold enrichment is very close to that of the six-fold one. The
dependence over the enrichment domain is clear, with geometrical enrichment being substantially
superior to topological enrichment. It should be stressed that, from an engineering point of view,
DOI: 10.1002/nme
1
topological enrichment
geometrical enrichment
topo std 4
geo std 4
no enrichment
0.1
(Kexact-Kfem)/Kexact
0.01
0.001
1e-04
1 10 100 1000
1/h
Figure 10. Convergence of K I for a semi-infinite crack.
1
topo std 4
geo std 4
no enrichment
0.1
(Kexact-Kfem)/Kexact)
0.01
0.001
1e-04
1 10 100 1000
1/h
Figure 11. Convergence of K II for a semi-infinite crack.
these results are in very good agreement with theory (less than 0.3% in relative error). The results
for K I and K II are similar, but those for K IV look somewhat wilder.
8.2. Griffith crack in an infinite domain

In this case, we consider a problem with boundary conditions given at the infinite. The boundary
conditions are ∞ ∞ ∞ ∞ ∞ −3
= 10 , = 0, = 0, D = 0 and D = 10 . As before, the computational
7
domain is finite; therefore, the exact solution is applied on the actual boundary of the domain. The
exact solution is detailed in Section 3.2. The exact values of the SIFs in this setting are determined
in Section 3.3 (see Table I). These are K I = 1.772454×107 and K IV = 1.772454×10−3 . The results
depicted in Figures 13 and 14 are similar to the previous case, with geometrical enrichment strategy
being substantially superior to topological enrichment.
DOI: 10.1002/nme
0.1
topo std 4
geo std 4
no enrichment
0.01
0.001
1e-04
1e-05
1 10 100 1000
1/h
Figure 12. Convergence of K IV for a semi-infinite crack.
1
topo std 4
geo std 4
no enrichment
0.1
0.01
0.001
1e-04
1 10 100 1000
1/h
Figure 13. Convergence of K I for a Griffith–Irwin crack.
9. CONCLUSIONS
We have devised an enrichment scheme based on the asymptotic expansion around a crack tip
in piezoelectric materials. This enrichment scheme uses the six eigenfunctions of the asymptotic
expansion. As such, it is optimal with regard to the problem to be solved (it contains the exact
solution for a semi-infinite crack in an infinite body). However, a comparison made with the
classical four-fold enrichment used for isotropic elasticity is interesting. It shows that, with an
accuracy higher than usually expected for engineering problems, the four-fold enrichment is almost
as efficient, concerning accuracy both in energy and in the SIFs. The good news is that it is simpler
to implement and involves less computational overhead, because it adds only four degrees of
freedom (dofs) per regular dof, instead of six. For piezoelectric problems, it is therefore advisable to
use the regular enrichment functions stemming from the isotropic elasticity. It should be mentioned
DOI: 10.1002/nme
0.1
topo std 4
geo std 4
no enrichment
0.01
0.001
1e-04
1 10 100 1000
1/h
Figure 14. Convergence of K IV for a Griffith–Irwin crack.
that similar results have been obtained independently [36] with alternative enrichment functions
for cracks in confined plasticity problems. Along with the appropriate crack propagation law for
piezoelectric materials, numerical crack propagation using level sets update algorithms will be
implemented and this will be the topic of a forthcoming publication.
APPENDIX A
A.1. Definition of the two-dimensional material characteristics

The two-dimensional material parameters introduced in Section 3.1 follow from the general three-
dimensional ones given in [8] by means of the so-called Voigt-notation-form for the properties ci j ,
ei j and i j .
1 1 −1 11
a11 = (c33 33 +e33
2
), a22 = (c11 33 +e31
2
), a12 = (c13 33 +e31 e33 ), a33 =
A A A B
1 1 e15
b21 = (c33 e31 −c13 e33 ), b22 = (c11 e33 −c13 e31 ), b13 =
A A B
c44 1
11 = , 22 = (c11 c33 −c13
2
)
B A
A = c11 (c33 33 +e33
2
)−c13
2
33 −e31 (2c13 e33 −c33 e31 ), B = c44 11 +e15
2
A.2. Griffith–Irwin crack solution

In order to explain the general solution procedure used in Section 3.2, we will consider the
Griffith–Irwin crack problem with polarization perpendicular to the crack. It will be shown that
DOI: 10.1002/nme
the algebraic manipulations are not very complicated, but one has to take care for the rigid body
motion and reference potential.
In the given case, the material axes 0x y coincide with 0 and is set to zero. The roots ak of
the characteristic equation, which is the same as (18), get the special form [15]:
(re) (im) (re) (im) (im) (re) (im) (im)
a1 = a1 +ia1 , a2 = −a1 +a1 , a3 = ia3 , (a1 > 0, a1 > 0, a3 > 0) (A1)
The other three roots represent the corresponding complex conjugate values. The constants k in
(26) are defined through
(b21 +b13 )ak2 +b22

k = − , k = 1, 2, 3 (A2)
11 ak2 +22
Applying the material equations (9) and (10), the deformation and electric field components can be
determined and after integration of these functions, one gets the representations of the displacements
u , u and the electric potential through the complex potentials k (k ):

3
u = 2RE (a11 ak +a12 −b12 k )k (k ) −∗ y +u ∗
2
k=1

3 1
u = 2RE (a12 ak +a22 −b22 k )k (k ) +∗ x +v∗
2
(A3)
k=1 a k

3
= −2RE (b13 +11 k )ak k (k ) +∗
k=1
The analytic functions k (k ) represent the antiderivatives of k (k ), that is,
dk (k )
k (k ) = , k = 1, 2, 3 (A4)
dk
∗ , u ∗ , v∗ and ∗ characterize the rigid body motion and a reference potential, which must be set
to zero (for instance). In order to solve the boundary value problem, one has to fulfill the boundary
conditions at the crack surfaces and at infinity.
Vanishing -stresses at the crack surfaces (||<l, = 0) result from (26) in
2x (re) (re) (im) (im) (im) (re) (re) (im)

((A1 + A2 )a1 +(A1 − A2 )a1 + A3 a3 )
2
l 2 −
(re) (re) (re) (im) (im) (im) (im) (im)
+2((B1 − B2 )a1 −(B1 + B2 )a1 − B3 a3 ) = 0 (A5)
(A5) is equivalent to the conditions

(re) (re) (im) (im) (im) (re) (re) (im)
(A1 + A2 )a1 +(A1 − A2 )a1 + A3 a3 =0
(A6)
(re) (re) (re) (im) (im) (im) (im) (im)
(B1 − B2 )a1 −(B1 + B2 )a1 − B3 a3 =0
DOI: 10.1002/nme
The corresponding equations for end up in

(im) (im) (im) (re) (re) (re)
A1 + A2 + A3 = 0, B1 + B2 + B3 =0 (A7)
and the equations for vanishing D at the crack surfaces have the form
(re) (re) (im) (im) (im)
yi1 (A1 − A2 )+ yr 1 (A1 + A2 )+ yr 3 A3 =0
(A8)
(im) (im) (re) (re) (re)
yi1 (B1 − B2 )− yr 1 (B1 + B2 )− yr 3 B3 =0
The definitions of the coefficients yi1 , yr 1 and yr 3 are given at the end of this section. Solving
(re) (re) (im) (re) (re) (re)
(A6), (A7) and (A8) with respect to A1 , A2 , A3 , B1 , B2 , B3 one gets
(re) 1 (im) (im) (im) (re)

A1 = (im)
(a11 A1 +a12 A2 −a3 A3 )
2a1
(re) −1 (im) (im) (im) (re)
A2 = (im)
(a12 A1 +a11 A2 +a3 A3 )
2a1
(re) 1 (im) (im) (im) (im)
B1 = (re)
(b11 B1 +b12 B2 +a3 B3 ) (A9)
2a1
(re) −1 (im) (im) (im) (im)
B2 = (re)
(b12 B1 +b11 B2 +a3 B3 )
2a1
(im) (im) (im) (re) yi1 (im) (im)
A3 = −(A1 + A2 ), B3 = (B − B1 )
yr 1 − yr 3 2
The parameters a11 , a12 , b11 and b11 depend on the material properties. The specific values
are given at the end. In accordance with the behavior at infinity shown in Figure 3, the following
relations are assumed:
∞
= lim = lim = lim = lim
= 0 (A10)
y→+∞ y→−∞ x→+∞ x→−∞
∞
= lim = lim = lim = lim
= 0 (A11)
y→+∞ y→−∞ x→+∞ x→−∞
D ∞ = lim D = lim D = lim D = lim D

= 0 (A12)
y→+∞ y→−∞ x→+∞ x→−∞
∞
= lim = lim = lim = lim (= 0) (A13)
y→+∞ y→−∞ x→+∞ x→−∞
D∞ = lim D = lim D = lim D = lim D (= 0) (A14)

y→+∞ y→−∞ x→+∞ x→−∞
Thus, the analytic functions k (k ) defined by means of (27) are prolongated continuously to
infinity with given constant values for the stresses and dielectric displacements. Connecting these
values with the coefficients of k (k ) in (27), it is necessary to assure the continuity of the applied
radical functions outside the crack and within the whole – -plane including infinity. Taking into
DOI: 10.1002/nme
account this circumstance, one gets the relations for ∞ ∞ ∞ ∞ ∞

, , D , and D by means of (26),
(27) and (A9)–(A14) in the form
(re) (im) (im)
2a1 (im) (im) a1 −a3 (re)
∞
= (im)
(A1 + A2 )+2 (im)
A3 (A15)
a1 a1
(re)

(im) (im) a (im) (im)
∞
= 2 a1 −a3 + 1 (yr 1 − yr 3 ) ·(A1 + A2 ) (A16)
yi1
(re)
(im)

a1 (im) (im) a3 (re)
D ∞ = 2 yi1 + y
(im) r 1
·(A1 − A2 )+2 y − yr 3
(im) r 1
A3 (A17)
a1 a1
(im) (im) (re) (im) (im)
∞
= a (A1 − A2 )+a3r A3 +b (B1 − B2 ) (A18)
D∞ = a (A(im) (im) (im)

1 + A2 )+b (B1 + B2(im) )+b3 B3(im) (A19)
The definition of the coefficients a , a3r , b , a , b and b3 is given at the end. The equation
(im) (im)
system (A15)–(A19) represents five equations for the determination of six unknowns: A1 , A2 ,
(re) (im) (im) (im)
A3 , B1 , B2 and B3 . Thus, the problem is under-determined so far. The meaning of this
fact will be explained in the following.
Now, we will restrict our considerations to symmetrical loading regarding the crack tip, i.e.
∞ ∞
= 0. Equating with zero, one finds from (A16) and from the fifth equation of (A9)
(im) (im) (im)
A2 = −A1 , A3 =0 (A20)
(im) (re)
With the aid of (A15), (A17) and (A20), A1 and A3 are defined through a linear combination
of ∞ ∞
and D y :
(im) (im) (im) (im)
(im) yr 3 a1 − yr 1 a3 a1 −a3
A1 = ∞
+ D ∞ (A21)
2as 2as
(re) (im) (re)
yr 1 a1 + yi1 a1
(re) a1
A3 = ∞
+ D∞ (A22)
as as
The value as is given at the end. Using the first relation of (A20) in connection with (A18) and
(im)
(A19) for ∞ ∞
= D = 0, which is valid in our case, B1 and B2(im) can be calculated by means of
(im) (im)
(im) (im) (re) x r 1 a3 − xi3 a1 (im)
B1 = −c A1 +cr A3 + (re) (im)
B3 (A23)
2(xi1 a1 − xr 1 a1 )
(im) (im)
(im) (im) (re) x r 1 a3 − xi3 a1 (im)
B2 = c A1 −cr A3 + (re) (im)
B3 (A24)
2(xi1 a1 − xr 1 a1 )
The constants c , cr , xi1 , xr 1 and xi3 are given at the end. Based on (A20)–(A24) and (A9), all
(im)
unknown coefficients (except B3 ) of the analytic functions in (27) are known linear combinations
DOI: 10.1002/nme
(im)
of the given parameters ∞ ∞
, D and the undefined constant B3 . Thus, the solution for this
(im) (im)
crack problem can be obtained, if B3 is determined. In order to identify B3 , it is necessary
to consider the corresponding displacements u , u and the electric potential in (A3) at infinity,
∞ ∗
that is, u ∞ ∞
, u and . Taking into account vanishing rigid body motion ( = u ∗ = v∗ = 0) and
∞
a vanishing reference potential (∗ = 0), u ∞ ∞
, u and have to be the form (based on (5), (6),
(9), (10) and (A10)–(A14)):
u∞
(, ) = a +a , u∞
(, ) = b +b , ∞ (, ) = c +c (A25)
The constants a , a , b , b , c and c are connected with ∞ ∞ ∞

, and D by means of (9) and
∞ ∞
(10) ( = D = 0!):
a = b21 D ∞ ∞
y +a12 (A26)
a +b = a33 ∞
(A27)
b = b22 D ∞ ∞
y +a22 (A28)
c = b13 ∞
(A29)
c = b22 ∞ ∞
−22 D y (A30)
On the other hand, a , a , b , b , c and c are related to the coefficients of the analytic functions
(k ) with respect to (A3)

2 3
a = lim RE (a11 ak +a12 −b12 k )k (k )
2
(A31)
x→∞ x k=1

2 3
a = lim RE (a11 ak +a12 −b12 k )k (k )
2
(A32)
y→∞ y k=1

2 3 1
b = lim RE (a12 ak +a22 −b22 k )k (k )
2
(A33)
x→∞ x k=1 ak

2 3 1
b = lim RE (a12 ak +a22 −b22 k )k (k )
2
(A34)
y→∞ y k=1 ak

−2
3
c = lim RE (b13 +11 k )ak k (k ) (A35)
x→∞ x k=1

−2
3
c = lim RE (b13 +11 k )ak k (k ) (A36)
y→∞ y k=1
For the considered case of ∞ = 0, Equations (A31) and (A34)–(A36) lead to the automatic
fulfillment of Equations (A26) and (A28)–(A30), respectively (with c = 0), based on (A20)–(A24)
DOI: 10.1002/nme
together with (A9). In the same manner, (A32) and (A33) define a and b by means of
(im) (im) (im)
2a11 11 a3 a f B3 2(a22 22 +(b22 )2 )a f B3
a = (im) 2
, b = (im) (re) (im) 2 2 (im)
11 (a3 ) −22 a3 ((a1 )2 +(a1 ) ) (22 −11 (a3 )2 ) (A37)
(im) (im) 2 (re) (im) (im) 2 (re)
a f = ((a1 +a3 ) +(a1 )2 )((a1 −a3 ) +(a1 )2 )>0
(im)
There is no doubt that a and b vanish only in the case of B3 = 0. The sum a +b follows
from (A37) with
(im)
2a f b f B3
a +b = (im) (re) (im) 2 2 (im)
≡0
a3 ((a1 )2 +(a1 ) ) (11 (a3 )2 −22 ) (A38)
(im) 2 (re) (im) 2 2
b f = a11 11 (a3 ) ((a1 )2 +(a1 ) ) −a22 22 −(b22 )2 ≡ 0
(re) (im) (im)
b f vanishes identically due to the fact that a1 , a1 and a3 fulfill the characteristic equation
(18). Thus, Equation (A27) is satisfied identically as well. Because of (A27) and ∞ = 0, non-
(im)
vanishing a and b are only imaginable in the case of a = −b
= 0 (B3
= 0), which represents
(im)
a rigid body rotation based on ∗ = b in (A3). In fact, one must set B3 equal zero to exclude
the rigid body rotation:
(im)
B3 =0 (A39)
Thus, all coefficients in (27) are known in form of a linear combination of ∞ ∞
and D based
on (A39), (A20)–(A24), and (A9), and the entire solution of the crack problem (∞ ∞

= 0, D y
= 0,
∞
= 0) is constructed and the problem regarding the under-definiteness in (A15)–(A19) is resolved.
Because of a = b = c = 0, this solution is completely symmetric with respect to the co-ordinate
(im)
axes. We have shown that B3 has the meaning of the rigid body rotation. In fact, the rigid body
rotation cannot be assigned arbitrarily to any of the constants Bi(im) or Bi(re) as one believes in
(im)
[15, 20, 21], but only to B3 .
Definitions of the parameters introduced above are as follows:
yr 1 = RE(1 ), yi1 = IM(1 ), yr 3 = RE(3 ), xr 1 = RE(1 a1 )
xi1 = IM(1 a1 ), xi3 = IM(3 a3 )
(In the considered case i follow from (A2)!)
−1 (re) (im) 1 (re) (im) 1 (re) (im)
a11 = (yi1 a1 +a1 D), a12 = (yi1 a1 −a1 D), b11 = (yi1 a1 +a1 D)
yi1 yi1 D
(re) (im)
1 (im) (re) a1 a3
b12 = (a D−yi1 a1 ), a = −2 (im)
A, a3r = 2 C
D 1 a1 a1(im)
(im)

−2 (im) (re) xr 1 a1
b = (yi1 B+2a1 a1 D), a = −2 xi1 −D , b = 2 xr 1 (re)
− xi1
D yi1 a1
DOI: 10.1002/nme
(im)
(re) (im)
a3 a1 a3
b3 = 2 xr 1 (re)
− xi3 , c = A·D, cr = C·D
a1 E 2E
(im) (im) (re) (re) (im) 2 (im) 2 (re) (im) 2

as = 2(yi1 (a1 −a3 )+a1 D), A = (a1 )2 +(a1 ) , B = (a1 ) −(a1 )2 −(a3 )
C = (a1(im) )2 −a1(im) a3(im) −(a1(re) )2 , D = yr 1 − yr 3 , E = a1(im) (yi1 B+2Da1(im) a1(re) )
A.3. Material constants appearing in the computations using the Voigt notation
The following constants are from Reference [37]:
Elastic constants [1010 N/m2 ]:
c11 = 14.02, c12 = 7.892

c33 = 11.58, c13 = 7.565
c44 = 2.527
Piezoelectric constants [C/m2 ]:
15 = 12.0
31 = −5.268
33 = 15.44
Dielectric constants [10−9 C/(Vm)]:
11 = 6.368
33 = 5.523
REFERENCES
1. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International
Journal for Numerical Methods in Engineering 1999; 46:131–150.
2. Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi
formulations. Journal of Computational Physics 1988; 79:12–49.
3. Bechet E, Minnebo H, Moës N, Burgardt B. Convergence and conditioning issues with x-fem in fracture
mechanics. International Journal for Numerical Methods in Engineering 2005; 64:1033–1056.
4. Ji H, Chopp D, Dolbow JE. A hybrid extended finite element/levelset method for modeling phase transformations.
International Journal for Numerical Methods in Engineering 2002; 54:1209–1233.
5. Song JH, Areias PMA, Belytschko T. A method for dynamic crack propagation and shear band propagation with
phantom nodes. International Journal for Numerical Methods in Engineering 2006; 67:868–893.
6. Areias P, Song JH, Belytschko T. Analysis of fracture in thin shells by overlapping paired elements. Computer
Methods in Applied Mechanics and Engineering 2006; 195:5343–5360.
7. Waisman H, Belytschko T. Parametric enrichment adaptivity by the extended finite element method. International
Journal for Numerical Methods in Engineering 2008; 73:1671–1692.
8. Pak YE. Linear electro-elastic fracture mechanics of piezoelectric materials. International Journal of Fracture
1992; 54(1):79–100.
DOI: 10.1002/nme
9. Qin QH. Fracture Mechanics of Piezoelectric Materials. WIT Press: Southampton, Boston, 2001.
10. Scherzer M, Kuna M. Combined analytical and numerical solution of 2d interface corner configurations between
dissimilar piezoelectric materials. International Journal of Fracture 2004; 127(1):61–99.
11. Zhang TY, Zhao M, Tong P. Fracture of piezoelectric ceramics. Advances in Applied Mechanics 2002; 38:147–289.
12. Kuna M. Finite element analyses of cracks in piezoelectric structures: a survey. Archives in Applied Mechanics
2006; 76:725–745.
13. Koy YL, Chiu KW, Marshall IH, Rajic N, Galea SC. Detection of disbonding in a repair patch by means of an
array of lead zirconate titanate and polyvinylidene fluoride sensors and actuators. Smart Materials and Structures
2001; 10:946–962.
14. Suo Z, Kuo CM, Barnett DM, Willis JR. Fracture mechanics for piezoelectric ceramics. Journal of the Mechanics
and Physics of Solids 1992; 40(4):739–765.
15. Sosa H. Plane problems in piezoelectric media with defects. International Journal of Solids and Structures 1991;
28(4):491–505.
16. Mußchelichwili NI. Einige Grundaufgaben zur mathematischen Elastizitaetstheorie. VEB Fachbuchverlag: Leipzig,
1971.
17. Savin GN. Distribution of Stresses at Holes (in Russian). Naukova Dumka: Kiew, 1968.
18. Landis CM. Energetically consistent boundary conditions for electromechanical fracture. International Journal
of Solids and Structures 2004; 41:6291–6315.
19. Enderlein M, Ricoeur A, Kuna M. Finite element techniques for dynamic crack analysis in piezoelectrics.
International Journal of Fracture 2005; 134:191–208.
20. Xu XL, Rajapakse RKND. Analytical solution for an arbitrarily oriented void/crack and fracture of piezoceramics.
Acta Materialia 1999; 47(6):1735–1747.
21. Xu XL, Rajapakse RKND. On plane cracks in piezoelectric solids. International Journal of Solids and Structures
2001; 38:7643–7658.
22. Lekhnitskii SG. Theory of Elasticity of an Anisotropic Body. Mir: Moscow, 1981.
23. Pan E. A BEM analysis of fracture mechanics in 2d anisotropic piezoelectric solids. Engineering Analysis with
Boundary Elements 1999; 23:67–76.
24. Babuska I, Melenk JM. The partition of unity method. International Journal for Numerical Methods in Engineering
1997; 40:727–758.
25. Strouboulis T, Babuska I, Copps K. The design and analysis of the generalized finite element method. Computer
Methods in Applied Mechanics and Engineering 2000; 181:43–69.
26. Strouboulis T, Copps K, Babuska I. The generalized finite element method: an example of its implementation and
illustration of its performance. International Journal for Numerical Methods in Engineering 2000; 47:1401–1417.
27. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal
for Numerical Methods in Engineering 1999; 44:601–620.
28. Moës N, Gravouil A, Belytschko T. Non-planar 3d crack growth by the extended finite element and level sets—
part i: Mechanical model. International Journal for Numerical Methods in Engineering 2002; 53:2549–2568.
29. Sih GC, Liebowitz H. Fracture, an Advanced Treatise, vol. 2. Academic Press: London, U.K., 1968.
30. Bui HD. Mécanique de la rupture fragile. Masson: Paris, 1978.
31. Laborde P, Pommier J, Renard Y, Salaün M. High order extended finite element method for cracked domains.
International Journal for Numerical Methods in Engineering 2005; 64:354–381.
32. Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed
body. International Journal of Fracture 1986; 30:79–102.
33. Abendroth M, Groh U, Kuna U, Ricoeur A. Finite element-computation of the electromechanical J-Integral for
2-D and 3-D crack analysis. International Journal of Fracture 2002; 114:359–378.
34. Yau J, Wang S, Corton H. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity.
Journal of Applied Mechanics 1980; 47:335–341.
35. Gosz M, Moran B. An interaction energy integral method for computation of mixed-mode stress intensity factors
along non-planar crack fronts in three dimensions. Engineering Fracture Mechanics 2002; 69:299–319.
36. Elguedj T, Gravouil A, Combescure A. Appropriate extended functions for x-fem simulation of plastic fracture
mechanics. Computer Methods in Applied Mechanics and Engineering 2005; 195(7–8):501–515.
37. Wang X, Shen YP. Exact solution for mixed boundary value problems at anisotropic piezoelectric bimaterial
interface and unification of various interface defects. International Journal of Solids and Structures 2002;
39:1591–1619.
DOI: 10.1002/nme

Application of The X-FEM To The Fracture of Piezoelectric Materials

Uploaded by

Copyright:

Available Formats

Application of The X-FEM To The Fracture of Piezoelectric Materials

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Application of The X-FEM To The Fracture of Piezoelectric Materials

Uploaded by

Copyright:

Available Formats

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING

Int. J. Numer. Meth. Engng 2009; 77:1535–1565

Application of the X-FEM to the fracture of

E. Béchet1, ∗, † , M. Scherzer2, ‡ and M. Kuna2

09596 Freiberg, Germany

Received 30 August 2007; Revised 8 July 2008; Accepted 22 July 2008

Copyright q 2008 John Wiley & Sons, Ltd.

3.1. Asymptotic expansion for a semi-infinite crack

*i (z i ) (b21 +b13 )ai2 +b22 *2 Ui (z i )

Figure 2. Definition of the material axes around the crack tip.

The complex variables di (

3.2. Analytical solutions for a Griffith–Irwin crack in an infinite medium

Figure 3. Griffith–Irwin crack inside an infinite piezoelectric solid.

Ak  (re) (im) (re) (im)

Figure 4. Griffith crack inclined to the load and poling directions.

Table I. K-factors of the inclined crack.

The values coincide for both tips. l = 1 m.

values at infinity can be calculated by means of

The K-factors presented in Table I are defined in the classical form

4. THE CLASSICAL X-FEM FORMULATION

4.1. Crack description

Figure 5. Cracked body with a level set representation of the crack.

The dependence of r and is the following:

where atan2 is the four-quadrant inverse tangent function.

4.2. Enrichment functions

5. NEW ENRICHMENT FUNCTIONS FOR CRACKS IN PIEZOELECTRIC MATERIALS

Figure 6. Comparison between ( ) (dashed line) and r ( ) = 1 (dotted–dashed line) appearing,

6. CONVERGENCE STUDY IN THE ENERGY NORM

Figure 8. Convergence for various enrichment strategies, with a mixed

topological enrichment (slope ~ 0.5)

Figure 9. Same as Figure 8, with a vanishing mechanical loading (D ∞ = 10−3 , ∞

6.1. Relative performance: geometrical enrichment vs topological enrichment

6.2. Relative performance: six-fold enrichment vs four-fold enrichment

7. COMPUTATION OF THE GENERALIZED SIFS

Q i j = 12 (kl εkl − Dk E k )i j −k j u k,i − D j ,i (44)

Expanding this derivation yields

appear. The variables εkl,i

{K Iaux1 = 1 K IIaux1 = 0 K III

{K Iaux2 = 0 K IIaux2 = 1 K III

I aux1,aux2 = Y M N K M K N = Y12 = Y21

8. CONVERGENCE OF THE SIFS

8.1. Semi-infinite crack in an infinite domain

Figure 10. Convergence of K I for a semi-infinite crack.

Figure 11. Convergence of K II for a semi-infinite crack.

8.2. Griffith crack in an infinite domain

Figure 12. Convergence of K IV for a semi-infinite crack.

Figure 13. Convergence of K I for a Griffith–Irwin crack.

Figure 14. Convergence of K IV for a Griffith–Irwin crack.

A.1. Definition of the two-dimensional material characteristics

A.2. Griffith–Irwin crack solution

(b21 +b13 )ak2 +b22

2x (re) (re) (im) (im) (im) (re) (re) (im)

(A5) is equivalent to the conditions

The corresponding equations for end up in

(re) 1 (im) (im) (im) (re)

D ∞ = lim D = lim D = lim D = lim D

D∞ = lim D = lim D = lim D = lim D (= 0) (A14)

account this circumstance, one gets the relations for ∞ ∞ ∞ ∞ ∞

D∞ = a (A(im) (im) (im)

i (z i ) (b21 +b13 )ai2 +b22 2 Ui (z i )

Ak (re) (im) (re) (im)

Figure 6. Comparison between () (dashed line) and r () = 1 (dotted–dashed line) appearing,

Q i j = 12 (kl εkl − Dk E k )i j −k j u k,i − D j ,i (44)

D∞ = lim D = lim D = lim D = lim D (= 0) (A14)

D∞ = a (A(im) (im) (im)

The constants a , a , b , b , c and c are connected with ∞ ∞ ∞