# Phase-field modeling of crack propagation in piezoelectric

## Transcript Of Phase-field modeling of crack propagation in piezoelectric

Phase-ﬁeld modeling of crack propagation in piezoelectric and ferroelectric materials with different electromechanical crack conditions

Amir Abdollahi, Irene Arias n

Laboratori de Calcul Numeric (LaCa N), Departament de Matematica Aplicada III, Universitat Politecnica de Catalunya (UPC), Campus Nord UPC-C2, E-08034 Barcelona, Spain

abstract

We present a family of phase ﬁeld models for fracture in piezoelectric and ferroelectric materials. These models couple a variational formulation of brittle fracture with, respectively, (1) the linear theory of piezoelectricity, and (2) a Ginzburg Landau model of the ferroelectric microstructure to address the full complexity of the fracture phenomenon in these materials. In these models, both the cracks and the ferroelectric domain walls are represented in a diffuse way by phase ﬁelds. The main challenge addressed here is encoding various electromechanical crack models (introduced as crack face boundary conditions in sharp models) into the phase ﬁeld framework. The proposed models are veriﬁed through comparisons with the corresponding sharp crack models. We also perform two dimensional ﬁnite element simulations to demonstrate the effect of the different crack face conditions, the electromechanical loading and the media ﬁlling the crack gap on the crack propagation and the microstructure evolution. Salient features of the results are compared with experiments.

1. Introduction

The design and implementation of electromechanical systems demand high performance materials regarding their coupling behavior and reliability. Piezoelectric ceramics are very prominent in this ﬁeld, exhibiting strong electromecha nical coupling with short response times. However, their inherent brittleness is a serious obstacle to their reliable operation in devices, which demands a deep understanding of the fracture behavior. Towards this goal, numerous theoretical and experimental investigations have been carried out during the past decades on the fracture of piezo ceramics. Comprehensive reviews of these works have been presented (Zhang and Gao, 2004; Schneider, 2007; Kuna, 2010). Most piezoelectrics also exhibit ferroelectric and ferroelastic switching behavior with macroscopic dielectric and butterﬂy hysteresis. For related modeling approaches, see the reviews by Kamlah (2001), Landis (2004b) and Huber (2005). Microstructural domains often nucleate and evolve under high electromechanical loadings and near load concentrations such as those produced at crack tips (Hackemann and Pfeiffer, 2003; Jones et al., 2007; Pojprapai et al., 2011). The interactions between the microstructure and the localized stress and electric ﬁelds near the crack tips are responsible for

n Corresponding author. Tel.: þ 34 934054181; fax: þ 34 934011825. E-mail address: [email protected] (I. Arias). URL: http://www.lacan.upc.edu/arias/ (I. Arias).

1

the complexity of fracture in ferroelectric materials. Furthermore, it is widely accepted that the crack face boundary conditions strongly affect the electromechanical ﬁelds and thus play an important role in these interactions.

To gain theoretical understanding of fracture in these materials, researchers have resorted to the linear theory of piezoelectricity, where microstructure effects are not taken into account. Due to their simplicity, these models are useful to study the basic concepts of the linear theory in the context of fracture mechanics and to evaluate the effects of individual and coupled electromechanical ﬁelds and different crack face boundary conditions (McMeeking, 1999, 2004; Landis, 2004a; Li et al., 2008). Such approaches have allowed researchers to identify the effect of remanent polarization on energy release rates in poled ferroelectrics (Haug and McMeeking, 2006). However, the intrinsic nature of most piezoelectric materials demands the consideration of the nonlinear microstructure effects. These include models inspired in plasticity theory and aimed at polycrystalline ferroelectric ceramics (McMeeking and Landis, 2002). A phenomenological constitutive theory has been introduced to describe implicitly the domain formation around the cracks (Landis, 2003; Wang and Landis, 2006; Sheng and Landis, 2007). Other models rely on an energy based switching criterion (Hwang et al., 1995), considering the local phase transformations near the crack tip under the assumption of small scale switching (Zhu and Yang, 1999, 1997; Yang and Zhu, 1998; Beom and Atlurib, 2003). Recently, phase ﬁeld or time dependent Devonshire Ginzburg Landau (TDGL) models are gaining a growing attention. These models explicitly describe the formation and evolution of individual ferroelectric domains in the framework of continuum mechanics (Shu and Bhattacharya, 2001; Zhang and Bhattacharya, 2005; Schrade et al., 2007; Su and Landis, 2007; Xu et al., 2009), and have been extended to accurately account for the stray ﬁelds (Yang and Dayal, 2011). Related models have also been developed in micromagnetics (DeSimone et al., 2001; DeSimone and James, 2002). The nucleation and growth of domains near crack tips have been studied using these microstructural models and the inﬂuence on the stress ﬁeld (Yang and Dayal, 2012) and the mechanical and electromechanical J integrals have been reported (Song et al., 2007; Wang and Zhang, 2007; Wang and Kamlah, 2009; Xu et al., 2010b; Li and Landis, 2011; Li and Kuna, 2012). For completeness, we mention that cohesive theories aimed at fracture in ferroelectric materials have also been proposed (Gao et al., 1997; Arias et al., 2006).

The above mentioned models of ferroelectric fracture consider stationary crack conﬁgurations. With this assumption it is not possible to capture the effects of the microstructure on the crack propagation behavior, specially when the crack tip ﬁeld interacts with obstacles such as defects or grain boundaries. Recently, a phase ﬁeld model (Miehe et al., 2010; Linder and Miehe, 2012) and a strong discontinuity approach (Linder et al., 2011) have been proposed to simulate propagating cracks in linear piezoelectric solids. Towards a more realistic approach for ferroelectric single crystals, we have introduced a phase ﬁeld model for the coupled microstructure and fracture evolution by tackling the full complexity of the phenomenon (Abdollahi and Arias, 2011a). With this model we have shown that the interaction of the microstructure and the crack leads to the slow fast crack propagation behavior observed in experiment (Fang et al., 2007). In another work, we have demonstrated that this interaction results in an anisotropic crack propagation in ferroelectric single crystals (Abdollahi and Arias, 2011b). We have recently extended the model to ferroelectric polycrystals (Abdollahi and Arias, 2012). We have considered so far, the simplest crack face boundary conditions in electromechanical fracture. A similar approach has also been proposed for the crack propagation and kinking in ferroelectrics, where the spontaneous polarization of the material is chosen as primary order parameter (Xu et al., 2010a). The phase ﬁeld models proposed by Abdollahi and Arias (2011a,b, 2012) are based on variational theories of ferroelectric ceramics (Shu and Bhattacharya, 2001; Zhang and Bhattacharya, 2005; Su and Landis, 2007; Schrade et al., 2007) and brittle fracture (Francfort and Marigo, 1998; Bourdin et al., 2000, 2008; Bourdin, 2007; Amor et al., 2009). The variational structure of these models allows us to naturally couple multiple physics. In addition, phase ﬁeld models are particularly interesting since a single partial differential equation governing the phase ﬁeld accomplishes at once (1) the tracking of the interfaces (cracks, domain walls) in a smeared way and (2) the modeling of the interfacial phenomena such as domain wall energies or crack face boundary conditions. This computational approach smears both the crack and the domain wall, and treats naturally crack and domain nucleation, crack branching, crack and domain wall merging, and interactions between multiple cracks and domains. This is in contrast with the sharp interface models of fracture such as cohesive zone models (Xu and Needleman, 1994; Camacho and Ortiz, 1996), extended ﬁnite element method (XFEM) (Moes et al., 1999), or the strong discontinuity approach (Oliver et al., 2002) and ferroelectric domain evolution (Loge and Suo, 1996; Li and Liu, 2004), which require the crack surfaces and domain walls to be tracked algorithmically. On the other hand, the peridynamic approach (Silling, 2000; Dayal and Bhattacharya, 2006) can be viewed as an alternative to the phase ﬁeld approach, smearing sharp discontinuities. However, the peridynamic theory is relatively new and it needs to be developed further to address electromechanical materials.

The ﬂexibility of phase ﬁeld models come at the expense of a high computational cost, since the width of the phase ﬁeld regularization of the domain wall and the crack must be resolved by the discretization. Another challenge, which is the main objective of this work, is to encode different sharp crack boundary conditions into the phase ﬁeld framework. We are particularly interested in different crack face boundary conditions because they have a strong effect on the fracture behavior of piezoelectrics and ferroelectrics, and ultimately on the reliability of the devices. For this purpose, we propose a general framework encompassing all the usual crack face boundary conditions in the context of phase ﬁeld models. The phase ﬁeld model of brittle fracture is viewed as a regularization of Grifﬁth’s sharp crack model. In this sense, the aim here is to converge to the corresponding sharp crack solutions of the speciﬁc fracture problem. We consider the different crack face boundary conditions proposed in the literature for sharp crack models in electromechanical materials and formulate them in the regularized phase ﬁeld model. Note that here the cracks are not boundaries of the computational domain but

2

rather features of the solution within the domain, and hence the different sharp crack conditions have to be modeled in the phase ﬁeld partial differential equations.

In the context of sharp crack models in electromechanical materials, the most common crack face boundary conditions in the literature can be classiﬁed as follows:

A. Uncoupled electrical/mechanical crack face conditions. Mechanical boundary conditions: These are mainly: (1) traction free crack faces and (2) cohesive zone models (Xu and Needleman, 1994; Camacho and Ortiz, 1996) introducing a traction separation law on the crack faces. Here, we consider only traction free crack boundary conditions. The encoding of the cohesive crack face conditions in the context of the phase ﬁeld models is beyond the scope of this paper. Electrical boundary conditions: These are mainly (1) permeable, (2) impermeable and (3) semi permeable crack models, each assuming different electrical properties of the crack gap. The permeable crack boundary conditions were the ﬁrst to be considered for electromechanical materials (Parton, 1976). These conditions assume that the crack does not exist electrically, i.e. crack faces are closed and electrical ﬁelds are not perturbed by the presence of the crack. In contrast, impermeable boundary conditions were proposed to deﬁne an open and electrically defective crack by assuming zero permittivity for the crack gap (Deeg, 1980). However, both the permeable and impermeable crack boundary conditions are not physically justiﬁable in many cases, since the effect of the medium ﬁlling the crack gap is neglected. As an improvement, semi permeable boundary conditions were introduced to treat the crack gap as a linear dielectric material with a ﬁnite permittivity (Hao and Shen, 1994). As a result, the electric charge can penetrate the crack gap. A physical inconsistency of the semi permeable conditions is that the stored electric charge in the crack gap induces a closing traction on the crack faces, which is not considered in these conditions (McMeeking, 2004). This theory has been shown to be variationally inconsistent. Furthermore, the electric ﬁeld inside the crack gap is limited by the dielectric strength of the crack gap ﬁlling medium. Above this level of electric ﬁeld, the medium experiences failure of its insulating properties, i.e. electrical breakdown.

