Modelling of Electro-Viscoelastic Materials through Rate Equations

Models of dielectric solids subject to large deformations are established by following a thermodynamic approach. The models are quite general in that they account for viscoelastic properties and allow electric and thermal conduction. A preliminary analysis is devoted to the selection of fields for the polarization and the electric field; the appropriate fields are required to comply with the balance of angular momentum and to enjoy the Euclidean invariance. Next, the thermodynamic restrictions for the constitutive equations are investigated using a wide set of variables allowing for the joint properties of viscoelastic solids, electric and heat conductors, dielectrics with memory, and hysteretic ferroelectrics. Particular attention is devoted to models for soft ferroelectrics, such as BTS ceramics. The advantage of this approach is that a few constitutive parameters provide a good fit of material behaviour. A dependence on the gradient of the electric field is also considered. The generality and the accuracy of the models are improved by means of two features. The entropy production is regarded as a constitutive property per se, while the consequences of the thermodynamic inequalities are made explicit by means of representation formulae.


Introduction
Electro-viscoelastic dielectrics are a subject of recent research in connection with electrically sensitive materials, where an external electric field results in a large deformation. The coupling of electromagnetic fields in deformable materials has been developed within continuum mechanics; both the non-instantaneous response and the rate-dependent behaviour indicate that a viscoelastic scheme is in order. Various approaches to the mathematical modelling of electro-mechanical materials have been developed in the literature (see, e.g., [1,2] and refs therein). The following are some ideas applied for the modelling.
As is carried out in other domains of mechanics, in [3] the physical deformation gradient F is given a multiplicative decomposition into elastic F e and viscous F v parts; the Cauchy-Green tensor B and the stretching tensor D are decomposed accordingly. A restriction to electroelastic solids is developed in [4]; systematic use is made of the Maxwell stress while the electric field and the electric displacement are decomposed additively into applied and self fields. A nonlinear treatment of the interaction between electric field and deformation is established in [5] using systematically Lagrangian fields and letting the stress be decomposed into mechanical and ponderomotive parts while the ponderomotive part is the sum of polarization stress and Maxwell stress. An approach to electro-rheological materials is developed in [6] by accounting for dissipation through the stretching tensor D in a viscous-like scheme.
To establish a thermodynamically consistent scheme for dielectrics undergoing large strains, we follow some ideas that have lately been applied to hysteretic phenomena [7][8][9]. First, the entropy production is viewed as a constitutive property which has to be determined using a constitutive equation in accordance with the second-law inequality (see the systematic procedure developed in [10]). Further, consistent with the objectivity principle, the constitutive equations are required to be invariant under the group of Euclidean transformations. This leads to a selection of appropriate fields for the electric field and the polarization. Furthermore, representation formulae for vectors and tensors allow the maximal generality of the thermodynamic restrictions on the constitutive equations.
Owing to the deformation of the body, the representatives of the electric field and the polarization are found by requiring the Euclidean invariance and the validity of the balance of angular momentum; this, in fact, leads to a pair of Lagrangian vectors. Next, the thermodynamic consistency is developed in a way that allows for rate-type equations in several models: hypo-electroelasticity, electro-viscoelasticity, heat conduction (and electric conduction) in dielectrics, electroelasticity with dielectric memory, and ferroelectric hysteresis.
The aim of this paper is to show that the thermodynamic consistency, in addition to ascertaining the physical admissibility of the constitutive functions, is also a guideline to the setting of the material model, thus providing a simple scheme for the selection of the parameters characterizing the material behaviour. To validate the resulting models, it is necessary to compare the results with the experimental observations. The work is mainly devoted to theoretical aspects. Yet, in § 9.3, particular attention is paid to models for soft ferroelectrics, such as BTS ceramics (BaTiO 3 doped with 7.5 mol% Sn). Virgin loops of Ba(Ti,Sn)O 3 ceramics in dependence on the tin content are devised in ( [11], Figure 6). Our model is able to very well describe materials of this type.

