Three-dimensional finite elements with embedded strong discontinuities to model failure in electromechanical coupled materials

Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Three-dimensional finite elements with embedded strong discontinuities to model failure in electromechanical coupled materials Christian Linder ⇑, Xiaoxuan Zhang Department of Civil and Environmental Engineering, Stanford University, Stanford, CA 94305, USA

a r t i c l e

i n f o

Article history: Received 7 November 2013 Received in revised form 15 January 2014 Accepted 25 January 2014 Available online 6 February 2014 Keywords: Failure in 3D solids Strong discontinuities Electromechanical coupling Piezoelectric ceramics Configurational forces Marching cubes algorithm

a b s t r a c t This paper presents new finite elements with embedded strong discontinuities to model failure in three dimensional electromechanical coupled materials. Following the strong discontinuity approach for plane electromechanical problems, the coupled boundary value problem is decomposed into a continuous global part and into a discontinuous local part where strong discontinuities in the displacement field and electric potential are introduced. Those are incorporated into general three-dimensional brick finite elements through nine mechanical separation modes and three new electrical separation modes. All the local enhanced parameters related to those modes can be statically condensed out on the element level, yielding a computationally efficient framework to model failure in electromechanical coupled materials. Impermeable electric boundary conditions are assumed along the strong discontinuities. Their initiation and orientation is detected through a configurational force driven failure criterion. A marching cubes based crack propagation concept is used to obtain smooth failure surfaces in the three dimensional problems of interest. Several representative numerical simulations are included and compared with experimental results of failure in piezoelectric ceramics to outline the performance of the new finite elements. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Electromechanical coupled materials play an important role as actuators, transducers or sensors in smart systems. Extensive research on their electromechanical properties has been performed over the years resulting in general frameworks and reviews of their phenomenological and micromechanical motivated constitutive response [1–5]. Their high ultimate strength and extremely low fracture energy makes those materials prone to develop cracks, affecting their mechanical properties like durability or strength. Numerous theoretical [6–11] and experimental [12–14] studies on the failure characteristics of, in particular, piezoelectric ceramics are available in the literature. In this work, a computational framework is proposed to model the onset and propagation of failure in three-dimensional electromechanical coupled materials. A large number of numerical techniques have been developed to model purely mechanical solids at failure such as cohesive finite element formulations [15,16], adaptive remeshing techniques [17–21], phase-field models [22–30], the extended

⇑ Corresponding author. Tel.: +1 6507232918. E-mail address: [email protected] (C. Linder). http://dx.doi.org/10.1016/j.cma.2014.01.021 0045-7825/Ó 2014 Elsevier B.V. All rights reserved.

144

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

finite element method [31–34], or the embedded finite element method [35–40]. Their extensions to model fracture in electromechanical coupled problems though is only in its infancy [41–50] and serves as the main motivation for this work. In particular, we consider the embedded finite element formulation, which introduces enhanced local kinematic parameters to describe failure zones such as cracks or shear bands through strong discontinuities in the displacement field for purely mechanical problems. This approach, originally developed in [35] is extended in [51,52] to three dimensions, in [53,54] to porous media, in [55–57] to structural problems such as beams and plates, in [58–61] to account for higher order kinematic approximations of the displacement field within the two dimensional (2D) and three dimensional (3D) continuum setting, in [62,63] to dynamic fracture, in [64] to multiple levels, and in [43–45] to 2D electromechanical coupled materials. Those latter works serve as the starting point for the extension to model failure in 3D electromechanical coupled problems in this work. Following [43], the overall electromechanical boundary value problem (BVP) is decomposed in this work into a continuous global problem, treated as the standard electromechanical BVP, and a discontinuous local problem, where jumps in the displacement field and electric potential are incorporated along the strong discontinuities. The decomposition of the method into a global and a local part results in an efficient numerical implementation by enhancing the standard displacement based, mixed, enhanced three-dimensional brick finite elements with local parameters carrying the information of jumps in the displacement field and electric potential. Those enhanced parameters, which all can be statically condensed out on the element level therefore leading to a computational efficient formulation, are associated with three newly introduced electrical separation modes and the nine mechanical separation modes proposed in [60] with a slight modification of the rotation mode around the normal of the failure surface by an in-plane shear mode in [61]. The smoothness and continuity of the final failure surfaces is ensured in this work through the incorporation of the marching cubes based crack propagation concept proposed in [61] for purely mechanical problems. One particular challenge for the modeling of failure in electromechanical coupled materials lies in the detection of the onset of failure and in the determination of the failure zone orientation. Existing failure criteria are based on the total energy release rate [65,8,66], the mechanical energy release rate [12], the local energy release rate [9,67], the strain energy density [68], the mechanical crack tip opening [69], or are based on a simplified form of the loss of ellipticity condition [43], among others. Extending the proposed treatment in [43] for plane problems to the 3D setting in this work has though resulted in an unsatisfactory performance with non-physical fracture zone orientations. It is for this reason, that we pursue a different strategy in this work through the application of a configurational force driven failure criterion based on the electromechanical coupled Eshelby stress tensor [70–73], where the numerical implementation proposed in [20,21,74–76] is adopted. As in [42], only the mechanical contribution of the electromechanical coupled configurational forces is considered as failure criterion in this work. The paper is organized as follows. Section 2 summarizes the classical electromechanical BVP within the continuum setting in Section 2.1 and the discrete setting in Section 2.2. The incorporation of strong discontinuities in 3D electromechanical coupled solids is outlined in Section 3. In particular, the incorporation of strong discontinuities in the electric potential is outlined in detail within the continuum setting in Section 3.1 and within the discrete setting in Section 3.2. The configurational force based failure criterion and crack propagation concept is discussed in Section 3.3. Section 4 describes the design of new finite elements through the development of three new electrical separation modes avoiding otherwise appearing locking phenomena in the 3D brick finite elements. The performance of the new finite elements is outlined in Section 5 based on three academic single element tests and two realistic tests in the form of a compact tension test and a three point bending test with different notch locations. The latter two problems allow us to validate the performance of the newly developed finite elements through a qualitative and quantitative comparison with experimental results available in the literature [12]. Several concluding remarks are given in Section 6. 2. The standard electromechanical BVP This section briefly summarizes the global governing equations of the classical electromechanical BVP within the infinitesimal range of interest in this work. We separately discuss the continuum and the discrete framework in Sections 2.1 and 2.2, respectively. Failure is not accounted for in this section with its treatment through the incorporation of strong discontinuities postponed to Section 3. 2.1. The electromechanical BVP in the continuum setting Let B  R3 be a electromechanical coupled solid with the primary unknowns of the mechanical displacement u and the electric potential u. The infinitesimal strain tensor e and the electric field e at material points x 2 B are then defined by

1 2

eðuÞ ¼ sym½ru ¼ ½ru þ ðruÞT  and eðuÞ ¼ ru

ð1Þ

in terms of the gradient operator r with respect to the coordinate x. For piezoelectric ceramics, the corresponding stress tensor r and the electric displacement field d follow as

r ¼ Ce  hT e and d ¼ he þ be

ð2Þ

145

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