B. Coupled electromechanical crack face conditions. The Energetically Consistent (EC) crack model was proposed to overcome the inconsistency of the semi permeable boundary conditions by considering not only the electric charge inside the crack gap, but also the corresponding induced closing traction on the crack faces (Landis, 2004a). In fact, the crack acts as a capacitor inside the material. EC conditions are believed to be the most physically realistic conditions on the crack faces in electromechanical materials.

C. Polarization boundary conditions. In contrast to piezoelectrics, the modeling of cracks in ferroelectric materials requires imposing boundary conditions on the polarization on the crack faces (Wang and Kamlah, 2010). There are two usual choices (Vendik and Zubko, 2000), namely (1) free polarization and (2) zero polarization crack conditions. The former is a homogeneous Neumann boundary condition for the polarization, leaving it unaffected by the presence of the crack, and hence dictated by the bulk material model. The latter is a homogeneous Dirichlet boundary condition for the polarization, thereby modeling an open crack ﬁlled with free space.

In the following, we ﬁrst summarize in Sections 2 and 3, the phase ﬁeld model of brittle fracture and the variational formulation of ferroelectric and piezoelectric materials, respectively. Next, in Section 4 we describe the formulation of the phase ﬁeld models for fracture in these materials. We focus on the phase ﬁeld formulation of the most relevant mechanical, electrical, and electromechanical crack face boundary conditions as described above. For each possible choice of crack face boundary conditions, a different phase ﬁeld formulation is proposed. The corresponding governing equations are presented in Section 4.4, along with a general solution procedure. The proposed models are exercised numerically in Section 5, where the effects of the different crack face boundary conditions, electromechanical loadings and crack gap media on the crack propagation are evaluated. Numerical simulation results are discussed and compared with experiments. The last section is the conclusion of this paper.

2. Phase-ﬁeld modeling of fracture in brittle materials

In the variational regularized formulation of Grifﬁth’s fracture theory ﬁrst proposed by Francfort and Marigo (1998), the

total energy of a body made of a brittle material and occupying a region O is obtained as the sum of the bulk and surface

energies (Bourdin et al., 2000):

Z

Z"

2

#

Ek½u,v ¼ Oðv2 þ ZkÞFðeðuÞÞ dO þ Gc O ð14kvÞ þk9rv92 dO, ð1Þ

where F is the elastic potential, u is the mechanical displacement, e is the strain deﬁned as eðuÞ ¼ 1=2ðru þruT Þ, and v is a

scalar phase ﬁeld describing a smooth transition in space between unbroken (v ¼1) and broken (v¼0) states of the

material. For a linear elastic material, FðeÞ ¼ 12e : C : e, where C is the forth order tensor of elastic moduli. By noting that v2

multiplies the elastic potential F in the bulk energy (ﬁrst integral), it is clear that the value v ¼0 effectively reduces the

stiffness of the material to zero in the fractured zone. Zk is a small residual stiffness to avoid the singularity of the bulk

energy in this zone. In the surface energy, Gc is the critical energy release rate or the surface energy in Grifﬁth’s theory

3

Fig. 1. Example of a smeared crack using the proﬁle in Eq. (2).

(Grifﬁth, 1921). k is a positive regularization constant, which regulates the size of the smeared fracture zone. It can be

shown mathematically, by way of G convergence, that when the regularization parameter k tends to zero, this regularized

theory converges to Grifﬁth’s brittle fracture model (Bourdin et al., 2008). In particular, the traction free conditions on the

crack faces of the sharp model are recovered in the limit of a vanishingly small regularization parameter. Working by

analogy, the different electrical and electromechanical conditions of the sharp crack model in piezoelectrics and

ferroelectrics are encoded into the phase ﬁeld framework in Section 4.

The crack propagation in this model results from the competition between the bulk and surface energy terms.

Deformation of an elastic body under load increases the elastic energy density F. When this value approaches a critical

value in a given region, it is energetically favorable for the system to decrease the value of v towards zero in that region in

order to release elastic energy. On the other hand, decreasing the value of v leads to an increase in the surface energy since

deviations from 1 are penalized. Furthermore, variations of v are also penalized in this second term, resulting in the

formation of smeared cracks whose width is governed by the regularization parameter k. By increasing the value of the

critical energy release rate Gc of the material, a higher value of the elastic energy density is required to nucleate or

propagate cracks. It can be shown that the second integral converges to the surface area of the crack when k tends to zero,

as expected in the sharp interface model.

The theory as outlined above may lead to crack healing, which makes it necessary to supplement it with an

irreversibility condition. Requiring the ﬁeld v to be a monotonically decreasing function of time is cumbersome, and in

practice, the ﬁeld v is frozen to 0 when and where it reaches a given small threshold g (Bourdin, 2007; Bourdin et al., 2008).

An illustration of the diffuse crack is presented in Fig. 1. This ﬁgure is obtained from a minimizer of Ek½u,v in two

dimensions. We show the proﬁle of v perpendicular to the crack in its wake. Denoting d(x) the distance to the crack of a

point x, it can be shown that this optimal proﬁle takes the form (Bourdin et al., 2008)

8

>< 0

if dðxÞ rak

vkðxÞ ¼ >: 1 exp

dðxÞ ak 2k

otherwise,

ð2Þ

where 2 ak indicates the width of the fully fractured region where v¼0. This width is given by the threshold g, and it is smaller than k, ak ¼ oðkÞ.

The total energy in Eq. (1) is quadratic and convex in v and u separately. Therefore, for a ﬁxed v or u, Ekðu,Þ and Ekð,vÞ can be efﬁciently minimized solving a linear system of equations, with the appropriate boundary conditions. As a

consequence, the numerical implementation of this theory is straightforward by means of an iterative algorithm minimizing

separately each ﬁeld in a staggered manner.

3. Variational formulation of electromechanical solids

The behavior and properties of electromechanical solids such as ferroelectric and piezoelectric materials can be deﬁned by a thermodynamical potential containing mechanical, electrical and electromechanical coupling energy terms. The form of this potential and number of parameters depend on the complexity of the material behavior. In ferroelectric ceramics, the polarization and strain state can evolve in a nonlinear fashion due to ferroelectric/ferroelastic domain switching. However, under the assumption of a small mechanical/electrical load, domain switching cannot occur and the material behavior can be expressed by the linear theory of piezoelectricity. In fact, a piezoelectric material can be viewed as a linear approximation of the ferroelectric material behavior near the spontaneous polarization and strain state.

Ignoring body loads and volume charges for simplicity, the total electromechanical enthalpy of a ferroelectric or piezoelectric body occupying a region O can be written in terms of the mechanical displacement u, the total polarization p

4

and the electric potential f, as

Z

Z

Z

H½u,p,f ¼ HðeðuÞ,p,rp,EðfÞÞ dO

t Á u dS þ of dS,

ð3Þ

O

GN,u

GN,f

where E is the electric ﬁeld deﬁned as E ¼ rf. t and o are the tractions and surface charge density, respectively. GN,u and GN,f are the parts of the boundary of the domain @O where mechanical and electrical Neumann boundary conditions are

applied.

To account for domain switching, we follow the phase ﬁeld model of ferroelectric domain evolution proposed by Zhang

and Bhattacharya (2005) and Su and Landis (2007). The electromechanical enthalpy density H of a ferroelectric material is formulated as

Hðe,p,rp,EÞ ¼ 1 e : C : e e : C : esðpÞ þUðrpÞþ wðpÞ E Á p e0 9E92,

ð4Þ

2

2

where the eigenstrain es depends on the polarization, U is the domain wall energy density penalizing sharp variations in the polarization, w is the phase separation or Landau potential, and e0 is the permittivity of free space. The ﬁrst term in Eq.

(4) is the elastic potential and the second term is the electromechanical coupling energy density. The combination of these

terms and energy function w is the total Landau Devonshire energy density penalizing deviations from the spontaneous

polarizations and strains of the material, hence introducing the anisotropy and nonlinearity of ferroelectric materials. The detailed formulation of these energy functions and the material constants are given in Appendix A.

The electromechanical enthalpy density in Eq. (4) contains all crystallographic and domain wall information of a ferroelectric material. However, only some of this information is required to model a piezoelectric material by assuming that the nonlinear switching of the polarization between the various crystallographic orientations does not occur. According to the linear theory of piezoelectricity, the electromechanical enthalpy density H of a piezoelectric material is written in terms of the strain and electric ﬁeld alone (Tieresten, 1969)

Hðe,

EÞ

¼

1 2

ðe

erÞ : C : ðe

er Þ

ðe

erÞ : eT Á E

E Á pr

12E Á KE,

ð5Þ

where e is the third order tensor of piezoelectric coupling constants, and eT is such that ½eT ðijÞk ¼ ½ekðijÞ, where the brackets

indicate the indices paired with strains. er is the remanent strain, pr is the remanent polarization, and K is the second

order dielectric tensor. The piezoelectric material constants can be obtained from the enthalpy density in Eq. (4) by

linearizing around the spontaneous polarization and strain state of the ferroelectric single crystal (Vo¨ lker et al., 2011).

The stresses and electric displacements are derived from the electromechanical enthalpy density H as

r ¼ @H , D ¼ @H :

ð6Þ

@e

@E

Then, given Eqs. (4) and (5), the stresses and electric displacements are obtained for the ferroelectric and piezoelectric models, respectively as

r ¼ C : ðe esÞ, D ¼ p þe0E,

ð7Þ

r ¼ C : ðe erÞ eT Á E, D ¼ pr þ KE þ e : ðe erÞ:

ð8Þ

The weak form of the mechanical and electrostatic equilibrium equations is then obtained as

Z

Z

Z

Z

0 ¼ dH½u,p,f; du,df ¼ sjideij dO DidEi dO

tidui dS þ odf dS:

ð9Þ

O

O

GN,u

GN,f

The polarization evolution in ferroelectric materials is generally obtained from a gradient ﬂow of the total electro

mechanical enthalpy with respect to the polarization (Zhang and Bhattacharya, 2005),

Z

Z

m p_ dp dO ¼ dH½u,p,f; dp ¼ @H dp dO,

ð10Þ

p

ii

O

O @pi i

where 1=mp 4 0 is the mobility of the process, solved together with Eq. (9). We take the remanent state of the piezoelectric material as the reference conﬁguration, therefore the remanent strain er