Notation
Throughout, we denote by Ω ⊂ E 3 the time-dependent region occupied by the body. The motion is described by the function χ on R × R with R the reference configuration and R the set of real numbers. Hence, Ω x = χ(X, t), X ∈ R, and we let ∇ = ∂ x , ∇R = ∂ X . The deformation gradient F is defined as F = ∇R χ or F iK = ∂ X K χ i . The function v = ∂ t χ provides the velocity. A superposed dot denotes the total time derivative; for any function f (x, t), we haveḟ = ∂ t f + (v · ∇) f . To avoid obvious ambiguities, the Green-Lagrange strain tensor is denoted by E (instead of E) and the stretching tensor by D (instead of D).

Balance Equations
We consider a ferroelectric, deformable body where dissipative properties are allowed to occur of a mechanical and electric character. To simplify the description of material properties, throughout, the electromagnetic fields are considered at the frame locally at rest with the body.
Let ρ and ρ R be the mass densities in Ω and R, respectively; the balance of mass leads to Let T be the mechanical Cauchy stress tensor and b the mechanical body force. The equation of motion can be written in the form where f P is the force per unit volume of an electric character. In stationary conditions or in the approximation of a negligible magnetic field, we have f P = (P · ∇)E, where E is the electric field and P the polarization (per unit volume). The balance of angular momentum results in skw(T + E ⊗ P) = 0, where skwA = 1 2 (A − A T ) denotes the skew part of a tensor A and ⊗ represents the dyadic product of two vectors, (E ⊗ P) ij = E i P j . Let ε be the internal energy (per unit mass), L the velocity gradient, L ij = ∂ x j v i , and p = P/ρ. Moreover, let q be the heat flux vector, r the energy supply (per unit mass), and J the electric current density. The balance of energy is expressed as Let θ be the absolute temperature and η the entropy density. The statement of the second law is: the inequality holds for any process compatible with the balance equations. The entropy production γ is assumed to be given by a constitutive function. Consequently, the process consists of η, q, r, γ and the other functions occurring in the balance equations.
Using the Helmholtz free energy

Representation Formulae
Let N be a given second-order tensor, |N| = 1. Then, for any second-order tensor Z, we can write If Z ⊥ is unknown, then we can represent Z ⊥ in the form where I is the fourth-order unit tensor and G is an arbitrary second-order tensor. Once g = Z · N is given, we can write the representation formula Likewise, if z is a vector and we know the inner product z · n = f with a unit vector n, then we can represent z in the form where w is an arbitrary vector.

Euclidean Invariance and Objectivity
A change in frame F → F * , given by a Euclidean transformation, maps x → x * in the form Under the transformation (7), F and E change as vectors, As shown in ( [10], Ch.15), invariant scalars, vectors, and tensors may involve F and E. The right Cauchy-Green tensor C and the Green-Lagrange strain tensor E, are invariant [12]. Instead, and then L is non-objective, namely, it is not a tensor relative to Euclidean transformations. By the standard decomposition where D is the stretching tensor and W is the spin, we have Vectors and tensors under Euclidean transformations are said to be objective; D is an objective tensor, W is not objective. Let be the second Piola stress. SinceĖ = F T DF then The referential heat flux and temperature gradient are invariant and so is the power In connection with the electric field E and the polarization P, we can consider the fields E E E = F T E, P P P = JF −1 P; E E E and P P P are the (Lagrangian) fields in reference configuration [13,14]. The fields E E E and P P P are invariant, P P P * = J * F * −1 P * = JFQ T QP = P P P; since J is invariant, so are the fields J p E E E and J q P P P for any p, q ∈ Z. Moreover, ∇R E E E is an invariant while the scalar |∇E E E| 2 is an invariant scalar. The time derivative of invariants is invariant too and, hence,Ė E E,Ṗ P P,Ṫ RR ,q R are invariants. For later purposes, we represent the power ρE ·ṗ in terms of E E E and P P P. Since ρ R = ρJ and P = J −1 F P P P then It then follows that ρE ·ṗ = (E ⊗ P) · L + J −1 E · FṖ P P.

