Finite element analysis of materially uniform dielectric elastomers

Abstract

In the definition of Noll, a body is uniform if all points are made of the same material. As shown by Noll himself and by Epstein and Maugin, uniformity makes the Helmholtz free energy depend on the material point exclusively through a tensor field, called uniformity tensor or implant tensor or material isomorphism. Uniformity is therefore a particular case of inhomogeneity. In turn, uniformity includes homogeneity as a particular case: indeed, homogeneity is attained when the uniformity tensor happens to be integrable. This work focuses on the non-linear large-deformation behaviour of uniform dielectric elastomers. Building on the foundational works of Toupin, Eringen and others, this work integrates continuum mechanics with electrostatics to develop a finite element framework for analysing uniform dielectric elastomers. This framework allows for considering the inherent inhomogeneity in materials exhibiting non-linear electromechanical coupling such as electro-active polymers. The inhomogeneity is assumed to be self-driven, i.e., not implied by the second law of thermodynamics: rather, it depends on the torsion of the connection (covariant derivative) induced by the uniformity tensor. A MATLAB^®-based finite element solver is developed and applied to the simulation of an electromechanical beam-type actuator. The solver is robust and capable of addressing various simulation scenarios. Numerical simulations demonstrate the significant impact of material uniformity on actuator performance. This research provides a tool for future applications in dielectric elastomers, particularly in sensors, actuators and bio-inspired robotics.

1. Introduction

Electro-elastic materials react to both mechanical forces and electric fields. When subjected to an electric field, they experience deformations and, when they are subjected to mechanical stress, their state of polarisation changes. A prime example of this phenomenon, in the small deformation range, is piezoelectricity, where electrical field and mechanical deformation are proportionally related. Other materials, like ferroelectrics and electro-active polymers (often abbreviated into EAPs), exhibit higher-order coupling effects between electrical and mechanical properties. Electro-active polymers are labelled as smart materials because of their versatility: they can be employed, e.g., in actuators, sensors and bio-inspired robotics [1]. Their lightweight nature, flexibility and efficiency in power usage give them an edge over traditional actuating materials [2]. With ongoing advancements in this area, electro-active polymers are anticipated to lead future material science innovations and engineering applications [3].

In this study, we focus on materially uniform dielectric elastomers and the simulation of their behaviour with the finite element method.

The theory of material uniformity was originally proposed in the seminal paper by Noll [4] and has been further refined through the influential contributions of Epstein and Maugin [5–7] and Epstein and Elżanowski [8]. This concept posits that all points in a uniform body are made of the same material: this implies that the Helmholtz free energy of the body depends on the material point exclusively through a tensor field, called uniformity tensor or implant tensor or material isomorphism, and representing anelastic deformations. Noll [4] also demonstrated that the torsion of the linear connection induced by the uniformity tensor $\mathit{K}$ indicates the presence of inhomogeneities in the body. To clarify, given a connection (also known as a covariant derivative) with Christoffel symbols ${\mathrm{\Gamma }}_{\mathrm{CB}}^A$ , the torsion is the third-order tensor $\mathbf{Tor}$ with components ${{\mathrm{Tor}}^A}_{\mathrm{BC}}={\mathrm{\Gamma }}_{\mathrm{BC}}^A-{\mathrm{\Gamma }}_{\mathrm{CB}}^A$ . For Noll’s uniformity tensor (material isomorphism) $\mathit{K}$ , the Christoffel symbols are expressed as ${{\mathrm{\Gamma }}^A}_{\mathrm{CB}}={{\mathrm{K}}^{\mathrm{A}}}_{\alpha }({{{\mathit{K}}^{-1})}^{\alpha }}_{\mathrm{C},\mathrm{B}}$ .

Epstein and Maugin [9] related the theory of material uniformity to the theories of plasticity, growth and anelastic behaviour in general, underscoring the importance of material inhomogeneity. The theory of uniformity has been related to the classical multiplicative decomposition of the deformation gradient [10–12], with applications to anisotropic yield criteria [13], or the problem of materials with Eshelby-like [14] inclusions [15]. Ciarletta and Maugin [16] contributed a biomechanics-focused perspective, introducing a finite strain second-gradient theory of uniformity for material growth and remodelling, integrating complex structural changes and mass transport through Eshelby-like stress tensors. Hamedzadeh et al. [17] employed the theory of uniformity to describe fibre recruitment and remodelling in biological tissues, with an example on the mechanics of the arterial wall.

The foundations of the theory of dielectric elastomers were laid by Toupin [18–20], who pioneered the integration of continuum mechanics with electrostatics and electrodynamics to develop a non-linear theory for polarisable solids. Later, Eringen [21] studied dielectric elastomers within the framework of the principle of virtual work. This trailblazing research set the stage for advances in large-deformation electro-magneto-mechanics, as evidenced by subsequent pivotal works from researchers like Lax and Nelson [22], Nelson [23] and Maugin [24]. The evolution of this field continued with contributions from Eringen and Maugin [25], Kovetz [26], Dorfmann and Ogden [27–29], Vu et al. [30], Steigmann [31], Trimarco [32,33], Ogden and Steigmann [34] and Maugin [35]. Bassiouny et al. [36,37] and Bassiouny and Maugin [38,39] introduced the first macroscopic model for the irreversible behaviour of ferroelectrics that is consistent with thermodynamics. These constitutive frameworks primarily rely on a phenomenological approach. However, taking a different route, researchers like Thylander et al. [40,41], Cohen et al. [42] and Thylander [43] have put forth electro-elastic models inspired by micro-mechanics. These large-deformation models consider the intricate dynamics of polymeric chain orientations and deformations under combined mechanical and electromechanical stresses and offer physically relevant material parameters.

Moreover, electro-active polymers have been studied in non-linear electro-elasticity or hyperelasticity contexts [28,44–46]. Bustamante et al. [47] introduced various electromechanical energy models and a comprehensive variational framework, comparing it with the findings by Toupin [18] and Ericksen [48]. They also examined the interaction of electric fields with deformable materials in quasi-static settings, focusing on stress measures, virtual work, variational principles and the suitability of total energy density function for electro-active elastomers.

Dorfmann and Ogden [49] focused on creating a non-linear framework that outlined the linearised behaviour of electro-elastic materials undergoing significant deformation in an electric field. Their findings emphasised that stability is closely tied to the values of the electromechanical coupling parameters within the constitutive equation. Building on the theories of electro-active polymers, Skatulla et al. [50] expanded these concepts to generalised continuum theories. This advanced continuum approach offers a method to consider scale effects. Dorfmann and Ogden [27,28], Bustamante [51], Sfyris et al. [52] and Dorfmann and Ogden [53] introduced models focused on isotropic electro-elasticity, largely built upon coupling invariants. Recognising the distinct inhomogeneous behaviour in particle-filled electro-active polymers, researchers such as Hossain and Steinmann [54] and Sharma and Joglekar [55] proposed mathematical frameworks of electro-elasticity that include transverse isotropy and dispersion-type anisotropy.

Variational methods have been developed and incorporated into the finite element method by Vu and Steinmann [56–58] to accurately model electro-elastostatics under large strain conditions. Saxena et al. [59,60] introduced an electro-viscoelastic model that elegantly separates the electric field into equilibrium and non-equilibrium components. Büschel et al. [61] developed a numerical method to simulate dielectric elastomers by taking into account both their electro-elastic properties and their time-dependent viscoelastic behaviour. In their model, they decomposed the deformation gradient into separate elastic and viscous components. However, they only allowed the mechanical component to evolve over time, while the electrical quantity remained time-independent. Ask et al. [62,63] focused on examining the electro-viscoelastic characteristics of polyurethane. They established a comprehensive framework for an electro-viscoelastic material and implemented it into finite elements.

Additionally, Ask et al. [64] presented a mixed finite element method specifically crafted to tackle the issue of volumetric locking observed in nearly incompressible materials. This innovative approach was effectively utilised to resolve the inverse motion problem. Similarly, Jabareen [65] proposed a mixed finite element method aimed at mitigating the volumetric locking challenges inherent in materials with near-incompressibility. Sharma and Joglekar [55] developed a finite element model of non-linear anisotropic dielectric elastomer actuators, building upon existing models and presenting a staggered solution for coupling non-linear equations. Kadapa and Hossain [66] introduced a finite element framework for computational electromechanics under finite strains, uniquely extending a two-field mixed displacement-pressure formulation to include coupled electromechanics.

In this work, we present a comprehensive study on the non-linear behaviour of dielectric elastomers, with a special focus on materially uniform dielectric elastomers. Within the scope of this research, the uniformity (which, we recall, is a particular case of inhomogeneity) of the material is modelled by a self-driven evolution law that accounts for non-linear electromechanical coupling, which is not implied by the second law of thermodynamics. Rather, uniformity is described by the torsion of the connection induced by the uniformity tensor $\mathit{K}$ , simulating the distortions present in real materials.

This view forms the cornerstone of our finite element framework, specifically developed to tackle the complexities arising from the inhomogeneity implied by the hypothesis of material uniformity and the non-linear electromechanical interactions observed in smart materials, such as electro-active polymers. The resulting MATLAB^®-based finite element solver, demonstrated through the analysis of a beam-type actuator, underscores the profound impact of the inhomogeneity brought about by the material uniformity, described by the uniformity tensor $\mathit{K}$ , on actuator performance. This insight into the inhomogeneous properties of smart materials provides a robust tool for future explorations and practical applications in the domain of dielectric elastomers, particularly in advanced fields like bio-inspired robotics, sensors and actuators.

Section 2 recalls basic continuum kinematics, the concept of connection (covariant derivative) and torsion of a connection, the definitions of fourth-order tensors that are relevant to constitutive theory and Maxwell’s equations for the case of electrostatics. Section 3 introduces the concept of material uniformity, in general, and for an electro-elastic body. Section 4 introduces the governing equations through the principle of virtual work and the dissipation inequality for uniform electro-elastic bodies, as well as the finite element formulation and the specific constitutive equations used in this work. The algorithm written in MATLAB is elucidated in section 5. Section 6 reports the finite element solution of our model for a benchmark problem: the electrically activated bending of a beam-type actuator. Finally, section 7 summarises our findings and their significance in advancing electro-elastic materials and their practical applications.

2. Theoretical background

This section starts by recalling the notation used in this work, which is based on that of Truesdell and Noll [67], Marsden and Hughes [68] and Epstein [69], as well as on that of previous studies [70–73]. Then, we introduce the concepts of connection and torsion, a set of fourth-order tensors that are useful in constitutive theory, and Maxwell’s equations for dielectric materials in the static case. Finally, we show the expression of the Helmholtz free energy for elastic dielectric materials.

2.1. Kinematics

The physical space is represented by the affine space $\mathcal{S}\equiv {\mathbb{E}}^3$ , constructed on ${\mathbb{R}}^3$ seen as both the point space and the modelling vector space [69]. The tangent space $T_x\mathcal{S}$ at point $x\in \mathcal{S}$ is the space of the vectors $\mathbf{v}$ tangent to all possible curves passing through $x$ . For a visual representation of these tangent vectors, see, e.g., [70]. The dual space of the tangent space $T_x\mathcal{S}$ is the cotangent space $T_x^{★}\mathcal{S}$ , the set of all covectors $\mathit{\phi }$ , taking vectors $\mathbf{v}\in T_x\mathcal{S}$ into scalars $〈\mathit{\phi }|\mathit{v}〉={\phi }_a{v}^a$ . We assume that the physical space $\mathcal{S}$ is endowed with a metric tensor field $\mathit{g}$ , called spatial metric tensor and defining the scalar product $\mathit{v}.\mathit{w}\equiv \mathit{s}(\mathit{v},\mathit{w})=v^a{g}_{\mathrm{ab}}{w}^b$ .

A continuum body is identified as one placement $\mathcal{B}$ in $\mathcal{S}$ , called reference configuration. Tangent and cotangent spaces at every point $X$ within $\mathcal{B}$ are defined similarly as in the physical space $\mathcal{S}$ , and so is the material metric tensor field $\mathit{G}$ . Following standard continuum mechanics notation [68], uppercase symbols indicate material quantities, while lowercase symbols denote spatial quantities. A motion of the body $\mathcal{B}$ during the time interval $\mathcal{I}$ is a mapping

$$\varphi :\mathcal{B}\times \mathcal{I}\to \mathcal{S}:(X,t)\mapsto x=\varphi (X,t),$$

(1)

that, at every instant of time $t$ , takes material points $X\in \mathcal{B}$ into spatial points $x\in \mathcal{S}$ . When we fix the time $t\in \mathcal{I}$ , we obtain the configuration (or deformation or placement) $\varphi (·,t)$ , which must be an embedding, meaning that its codomain restriction to the image $\varphi (\mathcal{B},t)$ is a diffeomorphism: a continuous and differentiable map with a continuous and differentiable inverse $\Phi (·,t)\equiv [\varphi (·,t){]}^{-1}$ .

The Lagrangian velocity is the spatial vector field $\stackrel{·}{\varphi }(·,t):\mathcal{B}\to T\mathcal{S}$ over the motion $\varphi$ , defined by

$$\stackrel{·}{\varphi }(X,t)\equiv \left({\partial }_t\varphi \right)(X,t).$$

(2)

The Eulerian velocity is the spatial vector field $\mathit{v}(·,t):\varphi (\mathcal{B},t)\to T\mathcal{S}$ defined by

$$\mathit{v}(x,t)=\mathit{v}(\varphi (X,t),t)=\stackrel{·}{\varphi }(X,t).$$

(3)

At every point $X\in \mathcal{B}$ , the tangent map of $\varphi$ is the differential (note that $x=\varphi (X,t)$ )

$$\mathit{F}(X,t)\equiv \left(T\varphi \right)(X,t):T_X\mathcal{B}\to T_x\mathcal{S}:\mathit{W}\mapsto \mathit{w}=[\mathit{F}(X,t)]\mathit{W},$$

(4)

called deformation gradient, such that, for every tangent vector $\mathit{W}$ in $T_X\mathcal{B}$ , the spatial vector

$$[\mathit{F}(X,t)]\mathit{W}=\left({\partial }_{\mathbf{W}}\varphi \right)(X,t)=\mathit{w},$$

(5)