and the remanent polarization pr are set to zero in the following. The poling direction of the piezoelectric material is implicitly encoded in the constitutive equations. It is noteworthy that this formulation, when applied to fracture mechanics as in Section 4, introduces the implicit assumption that the remanent polarization is perfectly balanced on the crack faces. This assumption is standard in fracture mechanics studies in piezoelectric materials and, for sharp crack models, it can be easily generalized to arbitrary levels of charge separation on the crack faces (Haug and McMeeking, 2006). Considering this assumption, the detailed formulation of Eq. (8) in Voigt form and the piezoelectric material constants are given in Appendix B.

4. Phase-ﬁeld modeling of fracture in piezoelectric and ferroelectric materials

The speciﬁc coupling between the phase ﬁeld model of brittle fracture in Section 2 and the models of piezoelectric and ferroelectric materials in Section 3 depends on the particular choice of the crack face conditions discussed in Section 1.

5

Based on these conditions, the presence of the crack can affect different terms of the electromechanical enthalpy density H in Eqs. (4) and (5) through the ﬁeld v. Therefore for each possible choice of crack face boundary conditions a different formulation of the enthalpy density H is required. Sections 4.1 4.3 deal with the uncoupled mechanical/electrical, the coupled electromechanical, and the polarization boundary conditions, respectively. Using the proposed formulations, the governing equations of the phase ﬁeld model of crack propagation in piezoelectric and ferroelectric materials are discussed in Section 4.4, along with the proposed solution algorithm.

4.1. Uncoupled mechanical/electrical crack face conditions

Most electromechanical crack models assume uncoupled mechanical and electrical boundary conditions on the crack faces. Although the uncoupled conditions are not physically realistic, they are good approximations of the coupled electromechanical conditions in some situations.

4.1.1. Traction free conditions Mathematically, these conditions are stated in the context of sharp crack models as

r þ Á n ¼ rÀ Á n ¼ 0,

ð11Þ

where the superscripts þ and denote the top and bottom crack faces and n is the unit normal. In the phase ﬁeld model

of brittle fracture, the traction free conditions are encoded by multiplying the elastic energy density F by the jump set

function ðv2 þ ZkÞ (see Eq. (1)). Working by analogy, in the case of the proposed electromechanical models, we multiply by the jump set function the ﬁrst two terms of the enthalpy density H associated with the elastic strains e. Therefore, the

elastic and electromechanical coupling potentials are annihilated in the fractured zone (v ¼0).

4.1.2. Traction free, electrically impermeable crack model In the case of impermeable cracks, the crack faces are treated as charge free surfaces, i.e. the normal components of the

electric displacement vanish on both crack faces as

D þ Á n ¼ DÀ Á n ¼ 0:

ð12Þ

The impermeable conditions imply that the crack cannot sustain any electric displacement inside the fractured zone. In the

phase ﬁeld model, we multiply all terms of H, involving the electric ﬁeld E, by the jump set function ðv2 þ ZkÞ. Therefore

the electric displacements vanish in the fractured zone (v ¼0), similar to the stresses. We illustrate in Appendix C (see Fig.

C2) that in the limit of small regularization parameter k, this formulation exhibits traction free and impermeable solutions

on the crack faces, as expected in the sharp crack model.

4.1.3. Traction free, electrically permeable crack model The permeable crack boundary conditions imply that the electric potential and the normal component of the electric

displacement are continuous across the crack faces, i.e. in the context of sharp crack models

f þ ¼ fÀ and D þ Á n ¼ DÀ Á n:

ð13Þ

In contrast to the impermeable crack model, encoding the permeable conditions requires that all terms of the enthalpy density H, involving the electric ﬁeld E, remain unmodiﬁed. As a way of example, the enthalpy density H of the traction free, electrically permeable crack model for the piezoelectric material is

H ¼ ðv2 þ ZkÞð12 e : C : e e : eT Á EÞ 12E Á KE:

ð14Þ

Accordingly, the electric displacements are

D ¼ ðv2 þ ZkÞe : e þKE:

ð15Þ

Since the traction free conditions result in a discontinuous displacement ﬁeld, Eq. (15) seems to be introducing a

discontinuity in the electric displacement through the coupling term e : e, in contradiction with Eq. (13). Nevertheless, we

show in Appendix C (see Fig. C1) that the effect of this discontinuity in the normal direction becomes negligible when the

regularization parameter k is small enough.

4.2. Electromechanical crack face conditions: energetically consistent crack model

As described in Section 1, the semi permeable crack conditions were introduced to model intermediate situations between the two over simplistic extremes of a completely closed crack (permeable) and an open crack ﬁlled with a medium of negligible conductivity (impermeable). This model assumes that the crack gap ﬁlling medium can sustain a certain degree of electric ﬁeld, thereby inducing an electric displacement inside the crack. In the context of sharp crack models, semi permeable crack face boundary conditions are formulated as

D þ Á n ¼ DÀ Á n ¼ Dc ¼ e0Ec,

ð16Þ

6

where Dc is the induced electric displacement and Ec is the electric ﬁeld in the crack gap. It is obtained as

Ec ¼

f þ fÀ

unþ

¼ uÀn

Df , Dun

ð17Þ

where the superscripts ‘‘ þ ’’ and ‘‘ ’’ denote the values on each crack face, Df is the potential drop across the crack gap and Dun is the crack opening displacement. The main deﬁciency stems from the fact that the electric ﬁeld inside the crack gap induces also tractions on the crack faces, which are not considered by the semi permeable conditions. It has been demonstrated that this model is not energetically consistent (Li et al., 2008). To overcome this inconsistency, energetically consistent (EC) crack face boundary conditions have been proposed in the context of sharp crack models (Landis, 2004a), in which

t þ ¼ scn, tÀ ¼ scn,

ð18Þ

where t þ and tÀ are the tractions acting on the crack faces and sc is the effective stress within the crack gap. According to

Landis (2004a) and assuming an inﬁnite breakdown strength of free space, the effective stress can be expressed as

sc ¼ 12e0E2c :

ð19Þ

In fact, EC conditions assume that the crack behaves electrically as a capacitor, storing electrical charge between the crack

faces. The associated electrical enthalpy is given by Z

Hc ¼ HcDun dS,

Sc

ð20Þ

where Sc denotes the surface (curve in 2D) of the sharp crack and the electrical enthalpy density Hc is

Hc ¼ 12e0E2c :

ð21Þ

In proposing the electrical enthalpy in Eq. (20), it is assumed that any electric ﬁeld and deformation components within

the gap parallel to the crack are negligible in comparison to the normal components. In the following, we develop the

general formulation of this enthalpy in the context of the phase ﬁeld models without any assumption regarding the

direction of the electric and deformation ﬁelds. Even though the present theory is geometrically linear, we need to

introduce geometric nonlinearity in selected terms to account for the emergence of a crack gap medium as a result of the

crack opening displacement. In fact, the enthalpy density in Eq. (4) can be obtained from the linearization of a general

enthalpy density of ferroelectrics in ﬁnite deformation electromechanics (Xiao and Bhattacharya, 2008). If this

linearization does not apply to the electrical enthalpy of free space (last term in Eq. (4)), we can account for the nonlinear enthalpy of the crack gap as follows. First, we express the electric ﬁeld in the deformed conﬁguration, Eb, in terms of the

reference gradient of the electric potential, Eb ¼ FÀT rf, where F is the deformation gradient tensor. Then, the electrical

enthalpy density of free space per unit deformed volume of material can be written as

Hb ¼ e0 9Eb92:

ð22Þ

2

This enthalpy can be decomposed into the enthalpy of the intact and fracture zones using the phase ﬁeld v

Hb ¼ v2Hb þ ð1 v2ÞHb :

ð23Þ

Ignoring all geometric nonlinearity in the intact zone, the total electrical enthalpy of free space becomes

e0 Z 2 2 e0 Z

2 ÀT

2

H0 ¼ 2 v 9rf9 dO 2 ð1 v Þ9F rf9 J dO,

ð24Þ

O

|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄOﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}

Hc

where J ¼ det F. The second term represents the phase ﬁeld counterpart of Eq. (20) and we can interpret Eb as Ec in the

fracture zone. Fig. 2 illustrates an opening crack where the electric ﬁelds are presented both for the sharp crack across the

faces and for the diffuse crack inside a partition of the fractured zone.

In Appendix C, we present numerical evidence substantiating the adequacy of this phase ﬁeld version of the EC

conditions. To provide a simple heuristic argument to support it, we consider a straight smeared crack along the x1

direction. The corresponding proﬁles of the ﬁelds u and f along the x2 direction are

8 >><

u

þ

x2 4 ak,

u þ uÀ

u þ þ uÀ

uðx2Þ ¼ >>: 2ak x2 þ 2

uÀ

ak r x2 r ak, x o ak,

ð25Þ

2

8>>>< ff þþ fÀ f þ þ fÀ x2 4 ak,

fðx2Þ ¼ >>> 2ak x2 þ 2

ak r x2 r ak,

ð26Þ

: fÀ

x2 o ak,

7

Fig. 2. Conﬁguration of an opening crack where the electric ﬁelds Ec and Ec are presented for the sharp and diffuse crack models, respectively. The crack faces are indicated by the dashed lines and a portion of the fractured zone is presented by the solid lines.

while the proﬁle of the ﬁeld v is given in Eq. (2) with dðxÞ ¼ 9x29. According to the second term of Eq. (24), the electrical

enthalpy of the smeared crack can then be written as

Hc ¼

e0 Z Z ak 9FÀT rf92J dx2 dS ¼

2 Sc Àak

e0 Z 2ak9FÀT rf92J dS ¼

2 Sc

e0 Z ðf þ fÀÞ2 2 S 2ak þ u þ uÀ dS,

c

ð27Þ

where we have used

21

03

FÀT ¼ 664 0

1 uþ

uÀ 775,

1 þ 2ak

2

3

0

rf ¼ 64 fþ fÀ 75,

2ak

u þ uÀ

J ¼ 1 þ 2ak :

When k tends to zero, the width ak also tends to zero and Eq. (27) converges to the enthalpy of the sharp crack in Eq. (20).

Of course, in practical computations the width of the smeared crack is not negligible and there is also a smooth transition of the ﬁeld v across the fracture zone. However, we illustrate in Appendix C (see Fig. C3) that for small enough values of the

regularization parameter k, the total electrical enthalpy of the diffuse crack accurately follows the sharp crack model.