Constraint on the Stress Tensor
The fields E E E and P P P enjoy Euclidean invariance. We now examine the validity of the constraint (1) which accounts for the balance of angular momentum. Consider the Clausius-Duhem inequality in the form (4) and observe Let θ, F, E, ∇θ be the set of variables for the functions φ, η, T, q, γ; of course, φ, η, and γ depend on F, E, and ∇θ through appropriate invariants. Substitution ofφ yields Owing to the arbitrariness of∇θ,θ and L,Ė, it follows that The constraint (1) results in and the requirement (13) To determine admissible relations, consider any objective fieldẼ of the form f (J)E E E. Hence, we let φ depend on F through E = (F T F − 1)/2 and jointly on F and E throughẼ; Since where Hence, using (14) and (15), it follows that the requirement (13) holds if and only if This relation holds identically for any electric field In light of expression (12) of the power, the pair P P P, E E E seems more suitable to describe the electric behaviour in deformable bodies. That is why we then proceed with the choice of E E E, i.e., f = 1, for the referential electric field. Consequently, we have We assume J is a vector and then we might consider J = F −1 J as an invariant current density, Incidentally, J J J equals the referential flux, say, J R , as for any vector field, such as q R = JF −1 q.

Thermodynamic Analysis
Motivated by the Euclidean invariance, we now investigate the Clausius-Duhem inequality (4) in the Lagrangian description. Hence, we consider J times inequality (4) and use the representations (8), (9), (10), and (16) of the powers T · L, q · ∇θ, ρE ·ṗ, and J · E to obtain Hereafter, we use the referential fields For later developments, it is convenient to consider the free energy Moreover, to save writing, we let Using (11) and the definition of T RR , we have Equation (17) is then written as

Constitutive Assumptions and Thermodynamic Restrictions
Viscoelasticity is a scheme that accounts for a persistent rate of the response under a constant action. This may suggest that we allow for rate equations of T, J, P, and q. Accounting for rate properties is also consistent with the modelling of hysteresis (in ferroelectrics). Thus, we might take (θ, F, E, P, ∇θ, T, q,Ḟ,Ė), as the set of independent variables. However, the dependence on the derivatives can occur only in an objective way. In particular, the function φ R can depend only on Euclidean invariants. Hence, we let and the like for η R and γ. The viscoelastic character is realized by lettingṪ RR ,q R ,Ṗ P P, and J J J be given by constitutive functions of The time derivative of φ R is computed and substituted in (20) where γ ≥ 0. The linearity and arbitrariness of ∇Rθ,Ë,Ë E E,θ, W implies that The symmetry condition in (22) is just the balance relation (1) of angular momentum. Thus, it follows that T RR ∈ Sym. Hence, (21) simplifies to Notice thatĖ,Ė E E, and ∇R θ can take independent values and, so far,Ṫ RR ,Ṗ P P, andq R are functions of the whole set of variables Γ. This allows the possibility of cross-coupling effects, which are usually negligible. We then consider models arising from independent entropy productions, namely, the entropy productions γ T , γ E , γ J , γ q being non-negative while It is worth remarking that the decomposition (24)-(27) of the entropy production does not hinder joint mechanical-electrical-thermal effects due to the dependence of the constitutive functions on Γ.

Electroelastic and Hypo-Electroelastic Models
As to (24), we start with the simple case ∂ T T T RR φ R = 0 and γ T = 0. Hence, we have If, instead, γ T > 0 then we may take where Ξ is a fourth-order positive definite tensor. Hence, thus showing a symmetric term ∂ E E E φ R , a non-symmetric dielectric term C −1 E E E ⊗ P P P, and a vis- a generalization of the Kelvin-Voigt viscoelastic model with E the fourth-order elasticity tensor. A different scenario follows from (24) when and, hence, we can determine the expression ofṪ RR using the representation formula (5) with for any second-order tensor function G of Γ. Let H be any non-singular fourth-order tensor function of Γ, deprived ofĖ, and choose G = HĖ. Then, we havė where E denotes a family of (possibly non-symmetric) fourth-order stiffness (or elastic) tensor functions parameterized by H. Equation (29) ascribes to T RR a hypoelastic character possibly parameterized by E E E and P P P. The corresponding stress then consists of the electroelastic dyadic product −C −1 E E E ⊗ P P P and the hypo-elettroelastic stress obeying the rate-type Equation (29). A special but significant class of hypo-electroelastic models is obtained by assuming that ∂ E E E φ R = 0. Consequently, (29) simplifies to We can look for H = 0 such that holds identically. If this is so, it follows from (31) that E = H, thus eliminating the dyadic term and thenṪ For definiteness, we now assume the free energy in a quadratic form with respect to T RR . Let G 0 = G T 0 be a non-singular fourth-order tensor and set Due to the symmetry of G 0 and the arbitrariness of T RR , we conclude that H = G 0 so that (33) becomesṪ RR = G 0Ė . If, further, we assume G 0 is a constant tensor, then an integration allows us to recover the linear electroelastic constitutive relation is an arbitrary initial value of the hypo-elettroelastic stress. Owing to (30), it follows