is the directional derivative of the motion $\varphi$ at $X$ in the direction $\mathit{W}$ . In coordinate charts $X^A$ in $\mathcal{B}$ and $x^a$ in $\mathcal{S}$ , the representation of $\mathit{F}$ takes the form FaA = ϕa,A. The determinant

$$J_{\mathit{F}}(X,t)=\det \mathit{F}(X,t)=\det \sqrt{[g\left(x\right)]}\det [\mathrm{F}(X,t)]\frac{1}{\det \sqrt{[G\left(X\right)]}},$$

(6)

called volume ratio measures volume changes [15,68,74]. In equation (6), $[g]$ , $[G]$ and $[\mathrm{F}]$ denote the representing matrices of $\mathit{g}$ , $\mathit{G}$ and $\mathit{F}$ , respectively. The right Cauchy–Green deformation tensor is the pull-back of the spatial metric $\mathit{g}$ , i.e.,

$$\mathit{C}(X,t)={\mathit{F}}^T(x,t)\mathit{g}\left(x\right)\mathit{F}(X,t),C_{\mathrm{AB}}(X,t)={{\mathrm{F}}^a}_A(X,t)g_{\mathrm{ab}}\left(x\right){{\mathrm{F}}^b}_B(X,t).$$

(7)

Occasionally, with abuse of notation, we shall omit indicating the arguments $(x,t)$ and $(X,t)$ in expressions such as equations (6) and (7).

Finally, we recall the volumetric–distortional decomposition of the deformation gradient $\mathbf{F}$ and the right Cauchy–Green deformation $\mathit{C}$ [75–79]:

$$\mathit{F}=J_{\mathit{F}}^{1/3}\overline{\mathit{F}},\mathit{C}=J_{\mathit{F}}^{2/3}\overline{\mathit{C}}.$$

(8)

2.2. Connection and torsion

Here, we briefly introduce the concept of connection, or covariant derivative, and torsion of a connection, in the same notation as in a previous work [80]. A connection can be defined on any differentiable manifold. For our purposes, we take the body $\mathcal{B}$ as the manifold.

A linear connection on $\mathcal{B}$ defines the vector field ${\nabla }_{\mathit{U}}\mathit{V}$ , called the covariant derivative of the vector field $\mathit{V}$ in the direction of the vector field $\mathit{U}$ . For any smooth vector fields $\mathit{U}$ , $\mathit{V}$ , $\mathit{W}$ and differentiable scalar fields $f$ and $g$ , the connection ∇ must satisfy

$${\nabla }_{\mathit{U}}(\mathit{V}+\mathit{W})={\nabla }_{\mathit{U}}\mathit{V}+{\nabla }_{\mathit{U}}\mathit{W},\mathrm{additivity}\mathrm{in}\mathrm{the}\mathrm{argument},$$

(9a)

$${\nabla }_{\mathit{U}}\left(f\mathit{V}\right)=f{\nabla }_{\mathit{U}}\mathit{V}+\left({\partial }_{\mathit{U}}f\right)\mathit{V},\mathrm{Leibniz}\prime \mathrm{rule}\mathrm{in}\mathrm{the}\mathrm{argument},$$

(9b)

$${\nabla }_{f\mathit{U}+g\mathit{V}}\mathit{W}=f{\nabla }_{\mathit{U}}\mathit{W}+g{\nabla }_{\mathit{V}}\mathit{W},\mathrm{linearity}\mathrm{in}\mathrm{the}\mathrm{direction}\mathrm{of}\mathrm{differentiation},$$

(9c)

where $\left({\partial }_{\mathit{U}}f\right)\left(X\right)$ is the directional derivative of $f$ with respect to the vector field $\mathit{U}$ and, for every $X\in \mathcal{B}$ , equals the differential $\left(\mathrm{D}f\right)\left(X\right)$ applied as a linear map onto the vector $\mathit{U}\left(X\right)$ , because of the hypothesis of differentiability of $f$ , i.e.,

$$\left({\partial }_{\mathit{U}}f\right)\left(X\right)=\underset{h\to 0}{\lim }\frac{f(X+h\mathit{U}\left(X\right))-f\left(X\right)}{h}=[\left(\mathrm{D}f\right)\left(X\right)]\mathit{U}\left(X\right).$$

(10)

Using the properties (9), in a system of coordinates $X^A$ with associated basis vectors ${\mathit{E}}_A$ , the components of the covariant derivative ${\nabla }_{\mathit{U}}\mathit{V}$ are obtained as

$$\begin{array}{c}\multicolumn{1}{c}{{\nabla }_{\mathit{U}}\mathit{V}=U^B{\nabla }_{{\mathit{E}}_B}\left(V^A{\mathit{E}}_A\right)=U^B(V_{,B}^A{\mathit{E}}_A+V^A{\mathrm{\Gamma }}_{\mathrm{AB}}^C{\mathit{E}}_C),}\\ \multicolumn{1}{c}{=U^B(V_{,B}^A+{\mathrm{\Gamma }}_{\mathrm{CB}}^A{V}^C){\mathit{E}}_A,}\\ \multicolumn{1}{c}{=U^B{V}_{;B}^A{\mathit{E}}_A,}\end{array}$$

(11)

where ${\mathrm{\Gamma }}_{\mathrm{AB}}^C$ are the Christoffel symbols of the connection, defined as

$${\nabla }_{{\mathit{E}}_B}{\mathit{E}}_A={\mathrm{\Gamma }}_{\mathrm{AB}}^C{\mathit{E}}_C,$$

(12)

and the derivatives

$$V_{;B}^A=V_{,B}^A+{\mathrm{\Gamma }}_{\mathrm{CB}}^A{V}^C,$$

(13)

are the covariant derivatives of the component $V^A$ .

The torsion of the connection ∇ is defined as the third-order tensor

$$\mathbf{Tor}(\mathit{U},\mathit{V})={\nabla }_{\mathit{U}}\mathit{V}-{\nabla }_{\mathit{V}}\mathit{U}-[\mathit{U},\mathit{V}],$$

(14)

where the Lie bracket $[\mathit{U},\mathit{V}]$ is

$$[\mathit{U},\mathit{V}{]}^B=U^A{V}_{,A}^B-V^A{U}_{,A}^B.$$

(15)

In components, the torsion reads

$${\mathrm{Tor}}_{\mathrm{BC}}^A={\mathrm{\Gamma }}_{\mathrm{BC}}^A-{\mathrm{\Gamma }}_{\mathrm{CB}}^A,$$

(16)

and, when the torsion is zero, the Christoffel symbols are symmetric, i.e., ${\mathrm{\Gamma }}_{\mathrm{BC}}^A={\mathrm{\Gamma }}_{\mathrm{CB}}^A$ .

2.3. Notable fourth-order tensors

This section presents a set of fourth-order tensors employed extensively within the theory of elasticity. These tensors are instrumental in the decomposition of symmetric second-order tensors into their spherical and deviatoric parts. Notably, certain forms of these tensors arise from the differentiation of the right Cauchy–Green deformation tensor $\mathit{C}$ , its inverse ${\mathit{C}}^{-1}$ , and its distortional part $\overline{\mathit{C}}$ with respect to $\mathit{C}$ . An exhaustive treatment of the various families of tensors of this type can be found in a prior publication [79].

The pulled-back mixed symmetric identity $\mathbb{I}$ , spherical operator ${\mathbb{K}}^{*}$ and deviatoric operator ${\mathbb{M}}^{*}$ are defined as

$$\mathbb{I}=\frac{1}{2}(\mathit{I}\underset{¯}{\otimes }\mathit{I}+\mathit{I}\overline{\otimes }\mathit{I}),{\mathrm{I}}_{\mathrm{CD}}^{\mathrm{AB}}=\frac{1}{2}({\delta }_C^A{\delta }_D^B+{\delta }_D^A{\delta }_C^B),$$

(17a)

$${\mathbb{K}}^{*}=\frac{1}{3}{\mathit{C}}^{-1}\otimes \mathit{C},{\left({\mathbb{K}}^{*}\right)}_{\mathrm{CD}}^{\mathrm{AB}}=\frac{1}{3}{\left({\mathit{C}}^{-1}\right)}^{\mathrm{AB}}{C}_{\mathrm{CD}},$$

(17b)

$${\mathbb{M}}^{*}=\mathbb{I}-{\mathbb{K}}^{*},{\left({\mathbb{M}}^{*}\right)}_{\mathrm{CD}}^{\mathrm{AB}}={\mathrm{I}}_{\mathrm{CD}}^{\mathrm{AB}}-{\left({\mathbb{K}}^{*}\right)}_{\mathrm{CD}}^{\mathrm{AB}},$$

(17c)

where $\mathit{I}$ is the second-order material identity tensor and we used the special tensor products ‘tensor-down’ $\underset{¯}{\otimes }$ and ‘tensor-up’ $\overline{\otimes }$ defined by Curnier et al. [81]. The tensor $\mathbb{I}$ acts as the identity on a material symmetric contravariant second-order tensor and as a symmetriser on a material non-symmetric contravariant second-order tensor, hence the terminology ‘symmetric identity’. The tensors ${\mathbb{K}}^{*}$ and ${\mathbb{M}}^{*}$ extract the spherical and deviatoric part of a symmetric contravariant second-order tensor, respectively, with respect to the pulled-back metric $\mathit{C}$ . Their transposes ${\mathbb{I}}^T$ , ${{\mathbb{K}}^{*}}^T$ and ${{\mathbb{M}}^{*}}^T$ do the same on covariant second-order tensors.

By raising all indices of the pulled-back operators (17) with the pull-back inverse metric tensor ${\mathit{C}}^{-1}$ , we obtain their contravariant counterparts

$${\mathbb{I}}^{♯*}=\frac{1}{2}({\mathit{C}}^{-1}\underset{¯}{\otimes }{\mathit{C}}^{-1}+{\mathit{C}}^{-1}\overline{\otimes }{\mathit{C}}^{-1}),{\left({\mathbb{I}}^{♯*}\right)}_{\mathrm{ABCD}}=\frac{1}{2}({\left({\mathit{C}}^{-1}\right)}^{\mathrm{AC}}{\left({\mathit{C}}^{-1}\right)}^{\mathrm{BD}}+{\left({\mathit{C}}^{-1}\right)}^{\mathrm{AD}}{\left({\mathit{C}}^{-1}\right)}^{\mathrm{BC}}),$$

(18a)

$${\mathbb{K}}^{♯*}=\frac{1}{3}({\mathit{C}}^{-1}\otimes {\mathit{C}}^{-1}),{\left({\mathbb{K}}^{♯*}\right)}_{\mathrm{ABCD}}=\frac{1}{3}\left({\left({\mathit{C}}^{-1}\right)}^{\mathrm{AB}}{\left({\mathit{C}}^{-1}\right)}^{\mathrm{CD}}\right),$$

(18b)

$${\mathbb{M}}^{♯*}={\mathbb{I}}^{♯*}-{\mathbb{K}}^{♯*},{\left({\mathbb{M}}^{♯*}\right)}_{\mathrm{ABCD}}={\left({\mathbb{I}}^{♯*}\right)}_{\mathrm{ABCD}}-{\left({\mathbb{K}}^{♯*}\right)}_{\mathrm{ABCD}},$$

(18c)

and by lowering all indices with the pull-back metric tensor $\mathit{C}$ , we obtain their covariant counterparts

$${\mathbb{I}}^{♭*}=\frac{1}{2}(\mathit{C}\underset{¯}{\otimes }\mathit{C}+\mathit{C}\overline{\otimes }\mathit{C}),{\left({\mathbb{I}}^{♭*}\right)}_{\mathrm{ABCD}}=\frac{1}{2}(C_{\mathrm{AC}}{C}_{\mathrm{BD}}+C_{\mathrm{AD}}{C}_{\mathrm{BC}}),$$

(19a)

$${\mathbb{K}}^{♭*}=\frac{1}{3}(\mathit{C}\otimes \mathit{C}),{\left({\mathbb{K}}^{♭*}\right)}_{\mathrm{ABCD}}=\frac{1}{3}\left(C_{\mathrm{AB}}{C}_{\mathrm{CD}}\right),$$

(19b)

$${\mathbb{M}}^{♭*}={\mathbb{I}}^{♭*}-{\mathbb{K}}^{♭*},{\left({\mathbb{M}}^{♭*}\right)}_{\mathrm{ABCD}}={\left({\mathbb{I}}^{♭*}\right)}_{\mathrm{ABCD}}-{\left({\mathbb{K}}^{♭*}\right)}_{\mathrm{ABCD}}.$$

(19c)

Contracting the contravariant tensors (18) and the covariant tensors (19), we obtain the identities

$${\mathbb{I}}^{♯*}:{\mathbb{I}}^{♭*}=\mathbb{I},{\mathbb{I}}^{♭*}:{\mathbb{I}}^{♯*}={\mathbb{I}}^T,$$

(20a)

$${\mathbb{K}}^{♯*}:{\mathbb{K}}^{♭*}={\mathbb{K}}^{*},{\mathbb{K}}^{♭*}:{\mathbb{K}}^{♯*}={{\mathbb{K}}^{*}}^T,$$

(20b)

$${\mathbb{M}}^{♯*}:{\mathbb{M}}^{♭*}={\mathbb{M}}^{*},{\mathbb{M}}^{♭*}:{\mathbb{M}}^{♯*}={{\mathbb{M}}^{*}}^T.$$

(20c)

Finally, we report the important results [77,79]

$$\frac{\partial {\mathit{C}}^{-1}}{\partial \mathit{C}}=-{\mathbb{I}}^{♯*},$$

(21a)

$$\frac{\partial \overline{\mathit{C}}}{\partial \mathit{C}}=J_{\mathit{F}}^{-2/3}{\mathbb{M}}^{*T},$$

(21b)

$$\frac{\partial J_{\mathit{F}}}{\partial \mathit{C}}=\frac{1}{2}{J}_{\mathit{F}}{\mathit{C}}^{-1},$$

(21c)

$$\frac{\partial J_{\mathit{F}}}{\partial \mathit{C}}\otimes \frac{\partial J_{\mathit{F}}}{\partial \mathit{C}}=\frac{1}{4}\left(J_{\mathit{F}}{\mathit{C}}^{-1}\right)\otimes \left(J_{\mathit{F}}{\mathit{C}}^{-1}\right)=\frac{3}{4}{J}_{\mathit{F}}^2{\mathbb{K}}^{♯*},$$