The ﬁrst term in Eq. (24) is the enthalpy term of the impermeable crack model. Therefore this equation indicates that the enthalpy density H of the EC crack model is obtained by adding the electrical enthalpy density of the diffuse crack to that of the impermeable crack model, which allows the displacement and electric potential ﬁelds to be discontinuous. Given the enthalpy density in the second term of Eqs. (24) and (6), the stresses and electric displacements are obtained in the fracture zone as

r ¼ J ð1 v2Þrc,

ð28Þ

2

D ¼ Jð1 v2ÞFÀ1Dc, with

Dc ¼ e0Ec ¼ e0FÀT E,

ð29Þ ð30Þ

rc ¼ FÀ1ðDc EcÞ þ ðDc EcÞFÀT

e0 9Ec92ðFÀ1 þ FÀT Þ:

2

ð31Þ

Note that the lower permittivity of the free space as compared to that of the material induces a high electric ﬁeld inside

the fracture zone, especially in a specimen under applied electrical loading. The free space can sustain a certain level of

electric ﬁeld, called dielectric strength, above which the free space experiences failure of its insulating properties, i.e.

electrical breakdown. This effect can be modeled mathematically as a bound constraint on the magnitude of the electric

ﬁeld inside the fracture gap

( 9FÀT rf9

9Ec9 ¼ Ed

if 9Ec9 r Ed, if 9Ec9 4 Ed,

ð32Þ

where Ed is the dielectric strength or the electrical discharge level of the free space. This constraint can be added to the governing equations. The corresponding Lagrange multiplier modiﬁes the electric displacement inside the fracture zone. See Appendix D for details. We evaluate the effect of the electrical breakdown constraint in Section 5.

8

In the above formulation, the free space is assumed to be the medium ﬁlling the crack gap. Other crack gap ﬁlling media

can be considered as well by replacing the free space permittivity e0 with the corresponding value em in Hc in Eq. (24).

4.3. Polarization boundary conditions

In addition to the crack face conditions discussed in the previous sections, crack face conditions on the polariza tion need to be supplied for ferroelectric materials. We discuss next the two usual choices and their phase ﬁeld implementation.

4.3.1. Free polarization crack face conditions The free polarization boundary conditions assume that the polarization distribution is unaffected by the presence of

the crack, and hence dictated by the bulk material model. Mathematically, they are written as

rp þ Á n ¼ rpÀ Á n ¼ 0:

ð33Þ

This condition can be encoded in the phase ﬁeld model by multiplying with the jump set function ðv2 þ ZkÞ the energy

terms involving the gradient of polarization, i.e. the domain wall energy density U.

4.3.2. Zero polarization crack face conditions In contrast to the bulk material, the polarization vanishes in the free space ﬁlling the crack gap. The zero polarization

boundary conditions assume that the polarization is continuous at the material free space interface, i.e.

p þ ¼ pÀ ¼ 0:

ð34Þ

This condition can also be formulated in a similar way, by multiplying with the jump set function the energy terms associated with the polarization. In Appendix C (see Fig. C4), we show that these methods indeed produce numerical solutions satisfying the free and zero polarization boundary conditions in a diffuse sense.

4.4. Phase ﬁeld model of crack propagation in piezoelectric and ferroelectric materials

Based on the discussion in Sections 4.1 4.3, the different electromechanical crack face boundary conditions along with

the corresponding electromechanical enthalpy density, stresses and electric displacements for piezoelectric and ferro

electric materials are summarized in Appendix E. We now can form the total electromechanical enthalpy of a possibly

fractured piezoelectric or ferroelectric body by adding the surface energy in Eq. (1) to the total electromechanical enthalpy

in Eq. (3) as

Z"

2

#

ð1 vÞ þ k9rv92 dO,

ð35Þ

Hs½u,p,f,v ¼ H½u,p,f,v þ Gc

O

4k

where H now depends also on the ﬁeld v through the modiﬁed electromechanical enthalpy density H for different crack

face boundary conditions.

To capture interactions between the fracture and the microstructure processes in ferroelectrics, the crack propagation

should not be overwhelmingly faster than the microstructure evolution, and vice versa. In practice, the relative kinetics of

the microstructure evolution and the crack propagation gives the two phenomena a chance to interact. In the absence of

detailed experimental or theoretical information on this kinetics, v is selected together with the polarization as primary

order parameters and ﬁnite mobilities are introduced for the fracture and microstructure processes. Then, the time

evolution of the system results from the gradient ﬂow of the primary order parameters, assuming that the displacement

and the electric ﬁeld adjust immediately to mechanical and electrostatic equilibrium. The gradient ﬂow of the polarization

and the weak form of the mechanical and electrostatic equilibrium have already been introduced in Eqs. (10) and (9),

respectively. Here, the gradient ﬂow of the ﬁeld v as an additional governing equation is obtained as

Z

Z

Z

@H dv dO Gc v 1 dv þ kv,idv,i dO,

ð36Þ

mv v_ dv dO ¼ dHs½u,p,f,v; dv ¼ O

O @v

O 4k

where 1=mv is the mobility of the fracture process. In fact this equation is a Ginzburg Landau type evolution equation,

where the left hand side term is the dissipation in the process zone and the mobility constant controls the rate of the energy dissipation (Hakim and Karma, 2009). A recent study, using a phase ﬁeld model of brittle fracture, indicated that the mobility parameter controls the crack velocity, particularly at the initial stages of an unstable crack propagation (Kuhn and Mu¨ ller, 2010). Therefore the proper selection of this parameter allows us to track crack propagation, even in mechanically unstable conditions.

The ﬁnite element equations are obtained from the weak forms in Eqs. (10), (9), and (36) with the enthalpy density H,

the stresses r, and the electric displacements D for each crack face condition given in Appendix E. Eqs. (10) and (36) are discretized in time with an implicit scheme from time tmÀ1 to tm ¼ tmÀ1 þ Dt. The Newton Raphson method is implemented for the implicit time integration of Eq. (10) due to the nonlinear terms of the phase separation energy (w).

Algorithm 1 presents a simple procedure to solve forward in time the coupled system in a staggered, iterative way. The

9

function g encodes the data for the applied electromechanical load as a function of the time step. Since the crack should

not be allowed to heal (irreversibility condition), the nodes reaching a value of v below a certain threshold g, are assigned a

ﬁxed value of v ¼0 for the rest of the calculation. Note that the piezoelectric models do not require the polarization data and the polarization evolution in step 6 of the algorithm is only computed for ferroelectrics. In this step the constraint for the polarization evolution is only used for the zero polarization conditions. In fact, with this constraint the mobility term of the gradient ﬂow is ﬁxed to zero in the fractured zone (v ¼0).

Algorithm 1. For crack propagation in piezoelectric and ferroelectric materials.

1:

Let m 0 and t0 0

2:

Set v0 vinit, p0 pinit , f0 0 and u0 0

3:

repeat

4:

m’mþ 1

5:

tm’tm 1 þ Dt

6:

Compute pm in Eq. (10) using pm 1, um 1, fm 1 and vm 1 under the constraint pm 0 for vm 1 0

7:

Compute um and fm in Eq. (9) using pm and vm 1 under the electromechanical load gðtmÞ

8:

Compute vm in Eq. (36) using pm, um, fm, and vm 1 under the constraint vm 0 for vm 1 r g and Hv 0 for Hv o 0

9:

until m n

We note that the EC conditions discussed in Section 4.2 introduce non linear terms into the ﬁnite element equations. Here,

we implement an internal loop in step 7 of Algorithm 1 to solve this non linearity. In each iteration of this loop, the non

linear terms are updated with the last values of the mechanical displacement and the electric potential, then Eq. (9) is solved.

This iteration continues until a steady state is reached for both the mechanical displacement and the electric potential.

An important aspect of electromechanical crack propagation deserves special attention. Generically, the enthalpy

functional can be written as

Z

Z"

2

#

O½ðv2 þ ZkÞHvðu,p,fÞ þ Hrestðu,p,fÞ dO þ Gc O ð14kvÞ þ k9rv92 dO: ð37Þ

Under high applied electrical load, Hv may become negative. Then, when minimizing the energy with respect to v, it may be favorable to localize in narrow regions high values of v, much above one. Such anomalous localization zones are not physically meaningful. In sharp crack models of electromechanical fracture, this issue manifests itself with the negative energy release rates found at high applied electric ﬁelds (Li et al., 2008). Numerically, we deal with this issue by limiting the maximum value of v to one, or alternatively, by setting negative values of Hv to zero in step 8 of Algorithm 1.

5. Numerical simulations and discussion

In this section, we examine the effects of the different crack face conditions on the crack propagation. The simulations of a propagating crack in a poled piezoelectric material under different applied electromechanical loadings are presented in Section 5.1. We also exercise the model with different crack gap ﬁlling dielectric ﬂuids and introducing the constraint of electrical breakdown as discussed in Section 4.2. Section 5.2 presents the results of crack propagation in a poled ferroelectric single crystal where the crack interacts with domain switching. The results of each section are discussed in detail.

5.1. Propagating cracks in piezoelectrics

To investigate the effects of the different crack face conditions on the crack propagation, a set of simulations is performed following Algorithm 1. A four point bending setup is considered. Fig. 3 presents the dimensions of the specimen along with the boundary conditions and the pre crack. Plane strain conditions and the material properties of PZT PIC 151 poled along the x2 axis (normal to the pre crack) are assumed (see Appendix B for the material parameters). The value of a is the length of the region with v ¼0 (i.e. the crack length). The pre crack with a length of a¼ 0.5 mm is introduced in the model using the proﬁle in Eq. (2). The model is discretized with a ﬁne mesh, with the smallest element size of h ¼10 4 mm in the vicinity of the pre crack. An adaptive mesh reﬁnement algorithm is employed to generate the mesh. A typical mesh

is presented in Fig. 4. The regularization parameter k is selected as four times the smallest element size, i.e. k ¼ 4 Â 10À4 mm, and the residual stiffness is set to Z ¼ 10À6. The inverse mobility of the fracture evolution is set to mv ¼ 0:1 Ns=m2. The point load P is applied incrementally using the load function P(t)¼At with a rate of A ¼100 N/s and a time step of Dt ¼ 10À2 s. A total number of 2.6 Â 104 time steps are performed for each simulation. The crack initiates and

becomes unstable when the load reaches P¼14 KN, at time t¼ 140 s. From this time on, the point load is ﬁxed and the crack length relative to the initial length (Da) is recorded at regular intervals, see Fig. 5 for the different crack face conditions. The ﬁnite mobility of the v ﬁeld allows us to follow a slow crack propagation in this period. During the crack propagation, external electric ﬁelds E¼ 71 MV/m (see Fig. 3 for the sign convention) are applied from time t ¼190 s. It is apparent that the impermeable crack is completely arrested by the application of both the positive and negative electrical

10

Amir Abdollahi, Irene Arias n

Laboratori de Calcul Numeric (LaCa N), Departament de Matematica Aplicada III, Universitat Politecnica de Catalunya (UPC), Campus Nord UPC-C2, E-08034 Barcelona, Spain

abstract

We present a family of phase ﬁeld models for fracture in piezoelectric and ferroelectric materials. These models couple a variational formulation of brittle fracture with, respectively, (1) the linear theory of piezoelectricity, and (2) a Ginzburg Landau model of the ferroelectric microstructure to address the full complexity of the fracture phenomenon in these materials. In these models, both the cracks and the ferroelectric domain walls are represented in a diffuse way by phase ﬁelds. The main challenge addressed here is encoding various electromechanical crack models (introduced as crack face boundary conditions in sharp models) into the phase ﬁeld framework. The proposed models are veriﬁed through comparisons with the corresponding sharp crack models. We also perform two dimensional ﬁnite element simulations to demonstrate the effect of the different crack face conditions, the electromechanical loading and the media ﬁlling the crack gap on the crack propagation and the microstructure evolution. Salient features of the results are compared with experiments.

1. Introduction

The design and implementation of electromechanical systems demand high performance materials regarding their coupling behavior and reliability. Piezoelectric ceramics are very prominent in this ﬁeld, exhibiting strong electromecha nical coupling with short response times. However, their inherent brittleness is a serious obstacle to their reliable operation in devices, which demands a deep understanding of the fracture behavior. Towards this goal, numerous theoretical and experimental investigations have been carried out during the past decades on the fracture of piezo ceramics. Comprehensive reviews of these works have been presented (Zhang and Gao, 2004; Schneider, 2007; Kuna, 2010). Most piezoelectrics also exhibit ferroelectric and ferroelastic switching behavior with macroscopic dielectric and butterﬂy hysteresis. For related modeling approaches, see the reviews by Kamlah (2001), Landis (2004b) and Huber (2005). Microstructural domains often nucleate and evolve under high electromechanical loadings and near load concentrations such as those produced at crack tips (Hackemann and Pfeiffer, 2003; Jones et al., 2007; Pojprapai et al., 2011). The interactions between the microstructure and the localized stress and electric ﬁelds near the crack tips are responsible for

n Corresponding author. Tel.: þ 34 934054181; fax: þ 34 934011825. E-mail address: [email protected] (I. Arias). URL: http://www.lacan.upc.edu/arias/ (I. Arias).

1

the complexity of fracture in ferroelectric materials. Furthermore, it is widely accepted that the crack face boundary conditions strongly affect the electromechanical ﬁelds and thus play an important role in these interactions.

To gain theoretical understanding of fracture in these materials, researchers have resorted to the linear theory of piezoelectricity, where microstructure effects are not taken into account. Due to their simplicity, these models are useful to study the basic concepts of the linear theory in the context of fracture mechanics and to evaluate the effects of individual and coupled electromechanical ﬁelds and different crack face boundary conditions (McMeeking, 1999, 2004; Landis, 2004a; Li et al., 2008). Such approaches have allowed researchers to identify the effect of remanent polarization on energy release rates in poled ferroelectrics (Haug and McMeeking, 2006). However, the intrinsic nature of most piezoelectric materials demands the consideration of the nonlinear microstructure effects. These include models inspired in plasticity theory and aimed at polycrystalline ferroelectric ceramics (McMeeking and Landis, 2002). A phenomenological constitutive theory has been introduced to describe implicitly the domain formation around the cracks (Landis, 2003; Wang and Landis, 2006; Sheng and Landis, 2007). Other models rely on an energy based switching criterion (Hwang et al., 1995), considering the local phase transformations near the crack tip under the assumption of small scale switching (Zhu and Yang, 1999, 1997; Yang and Zhu, 1998; Beom and Atlurib, 2003). Recently, phase ﬁeld or time dependent Devonshire Ginzburg Landau (TDGL) models are gaining a growing attention. These models explicitly describe the formation and evolution of individual ferroelectric domains in the framework of continuum mechanics (Shu and Bhattacharya, 2001; Zhang and Bhattacharya, 2005; Schrade et al., 2007; Su and Landis, 2007; Xu et al., 2009), and have been extended to accurately account for the stray ﬁelds (Yang and Dayal, 2011). Related models have also been developed in micromagnetics (DeSimone et al., 2001; DeSimone and James, 2002). The nucleation and growth of domains near crack tips have been studied using these microstructural models and the inﬂuence on the stress ﬁeld (Yang and Dayal, 2012) and the mechanical and electromechanical J integrals have been reported (Song et al., 2007; Wang and Zhang, 2007; Wang and Kamlah, 2009; Xu et al., 2010b; Li and Landis, 2011; Li and Kuna, 2012). For completeness, we mention that cohesive theories aimed at fracture in ferroelectric materials have also been proposed (Gao et al., 1997; Arias et al., 2006).

The above mentioned models of ferroelectric fracture consider stationary crack conﬁgurations. With this assumption it is not possible to capture the effects of the microstructure on the crack propagation behavior, specially when the crack tip ﬁeld interacts with obstacles such as defects or grain boundaries. Recently, a phase ﬁeld model (Miehe et al., 2010; Linder and Miehe, 2012) and a strong discontinuity approach (Linder et al., 2011) have been proposed to simulate propagating cracks in linear piezoelectric solids. Towards a more realistic approach for ferroelectric single crystals, we have introduced a phase ﬁeld model for the coupled microstructure and fracture evolution by tackling the full complexity of the phenomenon (Abdollahi and Arias, 2011a). With this model we have shown that the interaction of the microstructure and the crack leads to the slow fast crack propagation behavior observed in experiment (Fang et al., 2007). In another work, we have demonstrated that this interaction results in an anisotropic crack propagation in ferroelectric single crystals (Abdollahi and Arias, 2011b). We have recently extended the model to ferroelectric polycrystals (Abdollahi and Arias, 2012). We have considered so far, the simplest crack face boundary conditions in electromechanical fracture. A similar approach has also been proposed for the crack propagation and kinking in ferroelectrics, where the spontaneous polarization of the material is chosen as primary order parameter (Xu et al., 2010a). The phase ﬁeld models proposed by Abdollahi and Arias (2011a,b, 2012) are based on variational theories of ferroelectric ceramics (Shu and Bhattacharya, 2001; Zhang and Bhattacharya, 2005; Su and Landis, 2007; Schrade et al., 2007) and brittle fracture (Francfort and Marigo, 1998; Bourdin et al., 2000, 2008; Bourdin, 2007; Amor et al., 2009). The variational structure of these models allows us to naturally couple multiple physics. In addition, phase ﬁeld models are particularly interesting since a single partial differential equation governing the phase ﬁeld accomplishes at once (1) the tracking of the interfaces (cracks, domain walls) in a smeared way and (2) the modeling of the interfacial phenomena such as domain wall energies or crack face boundary conditions. This computational approach smears both the crack and the domain wall, and treats naturally crack and domain nucleation, crack branching, crack and domain wall merging, and interactions between multiple cracks and domains. This is in contrast with the sharp interface models of fracture such as cohesive zone models (Xu and Needleman, 1994; Camacho and Ortiz, 1996), extended ﬁnite element method (XFEM) (Moes et al., 1999), or the strong discontinuity approach (Oliver et al., 2002) and ferroelectric domain evolution (Loge and Suo, 1996; Li and Liu, 2004), which require the crack surfaces and domain walls to be tracked algorithmically. On the other hand, the peridynamic approach (Silling, 2000; Dayal and Bhattacharya, 2006) can be viewed as an alternative to the phase ﬁeld approach, smearing sharp discontinuities. However, the peridynamic theory is relatively new and it needs to be developed further to address electromechanical materials.

The ﬂexibility of phase ﬁeld models come at the expense of a high computational cost, since the width of the phase ﬁeld regularization of the domain wall and the crack must be resolved by the discretization. Another challenge, which is the main objective of this work, is to encode different sharp crack boundary conditions into the phase ﬁeld framework. We are particularly interested in different crack face boundary conditions because they have a strong effect on the fracture behavior of piezoelectrics and ferroelectrics, and ultimately on the reliability of the devices. For this purpose, we propose a general framework encompassing all the usual crack face boundary conditions in the context of phase ﬁeld models. The phase ﬁeld model of brittle fracture is viewed as a regularization of Grifﬁth’s sharp crack model. In this sense, the aim here is to converge to the corresponding sharp crack solutions of the speciﬁc fracture problem. We consider the different crack face boundary conditions proposed in the literature for sharp crack models in electromechanical materials and formulate them in the regularized phase ﬁeld model. Note that here the cracks are not boundaries of the computational domain but

2

rather features of the solution within the domain, and hence the different sharp crack conditions have to be modeled in the phase ﬁeld partial differential equations.

In the context of sharp crack models in electromechanical materials, the most common crack face boundary conditions in the literature can be classiﬁed as follows:

A. Uncoupled electrical/mechanical crack face conditions. Mechanical boundary conditions: These are mainly: (1) traction free crack faces and (2) cohesive zone models (Xu and Needleman, 1994; Camacho and Ortiz, 1996) introducing a traction separation law on the crack faces. Here, we consider only traction free crack boundary conditions. The encoding of the cohesive crack face conditions in the context of the phase ﬁeld models is beyond the scope of this paper. Electrical boundary conditions: These are mainly (1) permeable, (2) impermeable and (3) semi permeable crack models, each assuming different electrical properties of the crack gap. The permeable crack boundary conditions were the ﬁrst to be considered for electromechanical materials (Parton, 1976). These conditions assume that the crack does not exist electrically, i.e. crack faces are closed and electrical ﬁelds are not perturbed by the presence of the crack. In contrast, impermeable boundary conditions were proposed to deﬁne an open and electrically defective crack by assuming zero permittivity for the crack gap (Deeg, 1980). However, both the permeable and impermeable crack boundary conditions are not physically justiﬁable in many cases, since the effect of the medium ﬁlling the crack gap is neglected. As an improvement, semi permeable boundary conditions were introduced to treat the crack gap as a linear dielectric material with a ﬁnite permittivity (Hao and Shen, 1994). As a result, the electric charge can penetrate the crack gap. A physical inconsistency of the semi permeable conditions is that the stored electric charge in the crack gap induces a closing traction on the crack faces, which is not considered in these conditions (McMeeking, 2004). This theory has been shown to be variationally inconsistent. Furthermore, the electric ﬁeld inside the crack gap is limited by the dielectric strength of the crack gap ﬁlling medium. Above this level of electric ﬁeld, the medium experiences failure of its insulating properties, i.e. electrical breakdown.

B. Coupled electromechanical crack face conditions. The Energetically Consistent (EC) crack model was proposed to overcome the inconsistency of the semi permeable boundary conditions by considering not only the electric charge inside the crack gap, but also the corresponding induced closing traction on the crack faces (Landis, 2004a). In fact, the crack acts as a capacitor inside the material. EC conditions are believed to be the most physically realistic conditions on the crack faces in electromechanical materials.

C. Polarization boundary conditions. In contrast to piezoelectrics, the modeling of cracks in ferroelectric materials requires imposing boundary conditions on the polarization on the crack faces (Wang and Kamlah, 2010). There are two usual choices (Vendik and Zubko, 2000), namely (1) free polarization and (2) zero polarization crack conditions. The former is a homogeneous Neumann boundary condition for the polarization, leaving it unaffected by the presence of the crack, and hence dictated by the bulk material model. The latter is a homogeneous Dirichlet boundary condition for the polarization, thereby modeling an open crack ﬁlled with free space.

In the following, we ﬁrst summarize in Sections 2 and 3, the phase ﬁeld model of brittle fracture and the variational formulation of ferroelectric and piezoelectric materials, respectively. Next, in Section 4 we describe the formulation of the phase ﬁeld models for fracture in these materials. We focus on the phase ﬁeld formulation of the most relevant mechanical, electrical, and electromechanical crack face boundary conditions as described above. For each possible choice of crack face boundary conditions, a different phase ﬁeld formulation is proposed. The corresponding governing equations are presented in Section 4.4, along with a general solution procedure. The proposed models are exercised numerically in Section 5, where the effects of the different crack face boundary conditions, electromechanical loadings and crack gap media on the crack propagation are evaluated. Numerical simulation results are discussed and compared with experiments. The last section is the conclusion of this paper.

2. Phase-ﬁeld modeling of fracture in brittle materials

In the variational regularized formulation of Grifﬁth’s fracture theory ﬁrst proposed by Francfort and Marigo (1998), the

total energy of a body made of a brittle material and occupying a region O is obtained as the sum of the bulk and surface

energies (Bourdin et al., 2000):

Z

Z"

2

#

Ek½u,v ¼ Oðv2 þ ZkÞFðeðuÞÞ dO þ Gc O ð14kvÞ þk9rv92 dO, ð1Þ

where F is the elastic potential, u is the mechanical displacement, e is the strain deﬁned as eðuÞ ¼ 1=2ðru þruT Þ, and v is a

scalar phase ﬁeld describing a smooth transition in space between unbroken (v ¼1) and broken (v¼0) states of the

material. For a linear elastic material, FðeÞ ¼ 12e : C : e, where C is the forth order tensor of elastic moduli. By noting that v2

multiplies the elastic potential F in the bulk energy (ﬁrst integral), it is clear that the value v ¼0 effectively reduces the

stiffness of the material to zero in the fractured zone. Zk is a small residual stiffness to avoid the singularity of the bulk

energy in this zone. In the surface energy, Gc is the critical energy release rate or the surface energy in Grifﬁth’s theory

3

Fig. 1. Example of a smeared crack using the proﬁle in Eq. (2).

(Grifﬁth, 1921). k is a positive regularization constant, which regulates the size of the smeared fracture zone. It can be

shown mathematically, by way of G convergence, that when the regularization parameter k tends to zero, this regularized

theory converges to Grifﬁth’s brittle fracture model (Bourdin et al., 2008). In particular, the traction free conditions on the

crack faces of the sharp model are recovered in the limit of a vanishingly small regularization parameter. Working by

analogy, the different electrical and electromechanical conditions of the sharp crack model in piezoelectrics and

ferroelectrics are encoded into the phase ﬁeld framework in Section 4.

The crack propagation in this model results from the competition between the bulk and surface energy terms.

Deformation of an elastic body under load increases the elastic energy density F. When this value approaches a critical

value in a given region, it is energetically favorable for the system to decrease the value of v towards zero in that region in

order to release elastic energy. On the other hand, decreasing the value of v leads to an increase in the surface energy since

deviations from 1 are penalized. Furthermore, variations of v are also penalized in this second term, resulting in the

formation of smeared cracks whose width is governed by the regularization parameter k. By increasing the value of the

critical energy release rate Gc of the material, a higher value of the elastic energy density is required to nucleate or

propagate cracks. It can be shown that the second integral converges to the surface area of the crack when k tends to zero,

as expected in the sharp interface model.

The theory as outlined above may lead to crack healing, which makes it necessary to supplement it with an

irreversibility condition. Requiring the ﬁeld v to be a monotonically decreasing function of time is cumbersome, and in

practice, the ﬁeld v is frozen to 0 when and where it reaches a given small threshold g (Bourdin, 2007; Bourdin et al., 2008).

An illustration of the diffuse crack is presented in Fig. 1. This ﬁgure is obtained from a minimizer of Ek½u,v in two

dimensions. We show the proﬁle of v perpendicular to the crack in its wake. Denoting d(x) the distance to the crack of a

point x, it can be shown that this optimal proﬁle takes the form (Bourdin et al., 2008)

8

>< 0

if dðxÞ rak

vkðxÞ ¼ >: 1 exp

dðxÞ ak 2k

otherwise,

ð2Þ

where 2 ak indicates the width of the fully fractured region where v¼0. This width is given by the threshold g, and it is smaller than k, ak ¼ oðkÞ.

The total energy in Eq. (1) is quadratic and convex in v and u separately. Therefore, for a ﬁxed v or u, Ekðu,Þ and Ekð,vÞ can be efﬁciently minimized solving a linear system of equations, with the appropriate boundary conditions. As a

consequence, the numerical implementation of this theory is straightforward by means of an iterative algorithm minimizing

separately each ﬁeld in a staggered manner.

3. Variational formulation of electromechanical solids

The behavior and properties of electromechanical solids such as ferroelectric and piezoelectric materials can be deﬁned by a thermodynamical potential containing mechanical, electrical and electromechanical coupling energy terms. The form of this potential and number of parameters depend on the complexity of the material behavior. In ferroelectric ceramics, the polarization and strain state can evolve in a nonlinear fashion due to ferroelectric/ferroelastic domain switching. However, under the assumption of a small mechanical/electrical load, domain switching cannot occur and the material behavior can be expressed by the linear theory of piezoelectricity. In fact, a piezoelectric material can be viewed as a linear approximation of the ferroelectric material behavior near the spontaneous polarization and strain state.

Ignoring body loads and volume charges for simplicity, the total electromechanical enthalpy of a ferroelectric or piezoelectric body occupying a region O can be written in terms of the mechanical displacement u, the total polarization p

4

and the electric potential f, as

Z

Z

Z

H½u,p,f ¼ HðeðuÞ,p,rp,EðfÞÞ dO

t Á u dS þ of dS,

ð3Þ

O

GN,u

GN,f

where E is the electric ﬁeld deﬁned as E ¼ rf. t and o are the tractions and surface charge density, respectively. GN,u and GN,f are the parts of the boundary of the domain @O where mechanical and electrical Neumann boundary conditions are

applied.

To account for domain switching, we follow the phase ﬁeld model of ferroelectric domain evolution proposed by Zhang

and Bhattacharya (2005) and Su and Landis (2007). The electromechanical enthalpy density H of a ferroelectric material is formulated as

Hðe,p,rp,EÞ ¼ 1 e : C : e e : C : esðpÞ þUðrpÞþ wðpÞ E Á p e0 9E92,

ð4Þ

2

2

where the eigenstrain es depends on the polarization, U is the domain wall energy density penalizing sharp variations in the polarization, w is the phase separation or Landau potential, and e0 is the permittivity of free space. The ﬁrst term in Eq.

(4) is the elastic potential and the second term is the electromechanical coupling energy density. The combination of these

terms and energy function w is the total Landau Devonshire energy density penalizing deviations from the spontaneous

polarizations and strains of the material, hence introducing the anisotropy and nonlinearity of ferroelectric materials. The detailed formulation of these energy functions and the material constants are given in Appendix A.

The electromechanical enthalpy density in Eq. (4) contains all crystallographic and domain wall information of a ferroelectric material. However, only some of this information is required to model a piezoelectric material by assuming that the nonlinear switching of the polarization between the various crystallographic orientations does not occur. According to the linear theory of piezoelectricity, the electromechanical enthalpy density H of a piezoelectric material is written in terms of the strain and electric ﬁeld alone (Tieresten, 1969)

Hðe,

EÞ

¼

1 2

ðe

erÞ : C : ðe

er Þ

ðe

erÞ : eT Á E

E Á pr

12E Á KE,

ð5Þ

where e is the third order tensor of piezoelectric coupling constants, and eT is such that ½eT ðijÞk ¼ ½ekðijÞ, where the brackets

indicate the indices paired with strains. er is the remanent strain, pr is the remanent polarization, and K is the second

order dielectric tensor. The piezoelectric material constants can be obtained from the enthalpy density in Eq. (4) by

linearizing around the spontaneous polarization and strain state of the ferroelectric single crystal (Vo¨ lker et al., 2011).

The stresses and electric displacements are derived from the electromechanical enthalpy density H as

r ¼ @H , D ¼ @H :

ð6Þ

@e

@E

Then, given Eqs. (4) and (5), the stresses and electric displacements are obtained for the ferroelectric and piezoelectric models, respectively as

r ¼ C : ðe esÞ, D ¼ p þe0E,

ð7Þ

r ¼ C : ðe erÞ eT Á E, D ¼ pr þ KE þ e : ðe erÞ:

ð8Þ

The weak form of the mechanical and electrostatic equilibrium equations is then obtained as

Z

Z

Z

Z

0 ¼ dH½u,p,f; du,df ¼ sjideij dO DidEi dO

tidui dS þ odf dS:

ð9Þ

O

O

GN,u

GN,f

The polarization evolution in ferroelectric materials is generally obtained from a gradient ﬂow of the total electro

mechanical enthalpy with respect to the polarization (Zhang and Bhattacharya, 2005),

Z

Z

m p_ dp dO ¼ dH½u,p,f; dp ¼ @H dp dO,

ð10Þ

p

ii

O

O @pi i

where 1=mp 4 0 is the mobility of the process, solved together with Eq. (9). We take the remanent state of the piezoelectric material as the reference conﬁguration, therefore the remanent strain er

and the remanent polarization pr are set to zero in the following. The poling direction of the piezoelectric material is implicitly encoded in the constitutive equations. It is noteworthy that this formulation, when applied to fracture mechanics as in Section 4, introduces the implicit assumption that the remanent polarization is perfectly balanced on the crack faces. This assumption is standard in fracture mechanics studies in piezoelectric materials and, for sharp crack models, it can be easily generalized to arbitrary levels of charge separation on the crack faces (Haug and McMeeking, 2006). Considering this assumption, the detailed formulation of Eq. (8) in Voigt form and the piezoelectric material constants are given in Appendix B.

4. Phase-ﬁeld modeling of fracture in piezoelectric and ferroelectric materials

The speciﬁc coupling between the phase ﬁeld model of brittle fracture in Section 2 and the models of piezoelectric and ferroelectric materials in Section 3 depends on the particular choice of the crack face conditions discussed in Section 1.

5

Based on these conditions, the presence of the crack can affect different terms of the electromechanical enthalpy density H in Eqs. (4) and (5) through the ﬁeld v. Therefore for each possible choice of crack face boundary conditions a different formulation of the enthalpy density H is required. Sections 4.1 4.3 deal with the uncoupled mechanical/electrical, the coupled electromechanical, and the polarization boundary conditions, respectively. Using the proposed formulations, the governing equations of the phase ﬁeld model of crack propagation in piezoelectric and ferroelectric materials are discussed in Section 4.4, along with the proposed solution algorithm.

4.1. Uncoupled mechanical/electrical crack face conditions

Most electromechanical crack models assume uncoupled mechanical and electrical boundary conditions on the crack faces. Although the uncoupled conditions are not physically realistic, they are good approximations of the coupled electromechanical conditions in some situations.

4.1.1. Traction free conditions Mathematically, these conditions are stated in the context of sharp crack models as

r þ Á n ¼ rÀ Á n ¼ 0,

ð11Þ

where the superscripts þ and denote the top and bottom crack faces and n is the unit normal. In the phase ﬁeld model

of brittle fracture, the traction free conditions are encoded by multiplying the elastic energy density F by the jump set

function ðv2 þ ZkÞ (see Eq. (1)). Working by analogy, in the case of the proposed electromechanical models, we multiply by the jump set function the ﬁrst two terms of the enthalpy density H associated with the elastic strains e. Therefore, the

elastic and electromechanical coupling potentials are annihilated in the fractured zone (v ¼0).

4.1.2. Traction free, electrically impermeable crack model In the case of impermeable cracks, the crack faces are treated as charge free surfaces, i.e. the normal components of the

electric displacement vanish on both crack faces as

D þ Á n ¼ DÀ Á n ¼ 0:

ð12Þ

The impermeable conditions imply that the crack cannot sustain any electric displacement inside the fractured zone. In the

phase ﬁeld model, we multiply all terms of H, involving the electric ﬁeld E, by the jump set function ðv2 þ ZkÞ. Therefore

the electric displacements vanish in the fractured zone (v ¼0), similar to the stresses. We illustrate in Appendix C (see Fig.

C2) that in the limit of small regularization parameter k, this formulation exhibits traction free and impermeable solutions

on the crack faces, as expected in the sharp crack model.

4.1.3. Traction free, electrically permeable crack model The permeable crack boundary conditions imply that the electric potential and the normal component of the electric

displacement are continuous across the crack faces, i.e. in the context of sharp crack models

f þ ¼ fÀ and D þ Á n ¼ DÀ Á n:

ð13Þ

In contrast to the impermeable crack model, encoding the permeable conditions requires that all terms of the enthalpy density H, involving the electric ﬁeld E, remain unmodiﬁed. As a way of example, the enthalpy density H of the traction free, electrically permeable crack model for the piezoelectric material is

H ¼ ðv2 þ ZkÞð12 e : C : e e : eT Á EÞ 12E Á KE:

ð14Þ

Accordingly, the electric displacements are

D ¼ ðv2 þ ZkÞe : e þKE:

ð15Þ

Since the traction free conditions result in a discontinuous displacement ﬁeld, Eq. (15) seems to be introducing a

discontinuity in the electric displacement through the coupling term e : e, in contradiction with Eq. (13). Nevertheless, we

show in Appendix C (see Fig. C1) that the effect of this discontinuity in the normal direction becomes negligible when the

regularization parameter k is small enough.

4.2. Electromechanical crack face conditions: energetically consistent crack model

As described in Section 1, the semi permeable crack conditions were introduced to model intermediate situations between the two over simplistic extremes of a completely closed crack (permeable) and an open crack ﬁlled with a medium of negligible conductivity (impermeable). This model assumes that the crack gap ﬁlling medium can sustain a certain degree of electric ﬁeld, thereby inducing an electric displacement inside the crack. In the context of sharp crack models, semi permeable crack face boundary conditions are formulated as

D þ Á n ¼ DÀ Á n ¼ Dc ¼ e0Ec,

ð16Þ

6

where Dc is the induced electric displacement and Ec is the electric ﬁeld in the crack gap. It is obtained as

Ec ¼

f þ fÀ

unþ

¼ uÀn

Df , Dun

ð17Þ

where the superscripts ‘‘ þ ’’ and ‘‘ ’’ denote the values on each crack face, Df is the potential drop across the crack gap and Dun is the crack opening displacement. The main deﬁciency stems from the fact that the electric ﬁeld inside the crack gap induces also tractions on the crack faces, which are not considered by the semi permeable conditions. It has been demonstrated that this model is not energetically consistent (Li et al., 2008). To overcome this inconsistency, energetically consistent (EC) crack face boundary conditions have been proposed in the context of sharp crack models (Landis, 2004a), in which

t þ ¼ scn, tÀ ¼ scn,

ð18Þ

where t þ and tÀ are the tractions acting on the crack faces and sc is the effective stress within the crack gap. According to

Landis (2004a) and assuming an inﬁnite breakdown strength of free space, the effective stress can be expressed as

sc ¼ 12e0E2c :

ð19Þ

In fact, EC conditions assume that the crack behaves electrically as a capacitor, storing electrical charge between the crack

faces. The associated electrical enthalpy is given by Z

Hc ¼ HcDun dS,

Sc

ð20Þ

where Sc denotes the surface (curve in 2D) of the sharp crack and the electrical enthalpy density Hc is

Hc ¼ 12e0E2c :

ð21Þ

In proposing the electrical enthalpy in Eq. (20), it is assumed that any electric ﬁeld and deformation components within

the gap parallel to the crack are negligible in comparison to the normal components. In the following, we develop the

general formulation of this enthalpy in the context of the phase ﬁeld models without any assumption regarding the

direction of the electric and deformation ﬁelds. Even though the present theory is geometrically linear, we need to

introduce geometric nonlinearity in selected terms to account for the emergence of a crack gap medium as a result of the

crack opening displacement. In fact, the enthalpy density in Eq. (4) can be obtained from the linearization of a general

enthalpy density of ferroelectrics in ﬁnite deformation electromechanics (Xiao and Bhattacharya, 2008). If this

linearization does not apply to the electrical enthalpy of free space (last term in Eq. (4)), we can account for the nonlinear enthalpy of the crack gap as follows. First, we express the electric ﬁeld in the deformed conﬁguration, Eb, in terms of the

reference gradient of the electric potential, Eb ¼ FÀT rf, where F is the deformation gradient tensor. Then, the electrical

enthalpy density of free space per unit deformed volume of material can be written as

Hb ¼ e0 9Eb92:

ð22Þ

2

This enthalpy can be decomposed into the enthalpy of the intact and fracture zones using the phase ﬁeld v

Hb ¼ v2Hb þ ð1 v2ÞHb :

ð23Þ

Ignoring all geometric nonlinearity in the intact zone, the total electrical enthalpy of free space becomes

e0 Z 2 2 e0 Z

2 ÀT

2

H0 ¼ 2 v 9rf9 dO 2 ð1 v Þ9F rf9 J dO,

ð24Þ

O

|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄOﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}

