Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 44, 4, pp. 819-836, Warsaw 2006 INTERACTION OF POINT DEFECTS IN PIEZOELECTRIC MATERIALS – NUMERICAL SIMULATION IN THE CONTEXT OF ELECTRIC FATIGUE Oliver Goy Ralf Mueller Dietmar Gross Division of Mechanics, TU Darmstadt, Germany e-mail: goy@mechanik.tu-darmstadt.de The paper is concerned with the problem of electric fatigue in functional materials such as piezoelectric sensors and actuators. The fatigue degrades electromechanical properties with an increasing number of working cycles during which ionic and electronic charge carriers interact with each other within the bulk and on the interfaces of a material. This, in turn, influences local electric and mechanical fields coupled via the piezoelectric effect. It is assumed in the paper that the electric charges (vacancies) can be modelled as point defects, which tend to agglomerate and finally form clusters, espe- cially in the vicinity of electrodes. As a result, the piled up defects affect the distribution of polarisation. To solve the problem quantitatively, the analy- sed system is discretised with a grid of square cells and then transformed by FFT. Several numerical examples are formulated. Migration of a single defect, attraction of two and clustering ofmore point defects are thoroughly discussed and illustrated in the paper. Keywords: functionalmaterials, piezoelectricity, electric fatigue,pointdefect 1. Introduction Ferroelectric materials occur in a wide field of applications and can be found in actuators, sensors and in electronics. A general description of ferroelectrics can be found in Xu (1990). Especially in the field of actuators there exists a demand for high precision and sustainability. A material has to function during a high number of loading cycles, which can lead to fatigue phenome- na. Electric fatigue in functional materials encompasses a set of phenomena which result in degradation of material properties with an increasing number of electrical or mechanical cycles. Ionic and electronic charge carriers, later 820 O. Goy et al. modelled as point defects, interact with each other, with microstructural ele- ments in the bulk and on the interfaces. This causes a change of both, the local electric and mechanical fields, which results from the electromechanical coupling. Themain assumption at this stage is, that the point defects tend to agglomerate and to form clusters, especially in the vicinity of electrodes. As a result, the piled updefects could have an effect on the distribution of the pola- risation. Ferroelectrics change their polarisation and coupling behaviour by the do- main switchingmechanism, where a domain wall separating domains of diffe- rent polarisations is moving through the material. Defects like the mentioned defect clusters are able to slow down or even stop the domain wall movement, which leads to a reduced ability of domain switching, so that the coupling effect gradually decreases. Further investigations on the role of point defects during the electric fatigue of lead zirconate titanate (PZT) have been made in Lupascu (2004) and Lupascu and Roedel (2005). The influence of defects on the mobility of domain walls has been simulated in Schrade et al. (2006). There, gadolinium molybdate (GMO) has been used to model macroscopic defects on the surface (holes in an electrode, lateral defects), and defects in the bulk (polarisation defect). As a main result, it can be stated that defects can significantly slow downor even stop themotion of domainwalls, such that higher external electric fields are necessary to enforce further switching. The results are in good agreement with experimental observations. In addition to the investigation ofmacroscopic defects, effortswill bemade here to calculate the fields caused by point defects themselves, and to verify the parameters used.This can be regarded as the first step in getting a deeper insight in themechanismsof clusteringof pointdefects.Before consideringmo- re involved settings, an elementary investigation of a single defect and special configurations of two and three defects will be the topic of studies here. The analysis of those basic settings can provide a better understanding of complex situations. The calculation canbedonebydifferentmethods, numerical andanalytical ones. Here, a combination of Difference Methods and Fast Fourier Transform is used to obtain a numerical solution. This provides a fast and flexible tool for simulations concerning the interaction of defects. Ferroelectrics often exhi- bit transversely isotropic or orthotropic properties. For calculations, a linear theory with the coupling will be used. The interaction of pointdefectswill be treatedbymeans ofmaterial or con- figurational forces. More on that topic in non-coupled and coupled materials can be found e.g. in Gross et al. (2003), Gurtin (1995, 2000), Maugin (1993), Kienzler and Herrmann (2000). Configurational forces will be used in combi- nationwith thermodynamically admissible kinetic laws, to simulate the defect motion. The aim is to comprehend the defect agglomeration in ferroelectrics. Interaction of point defects in piezoelectric materials... 821 2. Basic equations In the following, we consider a continuous bodywhich occupies the domain B with the boundary ∂B.When neglecting inertia terms, the symmetric Cauchy stress tensor σ satisfies the equilibrium condition with volume forces f divσ+f =0 (2.1) The electric displacement D satisfies the electro static relation divD= q (2.2) where q represents electric charges in the bulk. Considering small deforma- tions, the infinitesimal strain tensor ε consists of the symmetric part of the displacement gradient ∇u, while the electric field E is the negative gradient of its potential ϕ ε= 1 2 ( ∇u+(∇u)⊤ ) and E=−∇ϕ (2.3) The constitutive equations for the electromechanical coupled system read σ= ICεe− Ib⊤E D=Ibεe+AE+P0 (2.4) where εe = ε−ε0 represents the elastic strain of thematerial. The elasticma- terial behaviour is described by the stiffness tensor ICwhich has the symmetry properties Cijkl = Cklij = Cjikl. The electromechanical coupling is characte- rised by the piezoelectric tensor Ib, where Ibijk = Ibjik and (Ib ⊤)ijk = Ibkij, while the dielectric properties are given by A with Aij = Aji. All material tensors will be assumed to be constant with respect to x. The irreversible quantities are the inelastic strain ε0 and the remanent polarisation P0. The constitutive equations can be obtained from the electric enthalpy function H(εe,E,x)= 1 2 εe : (ICεe)−E · (Ibεe)− 1 2 E · (AE)−P0 ·E (2.5) by differentiation σ= ∂H ∂ε and D=− ∂H ∂E (2.6) With (2.4) and (2.3), equations (2.1) and (2.2) can be rewritten in an index notation for Cartesian coordinates Cijkluk,jl+Ibkijϕ,kj = Cijklε 0 kl,j (2.7) Ibijkui,jk−Aijϕ,ij =Ibijkε 0 jk,i+q 822 O. Goy et al. In (2.7), we have taken advantage of the symmetries of the material tensors and neglected the influence of volume forces and polarisation. In this paper, the admitted boundary conditions on ∂B can be either uni- formor periodic. First, the uniformconditionswill be defined.Themechanical boundary conditions can be prescribed values for the displacement u or the traction vector t on the surface u=u0 on ∂Bu σn= t0 on ∂Bt (2.8) For the electric quantities, the electric potential ϕ or the surface charge Q0 can be set to ϕ=ϕ0 on ∂Bϕ Dn=−Q0 on ∂BQ (2.9) When considering periodic boundary conditions, the body B is a periodic cell, that is, it can be reproduced in two (plane problems), or three non-parallel directions.Thedisplacement andpotential fields can be continued periodically from one cell to another. When considering the defect interaction, it is necessary to define configu- rational forces, also known as driving forces acting (except for the definition of the sign) on each defect. The energy-momentum tensor Σ and thematerial forces g fulfil the following equation divΣ+g=0 (2.10) with Σ=H1− (∇u)⊤σ−∇ϕ⊗D (2.11) g=−(∇u)⊤f+∇ϕq− ∂H ∂x ∣∣∣ expl. +σ :∇ε0 where (∂(·)/∂x)|expl. denotes the explicit derivation with respect to x, when all other variables are kept constant, see Gross et al. (2003). 3. Point defects In a perfect crystal with an undisturbed lattice, all atoms are at perfect lattice sites. On the continuum level, this is modelled as a homogeneous medium. By creating vacancies or inserting foreign atoms, a disturbance in the strain Interaction of point defects in piezoelectric materials... 823 and charge state occurs. In continuum mechanics, this is modelled as a point defect with localised strain and charge concentration in its vicinity. For the mechanical defect, the assumption of an isotropic eigenstrain has beenmade ε0ij =αδ(x−η)δij (3.1) but this could be easily expanded to anisotropic eigenstrains. In (3.1), a point defect located at the position η is assumed. The only material parameter α depends on lattice parameters and can be either positive or negative. The positive sign represents a foreign atom too big for its lattice site, the negative one is a too small foreign atom or a vacancy. Similarly, free point charges at the point η are modelled with the parameter β q=βδ(x−η) (3.2) Here, the sign of β determines the sign of the charge of the defect. In the case of a homogeneous material, where H is not explicitly depending on x, with absence of volume forces f and by using (2.11)2, (3.1) and (3.2), thematerial forces on a point defect can be written as G= lim r→0 ∫ Vη(r) [ σ :1α∇δ(x−η)+∇ϕβδ(x−η) ] dV = (3.3) = [ −α∇(trσ)+β∇ϕ ]∣∣ x=η =−∇Λ ∣∣ x=η where Vη(r) represents a volume with the radius r and centre η, which is contracted to a point. The field Λ :=αtrσ−βϕ can be regarded as the force field of a defect. Putting a second defect into this field will cause migration along the direction of the field gradient if the kinetic law admits the diffu- sion. The second law of thermodynamics gives the following relation for the dissipation D D=−G ·v­ 0 (3.4) where v= ∂η/∂t is the velocity of the defect. With the simple law v=−cG c= const > 0 (3.5) the defect migration will be simulated later in Section 6.2. By considering (3.1) and (3.2), equations (3.6) read Cijkluk,jl+Ibkijϕ,kj = Cijklδklαδ,j(x−η) (3.6) Ibijkui,jk−Aijϕ,ij =Ibijkδjkαδ,i(x−η)+βδ(x−η) 824 O. Goy et al. In the case of periodic boundary conditions, the solution toEqs. (3.6) is uniqu- ely determined up to a constant value,which does not influence the periodicity of the results. The treatment of this constant will be described in Section 5. An analytical solution for general anisotropic materials of the coupled set of elliptic differential equations (3.6) can not be given in a closed form. Howe- ver, it is possible to derive integral representations of u and ϕ by means of Fourier- or Radon Transforms (also known as Plane Wave Transforms). For elastostatics, fundamental solutions and their derivatives can be found e.g. in Barnett (1972) or Vogel and Rizzo (1973). Application of this method to po- int defects in elastic materials can be found in Dobovs̆ek (1998, 2000). With respect to the implementation in the Boundary Element Method (BEM) the reader is referred to Schclar (1994). The advantage of integral or closed form solutions is that the considered body is infinite, so the boundaries have no influence on the solution. Another approach is a purely numerical one.The simplestmethod is aFiniteDifference (FD) scheme on a rectangular mesh. Boundary conditions can be prescribed in theDirichlet orNeumann form, or periodic conditions could be given. Ano- ther popular method is the Finite Element Method (FEM), where naturally the Dirichlet and Neumann boundary conditions are prescribed. For details the reader is referred to Mueller et al. (2006). With some efforts, periodic conditions can be implemented, but this is not further investigated here. The approach carried out in this paper is a combination of a FD scheme and the Discrete Fourier Transform (DFT). The major disadvantage of FD schemes is a large system matrix, which has to be solved. Although it is a sparsematrix,many efforts are necessary tomanage the allocating ofmemory for this operation, especially in 3D.Furthermore, it is very time consuming.To avoid these problems, the DFT is applied on the discretised set of equations. Afterwards, the derivatives can be arranged in a way, that for each grid point a small set of linear equations has to be solved in the Fourier transformation space, which has to be transformed back numerically. 4. Discrete Fourier Transformations TheFourier-Transform f̂(ξ) of a function f(x) is defined as F(f)= f̂(ξ)= ∫ Rn f(x)e−ix·ξ dξ (4.1) where x,ξ ∈ Rn. This transformation transfers the function f(x), which is often real valued (but can also be complex) into the complex function f̂(ξ). Interaction of point defects in piezoelectric materials... 825 In other words, the function f is transformed from the so called signal doma- in into the frequency domain. The reverse operation is the inverse Fourier- Transform F−1(f̂)= f(x)= 1 (2π)n ∫ Rn f̂(ξ)eix·ξ dξ (4.2) The Fourier-Transform has some useful properties • F ( ∂f ∂xj ) = iξjF(f) • F ( δ(x) ) =1 (4.3) • F(af+ bg)= aF(f)+ bF(g) When having only discrete values on a grid at hand, the continuous transform canbewritten in a discrete form, taking only the grid points into account.The integral over the infinite space has to be split up into sumsover the grid points, so finite intervals [0,a] and [0,b] can be computed in the x- and y-direction, respectively f̂(ξ)≈ f̂(p,q)= M−1∑ m=0 N−1∑ n=0 f(m,n)e−i2πpm/M e−i2πqn/N (4.4) with p=0,1, . . . ,M−1, q=0,1, . . . ,N−1. In this context, a transformation in 2D is considered, where p and q are integer numbers indicating grid positions, and M and N are the total numbers of points in each direction. So, with the distance of the grid points hx = a/(M − 1) and hy = b/(N − 1), the position vector transforms to x→ [mhx,nhy] ⊤. The inverse transformation reads as follows f(x)≈ f(m,n)= 1 MN M−1∑ p=0 N−1∑ q=0 f̂(p,q)ei2πpm/M ei2πqn/N (4.5) with m=0,1, . . . ,M−1, n=0,1, . . . ,N−1. It is in the nature of theDFT that the transformed function f is periodic, so the result of (4.5) will be periodic in m and n. Only quadratic domainswill be considered, so M = N and hx = hy = h. The properties of the discrete transformation used later are • F(fm±1,n)= f̂pqe ±i2πp/N (4.6) • F(fm,n±1)= f̂pqe ±i2πq/N 826 O. Goy et al. wheredefinitions fm±1,n := f(m±1,n) and f̂pq := f̂(p,q) havebeenused.The cost of a one-dimensional transform are O(N2) operations. This can be redu- ced to O(N log2N) operations by using the Fast-Fourier-Transform (FFT). For the examples in Section 6, theFFTW routines have been used, which are free available1. FFTW employs algorithms of the order O(N log2N) for all lengths of data vectors including prime numbers. Whendiscretising the domain, the space derivatives have to be substituted by difference formulas. The first derivatives of the real valued function f(x) with respect to the coordinates xi are ∇f :=    ∂f ∂x1 ∂f ∂x2    ≈    ∇ f hfm,n := 1 h [ fm+1,n−fm,n fm,n+1−fm,n ] (forward) ∇bhfm,n := 1 h [ fm,n−fm−1,n fm,n−fm,n−1 ] (backward) ∇chfm,n := 1 2h [ fm+1,n−fm−1,n fm,n+1−fm,n−1 ] (central) (4.7) The second derivatives become ∂2f ∂x21 ≈ fm+1,n−2fm,n+fm−1,n h2 ∂2f ∂x22 ≈ fm,n+1−2fm,n+fm,n−1 h2 (4.8) ∂2f ∂x1∂x2 ≈ fm+1,n−fm−1,n+fm,n+1−fm,n−1 2h2 In the following, only the central difference will be considered for the first derivatives for accuracy reasons. Finally, (4.7)3 and (4.8) can canbewritten as F ( ∇chfm,n ) = 1 2h   e+ω(p−1)− e−ω(p−1) e+ω(q−1)− e−ω(q−1)   f̂mn F (∂2f ∂x21 ) = 1 h2 ( e+ω(p−1)+e−ω(p−1)−2 ) f̂mn (4.9) ... where ω= i2π/N. By means of the Eulerian formulas, all exponential functions can also be expressed in trigonometric functions. It is obvious that thederivatives onagrid 1http://www.fftw.org Interaction of point defects in piezoelectric materials... 827 point (m,n) only depend on the Fourier transformed values of f̂ on the same grid points (m,n) andnot, as itwas the case in (4.7), on their neighbours.Due to that fact, each variable in theFourier space canbecalculated independently. Afterwards, an inverse transformation is performed. 5. Discretisation For calculation of the numerical solution for point defects, the domain B introduced in Section 2 has to be discretised. This is done by choosing a square cell of the edge length a with N grid points in each direction. The Fourier Transformgives periodic boundary conditions, so the cell is embedded into an infinite matrix of identical cells. For equations (3.6), the well known Voigt notation is used for σ, ε and the material tensors, taking advantage of the symmetries (ij → I, kl→J) (σ)I=̂    σ11 σ22 σ12    (ε)I=̂    ε11 ε22 2ε12    =    u1,1 u2,2 u1,2+u2,1    (5.1) (C)IJ=̂    C11 C12 0 C12 C22 0 0 0 C33    (b)kI=̂ ( 0 0 b13 b12 b22 0 ) where Cijkl →CIJ and Ibkij → bkI. Differential equations (3.6) become CIJεJ,j − bkIEk,j =CIJαIJδ(x−η) (5.2) biJεJ,i+AijEj,i = biJαIJδ(x−η)+βδ(x−η) with I = [1,1,0]⊤. When applying (4.7)3 and (4.8) to equations (4.9), the differential operators can be written as ∂ ∂xj : dj := 1 h isin (2π N (pj −1) ) ∂2 ∂xj∂xk : djk :=− 1 h2 sin (2π N (pj −1) ) sin (2π N (pk−1) ) (5.3) where pi=̂[p,q] ⊤. So the transformed equations can be expressed in a system of linear equations for fixed p and q Mûpq = r (5.4) 828 O. Goy et al. with M=    C11d11+C33d22 C12d12 (b21+ b13)d12 C33d11+C22d22 b13d11+ b22d22 sym. −(A11d11+A22d22)    (5.5) ûpq =    (û1)pq (û2)pq ϕ̂pq    r=    (C11+C12)d1αδ̂pq (C12+C22)d2αδ̂pq (b21+ b22)d2αδ̂+βδ̂pq    The discrete delta function δ with origin at the point (md,nd) can be easily obtained under consideration of property (4.3)2 δmn :=    1 h2 if m=md, n=nd 0 else (5.6) It should bementioned that ûpq, ϕ̂pq and δmn are representingmatrix entries of the field values at discrete grid points. The Discrete Fourier Transform of (5.6) becomes δ̂pq = 1 h2 e−i2π(mdp+ndq)/N (5.7) System (5.4) can be easily inverted and the field variables û can be trans- formed back by the numerical FFT. As described in Section 4, the FFTW routines are used for this purpose. It should be emphasised, that the index combination p=1, q=1 in (5.4) requires special treatment because coefficient matrix (5.5)1 becomes singular. In this special case, a constant term appears in the transformation, which can be chosen arbitrarily. It represents a rigid body translation in the mechanical displacement or a bias potential field in electrical terms. 6. Examples With the objective of getting a better physical understanding of the results, a two-dimensional representation has been chosen, where a plane strain state is assumed. Thematerial considered is a transversally isotropic PZT,which is poled along the y-axis. So, the material constants are Interaction of point defects in piezoelectric materials... 829 CIJ=̂    13.7 7.16 0 7.16 12.4 0 0 0 3.14    ·104 N mm2 bkI=̂ ( 0 0 10.4 −4.0 23.3 0 ) ·10−6 C mm2 (6.1) Aij=̂ ( 7.95 0 0 5.15 ) ·10−12 C Vmm Additionally, the parameters α and β will be chosen positive or negative to show the dependence of the fields on the charge and misfit of defects. A possibility to obtain these parameters is the adjusting of the continuous fields to atomistic simulations. Studies ofMgO resulted in the first estimate for the mechanical defect parameter α2. Here, a periodic cell is considered, where by means of Density Functional Theory (DFT) the displacement of atoms in the vicinity of the oxygen vacancy is determined. The continuous displacement field is adjusted to discrete values at the atom positions of the relaxed atom lattice. Hence, the value for α lies in the region of −10−18mm3. The electric parameter β will be roughly approximated by a multiple of an elementary charge e0 = 1.602 · 10 −19C. Assuming that we have an oxygen vacancy at hand, the parameter would be β=3.204 ·10−19mm3. The domain of interest is a square of the length a, where amesh of N×N gridpoints is applied.Thedistance h= a/(N−1) is determining thenumerical error. 6.1. A single point defect At first, the fields of a single point defect will be investigated. The material parameters are chosen as described in Section 6. The displacement fields in Fig.1a and Fig.1b exhibit plausible qualitative results. Material points in the domain are moved towards the defect, and the absolute values of the displacement show the singularity in the vicinity of the defect. Singular behaviour is also observed in the electric potential ϕ, shown inFig.1c.Due to theperiodicity, theperiodic variable u1 vanishes on the right and left boundary in Fig.1a, while u2 and ϕ vanish at the upper and lower boundary in Fig.1b and Fig.1c. All fields are dominated by the mechanical defect parameter α, and the weak anisotropy results in a slightly stiffer ma- terial response in the x1-direction. The isolines of the electric potential field caused only by β > 0 would be ellipses with vertical semi-minor axes. 2These results were obtained in cooperation with the department of Chemistry, University of TechnologyDarmstadt 830 O. Goy et al. Fig. 1. (a) u1, (b) u2 [mm] and (c) ϕ [V] with α=−1 ·10 −18mm−2, β=3.204 ·10−19C/mm2 The stress components and the electric displacement vector, shown in Fig.2 and Fig.3, respectively, are also dominated by the mechanical parame- ter α. The stress distribution looks very plausible: on the x1-axis the stress σ12 vanishes and there is tension (positive values for σ11) along lines with x2 = const andpressure(negative values for σ22) along lineswith x1 = const. The same behaviour can be seen on the x2-axis, where σ11 and σ22 change their sign. Due to the anisotropy, slightly higher stresses in the x1 direction can be observed. The electric displacement D1 is strongly influenced by the shear strains, which can be seen from the shear stresses σ12 in Fig.2c. The se- cond component D2 mostly results from the electromechanical coupling with the x2 strain component, see σ22 in Fig.2b. In the case of an electric defect only (β > 0), the electric displacement vectors would form a central vector field, where D points away from the centre. Another point of interest is the interaction of defects and to investigate how they migrate through the material. Due to the periodicity of the cell, there is no interaction of defects with the boundary. Hence, resulting driving force (3.3) on a single defect will vanish, and therefore it has no tendency to move. Mathematically spoken, the gradient on the Λ-field is vanishing at the point η in this case. On the basis of Fig.4a, the first estimation for defect interaction can be obtained.When inserting a second defect, the influence on Interaction of point defects in piezoelectric materials... 831 Fig. 2. (a) σ11, (b) σ22 and (c) σ [MPa] with α=−1 ·10 −18mm−2, β=3.204 ·10−19C/mm2 Fig. 3. (a)D1, (b)D2 [C/mm 2] with α=−1 ·10−18mm−2, β=3.204 ·10−19C/mm2 that defect is determined by the gradient of the Λ-field and vice versa. So the second defect will always take the path with the greatest ascend andmove in regionswithhighervalues.At thispoint, it canalreadybeanticipatedwhat the general trend of configurations withmore than one point defect will be. Thus, the second defect of the same properties would be attracted when inserted in the areas around the diagonals of the cell and would be repelled when placed along the direction of the poling axis or perpendicular to it. It should be emphasised here, that the defect fields are not reflected by boundaries due to the periodicity of the cell. But the periodic boundary conditions also imply that two repelled point defects come to an equilibrium position, where each 832 O. Goy et al. defect experiences the same force from the defect in the considered cell and the neighbouring cell. These predictions will be verified in Section 6.2. Fig. 4.Λ-field with (a) α dominant and (b) β dominant In Fig.4a, the field is dominated by the mechanical parameter α, and therefore the mechanical hydrostatic stress as well, see equation (3.3). Addi- tionally, the absolute values of the stress field caused by the parameter α are very large compared to those of thepotential field, seeFigs.2a,b andFig.1c. So the first term in (3.3) becomes dominant. When increasing the parameter β, the influence of the electric potential becomes obvious. First of all, the second term in (3.3) becomes more dominant. But also the topology of the electric potential field changes. So the region of negative Λ values widens, while the anisotropy in the dielectric constants becomes very limited. Furthermore, a gradient in each corner of the cell is observable. There are now repulsive and attractive forces caused bymechanical stresses and repulsive forces caused by the electric potential. Consequently, the area with the repulsive effect is get- ting widened. Defects will show the tendency to align diagonally due to local maxima of the Λ-field. 6.2. Selected configurations of multiple point defects In this section, several arrangements of two and three defectswill be investiga- ted. The interaction of two defects will be simulated bymeans of trajectories. At first, two defects arranged along horizontal (Fig.5a) and vertical lines (Fig.5b) are examined. The driving force is pointing outwards, so the defects will move apart from each other, as simulations show. The displayed force vector is a normalised vector and does not account for the strength of the driving force: G0 = −G/‖G‖. It should be emphasised that the two defects will takeupan equilibriumposition,where eachdefect has an equal distance to the other defect and to its counterpart in the neighbouring cell, as mentioned previously. Only slightly perturbed positions from the above mentioned horizontal and vertical arrangement, as in Fig.6a lead to completely different behavio- Interaction of point defects in piezoelectric materials... 833 Fig. 5.Λ-field for two defects with driving force G0 arranged (a) horizontally and (b) vertically ur. The direction of the gradient on the Λ-field does not coincide with the connecting line between the two defects anymore. When applying law (3.5), the trajectories of the two defects can be calculated by evaluating the confi- gurational forces at the defect positions andmoving them one grid point into the negative force direction. The trajectories of the two defects in Fig.6b can not be concluded from the initial direction of the configurational forces. The small disturbance from the vertical alignment drives the defects to a diagonal alignment, where both attract each other, so they finally collapse. Different behaviour is observed in Fig.7, where an equilibrium position is achieved by putting the defects far away from each other. The initial configurational forces predict the tendency of the two defects to push each other away. Here, the dependence on the distance of the two defects becomes obvious. Fig. 6. Two defects, (a) starting position with Λ-field and G0, (b) paths with © – start position,△ – end position At last, three defects are investigated with their initial driving forces. The first arrangement in Fig.8a shows that the upper defect is attracting the lower ones,while the latter showa repulsive interaction between each other.Because of the grid dependency, the symmetry can not be assured, and the force vec- tors exhibit slightly asymmetric directions. The driving forces have different 834 O. Goy et al. Fig. 7. Two defects, (a) starting position with Λ-field and G0, (b) paths with © – start position,△ – end position magnitudes in these scenarios. That would lead to different defect velocities and the motion is calculated by time integration of (3.5). So, simulations will not be as easy as with two defects (where both defects experienced the same force). Fig. 8. Three defects with Λ-field and G0 A similar result can be observed with the constellation in Fig.8b. This behaviour is plausible, because C11 andC22 do not differ much, so both con- figurations are nearly equivalent. 7. Conclusion A micromechanical approach for describing point defects in a piezoelectric material has been presented in this paper. The coupling and an anisotropic material have been taken into account. The energy momentum tensor and configurational forces onpoint defects have beenbriefly introduced.A solution for periodic boundary conditions has been obtained by using a FD scheme followedbyFFTof coupleddifferential equations.Therefore, a regular gridhas been required, which was chosen to be quadratic. After the transformation, a Interaction of point defects in piezoelectric materials... 835 small system of equation on each grid point have had to be solved, followed by a numerical back transformation to obtain the final solution. Several examples have been presented, where an oxygen vacancy has been simulated.Theresultingfieldshavebeenshownandthe significanceofdifferent defect parameters hasbeenpointedout.On thebasis of thefields of onedefect, predictions could bemade for defect migration. Another point of interest has been the simulation of defect movement by means of configurational forces, where the interaction of defects with simple kinetic laws has been studied. It has been shown that two vacancies in a piezoelectric material often have the tendency to unify when they get close to each other. There are lines in and perpendicular to the poling direction, where the defects repel each other, nearly all other alignments lead to a unification of both vacancies. Beside the dependence on the relative position of the defects, the initial distance is important. Two defects attract each other, if they are initially close enough, and they separate to an equilibrium position in all other cases. So it can be deduced, that defects, brought close enough together, will formdefect clusters. Then these clusterswill naturally interactwith other accumulated defects, and so on. This can slow down a domain wall and cause the fatigue phenomena mentioned at the beginning. References 1. Barnett D.M., 1972, The precise evaluation of derivatives of the anisotropic Green’s functions,Phys. Stat. Sol. (b), 49, 741-748 2. Dobovs̆ek I., 1998,Stressfields aroundananisotropicpointdefect in anelastic crystal, Z. Angew. Math. Mech., 78, 351-352 3. Dobovšek I., 2000, Stress fields around an anisotropic point defect in a quasi- continuum, Z. Angew. Math. Mech., 80, 379-380 4. Gross D., Kolling S., Mueller R., Schmidt I., 2003, Configurational forces and their application in solidmechanics,Europ. J.Mechanics, A/Solids, 22, 669-692 5. GurtinM.E., 1995,Thenature of configurational forces,Arch.RationalMech. Anal., 131, 67-100 6. Gurtin M.E., 2000, Configurational Forces as Basic Concept of Continuum Physics, Springer, Berlin, NewYork, Heidelberg 7. Kienzler R., Herrmann G., 2000, Mechanics in Material Space, Springer, Berlin, NewYork, Heidelberg 8. Lupascu D.C., 2004, Fatigue in Ferroelectric Ceramics and Related Issues, Springer, Heidelberg 836 O. Goy et al. 9. Lupascu D.C., Roedel J., 2005, Fatigue in bulk lead zirconate titanate ac- tuator materials,Advanced Engineering Materials, 7, 10, 882-898 10. MauginG.A., 1993,Material Inhomogeneities in Elasticity, Chapman&Hall, London, Glasgow, NewYork, Tokyo,Melbourne, Madras 11. Mueller R., Gross D., Lupascu D.C., 2006, Driving forces on domain walls in ferroelectric materials and interaction with defects, Comp. Mat. Sci., 35, 42-52 12. SchclarN.A., 1994,Anisotropic Analysis using Boundary Elements, Compu- tational Mechanics Publications, vol. 20, Topics in Engineering, Southampton UK and Boston USA 13. Schrade D., Mueller R., Gross D., Utschig T., Shur V.Y., Lupascu D.C., 2006, Interaction of domain walls with defects in ferroelectric materials, submitted for publication 14. Vogel S.M., Rizzo F.J., 1973, An integral equation formulation of three dimensional anisotropic elastostatic boundary value problems, Journal of Ela- sticity, 3, 203-216 15. Xu Y., 1990, Ferroelectric Materials and their Applications, Elsevier Science Publ., Amsterdam Oddziaływanie defektów punktowych w materiałach piezoelektrycznych – symulacje numeryczne w kontekście zmeczenia elektrycznego Streszczenie Praca dotyczy zagadnienia zmęczenia elektrycznego materiałów funkcjonalnych, takich jak piezoelektryki używane do wyrobu czujników i elementów wykonawczych. Zmęczenie elektryczne powoduje degradację właściwości elektromechanicznych z ro- snącą liczbą przebytych cykli roboczych,w czasie których jonowe ładunki elektryczne lub elektrony oddziaływują ze sobąwewnątrz i na skrajudanegomateriału.To z kolei wpływa na lokalne pola elektryczne i mechaniczne sprzężone ze sobą efektem piezo- elektrycznym. Założonow pracy, że ładunki elektryczne (dziury) będą zamodelowane defektami punktowymi, któremają skłonność do gromadzenia się i tworzenia skupisk, zwłaszczaw pobliżu elektrod.W rezultacie, nagromadzone defekty zmieniają rozkład polaryzacji układu. Do ilościowego rozwiązania zagadnienia zastosowano dyskretyza- cjębadanegoukładusieciąkwadratowychkomórek inastępnieprzeprowadzonoszybką transformatęFouriera.Obliczeniawykonano na kilku przykładach. Szczegółowoprze- dyskutowano problem migracji pojedynczego oraz wzajemnego przyciągania dwóch defektów, by wreszcie zilustrować proces tworzenia skupisk przez całe gromady de- fektów punktowych. Manuscript received March 24, 2006; accepted for print July 10, 2006