Electro-Viscoelastic Models
If ∂ T T T RR φ R = 0 and γ T > 0 then Formula (5) for any second-order function G. We can also writeṪ RR in the forṁ The corresponding stress then consists of the electroelastic dyadic product −C −1 E E E ⊗ P P P and the electro-mechanical symmetric stress obeying a rate-type equation in the form (35) for a given G. Several viscoelastic schemes follow from (34) depending on the choice of φ R and γ T . Solids are characterized by a stress dependence such that, asymptotically, T RR = G ∞ E with G ∞ , a positive-definite fourth-order tensor. Define where β > 0, while A and Λ are positive-definite fourth-and second-order tensors.
Notice that and let Hence, the representation formula (35) yieldṡ Consequently, upon choosing Equation (36) shows that T RR − G ∞ E evolves with a relaxation time τ that is influenced by the temperature and the electric field, At constant strain,Ė = 0, the solution T RR to (36) asymptotically is as expected for a solid model. Finally, we note that, upon letting G 0 = A −1 + G ∞ , we can write Equation (36) in the forṁ thus obtaining the model of the standard linear solid.

Heat Conduction in Dielectrics
Equations (26) and (27) have the same structure and we then restrict attention to (27). We look at the evolution of q R and observe that the sought functionq R is subject to (27). This means that the thermodynamic restriction is confined to the inner product ∂ q R φ R ·q R . Based on the representation formula (6), we consideṙ q R = (q R · n)n + (1 − n ⊗ n)w.
Letting ∂ q R φ R = 0 and choosing n = ∂ q R φ R /|∂ q R φ R | using (27), we havė Hence, the possible dependence of γ q and ∂ q R φ R on T RR and E E E allows the description of stress and electric field effects on heat conduction.
To show the flexibility of Equation (37) in the elaboration of (thermodynamically consistent) models, we show, e.g., how to obtain a Maxwell-Cattaneo equation. Let where the dots denote quantities independent of q R while λ, κ, µ may depend on θ, E E E and T RR . Since γ q ≥ 0, we let κ > 0. Substitution into (37) yieldṡ The two terms with q R · ∇R θ cancel by letting µ = −1/λθ. With this value of µ we find the equationq which is in the Maxwell-Cattaneo form with relaxation time τ = λθκ.
In stationary conditions we have the Fourier-like equation the positive value of κ required by γ q ≥ 0 implies that the conductivity is positive.

Electroelastic Materials with Dielectric Memory
Equations (24) and (25) have the same mathematical structure. Hence, we can show that (25) allows hypo-electroelastic models as well as descriptions with memory for the past.
For formal simplicity we ignore the dependence on the electric current J J J and the heat flux q R . Let φ R := ρ R φ(θ, E, T RR , E E E, P P P). Choosing γ E = 0, Equation (25) reduces to (∂ E E E φ R + P P P) ·Ė E E + ∂ P P P φ R ·Ṗ P P = 0.
Assume ∂ P P P φ R = 0. Then, the representation formula (6) can be applied by letting n = ∂ P P P φ R /|∂ P P P φ R | to obtaiṅ We now restrict our attention to constitutive equations with linear dependence onĖ E E. Hence, we let Ξ be any non-singular second-order tensor function of Γ, deprived of ∇ R θ,Ė E E, andĖ, and consider w = ΞĖ E E. It follows thaṫ Accordingly, can be viewed as a family of (possibly non-symmetric) permittivities in the form of second-order tensor-valued functions parameterized by Ξ. A particular model is obtained assuming that the Helmholtz free energy ψ is independent of P P P, and then ∂ P P P φ R = −E E E, so that We can choose Ξ, and, hence, the functionΞ(θ, E, E E E, P P P, J J J, q R ), such that holds identically. Consequently, the dyadic term vanishes and Equation (38) simplifies tȯ