Hc

where J ¼ det F. The second term represents the phase ﬁeld counterpart of Eq. (20) and we can interpret Eb as Ec in the

fracture zone. Fig. 2 illustrates an opening crack where the electric ﬁelds are presented both for the sharp crack across the

faces and for the diffuse crack inside a partition of the fractured zone.

In Appendix C, we present numerical evidence substantiating the adequacy of this phase ﬁeld version of the EC

conditions. To provide a simple heuristic argument to support it, we consider a straight smeared crack along the x1

direction. The corresponding proﬁles of the ﬁelds u and f along the x2 direction are

8 >><

u

þ

x2 4 ak,

u þ uÀ

u þ þ uÀ

uðx2Þ ¼ >>: 2ak x2 þ 2

uÀ

ak r x2 r ak, x o ak,

ð25Þ

2

8>>>< ff þþ fÀ f þ þ fÀ x2 4 ak,

fðx2Þ ¼ >>> 2ak x2 þ 2

ak r x2 r ak,

ð26Þ

: fÀ

x2 o ak,

7

Fig. 2. Conﬁguration of an opening crack where the electric ﬁelds Ec and Ec are presented for the sharp and diffuse crack models, respectively. The crack faces are indicated by the dashed lines and a portion of the fractured zone is presented by the solid lines.