(21d)

$$\frac{{\partial }^2{J}_{\mathit{F}}}{\partial {\mathit{C}}^2}=\frac{1}{4}{\mathit{C}}^{-1}\otimes \left(J_{\mathit{F}}{\mathit{C}}^{-1}\right)+\frac{1}{2}{J}_{\mathit{F}}\frac{\partial {\mathit{C}}^{-1}}{\partial \mathit{C}}=\frac{3}{4}{J}_{\mathit{F}}{\mathbb{K}}^{♯*}-\frac{1}{2}{J}_{\mathit{F}}{\mathbb{I}}^{♯*}.$$

(21e)

2.4. Equations of electrostatics

The formulation of electrostatics draws extensively from the established literature [27–29]. In the absence of volumetric densities of electric charges (no free charges in the interior of a dielectric) and in quasi-static conditions (negligible time derivatives of the fields), Maxwell’s equations reduce to the equations of electrostatics.

Let $\mathcal{R}$ be an arbitrary open subset of the body $\mathcal{B}$ with boundary $\partial \mathcal{R}$ and $\varphi (\mathcal{R},t)$ and $\varphi (\partial \mathcal{R},t)$ the current configuration of the region $\mathcal{R}$ and its boundary at time $t$ . Gauss’ law of conservation of charge states that the flux of the spatial electric induction field (or electric displacement field) $\mathit{d}$ over the arbitrary closed surface $\varphi (\partial \mathcal{R},·)$ must vanish, i.e.,

$$\underset{\varphi (\partial \mathcal{R},·)}{\int }〈\mathit{d}|\mathit{n}〉\mathrm{da}=0,$$

(22)

where $\mathit{n}$ is the spatial normal covector to the surface $\varphi (\partial \mathcal{R})$ .

Let $S$ be an arbitrary surface in the interior of the body $\mathcal{B}$ , the line $\partial S$ its closed boundary, and $\varphi (S,t)$ and $\varphi (\partial S,t)$ their current configuration at time $t$ . The quasi-static Faraday’s law states that the work of the spatial electric field $\mathbf{e}$ (which is a covector field [82,83]) on the arbitrary closed line $\varphi (\partial S)$ must identically vanish, i.e.,

$$\underset{\varphi (\partial S,·)}{\int }〈\mathit{e}|\mathit{t}〉d\ell =0,$$

(23)

where $\mathit{t}$ is the tangent vector to the line $\varphi (\partial S,t)$ .

The local forms of the integral equations (22) and (23) are obtained by applying Gauss’ theorem and Stokes’ theorem, respectively. Together with the constitutive equation for the dielectric induction, these read

$$\mathrm{div}\mathit{d}=0,$$

(24a)

$$\mathrm{curl}\mathit{e}=0,$$

(24b)

$$\mathit{d}={\epsilon }_0\mathbf{e}+\mathbf{\pi },$$

(24c)

where the divergence and the curl operators are taken with respect to the spatial coordinates $x^a$ . In the constitutive equation (24c), the tensor

$${\epsilon }_0={\epsilon }_0{\mathit{g}}^{-1},{\left({\epsilon }_0\right)}^{\mathrm{ab}}={\epsilon }_0{g}^{\mathrm{ab}}$$

(25)

is the constant dielectric permittivity tensor of vacuum, with ${\epsilon }_0$ (lightface symbol) being the (scalar) dielectric permittivity of vacuum, with SI units of farad per metre ( ${\epsilon }_0=8.854·{10}^{-12}\mathrm{F}/\mathrm{m}$ ), ${\epsilon }_0\mathbf{e}$ is the dielectric induction in vacuo, and $\mathbf{\pi }$ is the polarisation, representing the material’s response to the electric field $\mathbf{e}$ .

The material counterpart of equation (22) is obtained by performing the backward Piola transformation

$$\underset{\varphi (\partial \mathcal{R},·)}{\int }〈\mathit{d}|\mathit{n}〉\mathrm{d}a=\underset{\partial \mathcal{R}}{\int }〈J_{\mathit{F}}\mathit{d}{\mathit{F}}^{-T}|\mathit{N}〉\mathrm{d}A=\underset{\partial \mathcal{R}}{\int }〈\mathit{D}|\mathit{N}〉\mathrm{d}A=0,$$

(26)

where $\mathit{N}$ is the referential outward normal covector to the boundary of the referential region $\mathcal{R}\subseteq \mathcal{B}$ and $\mathit{D}$ is the material electric induction. The material counterpart of equation (23) is obtained as the pull-back

$$\underset{\varphi (\partial S,·)}{\int }〈\mathit{e}|\mathit{t}〉\mathrm{d}\ell =\underset{\partial S}{\int }〈{\mathit{F}}^T\mathit{e}|\mathcal{T}〉\mathrm{d}L=\underset{\partial S}{\int }〈\mathit{E}|\mathcal{T}〉\mathrm{d}L,$$

(27)

where $\mathcal{T}$ is the tangent vector to the line $\partial S$ and $\mathit{E}$ is the material electric field. The material counterparts of equation (24) are

$$\mathrm{Div}\mathit{D}=0,$$

(28a)

$$\mathrm{Curl}\mathit{E}=0,$$

(28b)

$$\mathit{D}={\mathcal{E}}_0\mathit{E}+\mathbf{\Pi },$$

(28c)

where the material fields are related to the spatial fields by [22,25,84]

$$\mathit{D}=J_{\mathit{F}}{\mathit{F}}^{-1}\mathit{d},D^A=J_{\mathit{F}}{\left({\mathit{F}}^{-1}\right)}_a^A{d}^a$$

(29a)

$$\mathit{E}={\mathit{F}}^T\mathbf{e},E_A={\mathrm{F}}_A^a{e}_a,$$

(29b)

$$\mathbf{\Pi }=J_{\mathit{F}}{\mathit{F}}^{-1}\mathbf{\pi },{\mathrm{\Pi }}^A=J_{\mathit{F}}{\left({\mathit{F}}^{-1}\right)}_a^A{\pi }^a,$$

(29c)

the divergence and the curl operators are taken with respect to the material coordinates $X^A$ , and the tensor

$${\mathcal{E}}_0=J_{\mathit{F}}{\mathit{F}}^{-1}{\epsilon }_0{\mathit{F}}^{-T}=J_{\mathit{F}}{\epsilon }_0{\mathit{C}}^{-1},{\left({\mathcal{E}}_0\right)}^{\mathrm{AB}}=J_{\mathbf{F}}{\left({\mathit{F}}^{-1}\right)}_a^A\left({\epsilon }_0{g}^{\mathrm{ab}}\right){\left({\mathit{F}}^{-T}\right)}_b^B=J_{\mathit{F}}{\epsilon }_0{\left({\mathit{C}}^{-1}\right)}^{\mathrm{AB}}$$

(30)

is the material dielectric permittivity tensor of vacuum. Note that, although the spatial dielectric permittivity tensor of vacuum ${\epsilon }_0={\epsilon }_0{\mathit{g}}^{-1}$ is constant, in the sense that it is given by the universal constant ${\epsilon }_0$ times the inverse spatial metric ${\mathit{g}}^{-1}$ , its material counterpart (30) depends on deformation.

The spatial equation (24b) and the material equation (28b) imply that the electric field is obtained from the gradient of a scalar field, called electric potential or often just scalar potential, i.e.,

$$\mathit{e}=-\mathrm{grad}\phi ,e_a=-{\phi }_{,a},$$

(31a)

$$\mathit{E}=-\mathrm{Grad}\mathcal{V},E_A=-{\mathcal{V}}_{,A},$$

(31b)

where $\phi$ is the spatial scalar potential and $\mathcal{V}$ is the material scalar potential, defined by

$$\mathcal{V}(X,t)=\phi (x,t).$$

(32)

The spatial gradient of $\phi$ and the material gradient of $\mathcal{V}$ are related via the chain rule:

$$\left(\mathrm{Grad}\mathcal{V}\right)(X,t)={\mathit{F}}^T(x,t)\left(\mathrm{grad}\phi \right)(x,t).$$

(33)

Finally, we note that, on the boundary of the body, the dielectric induction must satisfy the condition

$$〈\mathit{d}|\mathit{n}〉=-\omega ,d^a{n}_a=-\omega ,\mathrm{on}\varphi (\partial \mathcal{B}),$$

(34a)

$$〈\mathit{D}|\mathit{N}〉=-\mathrm{\Omega },D^A{\mathrm{N}}_A=-\mathrm{\Omega },\mathrm{on}\partial \mathcal{B},$$

(34b)

where $\mathit{n}$ is the normal covector to the boundary $\varphi (\partial \mathcal{B},·)$ , $\mathit{N}$ is the normal covector to the boundary $\partial \mathcal{B}$ , $\omega$ is the spatial surface charge density and $\mathrm{\Omega }$ is the material surface charge density, defined by

$$\mathrm{\Omega }(X,t)=J(X,t)\sqrt{{\mathit{C}}^{-1}(X,t):(\mathit{N}(X,t)\otimes \mathbf{N}(X,t))}\omega (x,t).$$

(35)

The factor of proportionality between $\mathrm{\Omega }$ and $\omega$ in equation (35) comes from the transformation of area elements, as discussed by Bonet and Wood [78, section 5.5.1].

3. Material uniformity

A body $\mathcal{B}$ whose material properties explicitly depend on the material point $X$ is called inhomogeneous. An inhomogeneous body $\mathcal{B}$ is called uniform if its material properties depend on the material point $X$ exclusively through a tensor field $\mathit{K}$ , called uniformity tensor or implant tensor. When the uniformity tensor $\mathit{K}$ happens to be integrable, the body is called homogeneous and there exists a reference configuration in which the material properties are independent of the material point $X$ . We will briefly present the theory of uniformity introduced by Noll [4] and further studied by Epstein and Maugin [5] and Epstein and Elżanowski [8], using the notation from Alhasadi et al. [80]. Subsequently, we will discuss the formulation of the energy density for an elastic dielectric material, drawing from the work of Toupin [18] and incorporating the perspective of the theory of uniformity.

3.1. Definition of materially uniform body

A material body $\mathcal{B}$ is called uniform if all of its points are made of the same material [4,5,8]. In the reference configuration $\mathcal{B}$ , two points $X$ and $Y$ may exhibit different material properties but, if the body is uniform, there exists a linear transformation that can ‘fix’ this discrepancy, making the properties at $Y$ identical to those at $X$ .

This ‘fix’ is a linear map called material isomorphism, denoted ${\mathit{K}}_{Y\to X}:T_Y\mathcal{B}\to T_X\mathcal{B}$ , mapping the tangent space of $Y$ into that at $X$ . We can think of the material isomorphism as the linear map that tells us how to transform the neighbourhood of $Y$ to make it match a neighbourhood of $X$ perfectly. Since this applies to any two points in the body, we can imagine a single ‘ideal’ point representing the entire material, called the archetype, whose tangent space is a vector space denoted $\mathcal{A}$ .

The archetype is like a blueprint for the material. For each point $X$ in the body, we have a map $\mathit{K}(X,t):\mathcal{A}\to T_X\mathcal{B}$ that transforms the archetype into the neighbourhood of that point. This map $\mathit{K}(X,t)$ is called the uniformity tensor or implant tensor. It is crucial to keep in mind that the uniformity tensor $\mathit{K}$ is generally non-integrable. This means it cannot be represented by the tangent map of a motion, unlike the deformation gradient $\mathit{F}$ . When $\mathit{K}$ happens to be integrable, implying that there exists a map $\wp$ such that $\mathit{K}=T\wp$ , then $\wp$ becomes a simple change of reference configuration. In this case, the body is referred to as homogeneous. Indeed, the inverse map ${\wp }^{-1}$ would bring the body into a configuration in which the material properties do not depend on the material point $X$ .

Due to the definition of uniformity tensor $\mathit{K}$ , the isomorphism at two points $Y$ and $X$ , ${\mathit{K}}_{Y\to X}:T_Y\mathcal{B}\to T_X\mathcal{B}$ between their tangent spaces can be obtained through transitivity. Given the values $\mathit{K}(Y,t)$ and $\mathit{K}(X,t)$ of the material isomorphism at these points, the isomorphism connecting $T_Y\mathcal{B}$ and $T_X\mathcal{B}$ is ${\mathit{K}}_{Y\to X}=\mathit{K}(X,t){\mathit{K}}^{-1}(Y,t)$ . Recall that $\mathit{K}(·,t)$ is a tensor field defined on $\mathcal{B}$ , i.e., it is a function of $X$ . The concept is illustrated in Figure 1.

Figure 1.

Representation of a materially uniform body $\mathcal{B}$ with the uniformity tensor (or implant tensor or material isomorphism) $\mathit{K}$ between the archetype (or reference crystal) $\mathcal{A}$ and the body $\mathcal{B}$ , and the isomorphism between two generic points $Y$ and $X$ in $\mathcal{B}$ .

A vector field $\mathit{V}(·,t):\mathcal{B}\to T\mathcal{B}$ is said to be $\mathit{K}$ -parallel [4] if, for every material point $X$ ,

$$\mathit{V}(X,t)=\mathit{K}(X,t)\mathfrak{v},$$

(36)

where $\mathfrak{v}$ is a fixed vector in the archetype $\mathcal{A}$ . The $\mathit{K}$ connection is defined as that connection such that the covariant derivative of a $\mathit{K}$ -parallel vector field vanishes identically. In components, we have

$$V_{;B}^A=V_{,B}^A+{\mathrm{\Gamma }}_{\mathrm{CB}}^A{V}^C=0\Rightarrow {\mathrm{\Gamma }}_{\mathrm{CB}}^A{V}^C=-V_{,B}^A,$$

(37)

where the $\mathit{K}$ -covariant derivative is denoted by a semicolon and the Christoffel symbols are ${\mathrm{\Gamma }}_{\mathrm{CB}}^A$ . Using definition (36) of $\mathit{K}$ -parallel vector field and the fact that $\mathfrak{v}$ is a fixed vector in $\mathcal{A}$ , we obtain

$$-V_{,B}^A=-{\mathrm{K}}_{\alpha ,B}^A{\mathfrak{v}}^{\alpha }=-{\mathrm{K}}_{\alpha ,B}^A{\left({\mathit{K}}^{-1}\right)}_C^{\alpha }{V}^C.$$