P P P =Ξ(θ, C, E E E)Ė E E.
Further, assume φ R in the quadratic form Since ∂ E E E φ R = ΣE E E − P P P then Equation (40) reads The validity of this relation for every vector E E E implies thatΞ = Σ. The paraelectric rate equation follows by letting Σ = 0 χ R , where 0 is the permittivity of free space and χ R = JF −1 χF −T is the electric susceptibility tensor in the material description. We then leṫ P P P = 0 χ RĖ E E.
As an example, consider a transversely isotropic material with easy axis m. The spatial electric susceptibility is then written in the form where χ and χ ⊥ are the electric susceptibilities in the direction parallel and perpendicular to m. Hence, it follows that If, instead, γ E > 0 then Equation (25) allows us to determine a large variety of dissipative electroelastic models. Assume ∂ P P P φ R = 0 and let n = ∂ P P P φ R /|∂ P P P φ R |. The representation formula (6) yieldṡ Applying Equation (41), we now establish two relevant classes of models, dielectrics with memory and hysteretic dielectrics.

Dielectrics with Memory
They are characterized by a polarization dependence that shows relaxation and that, asymptotically, is given by P P P = 0 χ ∞ E E E, with χ ∞ , where χ ∞ is a positive-definite secondorder tensor called relaxation susceptibility. To model this feature, we consider the free energy function where A and χ R are positive-definite second-order tensors while where α is a positive parameter possibly dependent on temperature and the scalar invariants of E and T RR . Consequently, and then For ease in writing, we define Hence, the representation formula (41) can be written in the forṁ A simple model arises by letting the dyadic term vanish. This happens by choosing (43)

then Equation (43) can be rewritten aṡ
which shows a time rate of P P P − 0 χ ∞ E E E with relaxation time 1/α. In essence, the models are characterized by the free energy φ and the entropy dissipation γ E . It is of interest to show that, while maintaining the same free energy, we can model hysteretic phenomena by letting γ E be proportional to |Ė E E| (or |Ṗ P P|).

Ferroelectric Hysteresis
Consider evolutions where only E E E and P P P are time-dependent while the remaining variables are constant. Restrict attention to cyclic processes in the time interval [t i , t f ] with (E E E, P P P)(t f ) = (E E E, P P P)(t i ), whencė Integration in time of (25) yields Hysteretic properties are now investigated by letting Hence, P P P and E E E are subject to Based on (41) and (45), we now determine a relation forṖ P P. Let Λ be a second-order tensor. By selecting w = ΛĖ E E, we havė A simple particular case emerges by letting Λ such that thus implying the vanishing of the dyadic term. Equation (47) shows that a linear relation is required between ∂ E E E φ R and ∂ P P P φ R . In terms of the free energy (42), the requirement (47) reads which implies Λ = 0 χ R . Hence, Equation (46) becomeṡ