while the proﬁle of the ﬁeld v is given in Eq. (2) with dðxÞ ¼ 9x29. According to the second term of Eq. (24), the electrical

enthalpy of the smeared crack can then be written as

Hc ¼

e0 Z Z ak 9FÀT rf92J dx2 dS ¼

2 Sc Àak

e0 Z 2ak9FÀT rf92J dS ¼

2 Sc

e0 Z ðf þ fÀÞ2 2 S 2ak þ u þ uÀ dS,

c

ð27Þ

where we have used

21

03

FÀT ¼ 664 0

1 uþ

uÀ 775,

1 þ 2ak

2

3

0

rf ¼ 64 fþ fÀ 75,

2ak

u þ uÀ

J ¼ 1 þ 2ak :

When k tends to zero, the width ak also tends to zero and Eq. (27) converges to the enthalpy of the sharp crack in Eq. (20).

Of course, in practical computations the width of the smeared crack is not negligible and there is also a smooth transition of the ﬁeld v across the fracture zone. However, we illustrate in Appendix C (see Fig. C3) that for small enough values of the

regularization parameter k, the total electrical enthalpy of the diffuse crack accurately follows the sharp crack model.

The ﬁrst term in Eq. (24) is the enthalpy term of the impermeable crack model. Therefore this equation indicates that the enthalpy density H of the EC crack model is obtained by adding the electrical enthalpy density of the diffuse crack to that of the impermeable crack model, which allows the displacement and electric potential ﬁelds to be discontinuous. Given the enthalpy density in the second term of Eqs. (24) and (6), the stresses and electric displacements are obtained in the fracture zone as