(38)

Comparison with equation (37) and the arbitrariness of $\mathit{V}$ yield the explicit expression of the Christoffel symbols of the $\mathit{K}$ -connection as

$${\mathrm{\Gamma }}_{\mathrm{CB}}^A=-{\left({\mathit{K}}^{-1}\right)}_C^{\alpha }{\mathrm{K}}_{\alpha ,B}^A={\mathrm{K}}_{\alpha }^A{\left({\mathit{K}}^{-1}\right)}_{C,B}^{\alpha }.$$

(39)

Note that, in the alternative expression ${\mathrm{\Gamma }}_{\mathrm{CB}}^A={\mathrm{K}}_{\alpha }^A{\left({\mathit{K}}^{-1}\right)}_{C,B}^{\alpha }$ , we used ${\mathrm{K}}_{\alpha }^A{\left({\mathit{K}}^{-1}\right)}_B^{\alpha }={\delta }_B^A$ . The torsion of the $\mathit{K}$ -connection is defined as in the general case described by equation (16) and is denoted by $\mathbf{Tor}$ .

3.2. Uniform elastic dielectric body

Let us consider the definition of material uniformity provided in section 3.1. An elastic dielectric body is uniform if its Helmholtz free energy per unit reference volume, evaluated at a point $X$ , can be derived from an archetypal energy density through a linear transformation of the material ‘legs’, represented by the arguments $\mathit{F}$ and $\mathit{E}$ . This transformation is precisely captured by the uniformity tensor $\mathit{K}$ and ensures that the archetypal energy density depends on the material point $X$ solely through $\mathit{K}$ . This concept was established generally by Noll [4] and thoroughly studied for elastic bodies by Epstein and Maugin [5]. We will follow and extend the approach presented by Epstein and Maugin [5] for the case of an elastic dielectric material and write

$$W=\hat{W}°(\mathit{F},\mathit{E},\mathcal{X})=J_{\mathit{K}}^{-1}\hat{\mathfrak{W}}°(\mathit{F}\mathit{K},{\mathit{K}}^T\mathit{E})=\stackrel{\vee }{W}°(\mathit{F},\mathit{E},\mathit{K}),$$

(40)

where $\mathcal{X}$ is the material identity map, i.e., $\mathcal{X}\left(X\right)=X$ , used to indicate the explicit dependence on the material point $X$ . The constitutive function $\hat{\mathfrak{W}}$ is the archetypal energy density and the factor $J_{\mathit{K}}^{-1}$ , in which $J_{\mathit{K}}=\det \mathit{K}$ , serves to account for the fact that $\hat{\mathfrak{W}}$ is per unit archetypal volume, while $\hat{W}$ is per unit reference volume. The constitutive function $\stackrel{\vee }{W}$ , instead, is used when we want to emphasise that the energy depends explicitly on $\mathit{F}$ , $\mathit{E}$ and $\mathit{K}$ .

By defining the elastic (or archetypal) deformation, mapping from the archetype to the current configuration, and the archetypal electric field, respectively, as

$$\mathfrak{F}=\mathit{F}\mathit{K},{\mathfrak{F}}_{\beta }^a={\mathrm{F}}_B^a{\mathrm{K}}_{\beta }^B,$$

(41a)

$$\mathfrak{E}={\mathit{K}}^T\mathit{E},{\mathfrak{E}}_{\alpha }={\left({\mathit{K}}^T\right)}_{\alpha }^B{E}_B,$$

(41b)

we can write the energy density (40) as

$$W=\hat{W}°(\mathit{F},\mathit{E},\mathcal{X})=J_{\mathbf{K}}^{-1}\hat{\mathfrak{W}}°(\mathfrak{F},\mathfrak{E})=\stackrel{\vee }{W}°(\mathit{F},\mathit{E},\mathit{K}).$$

(42)

Figure 2 depicts two key mappings:

• The uniformity tensor $\mathit{K}$ maps the archetype $\mathcal{A}$ to the reference configuration $\mathcal{B}$ ; the archetypal electric field $\mathfrak{E}$ is pushed-forward into the material electric field $\mathit{E}={\mathit{K}}^{-T}\mathfrak{E}$ (which is the inverse of equation (41b)).

• The elastic deformation $\mathfrak{F}$ maps the archetype $\mathcal{A}$ to the current configuration $\varphi (\mathcal{B},t)$ ; the archetypal electric field $\mathfrak{E}$ is pushed-forward into the spatial electric field $\mathit{e}={\mathfrak{F}}^{-T}\mathfrak{E}$ (obtained combining equations (29b) and (41b)).

Figure 2.

Mappings between the archetype $\mathcal{A}$ and the tangent spaces in the body and in the current configuration for a uniform body $\mathcal{B}$ . The red arrows represent the electric field in the archetype, the material electric field at the material point $X$ in the body $\mathcal{B}$ and the spatial electric field at the spatial point $x$ in the current configuration $\varphi (\mathcal{B},t)$ . Figure modified from Alhasadi et al. [72].

When we differentiate the free energy (42) with respect to $\mathit{F}$ and $\mathit{E}$ we obtain the identities

$$\frac{\partial \hat{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathcal{X})=J_{\mathit{K}}^{-1}\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{F}}°(\mathfrak{F},\mathfrak{E}){\mathit{K}}^T=\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K}),$$

(43a)

$$\frac{\partial \hat{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathcal{X})=J_{\mathit{K}}^{-1}\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{E}}°(\mathfrak{F},\mathfrak{E}){\mathit{K}}^T=\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K}).$$

(43b)

4. Governing equations

This section elucidates the governing equations for uniform non-linear electro-elastic solids. The formulation presented here treats the uniformity tensor $\mathit{K}$ as a non-standard kinematical variable, acted upon by non-standard (or configurational) forces. This picture is based on the works by Gurtin [85,86], Cermelli et al. [87], DiCarlo and Quiligotti [88], Grillo et al. [89,90] and Grillo and Di Stefano [91–94]. The dissipation inequality is enforced to determine the admissible form of the constitutive equations. Then, the finite element formulation is introduced and, finally, the specific constitutive equations used in the benchmark problem of section 6 are presented.

4.1. Principle of virtual work

Considering the motion $\varphi$ and the fields $\mathit{F}$ , $\mathcal{V}$ , $\mathit{E}$ and $\mathit{K}$ , their virtual variations are

$$\tilde{\varphi }=\varphi +\eta ,$$

(44a)

$$\tilde{\mathit{F}}=\mathit{F}+\mathrm{Grad}\eta ,$$

(44b)

$$\tilde{\mathcal{V}}=\mathcal{V}+\mathcal{Z},$$

(44c)

$$\tilde{\mathit{E}}=\mathit{E}-\mathrm{Grad}\mathcal{Z},$$

(44d)

$$\tilde{\mathit{K}}=\mathit{K}+\mathbf{\Xi },$$

(44e)

where $\eta$ denotes the Lagrangian virtual displacement field (the virtual variation of the motion $\varphi$ ), $\mathcal{Z}$ represents the virtual variation in the electric potential $\mathcal{V}$ , and $\mathbf{\Xi }$ is the virtual variation in the uniformity tensor $\mathit{K}$ . Recall that a virtual displacement is any displacement that is in harmony with the constraints of the system [95].

The principle of virtual work states that the internal virtual work ${\mathcal{W}}_{\mathrm{int}}$ equals the external virtual work ${\mathcal{W}}_{\mathrm{ext}}$ . Using the variations (44) and considering static conditions, the principle of virtual work reads

$$\underset{\mathrm{internal}\mathrm{virtual}\mathrm{work}{\mathcal{W}}_{\mathrm{int}}}{\underbrace{\underset{\mathcal{B}}{\int }(\mathit{T}:\mathrm{Grad}\eta +〈\mathit{D}|\mathrm{Grad}\mathcal{Z}〉+\left({\mathit{Y}}_i{\mathit{K}}^{-T}\right):\mathbf{\Xi })}}=\underset{\mathrm{external}\mathrm{virtual}\mathrm{work}{\mathcal{W}}_{\mathrm{ext}}}{\underbrace{\underset{\partial \mathcal{B}}{\int }(〈\mathfrak{t}|\eta 〉-\mathrm{\Omega }\mathcal{Z})+\underset{\mathcal{B}}{\int }(\mathit{f}|\eta +\left({\mathit{Y}}_e{\mathit{K}}^{-T}\right):\mathbf{\Xi })}},$$

(45)

where the completely material tensors ${\mathit{Y}}_i$ and ${\mathit{Y}}_e$ are the internal and external configurational body forces [85–87,90–94], respectively, and their linear compositions with ${\mathit{K}}^{-T}$ , i.e., the tensors ${\mathit{Y}}_i{\mathit{K}}^{-T}$ and ${\mathit{Y}}_e{\mathit{K}}^{-T}$ , exert work on the generalised virtual displacement $\mathbf{\Xi }$ . Furthermore, $\mathrm{\Omega }$ denotes the material surface charge density (see section 2.4), $\mathfrak{t}$ is the external surface traction force per unit referential area and $\mathit{f}$ is the body force per unit referential volume.

Applying the divergence theorem to the surface integrals in equation (45), substituting the identities

$$\mathit{T}:\mathrm{Grad}\eta =\mathrm{Div}\left(\mathit{T}\eta \right)-〈\mathrm{Div}\mathit{T})|\eta 〉,$$

(46a)

$$〈\mathit{D}|\mathrm{Grad}\mathcal{Z}〉=\mathrm{Div}\left(\mathcal{Z}\mathit{D}\right)-\left(\mathrm{Div}\mathit{D}\right)\mathcal{Z},$$

(46b)

in the first and second terms in the external virtual power ${\mathcal{W}}_{\mathrm{ext}}$ of equation (45), factorising the generalised virtual displacements $\eta$ , $\mathcal{Z}$ and $\mathrm{\Xi }$ and localising, we arrive at the governing equations and boundary conditions

$$\mathrm{Div}\mathit{T}+\mathit{f}=0,\mathrm{in}\mathcal{B},$$

(47a)

$$\mathrm{Div}\mathit{D}=0,\mathrm{in}\mathcal{B},$$

(47b)

$${\mathit{Y}}_i={\mathit{Y}}_e,\mathrm{in}\mathcal{B},$$

(47c)

$$\mathit{T}\mathit{N}=\mathfrak{t},\mathrm{on}\partial \mathcal{B},$$

(47d)

$$〈\mathit{D}|\mathit{N}〉=-\mathrm{\Omega },\mathrm{on}\partial \mathcal{B}.$$

(47e)

4.2. Dissipation inequality

In the isothermal case, the dissipation inequality reads [87,89,90]

$$-\stackrel{·}{W}+\mathit{T}:\stackrel{·}{\mathit{F}}-〈\mathit{D}|\stackrel{·}{\mathit{E}}〉+\left({\mathit{Y}}_i{\mathit{K}}^{-T}\right):\stackrel{·}{\mathit{K}}\ge 0.$$

(48)

Differentiating the referential energy density (42) with respect to time, we obtain

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{W}=[\stackrel{\vee }{W}°(\mathit{F},\mathit{E},\mathit{K})]·}\\ \multicolumn{1}{c}{=[\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K})]:\stackrel{·}{\mathit{F}}+〈\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K})|\stackrel{·}{\mathit{E}}〉+[\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K})]:\stackrel{·}{\mathit{K}},}\end{array}$$

(49a)

which, substituted in the dissipation inequality (48), gives

$$\begin{array}{c}\multicolumn{1}{c}{(\mathit{T}-\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K})):\stackrel{·}{\mathit{F}}+〈-\mathit{D}-\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K})|\stackrel{·}{\mathit{E}}〉+}\\ \multicolumn{1}{c}{+({\mathit{Y}}_i{\mathit{K}}^{-T}-\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K})):\stackrel{·}{\mathit{K}}\ge 0.}\end{array}$$

(50)

For the inequality (50) to be valid for any electromechanical process, it is necessary that, considering also equation (43),

$$\mathit{T}=\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K})=J_{\mathit{K}}^{-1}\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{F}}°(\mathfrak{F},\mathfrak{E}){\mathit{K}}^T,$$

(51a)

$$\mathit{D}=-\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K})=-J_{\mathit{K}}^{-1}\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{E}}°(\mathfrak{F},\mathfrak{E}){\mathit{K}}^T.$$

(51b)

The dissipation inequality (50) reduces to

$$({\mathit{Y}}_i{\mathit{K}}^{-T}-\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K})):\stackrel{·}{\mathit{K}}\ge 0,$$

(52)

and, using the definition

$${\mathit{L}}_{\mathit{K}}=\stackrel{·}{\mathit{K}}{\mathit{K}}^{-1}$$

(53)

of inhomogeneity rate tensor, it can be written as

$$({\mathit{Y}}_i-(\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K})){\mathit{K}}^T):{\mathit{L}}_{\mathit{K}}\ge 0.$$

(54)

Let us consider the Helmholtz free density (42) and calculate the derivative with respect to the uniformity tensor $\mathit{K}$ [5], considering the expressions (41a) and (41b), i.e.,

$$\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K})=-J_{\mathit{K}}^{-1}\hat{\mathfrak{W}}°(\mathfrak{F},\mathfrak{E}){\mathit{K}}^{-T}+{\mathbf{F}}^T{J}_v^{-1}\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{F}}°(\mathfrak{F},\mathfrak{E})+\mathit{E}\otimes J_{\mathit{K}}^{-1}\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{E}}°(\mathfrak{F},\mathfrak{E}).$$

(55)

By exploiting equation (42) and the identities

$$\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{F}}°(\mathfrak{F},\mathfrak{E})=J_{\mathit{K}}\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K}){\mathit{K}}^{-T}=J_{\mathit{K}}\mathit{T}{\mathit{K}}^{-T},$$

(56a)

$$\frac{\partial \hat{\mathfrak{W}}}{\partial \mathfrak{E}}°(\mathfrak{F},\mathfrak{E})=J_{\mathit{K}}\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K}){\mathit{K}}^{-T}=-J_{\mathit{K}}\mathit{D}{\mathit{K}}^{-T},$$

(56b)

arising from equation (51), we can put equation (55) in the form

$$\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K})=-[W{\mathit{I}}^T-{\mathit{F}}^T\mathit{T}+\mathit{E}\otimes \mathit{D}]{\mathit{K}}^{-T}=-\mathit{H}{\mathit{K}}^{-T},$$

