**INTRODUCTION**

The modern theories of Continuum Damage Mechanics are based on the thermodynamics of irreversible processes, the internal state variables theory and the relevant physical considerations (e.g. the assumption of distributed microcracks, homogenization concept, the definition of micromechanical damage variable and plasticity-damage coupling mechanism, etc.). For electromagnetic continuous media, most of the literature is devoted to concrete or one-dimensional models and no rigorous three-dimensional model for general, anisotropic, media is proposed, to the author's knowledge. This applies, in particular, to ferroelectrics materials, a subject of wide practical interest.

Based on the Principle of virtual power (for finite velocities and without thermodynamical restrictions)^{[1]} and on our model developed in Bassiouny *et al*.^{[2]} the present work aims at deducing a three-dimensional model relying on rigorous. Thermodynamics, for the description of damage evolution in ferroelectric ceramics. This includes the important, interaction effects of the electric field and the damage.

The essential ferroelectric property is the possibility of reversal, or change in orientation, of the polarization direction by an electric field. The hysteresis loop is a manifestation of this key property.

Many previous studies in damage mechanics were concerned with linear elastic damage mechanics for brittle materials. For non-linearly elastic and elastoplastic, solids, nonetheless, these methods are not applicable in general. A review of the main theories of Damage Mechanics may be found in Ju^{[3]}.

Damage and fracture in materials are generally induced by nucleation of microcavities and their growth and coalescence into macroscopic cracks in materials. Hence the effects of these internal material defects may be reasonably discussed by dividing them into a single or finite number of macroscopic cracks and a set of distributed microcavities preceding these macroscopic cracks.

The distributed defects in materials not only lead to the crack initiation and the final fracture, but also induce the material deterioration (material damage), such as decrease of strength, rigidity, stability and of residual life. Although the problems of the mechanical behavior of these microcavities are no less important than those of microscopic cracks, it is only recently that a new discipline has been developed to investigate the growth of microcavities and the mechanical behavior of damaged materials (i.e., materials with internal defects) by representing the effects of distributed cavities in terms of certain mechanical variables and is called Damage Mechanics or Continuum Damage Mechanics^{[3-8]}. Thus, continuum damage mechanics not only describes the behavior of damaged materials, but also discusses the process leading to the initiation of macroscopic cracks, or the process of crack growth at crack tips and is investigated from a combined view point of material science and continuum mechanics^{[9]}.

From a microscopic view point, it would be quite difficult to formulate a microscopic model of the material and attempt to calculate from it the response for specific dynamic conditions. In many situations of interest, the ferroelectric material is in polycrystalline form, a fact that further complicates the microscopic formulation because the single crystal properties must be averaged over the crystalline orientation. Alternately, one could formulate a macroscopic model of the electro**mechanical properties** of an arbitrary ferroelectric ceramic so that, the response to specific dynamic conditions could be determined from the general electromechanical equations.

It is salient, at this point, to recall some basic facts, definitions and hypotheses.

Damage can be defined as a collection of permanent, microstructural changes concerning material thermodynamical properties brought about in a material by a set of irreversible physical microcracking processes resulting from the application of thermodynamical loadings.

The selection of a damage variable should be based on proper microstructural considerations. For several types of material microstructure, microcracks develop in characteristic patterns and the microstructure can be assumed to be statistically homogeneous. These patterned damages are observed in fibrous composite laminates. concrete and ceramics^{[3]}.

Following^{[4-11]} let __M__ be a fourth-order tensor which characterizes
the state of damage and transforms the Cauchy stress tensor __σ__ into
the effective stress tensor

Different representations of __M__ in terms of __σ__ and __D__
are proposed in the literature^{[9-11]}. For isotropic damage, the mechanical
behavior of microcracks is independent of their orientation and depends only
on a scalar variable. Accordingly, M/ simply reduces to (1-d) __I__, where,
__I__ is the fourth rank unit tensor, thus equation (1.1)
reduces to

the coefficient (I-d) dividing the stress tensor __σ__ is a reduction factor associated with the amount of damage in the material first introduced by Kachanov^{[12]}.