r ¼ J ð1 v2Þrc,

ð28Þ

2

D ¼ Jð1 v2ÞFÀ1Dc, with

Dc ¼ e0Ec ¼ e0FÀT E,

ð29Þ ð30Þ

rc ¼ FÀ1ðDc EcÞ þ ðDc EcÞFÀT

e0 9Ec92ðFÀ1 þ FÀT Þ:

2

ð31Þ

Note that the lower permittivity of the free space as compared to that of the material induces a high electric ﬁeld inside

the fracture zone, especially in a specimen under applied electrical loading. The free space can sustain a certain level of

electric ﬁeld, called dielectric strength, above which the free space experiences failure of its insulating properties, i.e.

electrical breakdown. This effect can be modeled mathematically as a bound constraint on the magnitude of the electric

ﬁeld inside the fracture gap

( 9FÀT rf9

9Ec9 ¼ Ed

if 9Ec9 r Ed, if 9Ec9 4 Ed,

ð32Þ

where Ed is the dielectric strength or the electrical discharge level of the free space. This constraint can be added to the governing equations. The corresponding Lagrange multiplier modiﬁes the electric displacement inside the fracture zone. See Appendix D for details. We evaluate the effect of the electrical breakdown constraint in Section 5.

8

In the above formulation, the free space is assumed to be the medium ﬁlling the crack gap. Other crack gap ﬁlling media

can be considered as well by replacing the free space permittivity e0 with the corresponding value em in Hc in Eq. (24).

4.3. Polarization boundary conditions

In addition to the crack face conditions discussed in the previous sections, crack face conditions on the polariza tion need to be supplied for ferroelectric materials. We discuss next the two usual choices and their phase ﬁeld implementation.

4.3.1. Free polarization crack face conditions The free polarization boundary conditions assume that the polarization distribution is unaffected by the presence of

the crack, and hence dictated by the bulk material model. Mathematically, they are written as

rp þ Á n ¼ rpÀ Á n ¼ 0:

ð33Þ

This condition can be encoded in the phase ﬁeld model by multiplying with the jump set function ðv2 þ ZkÞ the energy

terms involving the gradient of polarization, i.e. the domain wall energy density U.

4.3.2. Zero polarization crack face conditions In contrast to the bulk material, the polarization vanishes in the free space ﬁlling the crack gap. The zero polarization

boundary conditions assume that the polarization is continuous at the material free space interface, i.e.

p þ ¼ pÀ ¼ 0:

ð34Þ

This condition can also be formulated in a similar way, by multiplying with the jump set function the energy terms associated with the polarization. In Appendix C (see Fig. C4), we show that these methods indeed produce numerical solutions satisfying the free and zero polarization boundary conditions in a diffuse sense.

4.4. Phase ﬁeld model of crack propagation in piezoelectric and ferroelectric materials

Based on the discussion in Sections 4.1 4.3, the different electromechanical crack face boundary conditions along with

the corresponding electromechanical enthalpy density, stresses and electric displacements for piezoelectric and ferro

electric materials are summarized in Appendix E. We now can form the total electromechanical enthalpy of a possibly

fractured piezoelectric or ferroelectric body by adding the surface energy in Eq. (1) to the total electromechanical enthalpy

in Eq. (3) as

Z"

2

#

ð1 vÞ þ k9rv92 dO,

ð35Þ

Hs½u,p,f,v ¼ H½u,p,f,v þ Gc

O

4k

where H now depends also on the ﬁeld v through the modiﬁed electromechanical enthalpy density H for different crack

face boundary conditions.

To capture interactions between the fracture and the microstructure processes in ferroelectrics, the crack propagation

should not be overwhelmingly faster than the microstructure evolution, and vice versa. In practice, the relative kinetics of

the microstructure evolution and the crack propagation gives the two phenomena a chance to interact. In the absence of

detailed experimental or theoretical information on this kinetics, v is selected together with the polarization as primary

order parameters and ﬁnite mobilities are introduced for the fracture and microstructure processes. Then, the time

evolution of the system results from the gradient ﬂow of the primary order parameters, assuming that the displacement

and the electric ﬁeld adjust immediately to mechanical and electrostatic equilibrium. The gradient ﬂow of the polarization

and the weak form of the mechanical and electrostatic equilibrium have already been introduced in Eqs. (10) and (9),

respectively. Here, the gradient ﬂow of the ﬁeld v as an additional governing equation is obtained as

Z

Z

Z

@H dv dO Gc v 1 dv þ kv,idv,i dO,

ð36Þ

mv v_ dv dO ¼ dHs½u,p,f,v; dv ¼ O

O @v

O 4k

where 1=mv is the mobility of the fracture process. In fact this equation is a Ginzburg Landau type evolution equation,

where the left hand side term is the dissipation in the process zone and the mobility constant controls the rate of the energy dissipation (Hakim and Karma, 2009). A recent study, using a phase ﬁeld model of brittle fracture, indicated that the mobility parameter controls the crack velocity, particularly at the initial stages of an unstable crack propagation (Kuhn and Mu¨ ller, 2010). Therefore the proper selection of this parameter allows us to track crack propagation, even in mechanically unstable conditions.

The ﬁnite element equations are obtained from the weak forms in Eqs. (10), (9), and (36) with the enthalpy density H,

the stresses r, and the electric displacements D for each crack face condition given in Appendix E. Eqs. (10) and (36) are discretized in time with an implicit scheme from time tmÀ1 to tm ¼ tmÀ1 þ Dt. The Newton Raphson method is implemented for the implicit time integration of Eq. (10) due to the nonlinear terms of the phase separation energy (w).

Algorithm 1 presents a simple procedure to solve forward in time the coupled system in a staggered, iterative way. The

9

function g encodes the data for the applied electromechanical load as a function of the time step. Since the crack should

not be allowed to heal (irreversibility condition), the nodes reaching a value of v below a certain threshold g, are assigned a

ﬁxed value of v ¼0 for the rest of the calculation. Note that the piezoelectric models do not require the polarization data and the polarization evolution in step 6 of the algorithm is only computed for ferroelectrics. In this step the constraint for the polarization evolution is only used for the zero polarization conditions. In fact, with this constraint the mobility term of the gradient ﬂow is ﬁxed to zero in the fractured zone (v ¼0).

Algorithm 1. For crack propagation in piezoelectric and ferroelectric materials.

1:

Let m 0 and t0 0

2:

Set v0 vinit, p0 pinit , f0 0 and u0 0

3:

repeat

4:

m’mþ 1

5:

tm’tm 1 þ Dt

6:

Compute pm in Eq. (10) using pm 1, um 1, fm 1 and vm 1 under the constraint pm 0 for vm 1 0

7:

Compute um and fm in Eq. (9) using pm and vm 1 under the electromechanical load gðtmÞ

8:

Compute vm in Eq. (36) using pm, um, fm, and vm 1 under the constraint vm 0 for vm 1 r g and Hv 0 for Hv o 0

9:

until m n

We note that the EC conditions discussed in Section 4.2 introduce non linear terms into the ﬁnite element equations. Here,

we implement an internal loop in step 7 of Algorithm 1 to solve this non linearity. In each iteration of this loop, the non

linear terms are updated with the last values of the mechanical displacement and the electric potential, then Eq. (9) is solved.

This iteration continues until a steady state is reached for both the mechanical displacement and the electric potential.

An important aspect of electromechanical crack propagation deserves special attention. Generically, the enthalpy

functional can be written as

Z

Z"

2

#

O½ðv2 þ ZkÞHvðu,p,fÞ þ Hrestðu,p,fÞ dO þ Gc O ð14kvÞ þ k9rv92 dO: ð37Þ

Under high applied electrical load, Hv may become negative. Then, when minimizing the energy with respect to v, it may be favorable to localize in narrow regions high values of v, much above one. Such anomalous localization zones are not physically meaningful. In sharp crack models of electromechanical fracture, this issue manifests itself with the negative energy release rates found at high applied electric ﬁelds (Li et al., 2008). Numerically, we deal with this issue by limiting the maximum value of v to one, or alternatively, by setting negative values of Hv to zero in step 8 of Algorithm 1.

5. Numerical simulations and discussion

In this section, we examine the effects of the different crack face conditions on the crack propagation. The simulations of a propagating crack in a poled piezoelectric material under different applied electromechanical loadings are presented in Section 5.1. We also exercise the model with different crack gap ﬁlling dielectric ﬂuids and introducing the constraint of electrical breakdown as discussed in Section 4.2. Section 5.2 presents the results of crack propagation in a poled ferroelectric single crystal where the crack interacts with domain switching. The results of each section are discussed in detail.

5.1. Propagating cracks in piezoelectrics

To investigate the effects of the different crack face conditions on the crack propagation, a set of simulations is performed following Algorithm 1. A four point bending setup is considered. Fig. 3 presents the dimensions of the specimen along with the boundary conditions and the pre crack. Plane strain conditions and the material properties of PZT PIC 151 poled along the x2 axis (normal to the pre crack) are assumed (see Appendix B for the material parameters). The value of a is the length of the region with v ¼0 (i.e. the crack length). The pre crack with a length of a¼ 0.5 mm is introduced in the model using the proﬁle in Eq. (2). The model is discretized with a ﬁne mesh, with the smallest element size of h ¼10 4 mm in the vicinity of the pre crack. An adaptive mesh reﬁnement algorithm is employed to generate the mesh. A typical mesh

is presented in Fig. 4. The regularization parameter k is selected as four times the smallest element size, i.e. k ¼ 4 Â 10À4 mm, and the residual stiffness is set to Z ¼ 10À6. The inverse mobility of the fracture evolution is set to mv ¼ 0:1 Ns=m2. The point load P is applied incrementally using the load function P(t)¼At with a rate of A ¼100 N/s and a time step of Dt ¼ 10À2 s. A total number of 2.6 Â 104 time steps are performed for each simulation. The crack initiates and

becomes unstable when the load reaches P¼14 KN, at time t¼ 140 s. From this time on, the point load is ﬁxed and the crack length relative to the initial length (Da) is recorded at regular intervals, see Fig. 5 for the different crack face conditions. The ﬁnite mobility of the v ﬁeld allows us to follow a slow crack propagation in this period. During the crack propagation, external electric ﬁelds E¼ 71 MV/m (see Fig. 3 for the sign convention) are applied from time t ¼190 s. It is apparent that the impermeable crack is completely arrested by the application of both the positive and negative electrical

10