(57)

${\mathit{I}}^T$ is the transpose of the material identity tensor, with components ${\delta }_A^B$ , and the expression within the square parentheses represents the negative of the electro-elastic Eshelby stress tensor $\mathit{H}$ [5,85,96,97], i.e.,

$$\mathit{H}=W{\mathit{I}}^T-{\mathit{F}}^T\mathit{T}+\mathit{E}\otimes \mathit{D},{H_A}^B=W{{\delta }_A}^B-{{\left({\mathit{F}}^T\right)}_A}^a{T_a}^B+E_A{D}^B.$$

(58)

By substituting equation (57) into the dissipation inequality (54), we arrive at

$$({\mathit{Y}}_i+\mathit{H}):{\mathit{L}}_{\mathit{K}}\ge 0.$$

(59)

The negative of the Eshelby stress $\mathit{H}$ is the non-dissipative part of the internal generalised force ${\mathit{Y}}_i$ . Indeed, in the particular case in which the internal generalised force ${\mathit{Y}}_i$ coincides with $-\mathit{H}$ , i.e.,

$${\mathit{Y}}_i=-\mathit{H}=\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K}){\mathit{K}}^T,$$

(60)

the generalised force ${\mathit{Y}}_i$ does not dissipate and equation (59) is satisfied with the equal sign. In this case, the evolution of the system is called self-driven [8]. For the case of self-driven evolution, equation (60) implies that the balance equation (47c) reduces to

$$\left({\mathit{Y}}_i{\mathit{K}}^{-T}\right)=\frac{\partial \hat{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathcal{X})=\left({\mathit{Y}}_e{\mathit{K}}^{-T}\right),$$

(61)

where we used also equation (57). Therefore, if we write the total virtual work (45) by expressing the Piola stress $\mathit{T}$ and the electric induction $\mathit{D}$ as in equation (51), and the energy $W$ as in equation (42), i.e.,

$$\begin{array}{c}\multicolumn{1}{c}{\underset{\mathrm{internal}\mathrm{virtual}\mathrm{work}{\mathcal{W}}_{\mathrm{int}}}{\underbrace{\underset{\mathcal{B}}{\int }(\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K}):\mathrm{Grad}\eta -〈\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K})|\mathrm{Grad}\mathcal{Z}〉+\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{K}}°(\mathit{F},\mathit{E},\mathit{K}):\mathbf{\Xi })}}=}\\ \multicolumn{1}{c}{=\underset{\mathrm{external}\mathrm{virtual}\mathrm{work}{\mathcal{W}}_{\mathrm{ext}}}{\underbrace{\underset{\partial \mathcal{B}}{\int }(〈\mathfrak{t}|\eta 〉-\mathrm{\Omega }\mathcal{Z})+\underset{\mathcal{B}}{\int }(〈\mathit{f}|\eta 〉+\left({\mathit{Y}}_e{\mathit{K}}^{-T}\right):\mathbf{\Xi })}},}\end{array}$$

(62)

and consider the condition (61) of self-driven evolution, the principle of virtual work reduces to

$$\underset{\mathcal{B}}{\int }(\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K}):\mathrm{Grad}\eta -〈\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K})|\mathrm{Grad}\mathcal{Z}〉)\equiv \underset{\partial \mathcal{B}}{\int }(〈\mathfrak{t}|\eta 〉-\mathrm{\Omega }\mathcal{Z})-\underset{\mathcal{B}}{\int }〈\mathit{f}|\eta 〉,$$

(63)

i.e., it coincides with the case in which the uniformity tensor $\mathit{K}$ is considered as an internal variable.

Following the work by Elżanowski and Preston [98], we express the self-driven evolution law as in terms of the torsion Tor of the connection induced by the uniformity tensor $\mathit{K}$ . They used the expression

$${\mathfrak{L}}_{\mathit{K}}=\mathbf{C}⋮({\mathit{K}}^{-1}\mathbf{Tor}:\mathit{K}\otimes \mathit{K}),{{\left({\mathfrak{L}}_{\mathit{K}}\right)}^{\rho }}_{\sigma }={{C^{\rho }}_{\alpha }}^{\beta \gamma }{{\left({\mathit{K}}^{-1}\right)}^{\alpha }}_I{\mathrm{Tor}}_{\mathrm{JL}}^I{\mathrm{K}}_{\beta }^J{\mathrm{K}}_{\gamma }^L,$$

(64)

where $\mathbf{C}$ is a completely archetypal fifth-order tensor with components ${C_{\sigma \alpha }^{\rho }}^{\beta \gamma }$ , the symbol ⋮ denotes triple contraction and

$${\mathfrak{L}}_{\mathit{K}}={\mathit{K}}^{-1}{\mathit{L}}_{\mathit{K}}\mathit{K}={\mathit{K}}^{-1}\stackrel{·}{\mathit{K}}$$

(65)

is the pull-back of the inhomogeneity rate tensor ${\mathit{L}}_{\mathit{K}}$ of equation (53) to the archetype [8]. By substituting equation (66) into equation (65), we obtain the component expression of the evolution law for $\mathit{K}$ as

$${\stackrel{·}{\mathrm{K}}}_{\sigma }^M={\mathrm{K}}_{\rho }^M{C_{\sigma \alpha }^{\rho }}^{\beta \gamma }{\left({\mathit{K}}^{-1}\right)}_I^{\alpha }{\mathrm{Tor}}_{\mathrm{JL}}^I{\mathrm{K}}_{\beta }^J{\mathrm{K}}_{\gamma }^L.$$

(66)

4.3. Self-driven evolution law

Building on the work of Elżanowski and Preston [98], we introduce self-driven evolution laws based on the torsion Tor of the connection induced by the uniformity tensor $\mathit{K}$ , as described in equation (64). Our analysis will focus on two specific forms of the uniformity tensor $\mathit{K}$ : the diagonal case and the nilpotent case [98]. In both cases, the uniformity tensor $\mathit{K}$ remains independent of the third material coordinate $X^3$ in the reference configuration, which thus serves as a ‘pivot axis’ for the uniformity tensor. Each case provides unique insights into the self-driven evolution of inhomogeneities.

4.3.1. Diagonal case

In the diagonal case, the representing matrix of the uniformity tensor is, unsurprisingly, diagonal:

$$[\mathit{K}(X^1,X^2,t)]=\left[\begin{array}{ccc}A(X^1,X^2,t)\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1/A(X^1,X^2,t)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right].$$

(67)

Substituting the uniformity tensor from equation (67) into the evolution law, which is formulated in terms of the torsion Tor, the matrix representation of ${\mathfrak{L}}_{\mathit{K}}$ is, as per equation (64),

$$\begin{array}{c}\multicolumn{1}{c}{[{\mathfrak{L}}_{\mathit{K}}]=\left[\begin{array}{ccc}\frac{{\partial }_{t}A}{A}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -\frac{{\partial }_{t}A}{A}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right],}\\ \multicolumn{1}{c}{=\left[\begin{array}{ccc}({C_{12}^1}^{12}-{C_{12}^1}^{21})A_{,1}-\frac{({C_{11}^1}^{21}-{C_{11}^1}^{12})A_{,2}}{A^2}\hfill & ({C_{22}^1}^{12}-{C_{22}^1}^{21})A_{,1}-\frac{({C_{21}^1}^{21}-{C_{21}^1}^{12})A_{,2}}{A^2}\hfill & 0\hfill \\ \multicolumn{1}{c}{({C_{12}^2}^{12}-{C_{12}^2}^{21})A_{,1}-\frac{({C_{11}^2}^{21}-{C_{11}^2}^{12})A_{,2}}{A^2}}& ({C_{22}^2}^{12}-{C_{22}^2}^{21})A_{,1}-\frac{({C_{21}^2}^{21}-{C_{21}^2}^{12})A_{,2}}{A^2}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right].}\end{array}$$

(68)

The archetypal velocity gradient ${\mathfrak{L}}_{\mathit{K}}$ shows a symmetric deformation where one axis is stretching (or shortening) at a rate proportional to ${\partial }_{t}A/A$ , while the perpendicular axis undergoes an equal and opposite rate of shortening (or stretching), satisfying the condition ${\left({\mathfrak{L}}_{\mathit{K}}\right)}_1^1=-{\left({\mathfrak{L}}_{\mathit{K}}\right)}_2^2$ . Since the eigenvalues are equal and opposite, ${\mathfrak{L}}_{\mathit{K}}$ describes pure shear. Equating the corresponding matrix entries in equation (68) and simplifying, we arrive at the conditions

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_1^1={\partial }_{t}A=({C_{12}^1}^{12}-{C_{12}^1}^{21})A{A}_{,1}-({C_{11}^1}^{21}-{C_{11}^1}^{12})A^{-1}{A}_{,2},$$

(69a)

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_2^1=0=({C_{22}^1}^{12}-{C_{22}^1}^{21})A{A}_{,1}-({C_{21}^1}^{21}-{C_{21}^1}^{12})A^{-1}{A}_{,2},$$

(69b)

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_1^2=0=({C_{12}^2}^{12}-{C_{12}^2}^{21})A{A}_{,1}-({C_{11}^2}^{21}-{C_{11}^2}^{12})A^{-1}{A}_{,2},$$

(69c)

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_2^2=-{\partial }_{t}A=({C_{22}^2}^{12}-{C_{22}^2}^{21})A{A}_{,1}-({C_{21}^2}^{21}-{C_{21}^2}^{12})A^{-1}{A}_{,2}.$$

(69d)

For equations (69b) and (69c) to be satisfied identically, it is necessary that

$${C_{22}^1}^{12}={C_{22}^1}^{21}\mathrm{and}{C_{21}^1}^{21}={C_{21}^1}^{12},$$

(70a)

$${C_{12}^2}^{12}={C_{12}^2}^{21}\mathrm{and}{C_{11}^2}^{21}={C_{11}^2}^{12}.$$

(70b)

Applying the condition ${\left({\mathfrak{L}}_{\mathit{K}}\right)}_1^1=-{\left({\mathfrak{L}}_{\mathit{K}}\right)}_2^2$ to equations (69a) and (69d), we obtain

$$\begin{array}{c}\multicolumn{1}{c}{({C_{12}^1}^{12}-{C_{12}^1}^{21})A{A}_{,1}-({C_{11}^1}^{21}-{C_{11}^1}^{12})A^{-1}{A}_{,2}=}\\ \multicolumn{1}{c}{[-({C_{22}^2}^{12}-{C_{22}^2}^{21})A{A}_{,1}+({C_{21}^2}^{21}-{C_{21}^2}^{12})A^{-1}{A}_{,2}].}\end{array}$$

(71)

For this condition to be satisfied identically, it is necessary that

$$({C_{12}^1}^{12}-{C_{12}^1}^{21})=({C_{22}^2}^{21}-{C_{22}^2}^{12})\mathrm{and}({C_{11}^1}^{21}-{C_{11}^1}^{12})=({C_{21}^2}^{12}-{C_{21}^2}^{21}).$$

(72)

Therefore, by applying the relationships between the material constants, we derive a single partial differential equation for the scalar function $A$ (which depends on the material coordinates $X^1$ and $X^2$ and on time $t$ ), which reads

$${\partial }_{t}A=C_{1}A{A}_{,1}-C_2{A}^{-1}{A}_{,2},$$

(73)

where $C_1=({C_{22}^2}^{21}-{C_{22}^2}^{12})$ and $C_2=({C_{21}^2}^{12}-{C_{21}^2}^{21})$ . Since we will assume a static electric load later in the numerical example, it is appropriate to consider steady-state conditions. Therefore, we set ${\partial }_{t}A=0$ in equation (73), which simplifies to

$$0=C_{1}A{A}_{,1}-C_2{A}^{-1}{A}_{,2}.$$

(74)

We solve the partial differential equation (74) using separation of variables and obtain the time-independent expression of $A$

$$A(X^1,X^2)=\pm \frac{\sqrt{2C_3{X}^1+2C_4}}{\sqrt{(C_5-2C_6{X}^2)}}.$$

(75)

To ensure that $A(X^1,X^2)$ is real, we must have $C_5-2C_6{X}^2>0$ , where $C_3=\frac{1}{C_1}$ , $C_4$ , $C_5$ and $C_6=\frac{1}{C_2}$ are material constants. This imposes constraints on $C_5$ and the term $2C_6{X}^2$ .

4.3.2. Nilpotent case

In the nilpotent case, the representing matrix of the uniformity tensor is upper-triangular:

$$[\mathit{K}(X^1,X^2,t)]=\left[\begin{array}{ccc}1\hfill & B(X^1,X^2,t)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right].$$

(76)

Substituting the uniformity tensor (76) into the evolution law (64) yields the matrix equation

$$[{\mathfrak{L}}_{\mathit{K}}]=\left[\begin{array}{ccc}0\hfill & {\partial }_{t}B\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right]=\left[\begin{array}{ccc}({C_{11}^1}^{21}-{C_{11}^1}^{12})B_{,1}\hfill & ({C_{11}^2}^{21}-{C_{11}^2}^{12})B_{,1}\hfill & 0\hfill \\ \multicolumn{1}{c}{({C_{21}^1}^{21}-{C_{21}^1}^{12})B_{,1}}& ({C_{21}^2}^{21}-{C_{21}^2}^{12})B_{,1}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right].$$

(77)

The archetypal velocity gradient ${\mathfrak{L}}_{\mathit{K}}$ has an off-diagonal element ${\partial }_{t}B$ , indicating shear along one axis. Equating the corresponding matrix elements in equation (77) and simplifying leads to the conditions

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_1^1=0=({C_{11}^1}^{21}-{C_{11}^1}^{12})B_{,1},$$

(78a)

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_2^1={\partial }_{t}B=({C_{11}^2}^{21}-{C_{11}^2}^{12})B_{,1},$$

(78b)

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_1^2=0=({C_{21}^1}^{21}-{C_{21}^1}^{12})B_{,1},$$

(78c)

$${\left({\mathfrak{L}}_{\mathit{K}}\right)}_2^2=0=({C_{21}^2}^{21}-{C_{21}^2}^{12})B_{,1}.$$

(78d)

Conditions (78a), (78c) and (78d) yield

$${C^1}_{11}^{21}={C^1}_{11}^{12},$$

(79a)