The value d = 0 corresponds to the undamaged state, d = d_{c} to the complete local rupture 0 <d_{c}< 1 and 0 < d_{c} < d, where 0 < d_{c} < 1 corresponds to a partially damaged state. Local stresses are redistributed to the undamaged material micro-bonds and therefore the effective stresses are higher than the nominal stresses. In addition, following hypothesis of strain equivalence is usually adopted.

“The strain associated with a damage state under the applied stress is equivalent to the strain associated with its undamaged state under the effective stress”.

Our model is motivated by assuming that, a volume element dV exists around an arbitrary point, __x__ in a body B containing distributed micro-voids. This volume element can be considered large enough to include a sufficient number of cavities or cracks, but is sufficiently small at the same time, so that the state of the stress, strain and void distributions are sufficiently uniform.

According to the notion of continuum damage mechanics, it is assumed that the effects of these microvoids at an arbitrary point __x__ can be described by an appropriate mechanical variable or damage field __D__(__x__) which we call damage variable or damage field^{[9,14]}. Since these damage variables are macroscopic variables, which describe mechanically the state of the internal structure, they are precisely internal state variables in the thermodynamic theories of constitutive equations.

**Three-dimensional evolution equations:** A phenomenological model capable of displaying electromechanical hysteresis effects and the evolution of the damage in ferroelectric materials can be constructed by assuming first specific dependence of the free energy 0 per unit volume on a set of state variables, including an internal vector variable __II__^{int} which describes an internal electric polarization^{[2]} and an internal second-rank tensor variable __D__(__x__) representing the damage field and including all damage effects^{[5,8]}. Thus, in the present model, the damage and the internal polarization” represent the only internal variables allowed. Several methods exist in the literature as to the way of choosing the damage variable^{[15]}. It is shown in^{[2]} that under some conditions, __II__^{int} is nothing but the residual polarization, up to a constant. Concerned with the poling process of ceramics, __II__^{int} represents the density of the “aligned” dipoles^{[15,16]} (3.23). As to the damage variable, it should be based on proper micromechanical considerations that take into account the distribution and evolution of different types of micro-defects^{[3]}. Our choice of the damage variable as a second rank tensor is motivated by the work of Murakami^{[8]} in which the anisotropic damage is shown to be represented by a second rank symmetric tensor.

We shall assume as in^{[2]} that the electric polarization __II__ is the sum of a reversible part __II__^{r} and a residual contribution __II__^{R}, that is:

Besides the assumption (2.1) of decomposition of the total polarization, we assume that the total strain tensor __ε__ is also the sum of a reversible, elastic-damage part __ε__^{e} and a plastic-damage part __ε__^{p}, assume that the deformation occurs in the frame of the theory of small perturbations), thus:

That is, __ε__^{e}(__ε__^{p}) includes not only the truly elastic (plastic) strain but also the added deformation due to active microcracks.

We shall consider the following dependence of the free energy ψ per unit volume on a set of (observable and internal) state variables and external parameters:

Where, θ is the absolute temperature measured from a reference temperature. We note that ψ should be concave in θ and convex in all other state variables __ε__^{e}, __II__^{r}, __II__^{int}, __D__.

We used the standard Cartesian tensor notation in a rectangular coordinate system. All considered fields are material fields written in the X_{K}- system and depending on the space-time variables (X_{K}, t). The superimposed dot denotes the material time derivative.

Accounting for the special form of the free energy (2.3) and the hypotheses (2.1) and (2.2), we may write the Clausins-Duhein inequality as^{[2]}:

wherein we have used the following definitions

Here, S is the entropy per unit volume, σ_{KL} are the components of the stress tensor, __E__^{L} is the local electric field^{[1,16,17]} representing the interaction between the electric dipoles and the mechanical structure of the solid and Q is the heat flux vector.

The Clausius-Duhem inequality (2.4) may be considered as a thermodynamical constraint imposed on the constitutive equations. Tensor Y_{KL} is the dual variable to the damage tensor and, therefore, is assumed to represent the damage strain energy release rate. In the absence of electric ordering effects and polarization inertia, __E__^{L} and the electromotive intensity ε are simply related by the intramolecular balance equation written here as^{[1,16,17,19-21]}:

This equation suggests a balance of forces acting on the charged particles forming the polarized material point of the medium.

We shall distinguish between dissipative processes that may depend on rates
and fluxes such as viscosity, electric relaxation and heat conduction and those
which may be referred to as instantaneous dissipative processes, in the sense
that they do not depend on the associated rates. Accounting for the thermodynamical
local-state axiom, the entropy S takes its usual thermostatic definition
and from (2.5) there remains the, following dissipative inequality:

where:

where, Φ_{1} and Φ_{2} refer to the two special classes
of dissipative processes referred to above. The state laws defined in (2.5)
must be supplemented by complementary laws for the fields that produce entropy
in Eq.(2.7). For these, we may postulate the existence of
dissipation potentials D_{1} and D_{2}, respectively homogeneous
of degree one and two in the rates, so that the complementary laws read

on one hand and

on the other.

For the first set (2.10), one introduces a dual dissipation potential D_{1}^{*}
via a Legendre Fenchel transformation^{[2]}, so that the evolution equations
will replace (2.10). Indeed, in analogy with rate-independent elastoplasticity
with hardening and magnetic hysteresis effects^{[20,21]} instead of
D_{1}^{*} one can directly introduce an electric and mechanical
loading function f (F_{1},.), F_{1} = {__σ__, __E__^{int},
__E__^{L}, __Y__} such that D_{1}^{*} becomes
the indicator of a convex set C = {F_{1} : f(F_{1}, .) <
0}, which defines the domain of reversible electric and mechanical processes
spanned by F_{1}. With this, Eq. (2.10) is replaced
by the following evolution equations:

where,
is called the “plastic multiplier”. It satisfies the following consistency
conditions

if we are on the

loading surface and

if we are in, or move to,

the elastic region.

Viscosity, electric relaxation and thermal conduction remain described by Eqs.
(2.11) with a dissipation function

If __σ__, __Y__, __E__^{int }and ε are functions of time at least derivable to the right, the Principle of maximal dissipation of Hill-Mandel yields the following orthogonality relation:

where, we have used the fact that __E__^{L} is the opposite of ε (the Maxwellian electric field). In the absence of plastic deformation and damage effects, an electric version of the celebrated. Drucker inequality^{[2,21,22]} can be obtained in the form:

on account of equation (2.8) and the convexity of the free
energy ψ. In its present form, equation (2.15) is useful
in determining the relation between the increments in the plastic strain, the
damage and the residual polarization. Careful measurements have to be carried
out as a verification of this relation.

The material coefficients are introduced by performing an expansion of the free energy in terms of its arguments. In doing so, a split into thermo-electro-elastic-damage contribution ψ_{d} and the energy stocked by the electrical internal variable ψ_{int} is achieved with for instance:

so that,

while the rest of equations (2.5) and (2.6)
remain unchanged. In the literature, it is some times taken that presence of
the damage variable in the free energy is just a multiplicative factor in the
form 1/1-d. The first of equations (2.18) clearly include
the effect of the reversible polarization on the strain energy release rate
by damage. An experimental verification is needed in order to evaluate the corresponding
coupling material coefficients.

When ψ_{int} is a homogeneous quadratic function of
, one obtains the following linear relations between the increments:

where:

Law (2.19) determines the increments in the internal electric field conditioned by incremental change in the internal polarization.

**Electromechanical couplings-loading surfaces:** Several examples of the loading functions have been studied in detail^{[2]}. The influence of stresses on these functions via various electromechanical couplings and both isotropic and kinematic hardenings were also studied.

Following^{[21-24]} where it was shown that if the loading surfaces act independently the normality rules (2.14) should be replaced by the following linear combinations, we write:

where, the undetermined factors
and
[i.e.; the “plastic” multipliers
(due to mechanical loadings) and
(due to electric loadings)] are non-negative. Thus they satisfy the consistency
conditions:

where, superscript q stands for the superscripts E or P.

We note that the normality rules (3.1)-(3.4) satisfy the orthogonality relation (2.15).

To study the influence of piezoelectric and electrostrictive couplings, let us assume the following expression for the actual fields:

**Piezoelectric couplings: **Piezoelectric materials have the property of becoming electrically polarized in response to an applied mechanical stress. This property known as piezoelectricity has an inverse, that is, when an electric field is applied, the material becomes strained. Thus the electric field created by the stress can be taken as:

where, P_{KLM} is a tensor of piezoelectric coefficients. These expressions
include the effect of damage within the effective quantities. Thus equations
(3.6) and (3.7) become

As seen from equations (3.10) and (3.11)
the damage will affect the actual fields through the effective stress tensor
and in the absence of the electric field the actual stress
in equation (3.11) reduces to the effective stress
.

**Electrostrictive couplings:** In the presence of centrosymmetry (no piezoelectricity), we may consider the following form for the coupled fields which account for the influence of electrostrictive couplings:

where,
is a tensor of electric susceptibilities and B_{KLMN} are the coefficients
of the electrostrictive tensor. Equations (3.6) and (3.7)
become

It may be noticed that we have combined electromechanical loadings surfaces, whereby it is proposed that there exist two loading functions f^{P} and f^{E} such that:

In analogy with elastoplasticity with hardening based on a Von-Mises yield condition, one may consider mechanical and electrical loading functions as follows:

where, E_{c} is a constant depending on the level of the electric loading (as also on the state of the stress)-it is called the coercive field (compare with the magnetic case^{[20,21]} and k_{0} is a positive constant (related to the material) specifying the level of stress at which plasticity occurs. It acquires its meaning clearly when A = 0.

The norms are defined as:

where,

Here (A_{KL}) is a positive-definite symmetric tensor which accounts for electric anisotropy: (B_{KLMN}) is a fourth order symmetric tensor and B_{KLMN} = B_{LKMN} = B_{KLMN} = B_{MNKL}.

Both the coercive electric field and the yield stress must be affected by the
actual fields given by equations (3.14) and (3.15).

Applying the normality rules (3.1)-(3.3), we obtain:

Integration of the later equation yields

This equation shows that the electric internal variable is nothing but the residual polarization up to a constant vector __b__ expressing a translation in the space of residual polarizations^{[2]}.

As a special case, one may choose the following form for the function A(__Y__, __D__) (compare with Maugin^{[24,25]}):

in such a way that (3.4) yields the evolution equation

In the case of pure electric loading function equation (3.22)
reduces to

where:

The plastic multipliers
and
are determined as usual by using the consistency conditions
and
We may assume that
where α is a scalar that represents the possible difference in the physical
dimension between
and
on account to the equation (2.19) we get

where:

where, h^{E} is the so-called “plastic hardening modulus” (due to electric loading) and < . > indicates the positive part.

On account of equations (3.29) and (3.30),
the evolution equation of the damage variable (3.26) becomes:

**Effect of damage on the material coefficients:** We search for an analytical expression for the thermodynamical potential, so that one may consider the simplest non trivial constitutive relations and rate laws sufficient to yield results comparable to experimental observations. For this purpose, we, used the following decomposition of the free energy:

We shall restrict our further considerations to a fairly thin specimen which is large in lateral extent and whose surfaces are stress free. We will assume the following dependence of the free energy on the various variables, indicating the influence of the damage on the material coefficients:

where, α and β take values 1, 2, 3,…,6 in Voigt’s notation.
Using the equations (2.5) we obtain the following constitutive
equations:

For the damage strain energy release rate due to an infinitesimal increase of damage at constant stress and temperature, we may write:

where, <<
θ_{0} < θ_{c}, with θ_{c}, is the
Curie temperature.
D_{α}, P_{kαβ}, D_{β}, τ_{αβ}D_{β},
C_{αβγ}D_{β}, Θk_{α}D_{α}
, are the changes in the inverse electric susceptibility, piezoelectricity,
thermoelastic coefficients, elastic coefficients and the change in the pyroelectric
vector and in the specific heat due to the effect of damage process and Γ
is the entropy per unit mass^{[26]}. The equation (4.6)
illustrates the effect of the electric polarization on the strain energy release
rate. In the absence of the electric polarization and temperature equation
(4.6) reduces to the following form^{[27]}:

As the constitutive equations (4.3)-(4.6)
include tensors of odd rank, they are valid for materials with no center of
symmetry. These equations constitute a fully coupled, three-dimensional thermodynamical
set of equations describing the effect of damage on ferroelectric ceramics.

**An application in elastoplasticity:** It is the purpose of this part to study the effects of damage and electromechanical couplings through a ferroelectric circular cylinder subjected to torsion and a longitudinal electric field.

The purely mechanical problem, that is the torsion problem, is a celebrated
one in elastoplasticity^{[22,23]}. We consider here a ferroelectric
cylinder of circular cross-section so that we have no warping. We also assume
that the Saint-Venant principle is verified. The traction on the lateral surface
are taken to be zero. The end of the cylinder is acted upon by a couple of moment,
Mz, along the OZ-axis.

The mechanical solution of this problem is well known and may be found in many places^{[24,25]}.

Consider now that the electric field is applied only in the direction of the axis of the cylinder so that only the electric behavior projected along this axis intervenes. Accounting for the fact that only one component, of the stress does not vanish, we can set the following definitions of physical quantities:

where, we have used the cylindrical system of coordinates (r,θ,z) and
n is a unit vector in the direction of the axis of the cylinder, εθ_{z}
is the only non-vanishing component of the stress and provided that. the symmetry
class of the material under consideration is the triclinic class^{[26]}
is the electric susceptibility constant.

For electric solids exhibiting electrostriction a coupling with the stress can be proposed so that, wc may consider the following form for the effective fields which take this type of coupling into account. In the one dimensional case we are led to:

In analogy with elastoplasticity with hardening based on a Von-Mises yield condition, one may consider mechanical and electrical loading functions f^{P}(σ, ε, A (Y, D)) and f^{E} (σ, ε, ε), respectively in one-dimensional form as follows:

where, we have replaced __E__^{int} by ε and used the definition (5.1) of the reversible part of the electric polarization.

The normality rules (3.1)-(3.3) assume the forms:

On account of equations (5.4) and (5.5),
the normality rules (5.6)-(5.9) become:

and the following one dimensional form of the energy release rate:

where, the plastic multiplier (due to electric loadings)
for instance, can be determined from the consistency condition
. One thus obtains:

while the electric hardening modulus h^{E} is given by

where:

Following Maugin^{[24,25]} one may assume that Q = Kε^{1/M} so that the electric hardening becomes:

the two coefficients K and M can be identified in a traction hardening test.

If α is the radius of the cylinder, M-the torsional moment, α-the torsion angle (i.e the twist per unit length of the cylinder) and c(α)-the radius of the plastic-elastic zone, then^{[15,16]}

where, and
is the value of the stress component σ_{θZ} in the pure mechanical
case (no electric field). Repeating the same steps as in^{[27]} but
for the damage case, we finally end up with the following formula for the torsional
moment, presently denoted M^{d}(α):

where
From equation (5.19) we note that the relative variation
in the estorsional moment due to damage is:

Equation (5.21) shows that Y should be negative. An experimental
measure of the torsional moments yields a numerical value for D. The sign of
the electrostrictive coefficient will play an important role in the value of
the moment M* and consequently the value of the relative variation.

Equations 5.10-5.13 clearly reflect the
effect of damage on the behavior of the cylinder. Comparison of these theoretical
considerations with experimental results should lead to the determination of
different parameters involved in the model.

**CONCLUSIONS**

A continuum approach was developed for the systematic modeling of a general, three dimensional anisotropic state of damage of ferroelectric ceramics.

The present work provides a consistent thermodynamical framework for analyzing the damage evolution and the electromechanical behavior of damage in ferroelectric materials. The paper focused on the influences of the damage process on the material coefficients of the considered ceramics, as well as on the influence of polarization on damage. Effects of different types of electromechanical couplings introduced through the effective quantities are studied. We obtain the additional polarization, stresses and entropy induced during the crack opening. Additional experimental verification of the model is still needed in order to evaluate the different coupling material coefficients.

As an application, the problem of a ferroelectric circular cylinder subjected to torsion and to a longitudinal electric field has been studied in one-dimensional case. We obtain the evolution equations governing the hysteresis behavior and the evolution of the damage of the ceramic material in one dimensional form.