One-Dimensional Models of Hysteresis
Suppose that the body is transversely isotropic and the spatial fields E and P are collinear in the direction m of easy polarization. Let (e 1 , e 2 , e 3 ) be an orthonormal basis where e 1 = m and assume E = Ee 1 , P = Pe 1 . The deformation gradient is taken in the form Hence, we have Consequently, E = (1 + ν)E and P = (1 − δ) 2 P are subject tȯ Owing to the transversely isotropic symmetry of the material, we let so that, in matrix form, χ = diag(χ , χ ⊥ , χ ⊥ ). Hence, it follows that Further, let Ae 1 = αe 1 .
The one-dimensional version of (48) iṡ where P(E) = ( 0 χ + 1/α)E. Except at inversion points (whereĖ = 0), we can divide bẏ E to obtain dP If ζ is independent ofĖ or depends on it at most through its sign, then (49) is invariant under the time change t → t * = ct, c > 0. We can, therefore, say that the equation is rate-independent. By letting ζ = 0, Equation (49) reduces to where the right-hand side represents the differential permittivity of a paraelectric/dielectric material. In general, letting we can write Equation (49) as a differential equation, for the unknown function P(E). Hysteretic effects are described by the second term, 2 sgnĖ, in that the slope of P(E) changes depending on the sign ofĖ. Since 2 is proportional to ζ, the vanishing of entropy production γ E implies 2 = 0 and the corresponding model describes non-dissipative hypo-dielectric materials. Consequently, 1 , which possibly depends on the values of P and E, represents the slope of the paraelectric curve P(E). When 2 = 0, we can view (51) as the differential electric permittivity. Since it is usually observed to be non-negative, we then require that according to the counterclockwise sense required by P dE ≤ 0. Summarizing the above analysis, we conclude that the model is characterized by three quantities: the paraelectric permittivity 1 = 0 χ , the hysteretic function ζ and the temperature-dependent function α. Since 1 is fully determined by the free energy φ R whereas 2 depends also on ζ, different models can be obtained starting from the same function φ R . As we will see in a while, the function 2 , which is connected with the entropy production through ζ, governs the hysteretic properties of the material.
It is worth considering the case where α 0 > 0 depends on the characteristic parameters of the material, such as θ C and χ . Since We conclude that, regardless of the form of ζ, as θ → θ C the curve P = P(E) represents the polarization curve of a paraelectric material.

Ferroelectric Soft Materials
We present a simple example of the theory that appears appropriate for materials with "ferroelectrically soft" behaviour, such as BTS ceramics (BaTiO 3 doped with a small amount of tin). They are distinguished by relatively high domain mobility and, thus, relatively easy polarization. The hysteresis loop of a soft material is, therefore, characterized by low coercive field strength and high spontaneous polarization.
The vanishing of χ as |E| approaches infinity is a way of modelling the saturation property of hysteresis. Starting from different initial states (E 0 , P 0 ), hysteresis cycles are obtained in the E-P plane by the systeṁ Since the model considered here is rate-independent, the hysteresis loops are independent of the frequency ω at which the alternating electric field varies. In Figure 1, we choose and α = ζ 0 = 2/3, so that τ h = 1 and 1 := 0 χ (E) = 3/ cosh 2 (2E). Different amplitudes A E are used in order to highlight the shape of switching (A E = 1.6) and nonswitching (A E = 0.4) ferroelectric cycles.  Figure 6). Their ferroelectric behavior changes from hysteretic (7.5 mol% Sn) to non-hysteretic (15 mol% Sn) polarization loops. Our model is able to very well describe materials of this type assuming that α is proportional to the molar content of tin. Indeed, as α increases the saturation limit is lowered and the hysteretic parameter τ h tends to vanish. In particular, Figure 1 corresponds to Ba(Ti,Sn)O 3 with 7.5 mol% Sn.

Dependence on the Gradient of the Electric Field
Spatial interaction in non-uniform electric fields are modelled by allowing for a dependence on the gradient of the electric (or the polarization) field [15,16]. For simplicity, we neglect heat conduction. We then let As we show in a moment, the dependence on ∇R E E E is allowed only if an extra-entropy flux occurs. The dependence on the second gradient ∇R ∇R E E E is considered so that all of the constitutive functions depend a priori on the same set of variables.
The Clausius-Duhem inequality (21) is modified to The arbitrariness ofË,Ë E E, ∇R ∇RĖ E E,θ, W implies that The possible dependence of k R onĖ E E allows us to write with the dots denoting terms independent of ∇RĖ E E. Hence, the arbitrariness of Consequently, to within terms independent ofĖ E E. For definiteness and simplicity, we let (55) hold exactly for k R . In view of (55), we have Hence, in light of (54), substitution of ∇R · k R in (53) yields where We notice that if, e.g., That is why we started with a possible dependence on ∇R ∇R E E E. As with (24)-(27), we notice that (56) holds if The major interest in the dependence on ∇R E E E is given by (57) in view of the variational derivative δ E E E φ R . Hence, we apply the representation (6) to (57). Letting n = ∂ P P P φ R /|∂ P P P φ R |, we have In the particular case of no dissipation, γ E = 0, and apart from the arbitrary contribution w, the representation shows a hypo-electric behavior, Here, though, the tensor Λ is also affected by ∇R E E E and ∇R ∇R E E E.