$${C^1}_{21}^{21}={C^1}_{21}^{12},$$

(79b)

$${C^2}_{21}^{21}={C^2}_{21}^{12}.$$

(79c)

The only non-trivial partial differential equation is (78b) and reads

$${\partial }_{t}B=C_7{B}_{,1},$$

(80)

where we set $C_7=({C_{11}^2}^{21}-{C_{11}^2}^{12})$ . At steady state, ${\partial }_{t}B=0$ and the governing equation become

$$B_{,1}=0,$$

(81)

i.e., $B$ is independent of the coordinate $X^1$ . A solution of equation (81) is simply

$$B(X^1,X^2)=B_0\left(X^2\right),$$

(82)

where $B_0$ is an arbitrary real-valued function of $X^2$ .

4.4. Finite element formulation

We consider a body $\mathcal{B}$ and subdivide it into elements, each with the same number $m$ of nodes. Within each element $e$ , we define the motion, $\varphi (X,t)$ , and the electric potential, $\mathcal{V}(X,t)$ , based on the nodal motion ${\overline{\varphi }}_{\mathrm{ei}}$ and the nodal electric potential ${\overline{\mathcal{V}}}_{\mathrm{ei}}$ at each node $i$ . The values of the motion and the electric potential are interpolated from the nodal values by means of the shape functions ${\mathrm{N}}_i\left(X\right)$ . Therefore, in the $\mathrm{eth}$ element ${\mathcal{B}}_e$ , we have the relationships

$${\varphi }^a(X,t)=\sum _{i=1}^m{\mathrm{N}}_i\left(X\right){\overline{\varphi }}_{\mathrm{ei}}^a\left(t\right),$$

(83a)

$$\mathcal{V}(X,t)=\sum _{i=1}^m{\mathrm{N}}_i\left(X\right){\overline{\mathcal{V}}}_{\mathrm{ei}}\left(t\right),$$

(83b)

$${\mathrm{F}}_A^a(X,t)=\sum _{i=1}^m\frac{\partial {\mathrm{N}}_i}{\partial X^A}\left(X\right){\overline{\varphi }}_{\mathrm{ei}}^a\left(t\right),$$

(83c)

$$E_A(X,t)=-\sum _{i=1}^m\frac{\partial {\mathrm{N}}_i}{\partial X^A}\left(X\right){\overline{\mathcal{V}}}_{\mathrm{ei}}\left(t\right),$$

(83d)

and, for the variations in motion ${\varphi }^a$ and electric potential $\mathcal{V}$ , and the corresponding variations in the deformation gradient ${\mathrm{F}}_A^a$ and electric field $E_A$ , we have

$${\eta }^a(X,t)=\sum _{i=1}^m{\mathrm{N}}_i\left(X\right){\overline{\eta }}_{\mathrm{ei}}^a\left(t\right),$$

(84a)

$$\mathcal{Z}(X,t)=\sum _{i=1}^m{\mathrm{N}}_i\left(X\right){\overline{\mathcal{Z}}}_{\mathrm{ei}}\left(t\right),$$

(84b)

$${\eta }_{|A}^a(X,t)=\sum _{i=1}^m\frac{\partial {\mathrm{N}}_i}{\partial X^A}\left(X\right){\overline{\eta }}_{\mathrm{ei}}^a\left(t\right),$$

(84c)

$$-{\mathcal{Z}}_{,A}(X,t)=-\sum _{i=1}^m\frac{\partial {\mathrm{N}}_i}{\partial X^A}\left(X\right){\overline{\mathcal{Z}}}_{\mathrm{ei}}\left(t\right).$$

(84d)

Incorporating these relationships into the principle of virtual work in the form equation (63), we arrive at a set of discrete equations, which can be linearised, by employing the Newton–Raphson iterative method, into

$$\left[\begin{array}{cc}\underset{3m\times 3m}{[{\mathrm{K}}_e^{\varphi \varphi }]}\hfill & \underset{3m\times m}{[{\mathrm{K}}_e^{\varphi \mathcal{V}}]}\hfill \\ \multicolumn{1}{c}{\underset{m\times 3m}{[{\mathrm{K}}_e^{\mathcal{V}\varphi }]}}& \underset{m\times m}{[{\mathrm{K}}_e^{\mathcal{V}\mathcal{V}}]}\hfill \end{array}\right]\left[\begin{array}{c}\underset{3m\times 1}{\{{\overline{\eta }}_e\}}\hfill \\ \multicolumn{1}{c}{\underset{m\times 1}{\{{\overline{\mathcal{Z}}}_e\}}}\end{array}\right]=-\left[\begin{array}{c}\underset{3m\times 1}{\{{\mathrm{F}}_e^{\varphi }\}}\hfill \\ \multicolumn{1}{c}{\underset{m\times 1}{\{{\mathrm{F}}_e^{\mathcal{V}}\}}}\end{array}\right],$$

(85)

where we recall that $m$ is the number of nodes in the element. These are linear equations in the nodal variables $\{{\overline{\eta }}_e\}$ and $\{{\overline{\mathcal{Z}}}_e\}$ . In compact form, we can write

$$\underset{4m\times 4m}{[{\mathrm{K}}_e]}\underset{4m\times 1}{\left\{{ū}_e\right\}}=-\underset{4m\times 1}{\left\{{\mathrm{F}}_e\right\}}.$$

(86)

The vector $\left\{{ū}_e\right\}$ is the local vector of the generalised displacements, given by

$$\underset{4m\times 1}{\left\{{ū}_e\right\}}=\left[\begin{array}{c}\underset{3m\times 1}{\left\{{\overline{\eta }}_e\right\}}\\ \underset{m\times 1}{\left\{{\overline{\mathcal{Z}}}_e\right\}}\end{array}\right].$$

(87)

The local stiffness matrix $[{\mathrm{K}}_e]$ is given by the blocks

$$\underset{3m\times 3m}{[{\mathrm{K}}_e^{\varphi \varphi }]}=\left[\left[\underset{{\mathcal{B}}_e}{\int }{{{A_a}^A}_b}^B\frac{\partial {\mathrm{N}}_i}{\partial X^A}\frac{\partial {\mathrm{N}}_j}{\partial X^B}\right]\right],$$

(88a)

$$\underset{3m\times m}{[{\mathrm{K}}_e^{\varphi \mathcal{V}}]}=-\left[\left[\underset{{\mathcal{B}}_e}{\int }{{\mathrm{D}}_a}^{\mathrm{AB}}\frac{\partial {\mathrm{N}}_i}{\partial X^A}\frac{\partial {\mathrm{N}}_j}{\partial X^B}\right]\right],$$

(88b)

$$\underset{m\times 3m}{[{\mathrm{K}}_e^{\mathcal{V}\varphi }]}=-\left[\left[\underset{{\mathcal{B}}_e}{\int }{{\mathrm{D}}_b}^{\mathrm{BA}}\frac{\partial {\mathrm{N}}_i}{\partial X^A}\frac{\partial {\mathrm{N}}_j}{\partial X^B}\right]\right],$$

(88c)

$$\underset{m\times m}{[{\mathrm{K}}_e^{\mathcal{V}\mathcal{V}}]}=\underset{{\mathcal{B}}_e}{\int }\left[\left[Y^{\mathrm{AB}}\frac{\partial {\mathrm{N}}_i}{\partial X^A}\frac{\partial {\mathrm{N}}_j}{\partial X^B}\right]\right],$$

(88d)

where

$$\mathbb{A}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \mathit{F}\partial \mathit{F}}°(\mathit{F},\mathit{E},\mathit{K}),{{A_a}_b^A}^B=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial {\mathrm{F}}_A^a\partial {\mathrm{F}}_B^b}°(\mathit{F},\mathit{E},\mathit{K}),$$

(89a)

$$\mathbf{D}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \mathit{F}\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K}),{{\mathrm{D}}_a}^{\mathrm{AB}}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial {\mathrm{F}}_A^a\partial E_B}°(\mathit{F},\mathit{E},\mathit{K}),$$

(89b)

$$\mathit{Y}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \mathit{E}\partial \mathit{E}}°(\mathit{F},\mathit{E},\mathit{K}),Y^{\mathrm{AB}}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial E_A\partial E_B}°(\mathit{F},\mathit{E},\mathit{K}),$$

(89c)

are the first elasticity tensor, the electro-elastic stiffness tensor and the dielectric permittivity tensor, respectively. The local vector of the nodal forces is given by the block vectors

$$\{{\mathrm{F}}_e^{\varphi }\}=\left[[\underset{{\mathcal{B}}_e}{\int }{T_a}^A\frac{\partial {\mathrm{N}}_i}{\partial X^A}-\underset{\partial {\mathcal{B}}_e}{\int }{\mathfrak{t}}_a{\mathrm{N}}_i-\underset{{\mathcal{B}}_e}{\int }{f}_a{\mathrm{N}}_i]\right],$$

(90a)

$$\{{\mathrm{F}}_e^{\mathcal{V}}\}=\left[[\underset{{\mathcal{B}}_e}{\int }{D}^A\frac{\partial {\mathrm{N}}_i}{\partial X^A}+\underset{\partial {\mathcal{B}}_e}{\int }\mathrm{\Omega }{\mathrm{N}}_i]\right],$$

(90b)

where we recall that the first Piola–Kirchhoff stress $\mathit{T}$ and dielectric induction $\mathit{D}$ are given by the constitutive equations (51), and the surface traction $\mathfrak{t}$ and surface charge density $\mathrm{\Omega }$ are given by equations (47d) and (34b), respectively.

In equation (88), the double brackets $\left[\left[·\right]\right]$ mean ‘matrix with components given by’ and, in equation (90), they mean ‘column matrix with components given by’. Care must be exercised in the use of the free indices $a$ and $b$ (spatial indices) and the indices $i$ and $j$ (nodal indices). The matrices and column vectors in equations (88) and (90) are built by mapping the pair of indices $(a,i)$ into an index $\alpha$ and the pair $(b,j)$ into an index $\beta$ , and both $\alpha$ and $\beta$ run from $1$ to $3m$ . When index $a$ is not present, as in equations (88c), (88d) and (90b), then we simply impose $\alpha =i$ and $\alpha$ runs from $1$ to $m$ . Similarly, when index $b$ is not present, as in equations (88b) and (88d), then we impose $\beta =j$ and $\beta$ runs from $1$ to $m$ . The mapping, for the case in which $a$ is present and the case in which and $b$ is present, is shown in Table 1.

Table 1.

Mapping the pairs of indices $(a,i)$ and $(b,j)$ .

$i$	$1$			$2$			…	$m$
$a$	1	2	3	1	2	3	…	1	2	3
$\alpha$	1	2	3	4	5	6	…	$3m-2$	$3m-1$	$3m$
$j$	$1$			$2$			…	$m$
$b$	1	2	3	1	2	3	…	1	2	3
$\beta$	1	2	3	4	5	6	…	$3m-2$	$3m-1$	$3m$

It is often convenient to express all quantities above in a local system of coordinates $({\gamma }_e^1,{\gamma }_e^2,{\gamma }_e^3)$ in the $\mathrm{eth}$ element. We define the transformations, each the inverse of the other,

$$X^A={\hat{X}}^A\left({\gamma }_e\right)\equiv {\hat{X}}^A({\gamma }_e^1,{\gamma }_e^2,{\gamma }_e^3),$$

(91a)

$${\gamma }_e^L={\hat{\gamma }}_e^L\left(X\right)\equiv {\hat{\gamma }}_e^L(X^1,X^2,X^3).$$

(91b)

The shape functions transform as

$${\tilde{\mathrm{N}}}_i\left({\gamma }_e\right)={\mathrm{N}}_i\left(X\right),$$

(92)

and their derivatives transform as

$$\frac{\partial {\tilde{\mathrm{N}}}_i}{\partial {\gamma }_e^L}\left({\gamma }_e\right)=\frac{\partial {\mathrm{N}}_i}{\partial X^B}\left(X\right)\frac{\partial {\hat{X}}^B}{\partial {\gamma }_e^L}\left({\gamma }_e\right),$$

(93a)

$$\frac{\partial {\mathrm{N}}_i}{\partial X^B}\left(X\right)=\frac{\partial {\tilde{\mathrm{N}}}_i}{\partial {\gamma }_e^L}\left({\gamma }_e\right)\frac{\partial {\hat{\gamma }}_e^L}{\partial X^B}\left(X\right),$$

(93b)

where the transformation matrices are clearly the inverse of each other:

$$\left[\left[\frac{\partial {\hat{X}}^B}{\partial {\gamma }_e^L}\left({\gamma }_e\right)\right]\right]={\left[\left[\frac{\partial {\hat{\gamma }}_e^L}{\partial X^B}\left(X\right)\right]\right]}^{-1}.$$

(94)

Finally, the element-level equation (85) must be expanded to the global dimension $4n$ , where $n$ is the total number of nodes in the finite element mesh of the body $\mathcal{B}$ and $4n$ is the total number of degrees of freedom of the problem [99,100]. We shall briefly outline the procedure later, in section 5. With the expansion, we obtain the global incremental equations

$$\left[\begin{array}{cc}\underset{3n\times 3n}{[{\mathrm{K}}^{\varphi \varphi }]}\hfill & \underset{3n\times n}{[{\mathrm{K}}^{\varphi \mathcal{V}}]}\hfill \\ \multicolumn{1}{c}{\underset{n\times 3n}{[{\mathrm{K}}^{\mathcal{V}\varphi }]}}& \underset{n\times n}{[{\mathrm{K}}^{\mathcal{V}\mathcal{V}}]}\hfill \end{array}\right]\left[\begin{array}{c}\underset{3n\times 1}{\{\overline{\eta }\}}\hfill \\ \multicolumn{1}{c}{\underset{n\times 1}{\{\overline{\mathcal{Z}}\}}}\end{array}\right]=-\left[\begin{array}{c}\underset{3n\times 1}{\{{\mathrm{F}}^{\varphi }\}}\hfill \\ \multicolumn{1}{c}{\underset{n\times 1}{\{{\mathrm{F}}^{\mathcal{V}}\}}}\end{array}\right].$$

(95)

In compact form, we have

$$\underset{4n\times 4n}{[\mathrm{K}]}\underset{4n\times 1}{\left\{ū\right\}}=-\underset{4n\times 1}{\left\{\mathrm{F}\right\}},$$