in terms of the fourth order elasticity moduli tensor C, the third order piezoelectric moduli tensor h and the second order dielectric moduli tensor b whose explicit forms are given in Remark 2.1. The electromechanical coupled solid is subjected to a volumetric mechanical loading qb and a volumetric electrical loading of charge density qe . As illustrated in the center of Fig. 1, the boundary of the body can be separated into @ B ¼ @ u B [ @ t B , with @ u B \ @ t B ¼ ;, where the displacements are  on @ u B and the traction rn ¼ t are imposed on @ t B . Similarly, @ B ¼ @ u B [ @ q B with @ u B \ @ q B ¼ ;, prescribed as u ¼ u  is imposed on @ q B .  is prescribed on @ u B and the surface charge density d  n ¼ q where the electric potential u ¼ u The two governing field equations, which arise from this electromechanical coupled problem are given in their strong form as

div½r þ qb ¼ 0 and div½d ¼ qe

ð3Þ

which can be converted through standard arguments into their weak form given as

Z B

r : sym½rdu dV 

Z

qb  du dV 

Z

B

t  du dA ¼ 0

ð4Þ

@t B

and

Z B

d  rdu dV þ

Z

qe du dV þ

Z

B

 du dA ¼ 0 q

ð5Þ

@q B

for all test functions du with du ¼ 0 on @ u B and for all test functions du with du ¼ 0 on @ u B . 2.2. The electromechanical BVP in the discrete setting The continuum representation in (4) and (5) of the above Section 2.1 is the starting point for the derivation of the discrete equations outlined in this section. To do so, the body B is discretized into nelem finite elements B he as illustrated in the center of Fig. 2. The displacement field uh and the electric potential uh are approximated by uh ¼ Nu d and uh ¼ Nu / in terms of standard shape function Nu ; Nu and the nodal values d and /, respectively. The strain field eh and the electric field eh are then approximated by

eh ¼ B u d and eh ¼ B u /

ð6Þ

 u and B  u of standard displacement based, mixed or enhanced finite element formuin terms of generic ‘‘B-bar’’ operators B lations. With those approximations, the continuum governing equations (4) and (5) can be reformulated as nelem

Ru ¼

A

e¼1 nelem

Ru ¼

A e¼1

Z Bhe

Z B he

 T r dV  B u  T d dV þ B u

Z B he

Z B he

NTu b

dV 

!

Z @ t B he

NTu qe dV þ

NTu t

Z @ q B he

dA

ð7Þ !

 dA NTu q

ð8Þ

and solved through an iterative method by bringing the residual equations (7) and (8) down to zero. Remark 2.1. The explicit expressions for the elasticity moduli tensor C, the piezoelectric moduli tensor h and the dielectric moduli tensor b used in (2) for an anisotropic linear piezoelectric material model are given in an invariant representation [77,44] as

9 Cijkl ¼ kdij dkl þ lðdik djl þ dil djk Þ þ 12 a1 ðdik aj al þ dil aj ak þ ai djl ak þ ai djk al Þ þ 2a2 ai aj ak al þ a3 ðdij ak al þ ai aj dkl Þ > = hikl ¼ b1 ai dkl  b2 ai ak al  12 b3 ðdil ak þ dik al Þ

bik ¼ 2ðc1 dik þ c2 ai ak Þ

> ;

ð9Þ

in terms of the Lamé constants fk; lg, the material parameters fa1 ; a2 ; a3 g, the poling direction a, the axial piezoelectric expansion coefficient b1 , the lateral piezoelectric expansion coefficient b2 , the piezoelectric shearing parameter b3 , and the dielectric constants fc1 ; c2 g. It is to be noted though that the newly developed 3D finite elements are not restricted to the constitutive response of piezoelectric ceramics. As outlined in [44], also electric non-linearities such as electric displacement saturation or dissipative poling processes of unpoled ferroelectric ceramics are possible and an account for polarization switching can be made in those new elements.

3. The strong discontinuity approach for electromechanical coupled materials Following [43,58,60], this section outlines how to model failure zones Cx  R2 characterized through strong discontinuities in the displacement field and the electric potential locally at a material point x in the continuous and discrete representation of the 3D electromechanical coupled body described in Section 2.

146

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

g together with the Fig. 1. The electromechanical BVP in the continuum setting with the mechanical loading fqb; tg and the electrical loading fqe ; q partition of the mechanical boundary @ B ¼ @ u B [ @ t B and the partition of the electrical boundary @ B ¼ @ u B [ @ q B are illustrated in the center. Strong discontinuities Cx characterized by jumps fsut; sutg in the displacement field u and electric potential u appearing locally at a material point x are shown on the left and on the right, respectively.

The main contribution of this work is to capture strong discontinuities in the electric potential within a 3D framework. This is described in detail within the continuum setting in Section 3.1 and within the discrete setting in Section 3.2. Only a short summary of the incorporation of strong discontinuities in the displacement field is provided in Remark 3.1 at the end of this section. Interested readers are referred to [60] for a more detailed information. Section 3.3 outlines the configurational force based failure criterion and crack propagation concept used in this work. 3.1. Strong discontinuities in the electric potential within the continuum setting To capture strong discontinuities in the electric potential within the continuum setting locally at a material point x 2 B x of a local neighborhood B x  B , where a failure zone Cx is detected by the failure criterion discussed in Section 3.3, the local response of the electric potential has to be captured. To do so, we follow the standard strong discontinuity approach for electromechanical problems given in [43] and decompose the total electric potential ul as

~ ðsul tÞ in B x ul ¼ u þ u

ð10Þ

~ depending on the electric potential jumps sul t along the into a global continuous part u and a local discontinuous part u strong discontinuity Cx as shown on the right of Fig. 1. Analogously, the total electric field can be written as

el ¼ eðuÞ þ ~eðsul tÞ in B x n Cx

ð11Þ

neglecting the resulting singular Dirac delta measures on Cx . To determine the newly introduced local unknown sul t, the local equilibrium

Z Cx

dsul tðd  n þ qC Þ dA ¼ 0

ð12Þ

in terms of the surface charge density q ¼ d  n and the ‘‘traction-like’’ qC has to be satisfied in a weak sense along the strong discontinuity Cx . An impermeable electric boundary condition inside the crack is considered in this work. Interested readers are referred to [6,78–80] for a detailed discussion about different types of electric boundary conditions. 3.2. Strong discontinuities in the electric potential within the discrete setting To capture the strong discontinuities in the electric potential locally in the discrete setting, we consider a single finite element B he  B h with a detected failure surface Che as illustrated on the right of Fig. 2. New local degrees of freedom nu , which carry the electric information of the strong discontinuity and can be statically condensed out on the element level, are introduced. The jumps suhl t in the local electric potential can then be approximated in terms of those local parameters by

suhl tðs; tÞ ¼ Ju ðs; tÞnu

ð13Þ

in terms of the electrical jump interpolation function Ju ðs; tÞ where s and t are local coordinates in an artificial planar failure surface obtained as an approximation of the actual, possibly non-planar failure surface illustrated in Fig. 2. We refer to [61] for more information on the construction of this artificial failure surface. Following [43], (11) is approximated as

 u /  Cu n ehl ¼ B u

ð14Þ

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

147

Fig. 2. The electromechanical BVP in the discrete setting with the body B h being decomposed into nelem finite elements Bhe . The mechanical and electrical g, the partition of the mechanical boundary @ B h ¼ @ u Bh [ @ t B h and the partition of the electrical boundary @ B h ¼ @ u B h [ @ q B h are loadings fqb; t; qe ; q illustrated in the center. Strong discontinuities Che characterized by jumps fsuhl t; suhl tg in the displacement field uh and the electric potential uh , appearing locally in finite element B he where failure is detected, are shown on the left and right, respectively.

~ðsul tÞ in (11) is where the approximation of the global term eðuÞ in (11) comes from (6), and the discontinuous term e approximated in terms of nu and the electrical ‘‘compatibility operator’’ Cu . The determination of this operator will be outlined in detail in Section 4 through the incorporation of three new electrical separation modes. Next, (12) is used to determine the newly introduced degrees of freedom nu . Their possible static condensation on the element level represents one of the major advantages of the strong discontinuity approach, which remains for its extension to electromechanical coupled 3D solids in this work. The discrete form of (12) is given as

reu ¼

Z B he

ETu d dV 

Z Che

JTu qC dA

ð15Þ

in terms of the electrical ‘‘equilibrium operator’’ Eu which ensures the equilibrium of electric components along the strong discontinuity Che and projects the normal component of the electric displacement field d from the Gauss points onto the failure surface Che . Following [43], the test functions of the electric potential jumps are approximated as polynomials of order q as

dsuhl t ¼

q X si tj dnhuiji

ð16Þ

i;j¼0 iþj6q

in terms of the local coordinates s and t. The electric ‘‘equilibrium operator’’ Eu is obtained as

h i h10i h01i h20i h11i h02i Eu ¼ Eh00i u ; Eu ; Eu ; Eu ; Eu ; Eu ; . . .

ð17Þ

1 hiji hiji with Ehiji whose closed form are identical to the u ¼  he e n in terms of the element size he and polynomial functions e expression used for capturing strong discontinuities in mechanical displacements given in [61].

3.3. Failure initiation and propagation The propagation criterion based on a simplified form of the loss of ellipticity condition used in [43] for electromechanical coupled materials in the plane setting is not suitable for the determination of failure initiation and the determination of the propagation direction in the 3D setting considered in this work due to the stress component arising in thickness direction caused by the substantial electrical concentration at the crack tip. Instead, a Griffith-type criterion constructed by the concept of configurational forces and the critical energy release rate is considered here in the form

^ g int Þ ¼ jg int j  g 6 0 with g int ¼ g =Le /ð c int

ð18Þ

where g int is the configurational force computed at the integration points within the individual finite element, Le is the length of the crack front in the associated element and g c is the critical energy release rate [20,21]. The configurational forces g int at the integration points are obtained by the interpolation from their nodal counterparts computed as nelem

g node ¼  A

e¼1

Z

 T Rh dV B Bhe

ð19Þ

 in terms of the B-bar matrix and the electromechanical Eshelby stress tensor R given as [71–73]

R ¼ H1  ðruÞT r  ðruÞ  d

ð20Þ

148

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

where H ¼ 12 e : Ce  e : hT e  12 e  be is the electrical enthalpy. Note that whereas g int and g node have units of a force due to the  int in (18) has the correct units of the configuration force as force per integration over the 3D volume element dV in (19), g length. Also, g node is computed from existing quantities obtained from the electromechanical coupled problem without solving an additional equation, as stated in [74]. In line with the experimental investigation in [12], we assume that the evolution of the electromechanical damage is mechanically driven so that only the mechanical part of the electromechanical Eshelby stress tensor is considered. In particular, the electric contributions in both the enthalpy H and the fields are neglected in this work. The interested reader is referred to the literature for a more elaborate discussion and for alternative strategies on the challenging topic on how to find more physical rooted failure initiation criteria for electromechanical coupled materials [8,9,12,42–45,75,81–86]. It is to be noted, that the proposed numerical technique, and the newly developed 3D finite elements in this work do not rely on a particular failure initiation criterion and may easily accommodate alternative approaches. In addition to the detection of failure onset through (18), the orientation of a smooth failure surface within the solid must be determined. To do so, we make use of the recently developed marching cubes based crack propagation concept [61], which requires two tangential directions of the potential failure surface within each finite element as input. In addition to the tangential direction obtained from

s ¼ g int =jg int j

ð21Þ

we choose the unit vector in thickness direction as the second tangential direction t. This assumption is reasonable for the examples shown in Section 5, where only slight variations of the orientation of the failure surface in thickness direction are observed. Those tangential directions are used as input for the global tracking algorithm proposed in [51,87] to solve a heat conduction like problem, resulting in a heat-flux like level-set quantity. A unique and globally continuous isosurface corresponding to this level-set value is ensured within the individual brick finite elements by our recent application of the marching cubes algorithm to fracture problems [61]. The actual failure surface propagation subsequently takes place along that level set surface for those elements, where the failure criterion (18) is meet. Since needed for the strong discontinuity formulation, the normal of the failure surface can then be recomputed based on the actual geometry of the obtained failure surface. Remark 3.1. The strong discontinuities in the displacement field illustrated on the left of Fig. 1 in the continuous setting and on the left of Fig. 2 in the discrete setting are incorporated in this paper by following the recent works [60,61] on purely mechanical problems. There, the total mechanical displacement field ul in B x and the total strain field el in B xnCx are decomposed into a global continuous part and a local discontinuous part as

~ ðsul tÞ and ul ¼ u þ u

el ¼ eðuÞ þ ~eðsul tÞ

ð22Þ

in terms of displacement jumps sul t along the strong discontinuity Cx . The newly introduced unknown sul t must fulfill the local equilibrium equation in a weak sense as

Z

dsul t  ðrn  t C Þ dA ¼ 0

ð23Þ

Cx

between the traction t ¼ rn coming from the bulk of the materials and the traction t C resulting from the constitutive relation along the strong discontinuity Cx . In the discrete setting, new degrees of freedom nu are introduced in each element B he  B h , where a failure surface Che is detected, to approximate the displacement jumps and the local strain fields as

suhl tðs; tÞ ¼ Ju ðs; tÞnu

and

ehl ¼ B u d þ Cu nu

ð24Þ

in terms of the mechanical jump interpolation function Ju ðs; tÞ and the mechanical ‘‘compatibility operator’’ Cu with their close form given in Remark 4.1. The newly introduced degrees of freedom nu along the failure surface Cx , which again can be statically condensed out on the element level, are determined by the discrete form of (23), given as

reu ¼ 

Z B he

ETu r dV 

Z Che

JTu t C dA

ð25Þ

in terms of the mechanical ‘‘equilibrium operator’’ Eu . This operator has the form

  h10i h01i h20i h11i h02i Eu ¼ Eh00i u ; Eu ; Eu ; Eu ; Eu ; Eu ; . . .

ð26Þ

 hijni hijsi hijti  with Ehiji ¼  h1e ehiji ðn  kÞs with k ¼ n; s or t, ð  Þs is the symmetric part of ð  Þ and he is where Ehijki u ¼ Eu ; Eu ; Eu u the element size. The explicit form for ehiji can again be found in [61]. Remark 3.2. We refer to [43–45] for the general framework of the numerical implementation of electromechanical coupled solids with an account for failure through the strong discontinuity approach in the plane setting, which conceptually is identical to the implementation performed in this work, and to [60,61] for the general framework of the numerical implementation of purely mechanical problems in the 3D context.

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

149

4. Finite element design What remains to be shown is the closed form construction of the jump interpolation matrices fJu ; Ju g and the compatibility operators fCu ; Cu g. The derivation of the mechanical parts fJu ; Cu g is achieved by the incorporation of nine different mechanical separation modes initially proposed in [60] with a slight modification of the rotation mode around the normal of failure surface by an in-plane shear mode proposed in [61]. A short summary of those results is provided in Remark 4.1. The focus of the discussion here will though again be laid towards the derivation of the electrical parts fJu ; Cu g. Those have been initially derived in [43] for electromechanical coupled materials in the plane setting and will be extended to the fully 3D setting in this section. Following the general strategy of identifying certain separation modes as being critical for the accuracy of the strong discontinuity approach [55–63], we propose three new electrical separation modes in this work, incorporate those directly into the finite element formulation, and thereby are in a position to obtain the desired closed form solutions of fJu ; Cu g. 4.1. One constant electrical separation mode þ



Consider a single finite element B he which is fully separated into two parts B he and B he by a discontinuity Che . For the constant electrical separation mode illustrated on the left of Fig. 3, the electric potential is chosen as

( h00i /A

¼

h00i

ju ¼ nh00i u

for node A 2 B he

0

otherwise

þ

ð27Þ

based on which the constant electric potential jumps can be obtained directly as

 h00i  suhl th00i ¼ ju 

Che

¼ nh00i u :

ð28Þ

Insertion of (28) into (13) results in a direct solution for the jump interpolation Jh00i u for the constant electrical separation mode as

Jh00i u ¼ 1:

ð29Þ

For the constant electrical separation mode, the total electric field can be rewritten as h00i ehl;h00i ¼ ehh00i  Ch00i u nu

ð30Þ

based on (14), with the total electric field and the global electric field given as

ehl;h00i ¼ 0 and ehh00i ¼ 

X

 A nh00i : B u u

ð31Þ

þ A2B eh

Insertion of (31) into (30) yields the explicit form for the compatibility operator Ch00i u corresponding to the constant electrical separation mode in the form

Ch00i u ¼

X

A B u

ð32Þ

þ A2B he þ

as the summation of the electrical ‘‘B-bar’’ operators over the nodes A residing on B he .

h Fig. 3. Illustration of the electrical separation modes. The constant electrical separation mode is shown on the left in terms of a constant jump nh00i u at Ce with a constant but different electric potential (distinguished by its color) on both parts of the element. The linear separation mode along the s-direction is h h h01i shown in the center in terms of jumps nh10i u at Ce and the linear separation mode along the t-direction is shown on the right in terms of jumps nu at Ce . The  two latter modes are characterized by a constant electric potential in Bhe and a linearly distributed electric potential in the tangential directions s and t in þ Bhe (the varying color shows the linearly changing values of the electric potential).

150

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

4.2. Two linear electrical separation modes In addition to the constant electrical separation mode introduced in Section 4.1, two linear electrical separation modes are proposed to improve the behavior of the finite element formulation. We first consider the linear separation mode along the s-direction, as illustrated in the center of Fig. 3. The electrical potential corresponding to this mode is given as

( h10i /A

¼

h10i  ju ¼ nh10i u s  xA

for node A 2 B he

0

otherwise

þ

ð33Þ

A ¼ sA , where x A ¼ xA  xC with xA as the location of node A and xC as the center of failure surface and sA is the local with s  x coordinate of node A along the tangential direction s of the failure surface. Based on (33), the linear electrical jumps along the s-direction can be obtained as

 h10i  suhl th10i ¼ ju 

¼ s nh10i u :

Che

ð34Þ

Insertion of (34) into (13) yields the explicit form for the jump interpolation Jh10i u corresponding to this linear electrical separation mode in the form

Jh10i u ¼ s:

ð35Þ

For this linear electrical separation mode, the total electric field can be rewritten as h10i ehl;h10i ¼ ehh10i  Ch10i u nu

ð36Þ

based on (14), with the total electric field and the global electric field given as

ehl;h10i ¼ nh10i u sH Ch e

and ehh10i ¼ 

X

 A nh10i sA B u u

ð37Þ

þ A2B he þ



in terms of the Heaviside function H Ch which has a value 1 in B he and a value 0 in B he . Insertion of (37) into (36) yields the e explicit form for the compatibility operator Ch10i u corresponding to this linear electrical separation mode along the s-direction in the form

Ch10i u ¼ sH Ch  e

X

A sA B u

ð38Þ

þ A2Beh þ

as the summation of the electrical ‘‘B-bar’’ operators over the nodes A residing on B he and an additional part in terms of the Heaviside function H Ch accounting for the nonzero total electric field of this linear electrical separation mode. e Similarly, the linear electrical separation mode along the t-direction, as illustrated on the right of Fig. 3, yields expressions for the jump interpolation function and the compatibility operator as

Jh01i u ¼t

and Ch01i u ¼ tH Ch  e

X

A tA B u

ð39Þ

þ A2B he

A as the local coordinate of node A along the tangential direction t of the failure surface. with tA ¼ t  x þ As in the plane setting [43], the choice of B he does not affect the design of the element for the electrical separation modes in the 3D setting. Also, as in the plane setting [43], the columns of Cu corresponding to the different electric separation modes are linearly independent. Remark 4.1. The nine different mechanical separation modes used in [60,61] yield the closed form expressions of fJu ; Cu g. Those are summarized here as follows. h i h00i h00i h00ni h00si h00ti for the three constant The jump interpolation matrix Ju and the compatibility operator Cu ¼ Cu ; Cu ; Cu separation modes is given as

Jh00i ¼ ½n; s; t ; u

Cuh00ki ¼ 

X

Ak B u

ð40Þ

þ A2B he

where k stands for n; s or t. h i h10i h10i h10ni h10si h10ti for the three linear The jump interpolation matrix Ju and the compatibility operator Cu ¼ Cu ; Cu ; Cu separation modes along the s-direction are given as

Jh10i ¼ ½sn; ss; st ; u

Cuh10ni ¼ 

X

þ A2B he

 A ðn  sÞa x A ; Cuh10ki ¼ ðk  sÞs H Ch  B u e

X

 A ðk  sÞx A B u

ð41Þ

þ A2B he

where k stands for s or t; ð  Þs is the symmetric part of ð  Þ, and ð  Þa ¼ ð  Þ h  ð  Þ. i h01i h01i h01ni h01si h01ti for the three linear Similarly, the jump interpolation matrix Ju and the compatibility operator Cu ¼ Cu ; Cu ; Cu separation modes along the t-direction are given as

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

Jh01i ¼ ½tn; ts; tt ; u

Cuh01ni ¼ 

X

 A ðn  tÞa x A ; Ch01ki B ¼ ðk  tÞs H Ch  u u e

þ A2B he

X

 A ðk  tÞx A B u

151

ð42Þ

þ A2Beh

where k stands again for s or t.

5. Representative numerical simulations The numerical simulations presented in this section will outline the performance of the newly developed finite elements. Three academic single element tests are performed in Section 5.1 to outline their locking free property. A compact tension test in Section 5.2 and a three point bending test with three different locations of a pre-existing notch in Section 5.3 are investigated numerically and compared with experimental results obtained in [12]. The smallest finite element size used in those latter two simulations in any dimension is 0.14 mm. Qualitative comparisons of the predicted crack path for the electromechanically loaded PZT-4 piezoelectric ceramic specimens in the latter two tests, as well as quantitative results capturing the ultimate load sensitivity reported in [12] with regard to the amount of the applied electric field are given. 5.1. Single element tests Three academic numerical simulations are chosen specifically to test the newly proposed separation modes for elements experiencing strong discontinuities in the electric field. A block with dimension of a a a where a = 1 is discretized by a single bilinear brick finite element. The element is separated into two parts B  and B þ by a pre-existing mechanically fully softened and electrically impermeable strong discontinuity, which ensures a completely decoupled response of these two parts of the element. All the eight nodes of the element are mechanically restrained, whereas the applied electrical boundary conditions vary in each test. In the first element test, a nonzero constant electric potential d is applied at all nodes in B þ , as illustrated on the left of Fig. 4. The electric potential remains zero in the part of B  in the element due to the pre-existing electrical impermeable strong discontinuity C. The analytical solution of the electric field can be easily obtained from (1) as

e ¼ 0 in B 

and e ¼ 0 in B þ

ð43Þ

for both parts of the element. The right part of Fig. 4 shows the obtained numerical results which is in complete agreement with the analytical solution (43), which confirms the correctness of the constant electrical separation mode proposed in Section 4.1. Next, two different sets of nodal values for the electric potential, as shown in the top left illustration of Fig. 5, are applied on the nodes of part B þ of the element. This results in a nonzero constant electric field in m direction, which either corresponds to the tangential s- or t-direction. Similarly as in the previous test, the analytical solution of the electric fields is obtained as

e ¼ 0 in B 

and en ¼ 0;

em ¼

2d ; a

ek ¼ 0 in B þ

ð44Þ

for both parts of the element. Here k ¼ t when m ¼ s, and k ¼ s when m ¼ t with en ; es and et as the electric fields in normal and the two tangential directions of the strong discontinuity C. The numerical results using both, constant and linear electrical separation modes, are illustrated in Fig. 5, which again are in exact agreement with the analytical solution (44). One

Fig. 4. Single element test: illustration of the element test setting and the obtained numerical results. A single element has a mechanical fully softened and electrical impermeable pre-existing crack and is mechanically restrained at all nodes. An electric potential representing a constant electrical separation mode is applied at the nodes in B to generate a constant non-zero electric potential distribution in Bþ . The computed resulting electric field in both parts of the element is illustrated on the right, which exactly captures the analytical solution (43), being zero in all three directions n; s, and t.

152

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

Fig. 5. Single element test: illustration of the element test setting and the obtained numerical results. A single element has a mechanical fully softened and electrical impermeable pre-existing crack and is mechanically restrained at all nodes. Two sets of electric potential representing the two linear electric separation modes in m ¼ s (and k ¼ t) or m ¼ t (and k ¼ s) direction are applied at the nodes in B to generate a linear electric potential distribution in B þ . The computed resulting electric fields for these two tests have the same magnitude illustrated on the right and differ only by their direction in the part B þ . Locking is observed when only the constant electrical separation mode is used. Contrary, the analytical solution (44) is captured exactly with the incorporation of the two linear electrical separation modes.

can also observe that using only the constant electric separation mode is insufficient to capture the exact solution. This observed locking phenomenon is overcome when using all the three newly developed electrical separation modes. 5.2. Compact tension test The second example is the compact tension test of a PZT-4 piezoelectric ceramics block with dimensions of 25:5 19:1 5:1 mm3. The material parameters of PZT-4 piezoelectric ceramics used in the simulation are summarized in Table 1. The geometry, loading and boundary conditions of the problem are depicted on the left of Fig. 6. The specimen with a horizontal pre-existing notch with a width of 0.46 mm and a length of 11.5 mm is loaded mechanically by the application of displacements through two rigid cylindrical bars in opposite directions and electrically by the application of a zero electric potential ubottom on the bottom surface and a non-zero electric potential utop at the top surface. A horizontal failure

Table 1 Material parameters for PZT-4 piezoelectric ceramics in the bulk and along the strong discontinuity used in the compact tension test in Section 5.2 and the three point bending test in Section 5.3. Name

Parameter

Value

Lamé parameter

k

7:78 104

N/mm2

Lamé parameter

l a1 a2 a3

4

3:06 10

N/mm2

1:0 104

N/mm2

Material parameter Material parameter Material parameter Axial piezoelectric expansion Lateral piezoelectric expansion Piezoelectric shearing Dielectric constant

b1 b2 b3

Dielectric constant

c1 c2

Fracture energy release rate

Gf

Unit

2

5:0 10

3:5 103 6:98 6:06 26.88, 16.0

N/mm2 N/mm2

3:0 103

N/(kV mm) N/(kV mm) N/(kV mm) mC/(kV m)

2:65 104 2.34, 3.08

N/m

mC/(kV m)

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

153

Fig. 6. Compact tension test: illustration of geometry, loading, boundary conditions and the finite element discretization of the PZT-4 piezoelectric ceramic block with a pre-existing horizontal notch. The specimen is loaded mechanically by the application of opposite displacements through two rigid cylindrical bars close to the left top and bottom corners. The electric loading is achieved through a zero electric potential ubottom at the bottom surface and a non-zero electric potential utop at the top surface. The poling direction is vertical from the bottom towards the top. The finite element meshes used for the simulations consist of a coarse mesh with 3315 Q1 brick elements and a fine mesh with 5435 Q1 brick elements with five elements in the thickness direction for both.

Fig. 7. Compact tension test: evolution of crack propagation and illustration of the electric potential using the meshes in Fig. 6 with an applied electric potential utop ¼ 9:95 kV. The material parameters b3 ¼ 26:88 N/(kV mm) and Gf ¼ 2:34 N/mm from [12] are used.

surface propagating from the tip of the pre-existing notch towards the right surface of the specimen is expected based on experimental results reported in [12] and based on numerically results reported in [42,43] for the 2D setting. The finite element discretization considers a coarse mesh with 3315 Q1 brick elements and a fine mesh with 5435 Q1 brick elements with five elements in thickness direction for both meshes, as shown in the center and on the right of

154

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

20 10

Reaction force [N]

0 −10 −20 −30 −40

constant speration modes linear seperation modes

−50 0

0.5

1 1.5 Applied displacement u[mm]

2 −3

x 10

Fig. 8. Compact tension test: reaction force versus applied displacement relation for the coarse mesh in Fig. 6 and different separation modes using the material parameter in Table 1 with b3 ¼ 26:88 N/kV mm and Gf ¼ 2:34 N/mm for an applied electric potential utop ¼ 9:95 kV. Significant improvements of the numerical results can be observed by incorporating the linear separation modes.

Fig. 6. The failure surfaces and electric potential at different loading stages for these two meshes with an applied electric potential utop ¼ 9:95 kV are shown in Fig. 7 with the material parameters of b3 ¼ 26:88 N/(kV mm) for the piezoelectric shearing parameter and Gf ¼ 2:34 N/m for the fracture energy release rate reported in [12]. To obtain the smooth failure surface shown in Fig. 7, it was essential for us to use the marching cubes based propagation concept proposed in [61] and the configurational force based failure criterion (18). Next, the influence of the separation modes is investigated. Simulations are carried out for the coarse mesh in Fig. 6 using constant and linear separation modes and the same material parameters as above with an applied electric potential of utop ¼ 9:95 kV. The reaction force versus applied displacement relation is reported in Fig. 8 where locking is observed when using constant separation modes, which is significantly improved by the incorporation of the linear separation modes. When performing a quantitative evaluation of the simulation results, it turned out that the peak value of the reaction force in Fig. 8 underestimates the experimentally observed fracture load given in [12]. Following [42], we therefore investigated the effect of the piezoelectric shearing constant b3 and the critical energy release rate Gf on the resultant fracture load. In particular, the value b3 ¼ 16:00 N/(kV mm) reported in [42] is used together with a fitted value for the fracture energy Gf ¼ 3:08 N/m. Using those parameters, the reaction force versus applied displacement relation for an applied electric potential utop ¼ 9:95 kV using the coarse mesh in Fig. 6 is reported in Fig. 9. The fracture loads corresponding to the ultimate strength of the specimen using the two meshes in Fig. 6, for different electric fields (varying from 5 kV/cm to +10 kV/cm), with those four combinations of b3 and Gf are reported in Fig. 10. It can be observed that the piezoelectric shearing parameter b3 and the critical energy release rate Gf show a large influence on the numerical results and that an excellent match of the fracture loads between the numerical results and the experimental data is obtained for different applied electric fields when using b3 ¼ 16:00 N/(kV mm) and Gf ¼ 3:08 N/m. β =−26.88, G =2.34 3

Reaction force [N]

f

β =−26.88, G =3.08

60

3

f

β3=−16.00, Gf=2.34

40

β =−16.00, G =3.08 3

f

20 0 −20 −40 0

1

2 3 Applied displacement u[mm]

4

5 −3

x 10

Fig. 9. Compact tension test: reaction force versus applied displacement relation for the coarse finite element mesh in Fig. 6 with different combinations of b3 (unit: N/(kV mm)) and Gf (unit: N/m) for an applied electric potential utop ¼ 9:95 kV.

155

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160 200

200 exprimental data coarse mesh fine mesh

exprimental data coarse mesh fine mesh 150

Fracture load [N]

Fracture load [N]

150

100

50

0

−50

100

50

0

−5

−2.5

0

2.5 5 electric field [kV/cm]

7.5

−50

10

200

−5

−2.5

0

2.5 5 electric field [kV/cm]

exprimental data coarse mesh fine mesh 150

Fracture load [N]

150

Fracture load [N]

10

200 exprimental data coarse mesh fine mesh

100

50

0

−50

7.5

100

50

0

−5

−2.5

0

2.5 5 electric field [kV/cm]

7.5

10

−50

−5

−2.5

0

2.5 5 electric field [kV/cm]

7.5

10

Fig. 10. Compact tension test: influence of the piezoelectric shearing parameter b3 and the critical energy release rate Gf on the fracture loads using the meshes in Fig. 6 with different applied electric fields.

Fig. 11. Three point bending test: illustration of geometry loading and boundary conditions of the PZT-4 piezoelectric ceramic block with a pre-existing vertical notch at three different locations (A, B, and C). The specimen is mechanically supported at two locations on the bottom surface and loaded mechanically by the application of a vertical downward displacement at the center of the top surface. The electric loading is applied through a zero electric potential uleft at the left surface and a non-zero electric potential uright at the right surface. The poling direction is horizontally oriented from the left towards the right.

156

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

coarse meshes

fine meshes

Position A

Position B

Position C

Fig. 12. Three point bending test: illustration of the finite element meshes used for the simulation. Coarse meshes consisting of 2735, 2975, 2905 Q1 brick elements and fine meshes consisting of 4845, 4895, 5055 Q1 brick elements for notch locations A, B, and C are considered with five elements in thickness direction for each mesh.

5.3. Three point bending test Finally, a three point bending test of a PZT-4 piezoelectric ceramics beam with three different notch locations is carried out. The geometry, loading and boundary conditions of the problem are illustrated in Fig. 11 where a rectangular block with dimensions of 17:1 9 5:1 mm3 and with a pre-existing vertical notch of a width of 0.46 mm and a length of 4 mm at three different locations is restrained mechanically at two positions on the bottom surface and loaded mechanically by applying a vertical downward displacement at the center of the top surface and electrically by applying a zero electric potential uleft at the left surface and a non-zero electric potential uright at the right surface. The same material parameters as before are used. Different failure surfaces propagating from the tip of the pre-existing notch locations toward the center of the top surface are expected as reported experimentally [12] and numerically within the plane setting [42,43].

coarse meshes

fine meshes

Position A

Position B

Position C

Fig. 13. Three point bending test: illustration of the final failure surface together with the experimental results for the different notch locations A, B, and C using coarse and fine meshes with uright ¼ 9:95 kV, b3 ¼ 26:88 N/(kV mm) and Gf ¼ 2:34 N/m. The front view of the final failure surfaces obtained from experiments reported in [12] are given as the shaded region on the front surfaces of the specimen.

157

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

The finite element discretization considers both, coarse meshes of 2735, 2975, 2905 Q1 brick elements, and fine meshes of 4845, 4895, 5055 Q1 brick elements for the notch positions A, B, and C, respectively. All meshes have five elements in thickness direction, as shown in Fig. 12. The final failure surfaces are shown in Fig. 13 for an applied electric potential uright ¼ 9:95 kV using the original material parameter used in [12] with the piezoelectric shearing parameter b3 ¼ 26:88 N/(kV mm) and the fracture energy release rate Gf ¼ 2:34 N/m where the experimental results reported in [12] are given as the shaded region on the front surface of the specimen in Fig. 13. A very good match between the numerical results and the experimental results can be observed. A similar numerical study as in Section 5.2 concerning the dependency of the ultimate fracture load on the piezoelectric shearing constant b3 and the critical energy release rate Gf is performed for the three point bending test. The reaction force versus applied displacement relation is shown in Fig. 14 for the four combinations of the different values of b3 ¼ 26:88 N/ (kV mm), b3 ¼ 16:00 N/(kV mm) and Gf ¼ 2:34 N/m, Gf ¼ 3:08 N/m when using the coarse mesh and notch location A in Fig. 12 for an applied electric potential of uright ¼ 9:95 kV. Similar as for the compact tension test, a large influence of b3 and Gf on the reaction force curve is observed. Finally, the influence of the applied electric field on the fracture load is investigated quantitatively for the different meshes in Fig. 12 with the best combination again obtained for the piezoelectric shearing constant b3 ¼ 16:00 N/(kV mm) and the critical energy release rate with Gf ¼ 3:08 N/m for different applied electric potentials uright (the electric field varies from 5 kV/cm to +10 kV/cm) applied at the right surface. The results are reported in Fig. 15 and an excellent match can be observed between the numerical and the experimental results for both, the coarse and fine meshes.

100

β =−26.88, G =2.34 3 3

80 Reaction force [N]

f

β =−26.88, G =3.08 f

β3=−16.00, Gf=2.34 β =−16.00, G =3.08 3

60

f

40 20 0 −20

0

0.5

1 1.5 2 2.5 Applied displacement u[mm]

3

3.5 −3

x 10

250

250

200

200

Fracture load [N]

Fracture load [N]

Fig. 14. Three point bending test: reaction force versus applied displacement relation for the coarse finite element mesh in Fig. 12 and notch position A with different combinations of b3 (unit: N/(kV mm)) and Gf (unit: N/m) for an applied electric potential uright ¼ 9:95 kV.

150

100

exprimental data A exprimental data B exprimental data C Notch A coarse mesh Notch B coarse mesh Notch C coarse mesh

50

0

−5

−2.5

0

2.5 5 electric field [kV/cm]

150

100

exprimental data A exprimental data B exprimental data C Notch A fine mesh Notch B fine mesh Notch C fine mesh

50

0

7.5

10

−5

−2.5

0

2.5 5 electric field [kV/cm]

7.5

10

Fig. 15. Three point bending test: influence of applied electric field on the fracture loads of the three point bending test for the three different notch locations A, B, and C using the coarse (left) and fine (right) meshes in Fig. 12 with b3 ¼ 16:00 N/(kV mm) and Gf ¼ 3:08 N/m.

158

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

6. Conclusion We have presented in this work new finite elements with embedded strong discontinuities to model failure in 3D electromechanical coupled materials. This has been achieved based on an extension of the framework presented in [43] for the plane setting to the fully 3D setting in this work. In addition to nine purely mechanical separation modes, three new electrical separation modes are developed and incorporated into the new finite element formulation to accurately capture jumps in the electric potential without the appearance of locking. Those strong discontinuities are detected by a configurational force based failure criterion, which only considers the mechanical contribution in the electromechanical Eshelby stress tensor. The smoothness and continuity of the 3D failure surfaces is ensured by the incorporation of the marching cubes based failure propagation concept previously developed in [61]. The performance of the newly proposed finite elements have been evaluated qualitatively and quantitatively by three academic single element tests and two realistic simulations. In particular, the dependency of the fracture load on the applied electric field is accurately captured, after performing a parameter study to investigate the effect of the piezoelectric shearing parameter and the fracture energy release rate. Future extensions are e.g., possible when dealing with electric boundary conditions inside the crack being different from the impermeable assumption used in this work [6,78–80]. Furthermore, the account of domain switching in the vicinity of the crack tip [88,89] will improve the overall performance of the proposed numerical technique. Also, extensions to model failure in electroactive polymers are planned. Those will rely on a combination of this work with our recently developed constitutive models for soft matter materials [90–92]. Acknowledgments Financial support was provided by the Haythornthwaite Research Initiation Grant of the ASME Applied Mechanics Division, and the UPS Endowment Fund at Stanford University. This support is gratefully acknowledged. References [1] C. Truesdell, R. Toupin, The classical field theories, in: S. Flügge (Ed.), Principles of Classical Mechanics and Field Theory, Encyclopedia of Physics, vol. III/1, Springer-Verlag, Berlin, 1960, pp. 226–795. [2] L.D. Landau, E.M. Lifshitz, The classical theory of fields, Course of Theoretical Physics, vol. 2, Pergamon Press, Oxford, 1975. [3] M. Kamlah, Ferroelectric and ferroelastic piezoceramics – modeling of electromechanical hysteresis phenomena, Continuum Mech. Thermodyn. 13 (2001) 219–268. [4] C.M. Landis, Non-linear constitutive modeling of ferroelectrics, Curr. Opin. Solid State Mater. Sci. 8 (2004) 59–69. [5] J.E. Huber, Micromechanical modelling of ferroelectrics, Curr. Opin. Solid State Mater. Sci. 9 (2005) 100–106. [6] V.Z. Parton, Fracture mechanics of piezoelectric materials, Acta Astronaut. 3 (1976) 671–683. [7] Y.E. Pak, Linear electro-elastic fracture mechanics of piezoelectric materials, Int. J. Fract. 54 (1) (1992) 79–100. [8] Z. Suo, C. Kuo, D. Barnett, J. Willis, Fracture mechanics for piezoelectric ceramics, J. Mech. Phys. Solids 40 (4) (1992) 739–765. [9] H. Gao, T. Zhang, P. Tong, Local and global energy release rates for an electrically yielded crack in a piezoelectric ceramic, J. Mech. Phys. Solids 45 (1997) 491–510. [10] R.M. McMeeking, Towards a fracture mechanics for brittle piezoelectric and dielectric materials, Int. J. Fract. 108 (2001) 25–41. [11] T.-Y. Zhang, M. Zhao, P. Tong, Fracture of piezoelectric ceramics, Adv. Appl. Mech. 38 (2002) 147–289. [12] S. Park, C.T. Sun, Fracture criteria for piezoelectric ceramics, J. Am. Ceram. Soc. 78 (6) (1995) 1475–1480. [13] H. Wang, R.N. Singh, Crack propagation in piezoelectric ceramics: effects of applied electric fields, J. Appl. Phys. 81 (1997) 7471–7479. [14] R. Fu, T.Y. Zhang, Effects of an electric field on the fracture toughness of poled lead zirconate titanate ceramics, J. Am. Ceram. Soc. 83 (2000) 1215–1218. [15] X. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, J. Mech. Phys. Solids 42 (1994) 1397–1434. [16] M. Ortiz, A. Pandolfi, Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis, Int. J. Numer. Methods Eng. 44 (1999) 1267–1282. [17] M. Ortiz, J. Quigley, Adaptive mesh refinement in strain localization problems, Comput. Methods Appl. Mech. Eng. 90 (1991) 781–804. [18] A. Pandolfi, M. Ortiz, An efficient adaptive procedure for three-dimensional fragmentation simulations, Eng. Comput. 18 (2002) 148–159. [19] P. Heintz, F. Larsson, P. Hansbo, K. Runesson, Adaptive strategies and error control for computing material forces in fracture mechanics, Int. J. Numer. Methods Eng. 60 (2004) 1287–1299. [20] C. Miehe, E. Gürses, A robust algorithm for configurational-force-driven brittle crack propagation with r-adaptive mesh alignment, Int. J. Numer. Methods Eng. 72 (2007) 127–155. [21] E. Gürses, C. Miehe, A computational framework of three-dimensional configurational-force-driven brittle crack propagation, Comput. Methods Appl. Mech. Eng. 198 (2009) 1413–1428. [22] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Eng. 199 (2010) 2765–2778. [23] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically-consistent phase field models of fracture: variational principles and multifield FE implementations, Int. J. Numer. Methods Eng. 83 (2010) 1273–1311. [24] C. Kuhn, R. Müller, A new finite element technique for a phase field model of brittle fracture, J. Theor. Appl. Mech. 49 (4) (2011) 1115–1133. [25] H. Gomez, T.J. Hughes, Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models, J. Chem. Phys. 230 (13) (2011) 5310–5327. [26] M. Hofacker, C. Miehe, Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation, Int. J. Fract. 178 (1-2) (2012) 113–129. [27] M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J. Hughes, C.M. Landis, A phase-field description of dynamic brittle fracture, Comput. Methods Appl. Mech. Eng. 217 (2012) 77–95. [28] M. Hofacker, C. Miehe, A phase field model of dynamic fracture: robust field updates for the analysis of complex crack patterns, Int. J. Numer. Methods Eng. 93 (3) (2013) 276–301. [29] C.V. Verhoosel, R. Borst, A phase-field model for cohesive fracture, Int. J. Numer. Methods Eng. 96 (1) (2013) 43–62. [30] M.N. da Silva Jr., F.P. Duda, E. Fried, Sharp-crack limit of a phase-field model for brittle fracture, J. Mech. Phys. Solids 61 (11) (2013) 2178–2195. [31] T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Methods Eng. 45 (1999) 601–620. [32] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Methods Eng. 46 (1999) 131–150.

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

159

[33] G. Wells, L. Sluys, A new method for modelling cohesive cracks using finite elements, Int. J. Numer. Methods Eng. 50 (2001) 2667–2682. [34] A. Hansbo, P. Hansbo, A finite element method for the simulation of strong and weak discontinuities in solid mechanics, Comput. Methods Appl. Mech. Eng. 193 (2004) 3523–3540. [35] J.C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comput. Mech. 12 (1993) 277–296. [36] F. Armero, K. Garikipati, An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids, Int. J. Solids Struct. 33 (1996) 2863–2885. [37] R. Regueiro, R. Borja, A finite element model of localized deformation in frictional materials taking a strong discontinuity approach, Finite Elem. Anal. Des. 33 (1999) 283–315. [38] R. Borja, A finite element model for strain localization analysis of strongly discontinuous fields based on standard Galerkin approximation, Comput. Methods Appl. Mech. Eng. 190 (2000) 1529–1549. [39] J. Oliver, A.E. Huespe, M.D.G. Pulido, E. Samaniego, On the strong discontinuity approach in finite deformation settings, Int. J. Numer. Methods Eng. 56 (2003) 1051–1082. [40] R. Radulovic, O.T. Bruhns, J. Mosler, Effective 3d failure simulations by combining the advantages of embedded strong discontinuity approaches and classical interface elements, Eng. Fract. Mech. 78 (12) (2011) 2470–2485. [41] E. Béchet, M. Scherzer, M. Kuna, Application of the X-FEM to the fracture of piezoelectric materials, Int. J. Numer. Methods Eng. 77 (2009) 1535–1565. [42] C. Miehe, F. Welschinger, M. Hofacker, A phase field model of electromechanical fracture, J. Mech. Phys. Solids 58 (2010) 1716–1740. [43] C. Linder, D. Rosato, C. Miehe, New finite elements with embedded strong discontinuities for the modeling of failure in electromechanical coupled solids, Comput. Methods Appl. Mech. Eng. 200 (2011) 141–161. [44] C. Linder, C. Miehe, Effect of electric displacement saturation on the hysteretic behavior of ferroelectric ceramics and the initiation and propagation of cracks in piezoelectric ceramics, J. Mech. Phys. Solids 60 (2012) 882–903. [45] C. Linder, An analysis of the exponential electric displacement saturation model in fracturing piezoelectric ceramics, Tech. Mech. 32 (2012) 53–69. [46] R. Bhargava, K. Sharma, A study of finite size effects on cracked 2-d piezoelectric media using extended finite element method, Comput. Mater. Sci. 50 (6) (2011) 1834–1845. [47] A. Abdollahi, I. Arias, Phase-field modeling of crack propagation in piezoelectric and ferroelectric materials with different electromechanical crack conditions, J. Mech. Phys. Solids 60 (12) (2012) 2100–2126. [48] H. Nguyen-Vinh, I. Bakar, M. Msekh, J.-H. Song, J. Muthu, G. Zi, P. Le, S. Bordas, R. Simpson, S. Natarajan, T. Lahmer, T. Rabczuk, Extended finite element method for dynamic fracture of piezo-electric materials, Eng. Fract. Mech. 92 (2012) 19–31. [49] Z. Wilson, M. Borden, C.M. Landis, A phase-field model for fracture in piezoelectric ceramics, Int. J. Fract. 183 (2013) 135–153. [50] S. Nanthakumar, T. Lahmer, T. Rabczuk, Detection of flaws in piezoelectric structures using extended FEM, Int. J. Numer. Methods Eng. 96 (6) (2013) 373–389. [51] J. Oliver, A.E. Huespe, E. Samaniego, E.W.V. Chaves, Continuum approach to the numerical simulation of material failure in concrete, Int. J. Numer. Anal. Methods Geomech. 28 (7-8) (2004) 609–632. [52] J. Mosler, G. Meschke, 3D modelling of strong discontinuities in elastoplastic solids: fixed and rotating localization formulations, Int. J. Numer. Methods Eng. 57 (2003) 1553–1576. [53] P. Steinmann, A finite element formulation for strong discontinuities in fluid-saturated porous media, Mech. Cohes. Frict. Mater. 4 (1999) 133–152. [54] C. Callari, F. Armero, Finite element methods for the analysis of strong discontinuities in coupled poro-plastic media, Comput. Methods Appl. Mech. Eng. 191 (2002) 4371–4400. [55] D. Ehrlich, F. Armero, Finite element methods for the analysis of softening plastic hinges in beams and frames, Comput. Mech. 35 (2005) 237–264. [56] F. Armero, D. Ehrlich, Finite element methods for the multi-scale modeling of softening hinge lines in plates at failure, Comput. Methods Appl. Mech. Eng. 195 (2006) 1283–1324. [57] F. Armero, Strong discontinuities in antiplane/torsional problems of computational failure mechanics, Int. J. Fract. 178 (2012) 3–32. [58] C. Linder, F. Armero, Finite elements with embedded strong discontinuities for the modeling of failure in solids, Int. J. Numer. Methods Eng. 72 (2007) 1391–1433. [59] F. Armero, C. Linder, New finite elements with embedded strong discontinuities in the finite deformation range, Comput. Methods Appl. Mech. Eng. 197 (2008) 3138–3170. [60] F. Armero, J. Kim, Three-dimensional finite elements with embedded strong discontinuities to model material failure in the infinitesimal range, Int. J. Numer. Methods Eng. 91 (2012) 1291–1330. [61] C. Linder, X. Zhang, A marching cubes based failure surface propagation concept for three-dimensional finite elements with non-planar embedded strong discontinuities of higher-order kinematics, Int. J. Numer. Methods Eng. 96 (2013) 339–372. [62] F. Armero, C. Linder, Numerical simulation of dynamic fracture using finite elements with embedded discontinuities, Int. J. Fract. 160 (2009) 119–141. [63] C. Linder, F. Armero, Finite elements with embedded branching, Finite Elem. Anal. Des. 45 (2009) 280–293. [64] C. Linder, A. Raina, A strong discontinuity approach on multiple levels to model solids at failure, Comput. Methods Appl. Mech. Eng. 253 (2013) 558– 583. [65] Y.E. Pak, Crack extension force in a piezoelectric material, J. Appl. Mech. 57 (3) (1990) 647–653. [66] Q. Li, M. Kuna, Evaluation of electromechanical fracture behavior by configurational forces in cracked ferroelectric polycrystals, Comput. Mater. Sci. 57 (2012) 94–101. [67] D. Fang, B. Liu, K. Hwang, Energy analysis on fracture of ferroelectric ceramics, Int. J. Fract. 100 (4) (2000) 401–408. [68] S. Shen, T. Nishioka, Fracture of piezoelectric materials: energy density criterion, Theor. Appl. Fract. Mech. 33 (1) (2000) 57–65. [69] C. Xue, H. Yong, Y. Zhou, The mechanical crack tip opening displacement fracture criterion in piezoelectric ceramics, Eng. Fract. Mech. 96 (2012) 606– 614. [70] Y.E. Pak, G. Herrmann, Conservation laws and the material momentum tensor for the elastic dielectric, Int. J. Eng. Sci. 24 (8) (1986) 1365–1374. [71] R. Kienzler, G. Herrmann, Mechanics in Material Space: with Applications in Defect and Fracture Mechanics, Springer-Verlag, New York, 2000. [72] G.A. Maugin, Configurational Forces: Thermomechanics, Physics, Mathematics, and Numerics, Taylor & Francis, USA, 2010. [73] S.R.M. Khalaquzzaman, B. Xu, R. Müller, Computational homogenization of piezoelectric materials using FE2 to determine configurational forces, Tech. Mech. 32 (2012) 21–37. [74] R. Müller, G. Maugin, On material forces and finite element discretizations, Comput. Mech. 29 (2002) 52–60. [75] R. Mueller, S. Kolling, D. Gross, On configurational forces in the context of the finite element method, Int. J. Numer. Methods Eng. 53 (7) (2002) 1557– 1574. [76] C. Miehe, E. Gürses, M. Birkle, A computational framework of configurational-force-driven brittle fracture based on incremental energy minimization, Int. J. Fract. 145 (4) (2007) 245–259. [77] J. Schröder, D. Gross, Invariant formulation of the electromechanical enthalpy function of transversely isotropic piezoelectric materials, Arch. Appl. Mech. 73 (2004) 533–552. [78] W.F. Deeg, The analysis of dislocation, crack and inclusion problems in piezoelectric solids (Ph.D. thesis), Stanford University, 1980. [79] T.H. Hao, Z.Y. Shen, A new electric boundary condition of electric fracture mechanics and its application, Eng. Fract. Mech. 47 (1994) 793–802. [80] C.M. Landis, Energetically consistent boundary conditions for electromechanical fracture, Int. J. Solids Struct. 4 (2004) 6291–6315. [81] A.G. Tobin, Y.E. Pak, Effect of electric fields on fracture behavior of PZT ceramics, in: V.K. Varadan (Ed.), Proceedings of SPIE, Smart Structure and Materials, vol. 78, International Society for Optics and Photonics, 1993, pp. 78–86. [82] R. McMeeking, Crack tip energy release rate for a piezoelectric compact tension specimen, Eng. Fract. Mech. 64 (1999) 217–244.

160

C. Linder, X. Zhang / Comput. Methods Appl. Mech. Engrg. 273 (2014) 143–160

[83] H. Jelitto, H. Keßler, G. Schneider, H. Balke, Fracture behavior of poled piezoelectric PZT under mechanical and electrical loads, J. Eur. Ceram. Soc. 25 (5) (2005) 749–757. [84] H. Jelitto, F. Felten, C. Häusler, H. Kessler, H. Balke, G. Schneider, Measurement of energy release rates for cracks in PZT under electromechanical loads, J. Eur. Ceram. Soc. 25 (12) (2005) 2817–2820. [85] R. Mueller, D. Gross, D. Lupascu, Driving forces on domain walls in ferroelectric materials and interaction with defects, Comput. Mater. Sci. 35 (1) (2006) 42–52. [86] M. Kuna, Fracture mechanics of piezoelectric materials – where are we right now?, Eng Fract. Mech. 77 (2010) 309–326. [87] J. Oliver, A.E. Huespe, E. Samaniego, E.W.V. Chaves, On strategies for tracking strong discontinuities in computational failure mechanics, in: H. Mang, F. Rammerstorfer, J. Eberhardsteiner (Eds.), Fifth World Congress on Computational Mechanics, Vienna University of Technology, Austria, December 7, 2002, pp. 1–15. [88] B.-X. Xu, D. Schrade, D. Gross, R. Mueller, Phase field simulation of domain structures in cracked ferroelectrics, Int. J. Fract. 165 (2010) 163–173. [89] B.-X. Xu, D. Schrade, D. Gross, R. Mueller, Fracture simulation of ferroelectrics based on the phase field continuum and a damage variable, Int. J. Fract. 166 (2010) 163–172. [90] C. Linder, M. Tkachuk, C. Miehe, A micromechanically motivated diffusion-based transient network model and its incorporation into finite rubber viscoelasticity, J. Mech. Phys. Solids 59 (2011) 2134–2156. [91] M. Tkachuk, C. Linder, The maximal advance path constraint for the homogenization of materials with random network microstructure, Philos. Mag. 92 (2012) 2779–2808. [92] A. Raina, C. Linder, A homogenization approach for nonwoven materials based on fiber undulations and reorientation, J. Mech. Phys. Solids 65 (2014) 12–34.