Relation to Other Approaches
Hysteresis models of ferroelectric materials have been established through various approaches. Some models are based on hysteresis operators. This is, e.g., the case in [17] (see also [18,19]) where a Preisach operator P is considered in a one-dimensional setting, with a potential U dependent on an auxiliary state function q(ε, E), ε representing the strain and E the electric field. It is of interest that the assumed free energy involves a joint dependence on ε and E via U (E/ f (ε)), which is quite analogous with the dependence of ζ on [P − P(E)] 2 .
Thermodynamically consistent models are based on the compatibility with the second law of thermodynamics expressed by the Clausius-Duhem inequality. Yet the schemes adopted, and the mathematical procedure, are quite different to the present one in that they involve additional internal variables [20][21][22][23]. In [20], the thermodynamic potential, the enthalpy H, is a function of the strain ε, the electric field E, and an internal variable ξ. The strain ε and the electric displacement D are decomposed additively in reversible parts ε e , D e and remanent parts, ε r , D r , with D r being equal to the polarization P r . The Clausius-Duhem inequality is written in the form where σ is the Cauchy stress. Next, ε, E, ε r , P r are taken as the variables and the inequality is splitted in σ = ∂ ε H, D = −∂ E H, ∂ ε r H ·ε r + ∂ P r H ·Ṗ r ≤ 0, as though ε, E, ε r , P r were independent. The internal variable ξ is identified with the pair ε r , P r and, hence, the evolution of ξ is subject tõ σ ·ε r +Ẽ ·Ṗ r ≥ 0, whereσ = −∂ ε r H,Ẽ = −∂ P r H. By appealing to the principle of maximum remanent dissipation and introducing a non-positive switching function Φ(σ,Ẽ), it is found thaṫ ε r = λ∂σ Φ(σ,Ẽ),Ṗ r = λ∂ẼΦ(σ,Ẽ); the loading/unloading conditions are λ ≥ 0, Φ(σ,Ẽ) ≤ 0 and λΦ(σ,Ẽ) = 0. Ref. [21] models ferroelectric materials via the deformation gradient F, the electric field E, the temperature θ and an internal variable Q. In terms of the free energy Ψ it is found that where T R is the first Piola stress, and β ·Q ≥ 0, β := −ρ R ∂ Q Q Q Ψ.
The evolution of Q is then assumed to be governed by a dissipation potential function Φ(Q) so that β = ∂Q Q Q Φ and the thermodynamic consistency is claimed to be achieved if Φ is a convex function.
The recourse to the function Φ, to characterize the evolution of the internal variable, may be viewed as the analogue of the present modelling directly via the entropy production γ. The advantage of our procedure is the consistent thermodynamic analysis via the Clausius-Duhem inequality. Furthermore, we avoid any decomposition of the pertinent fields, thus making the connection with experimental results more immediate.

Conclusions
This paper develops models of deformable dielectric solids. The models are quite general in that they account for viscoelastic properties and allow electric and thermal conduction. The procedure involves some conceptual points. First, the appropriate electric field and polarization are selected on the basis of two requirements: the fields are required to comply with the balance of angular momentum, skw(T + E ⊗ P) = 0, and enjoy the Euclidean invariance. Among the possible choices we considered E E E = F T E, P P P = JF −1 P as the invariant electric field and electric polarization. Next, we investigated the thermodynamic requirements for constitutive equations involving θ, E, T RR , E E E, P P P, J J J, q R , ∇R θ,Ė,Ė E E, where E is the Green-Lagrange strain, T RR = T RR + C −1 E E E ⊗ P P P, C = F T F, and J J J is the referential current density. The resulting equations (24)-(27) are examined using the representation formulae and this allows us to find the hypo-elastic and hypoelectric behaviour. Further, the entropy production γ is considered as a constitutive function, thus leading to models of hypo-electroelasticity, electro-viscoelasticity, heat conduction (and electric conduction) in dielectrics, electroelasticity with dielectric memory, and hysteretic dielectrics. As a sigificant generalization, often considered in the literature [16], we show the properties induced by considering the gradients ∇R E E E, ∇R ∇R E E E among the independent variables.