(96)

where $[\mathrm{K}]$ , $\left\{ū\right\}$ and $\{\mathrm{F}\}$ are the global stiffness matrix, the global vector of the generalised nodal displacements and the global vector of the generalised nodal forces, respectively.

The system of equation (95) is iteratively solved for the increments in Lagrangian displacement field $\overline{\eta }$ and the electric field $\overline{\mathcal{Z}}$ . At each iteration $z$ , the nodal values of the configuration $\varphi$ and the electric field $\mathcal{V}$ are updated from the increments, i.e.,

$${\underset{3n\times 1}{\{\overline{\varphi }\}}}^{z+1}={\underset{3n\times 1}{\{\overline{\varphi }\}}}^z+\underset{3n\times 1}{\{\overline{\eta }\}},{\underset{n\times 1}{\{\overline{\mathcal{V}}\}}}^{z+1}={\underset{n\times 1}{\{\overline{\mathcal{V}}\}}}^z+\underset{n\times 1}{\{\overline{\mathcal{Z}}\}}.$$

(97)

This iterative process is repeated until the set convergence criterion is met.

4.5. Constitutive equations

We confine our attention to a nearly incompressible electrostrictive polyurethane with rubber-like properties [62], which undergoes large deformations.

The referential energy (42), which we report here for convenience,

$$W=\hat{W}°(\mathit{F},\mathit{E},\mathcal{X})=J_{\mathit{K}}^{-1}\hat{\mathfrak{W}}°(\mathfrak{F},\mathfrak{E})=\stackrel{\vee }{W}°(\mathit{F},\mathit{E},\mathit{K}),$$

must be frame-indifferent [67,101,102]. Frame-indifference requires that the referential constitutive functions $\hat{W}$ and $\stackrel{\vee }{W}$ must depend on the deformation gradient $\mathit{F}$ through the right Cauchy–Green deformation tensor $\mathit{C}={\mathit{F}}^T\mathit{g}\mathit{F}$ defined in equation (7) and that the archetypal constitutive function $\hat{\mathfrak{W}}$ must depend on the elastic deformation $\mathfrak{F}$ through the elastic right Cauchy–Green tensor [103,104]

$$\mathfrak{C}={\mathfrak{F}}^T\mathit{g}\mathfrak{F}={\mathit{K}}^T{\mathit{F}}^T\mathit{g}\mathit{F}\mathit{K}={\mathit{K}}^T\mathit{C}\mathit{K}.$$

(98)

Indeed, the entirely material right Cauchy–Green tensor $C_{\mathrm{AB}}={\mathrm{F}}_A^a{g}_{\mathrm{ab}}{\mathrm{F}}_B^b$ and the completely archetypal elastic right Cauchy–Green tensor ${\mathfrak{C}}_{\alpha \beta }={\mathfrak{F}}_{\alpha }^a{g}_{\mathrm{ab}}{\mathfrak{F}}_{\beta }^b$ are unaffected by spatial rotations, which are performed on the spatial legs of the tensors involved (lowercase Latin indices). Finally, the constitutive function $\stackrel{\vee }{W}$ must depend explicitly on the uniformity tensor $\mathit{K}$ . Thus, we redefine the energy (42) as

$$W=\hat{W}°(\mathit{C},\mathit{E},\mathcal{X})=J_{\mathit{K}}^{-1}\hat{\mathfrak{W}}°(\mathfrak{C},\mathfrak{E})=\stackrel{\vee }{W}°(\mathit{C},\mathit{E},\mathit{K}).$$

(99)

As usually assumed in plasticity, we impose that the anelastic deformation described by the uniformity tensor $\mathit{K}$ be isochoric (i.e., volume-preserving), which implies, by virtue of the the multiplicative decomposition (41a),

$$\det \mathit{K}=1\Rightarrow J_{\mathfrak{F}}=\det \mathfrak{F}=\det \left(\mathit{F}\mathit{K}\right)=\det \mathit{F}=J_{\mathit{F}}.$$

(100)

We adopt an energy similar to that in the paper by Ask et al. [62], which is written with the volumetric–distortional decomposition of $\mathit{F}$ and $\mathit{C}$ , seen in equation (8). As highlighted by Simo et al. [105], Lubliner [106] and Bonet [107], since the plastic deformation (the inverse of the uniformity tensor $\mathit{K}$ ) is isochoric, multiplying both sides of the multiplicative decomposition (41a) by $J_{\mathit{F}}^{-1/3}=J_{\mathfrak{F}}^{-1/3}$ gives

$$\overline{\mathfrak{F}}=J_{\mathfrak{F}}^{-1/3}\mathfrak{F}=J_{\mathit{F}}^{-1/3}\mathfrak{F}=J_{\mathit{F}}^{-1/3}\mathit{F}\mathit{K}=\overline{\mathit{F}}\mathit{K},$$

(101)

where $\overline{\mathfrak{F}}$ is the isochoric part of the elastic deformation. Thus, we rewrite equation (99) as

$$W=\hat{W}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathcal{X})=\hat{\mathfrak{W}}°(J_{\mathit{F}},\overline{\mathfrak{C}},\mathfrak{E})=\stackrel{\vee }{W}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K}),$$

(102)

where $\overline{\mathfrak{C}}={\mathit{K}}^T\overline{\mathit{C}}\mathit{K}$ . The archetypal energy is

$$\begin{array}{c}\multicolumn{1}{c}{\hat{\mathfrak{W}}°(J_{\mathit{F}},\overline{\mathfrak{C}},\mathfrak{E})=\frac{k}{2}{[J_{\mathit{F}}-1]}^2+\frac{\mu }{2}[\overline{\mathfrak{C}}:{\mathfrak{G}}^{-1}-3]+}\\ \multicolumn{1}{c}{+c_m(\left({\mathfrak{G}}^{-1}\overline{\mathfrak{C}}{\mathfrak{G}}^{-1}\right):\mathfrak{E}\otimes \mathfrak{E}-1)+c_e({\mathfrak{G}}^{-1}:\mathfrak{E}\otimes \mathfrak{E}),}\end{array}$$

(103)

where the parameters $k$ , $\mu$ , $c_m$ and $c_e$ denote the bulk modulus, shear modulus, coefficient of electrostriction, a parameter related to the dielectric permeability of the material and a plastic modulus, respectively, and ${\mathfrak{G}}^{-1}$ is the inverse archetypal metric tensor. Now, using the definition $\overline{\mathfrak{C}}={\mathit{K}}^T\overline{\mathit{C}}\mathit{K}$ and $\mathfrak{E}={\mathit{K}}^T\mathit{E}$ and the inverse plastic right Cauchy–Green tensor ${\mathit{C}}_{\mathit{K}}^{-1}=\mathit{K}{\mathfrak{G}}^{-1}{\mathit{K}}^T$ , we can express the energy equation (103) with the help equation (102) in the form

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{\vee }{W}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})=\frac{k}{2}{[J_{\mathit{F}}-1]}^2+\frac{\mu }{2}[\overline{\mathit{C}}:{\mathit{C}}_{\mathit{K}}^{-1}-3]+}\\ \multicolumn{1}{c}{+c_m(\left({\mathit{C}}_{\mathit{K}}^{-1}\overline{\mathit{C}}{\mathit{C}}_{\mathit{K}}^{-1}\right):\mathit{E}\otimes \mathit{E}-1)+c_e({\mathit{C}}_{\mathit{K}}^{-1}:\mathit{E}\otimes \mathit{E}).}\end{array}$$

(104)

Following Federico [79], the second Piola–Kirchhoff stress $\mathit{S}$ is expressed in terms of its pulled-back spherical and pulled-back deviatoric parts as

$$\mathit{S}=-p{J}_{\mathit{F}}{\mathit{C}}^{-1}+J_{\mathit{F}}^{-2/3}{\mathbb{M}}^{*}:\tilde{\mathit{S}}={\mathit{S}}_{\mathrm{Sp}{\mathrm{h}}^{*}}+{\mathit{S}}_{\mathrm{De}{\mathrm{v}}^{*}},$$

(105)

where ${\mathbb{M}}^{*}$ is the pull-back mixed deviatoric operator defined in equation (17) and the hydrostatic pressure $p$ and the second Piola–Kirchhoff pseudo-stress $\tilde{\mathit{S}}$ are defined by

$$-p=\frac{\partial \stackrel{\vee }{W}}{\partial J_{\mathit{F}}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K}),\tilde{\mathit{S}}=2\frac{\partial \stackrel{\vee }{W}}{\partial \overline{\mathit{C}}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K}),$$

(106)

and, for the energy (103), they take the values

$$-p=k[J_{\mathit{F}}-1],\tilde{\mathit{S}}=\mu {\mathit{C}}_{\mathit{K}}^{-1}+2c_m\left({\mathit{C}}_{\mathit{K}}^{-1}\mathit{E}\right)\otimes \left({\mathit{C}}_{\mathit{K}}^{-1}\mathit{E}\right).$$

(107)

The first Piola–Kirchhoff stress $\mathit{T}$ is obtained from the second Piola–Kirchhoff stress $\mathit{S}$ by pushing the first leg of $\mathit{S}$ forward and lowering the first index of the resulting tensor by means of the the spatial metric tensor $\mathit{g}$ , i.e.,

$$\mathit{T}=\mathit{g}\mathit{F}\mathit{S}.$$

(108)

Thus, by substituting equation (105) into equation (108), the first Piola–Kirchhoff stress tensor $\mathit{T}$ becomes

$$\mathit{T}=-p{J}_{\mathit{F}}{\mathit{F}}^{-T}+J_{\mathit{F}}^{-2/3}\mathit{g}\mathit{F}({\mathbb{M}}^{*}:\tilde{\mathit{S}})={\mathit{T}}_{\mathrm{Sp}{\mathrm{h}}^{*}}+{\mathit{T}}_{\mathrm{De}{\mathrm{v}}^{*}}.$$

(109)

The material electric induction vector $\mathit{D}$ (equation (51b)) is obtained as

$$\mathit{D}=-\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{E}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})=-2[c_m{\mathit{C}}_{\mathit{K}}^{-1}\overline{\mathit{C}}{\mathit{C}}_{\mathit{K}}^{-1}\mathit{E}+c_e{\mathit{C}}_{\mathit{K}}^{-1}\mathit{E}].$$

(110)

The electromechanical tangent stiffness tensors, in the form used in the incremental equations (85) and (88) are written as a function of $\mathit{F}$ , and must be related to the tensors obtained using the frame-indifferent energy (102). We shall forgive ourselves for the abuse of notation of using $\stackrel{\vee }{W}$ for both the function of $\mathit{F}$ and the function of $\mathit{C}$ .

We start from the second elasticity tensor

$$\mathbb{C}=4\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \mathit{C}\partial \mathbf{C}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K}).$$

(111)

Since the energy (103) is decoupled in $J_{\mathit{F}}$ and $\overline{\mathit{C}}$ , we have the simplified expression [79]

$$\begin{array}{c}\multicolumn{1}{c}{\mathbb{C}=3J(J\mathrm{K}-p){\mathbb{K}}^{♯*}+2\mathrm{Jp}{\mathbb{I}}^{♯*}+}\\ \multicolumn{1}{c}{+J^{-4/3}{\mathbb{M}}^{*}:\tilde{\mathbb{C}}:{\mathbb{M}}^{*T}+\frac{2}{3}{J}^{-2/3}(\mathit{C}:\tilde{\mathit{S}}){\mathbb{M}}^{♯*}-\frac{2}{3}[{\mathit{C}}^{-1}\otimes {\mathit{S}}_{\mathrm{De}{\mathrm{v}}^{*}}+{\mathit{S}}_{\mathrm{De}{\mathrm{v}}^{*}}\otimes {\mathit{C}}^{-1}],}\end{array}$$

(112)

where the fourth-order tensors ${\mathbb{I}}^{♯*}$ , ${\mathbb{K}}^{♯*}$ and ${\mathbb{M}}^{♯*}$ are defined in Equation (18), and the large-deformation bulk modulus and the fictitious elasticity tensor are

$$\mathrm{K}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial J_{\mathit{F}}^2}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})\equiv k,\tilde{\mathbb{C}}=4\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \overline{\mathit{C}}\partial \overline{\mathit{C}}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})=\mathbb{O}.$$

(113)

For the energy (103), the former reduces to the bulk modulus $k$ , and the latter vanishes. The second elasticity tensor (112) then becomes

$$\mathbb{C}=3J(Jk-p){\mathbb{K}}^{♯*}+2\mathrm{Jp}{\mathbb{I}}^{♯*}+\frac{2}{3}{J}^{-2/3}(\mathit{C}:\tilde{\mathit{S}}){\mathbb{M}}^{♯*}-\frac{2}{3}[{\mathit{C}}^{-1}\otimes {\mathit{S}}_{\mathrm{De}{\mathrm{v}}^{*}}+{\mathit{S}}_{\mathrm{De}{\mathrm{v}}^{*}}\otimes {\mathit{C}}^{-1}].$$

(114)

The first elasticity tensor [68] is obtained from the second as

$$\mathbb{A}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \mathit{F}\partial \mathit{F}}°(J_{\mathit{F}},\overline{\mathit{F}},\mathit{E},\mathit{K})=\mathit{g}\underset{¯}{\otimes }\mathit{S}+\left(\mathit{g}\mathit{F}\underset{¯}{\otimes }\mathit{I}\right):\mathbb{C}:\left({\mathit{F}}^T\mathit{g}\underset{¯}{\otimes }{\mathit{I}}^T\right),$$

(115a)

$${A_a}_c^D=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial {\mathrm{F}}_B^a\partial {\mathrm{F}}_D^c}°(J_{\mathit{F}},\overline{\mathit{F}},\mathit{E},\mathit{K})=g_{\mathrm{ab}}{S}^{\mathrm{AB}}+g_{\mathrm{ab}}{\mathrm{F}}_A^b{C}^{\mathrm{ABCD}}{g}_{\mathrm{cd}}{\mathrm{F}}_C^d.$$

(115b)

The electro-elastic stiffness tensor is obtained most easily by differentiating the expression for $-\mathit{D}$ , obtained by taking the negative of equation (110), with respect $\mathit{F}$ and considering the derivative of $\mathit{C}$ with respect to $\mathit{F}$ , i.e.,

$$\mathbf{D}=\frac{\partial \stackrel{\vee }{W}}{\partial \mathit{F}\partial \mathit{E}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})=4c_m{J}_{\mathit{F}}^{-2/3}[\mathit{g}\mathit{F}{\mathbb{M}}^{*}:\left({\mathit{C}}_{\mathit{K}}^{-1}\underset{¯}{\otimes }{\mathit{C}}_{\mathit{K}}^{-1}\right)]\mathit{E},$$

(116a)

$${{\mathrm{D}}_a}^{\mathrm{BC}}=\frac{\partial \stackrel{\vee }{W}}{\partial {\mathrm{F}}_B^a\partial E_C}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})=4c_m{J}_{\mathit{F}}^{-2/3}\left[g_{\mathrm{ab}}{\mathrm{F}}_A^b{\left({\mathbb{M}}^{*}\right)}_{\mathrm{DP}}^{\mathrm{AB}}{\left({\mathit{C}}_{\mathit{K}}^{-1}\right)}^{\mathrm{CD}}{\left({\mathit{C}}_{\mathit{K}}^{-1}\right)}^{\mathrm{PQ}}\right]E_Q,$$

(116b)

where we considered the symmetry of ${\mathit{C}}_{\mathit{K}}^{-1}$ in ${\left({\mathit{C}}_{\mathit{K}}^{-1}\right)}^{\mathrm{CD}}={\left({\mathit{C}}_{\mathit{K}}^{-1}\right)}^{\mathrm{DC}}$ . Finally, the dielectric permittivity tensor is

$$\mathit{Y}=\frac{{\partial }^2\stackrel{\vee }{W}}{\partial \mathit{E}\partial \mathit{E}}°(J_{\mathit{F}},\overline{\mathit{C}},\mathit{E},\mathit{K})=-2[c_m{\mathit{C}}_{\mathit{K}}^{-1}\overline{\mathit{C}}{\mathit{C}}_{\mathit{K}}^{-1}+c_e{\mathit{C}}_{\mathit{K}}^{-1}].$$

(117)

5. Algorithm: finite element implementation in MATLAB

In this section, we describe a step-by-step guide for analysing dielectric elastomers with the finite element method on MATLAB. This tailored walkthrough focuses on the core computations for displacements and electric potential, and the iterative equation-solving with the Newton–Raphson method.

In the example problem that we study, we consider self-driven evolution, as described in section 4.3. Specifically, we examine two self-driven evolution laws: the diagonal case and the nilpotent case, both at steady state.

5.1. Initialisation

• Set maximum Newton–Raphson iterations $\mathrm{ite}{\mathrm{r}}_{\max }=10,000$ and error tolerance $\mathrm{to}{\mathrm{l}}_{\mathrm{error}}=1.0\times {10}^{-9}$ .

• Set the number of Gauss integration points ${\mathrm{N}}_{\mathrm{Gauss}}=27$ .

• Initialise matrices for referential nodal coordinates.

• Initialise assembly matrices (the matrices expanding from element dimension to global dimension).

5.2. Mesh construction and Gauss integration setup

• Construct the connectivity array for hexahedral (brick-like) elements.

• Define trilinear shape function ${\mathrm{N}}_i$ for each node $i$ in the $m$ -node element $e$ .

• Construct the $3m\times 1$ vector of the nodal coordinates ( $3$ coordinates for each of the $m$ nodes).

• Define the local coordinates $({\gamma }_e^1,{\gamma }_e^2,{\gamma }_e^3)$ of the Gauss integration points.

• Define the weights $(W_{{\gamma }_e^1},W_{{\gamma }_e^2},W_{{\gamma }_e^3})$ of the local coordinates $({\gamma }_e^1,{\gamma }_e^2,{\gamma }_e^3)$ .

• Set up boundary conditions.

• Initialise $3n\times 1$ global nodal motion $\{\overline{\varphi }\}$ and $n\times 1$ global nodal electric potential $\{\overline{\mathcal{V}}\}$ of equation (97).

5.3. Loop for each Newton–Raphson iteration

5.3.1. Newton–Raphson loop

For each Newton–Raphson iteration ( $\mathrm{ite}{\mathrm{r}}_{\mathrm{Newton}}=1$ to $\mathrm{ite}{\mathrm{r}}_{\max }$ ):

• Initialise $4n\times 4n$ global stiffness matrix $[\mathrm{K}]$ and $4n\times 1$ global nodal force vector $\{\mathrm{F}\}$ .

• Loop over all elements ( $\mathrm{ite}{\mathrm{r}}_{\mathrm{el}}=1$ to ${\mathrm{N}}_{\mathrm{el}}$ ), where ${\mathrm{N}}_{\mathrm{el}}$ is the total number of elements.

• Assemble $4n\times 4n$ global stiffness matrix $[\mathrm{K}]$ and $4n\times 1$ global generalised force vector $\{\mathrm{F}\}$ from local matrices and vectors.

• Apply boundary conditions to $[\mathrm{K}]$ and $\{\mathrm{F}\}$ .

• Solve for $3n\times 1$ global motion $\{\overline{\varphi }\}$ and $n\times 1$ global electric potential $\{\overline{\mathcal{V}}\}$ .

• Check convergence: if error is below $\mathrm{to}{\mathrm{l}}_{\mathrm{error}}$ , exit the loop.

5.3.2. Post-processing

• Output the solution to a Tecplot ASCII file for visualisation

• Display current iteration number and RMS (root mean square) error, calculated as part of the Newton–Raphson iterative process, as

$$RMS=\sqrt{\frac{1}{4n}{\underset{4n\times 1}{\left\{ū\right\}}}^T\underset{4n\times 1}{\left\{ū\right\}}},$$

where $n$ is the total number of nodes, so that $4n$ is the total number of degrees of freedom

• Write the solution to separate files for visualisation in a software package such as VisIt (open source) or Tecplot (Bellevue, WA, USA)

5.3.3. Plot convergence history

• Plot the RMS error against Newton–Raphson iterations $\mathrm{ite}{\mathrm{r}}_{\mathrm{Newton}}$ to visualise the convergence of the solution.

6. Numerical simulations

To evaluate the finite element implementation, we consider a beam-type actuator capable of bending in response to an applied electric field. The actuator is shown in Figure 3 and it comprises two perfectly bonded layers made of the same dielectric elastomer with identical material properties. The top layer is sandwiched between two (infinitely) compliant electrodes. Each layer has a thickness of $0.5\mathrm{mm}$ , and the overall dimensions are $50\mathrm{mm}$ in length ( $X^1$ -direction), $4\mathrm{mm}$ in width ( $X^2$ -direction) and $1\mathrm{mm}$ in thickness ( $X^3$ -direction). A potential difference $\Delta \mathcal{V}$ is applied between the two electrodes sandwiching the layer and Dirichlet boundary conditions are prescribed:

• At the fixed end $(X^1=0.0)$ : The motion $\varphi$ is constrained to zero in all three directions;

• At the bottom layer $(X^3=0.0)$ : The electric potential is grounded, i.e, $\mathcal{V}=0$ , which serves as the reference potential level;

• At the midpoint in the $X^3$ direction: The electric potential $\mathcal{V}$ is also set to zero;

• Between the bottom layer and midpoint in the $X^3$ direction: The electric potential $\mathcal{V}$ varies linearly from 5000 units just above the top layer to zero at the midpoint, creating a potential gradient.

Figure 3.

The bilayer dielectric elastomer beam of our example.

The mesh parameters, domain limits and material parameters (the units are arbitrary, but assumed consistent) are defined as in Table 2.

Table 2.

Simulation parameters; in the material parameters, the units are arbitrary but assumed consistent.

Domain limits		Mesh parameters		Material parameters
Direction	Range (mm)	Axis	Elements	Parameter	Value
$X^1$	$0.0$ to $50.0$	$X^1$	$50$	$k$	$1.830\times {10}^8$
$X^2$	$0.0$ to $4.0$	$X^2$	$8$	$\mu$	$1.307\times {10}^6$
$X^3$	$0.0$ to $1.0$	$X^3$	$6$	$c_e$	$1.000\times {10}^0$
				$c_m$	$1.763\times {10}^{-9}$
				$\beta$	$1.307\times {10}^6$
				$C_3$	$0.22$
				$C_4$	$0.014$
				$C_5$	$0.035$
				$C_6$	$0.22$

6.1. Case of homogeneous material

Let us first analyse the scenario in which $\mathit{K}$ is an orthogonal tensor independent of the point $X$ and therefore trivially integrable: thus, the body is homogeneous. The left panel in Figure 4 shows the homogeneous beam deformation under a strong electric field. This behaviour is desirable for applications requiring consistent and predictable actuation upon electrical stimulation. The colour gradient represents the potential difference across the beam, with the deformation pattern offering a visual representation of the actuator response. The RMS error analysis for the homogeneous case is presented in the right panel of Figure 4. This graph offers insights into the numerical model convergence, particularly emphasising the stability of the simulation. The RMS error plot is crucial for assessing the computational model performance and ensuring the reliability of the simulation results.

Figure 4.

Deformation of a homogeneous beam under an applied electric field (left) and corresponding RMS error analysis (right).

6.2. Case of uniform material with self-driven evolution

Two distinct forms of the uniformity tensor $\mathit{K}$ are assumed based on evolution laws driven by the torsion $\mathbf{Tor}$ , as outlined in section 4.3: the diagonal case and the nilpotent case. The deformations are profoundly different from the case of homogeneous material seen in section 6.1.

The left panel of Figure 5 shows the deformation in the case in which the uniformity tensor $\mathit{K}$ is a diagonal tensor featuring the function $A(X^1,X^2)$ , as presented in equation (67). The function $A$ , generated based on the self-driven evolution law in equation (64), takes the form shown in equation (75), and causes the beam to warp towards the left, as shown in Figure 5. The right panel of Figure 5 shows the corresponding RMS plot.

Figure 5.

Warping caused by a diagonal uniformity tensor $\mathit{K}$ (left) and corresponding RMS error analysis (right).

The left panel of Figure 6 illustrates the case of nilpotent uniformity tensor $\mathit{K}$ , featuring the function $\mathit{K}$ , as presented in equation (76). This function was obtained from the torsion-induced evolution law as an arbitrary function of $X^2$ , in equation (82). For simplicity, we studied case of constant function $B$ , i.e., $B(X^1,X^2)=\pm B_0$ . The left panel of Figure 6 shows the case $B(X^1,X^2)=+B_0$ . The case $B(X^1,X^2)=-B_0$ is identical, except that it is mirrored with respect to the $X^1$ - $X^3$ plane. The RMS plot, identical for both cases $B(X^1,X^2)=\pm B_0$ , is reported in the right panel of Figure 6.

Figure 6.

Warping caused by a nilpotent uniformity tensor $\mathit{K}$ with $B=+B_0$ (left) and corresponding RMS error analysis (right).

Figure 7 is a comparison of all cases studied: homogeneous (uncoloured meshed beam), diagonal and nilpotent. Both cases of nilpotent case, $B(X^1,X^2)=\pm B_0$ are shown: note how the two cases mirror each other.

Figure 7.

Comparison of the deformation of a homogeneous beam actuator (uncoloured mesh) with the deformation of beam actuators influenced by a diagonal or nilpotent uniformity tensor $\mathit{K}$ , in response to an applied electric field.

7. Summary and outlook

In this study, we presented a model for uniform electro-elastic materials, in which we combined the theory of dielectric elastomers, rooted in the foundational works by Toupin [18–20] and others, the theory of material uniformity of Noll [4] in the formulation by Epstein and Maugin [5], and the theory of configurational mechanics of Eshelby [96,97], in the picture by Gurtin [85,86], Cermelli et al. [87] and DiCarlo and Quiligotti [88]. In this picture, the uniformity tensor $\mathit{K}$ is seen as an additional non-standard kinematical variable, acted upon by non-standard configurational forces. For this preliminary study, we chose the particular case of what Epstein and Elżanowski [8] call self-driven evolution, in which the internal configurational force ${\mathit{Y}}_i$ coincides with the negative of the Eshelby stress $\mathit{H}$ (equation (60)) and therefore does not dissipate (equation (59)).

In our finite element implementation, we specifically selected a case where the uniformity tensor $\mathit{K}$ was predefined from the outset, with a torsion-driven evolution law, equation (64), which may not necessarily descend from the second law of thermodynamics. We developed a MATLAB-based finite element solver, capable of handling various electromechanical simulations and studied, as a benchmark problem, a beam-type actuator. The example actuator studied here consists of two dielectric elastomer layers, with the top layer sandwiched between two infinitely compliant and thin electrodes, one set to a certain potential, and the other grounded to zero. Other configurations, which are representative of actuators commonly found in industrial practice, could include three electrodes. For instance, the top electrode could be set to a positive potential ${\mathcal{V}}_0$ , the middle one grounded to zero and the bottom electrode at a negative potential $-{\mathcal{V}}_0$ . Three cases were studied:

• A perfectly homogeneous beam;

• A beam with a diagonal uniformity tensor, characterised by the function $A(X^1,X^2)$ ;

• A beam with a nilpotent uniformity tensor, characterised by the shear function $B(X^1,X^2)$ .

The diagonal case resulted in the beam bending and warping towards the left. The nilpotent case resulted in a different warping behaviour depending on the sign of the shear function $B$ : positive values led to warping to the right, while negative values caused warping to the left.

The analysis confirmed the convergence of the numerical solution in all scenarios and demonstrated the capability of our in-house finite element solver to accurately capture the complex interactions between electrical stimulation and material inhomogeneity in this type of actuator.

Our next steps shall be the study of two cases. The first consists in treating the uniformity tensor $\mathit{K}$ as an internal variable, obtained via an evolution law abiding the second law of thermodynamics, in the approach of Epstein and Maugin [9] and Epstein and Elżanowski [8]. The second consists in treating the uniformity tensor $\mathit{K}$ as an additional kinematical variable, as in section 4 of this work, in the approach of Gurtin [85,86]. However, rather than assuming that the internal configurational force ${\mathit{Y}}_i$ equals the negative of the Eshelby stress $\mathit{H}$ , we shall define the dissipative part of ${\mathit{Y}}_i$ constitutively [87,89–94]. This will require the modification of our finite element solver, to include the additional degrees of freedom associated with the uniformity tensor $\mathit{K}$ .