Stability analysis of charge-controlled soft dielectric plates

We examine the stability of a soft dielectric plate deformed by the coupled effects of a mechanical pre-stress applied on its lateral faces and an electric field applied through its thickness under charge control. The electric field is created by spraying charges on the major faces of the plate: although in practice this mode of actuation is harder to achieve than a voltage-driven deformation, here we find that it turns out to be much more stable in theory and in simulations. First we show that the electromechanical instability based on the Hessian criterion associated with the free energy of the system does not occur at all for charge-driven dielectrics for which the electric displacement is linear in the electric field. Then we show that the geometric instability associated with the formation of small-amplitude wrinkles on the faces of the plate that arises under voltage control does not occur either under charge control. This is in complete contrast to voltage-control actuation, where Hessian and wrinkling instabilities can occur once certain critical voltages are reached. For the mechanical pre-stresses, two modes that can be implemented in practice are used: equi-biaxial and uni-axial. We confirm the analytical and numerical stability results of homogeneous deformation modes with Finite Element simulations of real actuations, where inhomogeneous fields may develop. We find complete agreement in the equi-biaxial case, and very close agreement in the uni-axial case, when the pre-stress is due to a dead-load weight. In the latter case, the simulations show that small inhomogeneous effects develop near the clamps, and eventually a compressive lateral stress emerges, leading to a breakdown of the numerics.


Introduction
Soft dielectric materials can undergo large actuation stretches when a potential difference is induced in the material. Typically, compliant electrodes such as carbon grease are smeared onto the faces of a soft dielectric elastomer plate and a voltage is applied across the thickness of the material. As the voltage increases the material gradually expands in area until a maximum voltage is reached, at which point a rapid large deformation known as snap-through occurs [29]. The large actuation achieved due to the snap-through behaviour is desirable for many applications but is difficult to achieve in practice. Snapthrough is often prevented by electric breakdown [29] or by instabilities such as inhomogeneities [1], compression failure [5], band localisation [10], wrinkles [17,20,21,28,24,3], membrane wrinkling [12], etc.
Various methods have been proposed for avoiding electric breakdown without sacrificing the large actuation. For example, if the material is pre-stretched before the voltage is applied, electric breakdown may be avoided, but the stretch gain achieved might be reduced [23]. Another method proposed is charge-controlled actuation, as shown experimentally by Keplinger et al. [15] and theoretically by Li et al. [16]. In charge-controlled actuation, charges of opposite signs are sprayed on opposite planar surfaces of a dielectric plate, inducing a potential difference, and hence an electric field in the dielectric, thereby inducing a deformation. In principle this method of actuation annihilates the possibility of snap-through because the theoretical charge-stretch loading curves are monotonic [16]. In this paper, we investigate the stability of a charge-driven dielectric plate, which has not been considered previously and the results of which are significantly different from those for voltage control.
We first focus on equi-biaxial loading and show that charge-controlled actuation is stable since the Hessian criterion-or rather, its version for this problem-for onset of instability is never met (Section 2.2). This result is far from straightforward to obtain, because the Hessian determinant of the energy density is always negative, from which it could erroneously be concluded that the actuation is unstable. In fact, we show that the second variation of the free energy of the whole system is always positive, which ensures stability throughout. This is in sharp contrast to the corresponding situation for voltagecontrolled actuation, which is well-known [29] to become unstable once a critical voltage is reached.
We then highlight another new, and complementary, feature of charge control by showing that charge-controlled actuation is also stable with respect to geometric instability because, provided the material is pre-stretched, small-amplitude inhomogeneous wrinkled solutions superposed on the large homogeneous actuation do not develop (Section 3). Again, this contrasts with the situation for voltage-controlled actuation, for which dielectric plates eventually wrinkle under sufficiently large voltages [7,8,23].
In Section 4 we model the experiments of Keplinger et al. [15] where a plate was prestretched by a weight prior to charge-controlled actuation. We thus study the stability of a homogeneously deforming plate under uni-axial tension and charge-actuation and again we find Hessian-based and geometric stability in this case, again contrary to the corresponding situation for voltage-controlled actuation.
Finally in Section 5 we use Finite Element simulations to account for the finite dimensions of plates. We find that in the equi-biaxial case there are no differences between the results of the homogeneous loading analytical modelling and those of the Finite Element method, because the plate is free to stretch laterally and the loading curves are indeed monotonic (no snap-through). However, for the uni-axial case we find that the clamping of the plate required to apply the weight leads to non-homogeneous deformations with local variations of stresses and strains compared to the homogeneous solution, and that these effects build up and eventually lead to a breakdown of the simulation. We identify the point of breakdown as corresponding to the appearance of compressive stresses in the plate.

Equations of electroelasticity
Consider the stress-free reference configuration B r of an electroelastic material in the absence of an electric field and applied mechanical loads. Points in B r are labelled by the position vector X. When subject to loads and an electric field under static conditions the material occupies the configuration B, with the material point X now at x. Let F = Gradx denote the deformation gradient from B r to B, where Grad is the gradient operator with respect to X. We denote by E and D, respectively, the electric field and electric displacement vectors in B, and by τ the Cauchy stress tensor (which in general depends on F and either E or D).
It has been found advantageous [6,29] to formulate constitutive equations in terms of the Lagrangian field variables, denoted E L , D L , and the nominal stress tensor T , which are related to E, D and τ by the following pull-back operations (from B to B r ) where J = det F . The constitutive equations are based on the use of so-called 'total' energy functions, depending either on F and E L , denoted Ω, or on F and D L , denoted Ω * , with the (partial) Legendre transform connection Henceforth, we confine attention to incompressible materials, so that the constraint J ≡ 1 is in force. Then we have the constitutive relations (where p and p * are Lagrange multipliers associated with the constraint, in general with p * = p), and The governing equations are Div T = 0, where Div and Curl are the divergence and curl operators with respect to X. We shall consider the situation in which there is no external field, so that on the boundary ∂B r of B r the standard electric boundary conditions associated with the equations (5) are simply where N is the unit outward normal on ∂B r , t A is the applied mechanical traction per unit area of ∂B r and σ F is the surface charge density per unit area of ∂B r .
In considering applications to dielectric elastomers, which are isotropic electroelastic materials, the functional dependence of Ω and Ω * can be expressed in terms of five invariants. First of all, the isotropic purely kinematic invariants defined by where c = F T F is the right Cauchy-Green deformation tensor. Secondly, invariants associated with E L , which typically are taken to be as in [6], and, thirdly, invariants associated with D L , here defined by as used in [6] in different notation. The expanded forms of the constitutive relations (3), when converted to Eulerian form using (1) with J = 1, are where I is the identity tensor, b = F F T is the left Cauchy-Green deformation tensor, and Ω * i = ∂Ω * /∂I * i , i = 4, 5, 6.

Specialization to biaxial deformations of a plate
We now consider the application of the above theory to the biaxial deformation of a rectangular plate.
The plate has sides of lengths L 1 , L 2 , L 3 in the reference configuration B r , where L 2 = H is the thickness of the plate, which is small compared with its lateral dimensions. Mechanical loads are applied in the 1 and 3 directions; also, a potential difference, say V , exists between the major surfaces of the plate and the associated charges on the surfaces are denoted ±Q. As a result, the plate is stretched homogeneously with stretches λ 1 and λ 3 parallel to the major surfaces, and, by incompressibility, a stretch λ 2 = λ −1 1 λ −1 3 normal to the major surfaces. The potential difference generates an electric field with a single component E = E 2 , associated with an electric displacement component D = D 2 . The corresponding components of the Lagrangian fields are E L = λ 2 E and D L = λ −1 2 D. In terms of the potential difference V and the associated charges ±Q on the surfaces, we have the simple connections Thus, for a fixed potential, E L is fixed, while fixed charge Q corresponds to fixed D L . For this combination of deformation and electric field, Ω and Ω * specialize accordingly. The invariants are now given in terms of the independent stretches λ 1 , λ 3 and E L and D L by We denote the specializations of Ω and Ω * by ω and ω * , respectively, and the independent variables by (λ 1 , λ 3 , E L ) and (λ 1 , λ 3 , D L ), respectively, with, from the connection (2), Since the resulting deformation is purely biaxial the corresponding nominal stress is coaxial with the edges of the plate; we denote its components by t 1 , t 2 , t 3 .
We now assume that there is no mechanical traction on the major faces of the plate so that the boundary condition (6) 1 yields t 2 = 0. Then, on elimination of the hydrostatic stress from (10) and (11), we obtain the simple formulas and from (13) The particular case of equi-biaxial deformations is of special interest, for then, with λ 1 = λ 3 = λ, and incompressibility giving λ 2 = λ −2 , we may introduce the following further specialisations of the total energy functions, We also have t 1 = t 3 = t, say, so that For illustration, we now consider models for which D = εE, where ε, the material permittivity, is taken to be a constant. These are "ideal" dielectrics in the terminology of Suo [25]. Note that this linear relationship has recently been verified [31] using experimental data for low to moderate values of the electric field for the acrylic dielectric elastomer VHB 4905. In general, ε may depend on the deformation, as has been shown in [27], for example, for the acrylic dielectric elastomer VHB 4910, but for our present purposes we consider it to be a material constant.
Then, Ω and Ω * have the forms and, for the biaxial deformations of a plate considered above, where w(λ 1 , λ 3 ) = W (I 1 , I 2 ) with I 1 and I 2 given by (15) and D L = λ 2 1 λ 2 3 E L . For equibiaxial deformations, we havẽ wherew(λ) = w(λ, λ) and D L = ελ 4 E L . For our subsequent applications we consider two representative energy density functions, a neo-Hookean dielectric model and a Gent dielectric model, defined, in the two representations, by where µ is the shear modulus in the absence of an electric field and J m is a stiffening parameter. Notice that Ω G recovers Ω nH and Ω * G recovers Ω * nH in the limit J m → ∞. We now express the equations in dimensionless form by defining the following quantitiesω so that D 0 = λ 4 E 0 in the equi-biaxial case, and Based on eitherω orω * , we now obtain the expression for D 0 in terms of λ and s for the neo-Hookean and Gent dielectric models as respectively, and note that the latter reduces to the former when J m → ∞.  [11] has been used here. Also shown in Figure 1(a) is the curve D 0 = (λ 6 + 5)/3, which cuts the fixed s curves at points where E 0 is a maximum in Figure 1(b), which also shows the corresponding dashed curve. Similarly for the Gent dielectric in Figures 2(a) and 2(b), although, for the larger values of s, there is no maximum in (b) and no corresponding intersection.
We now turn to the analysis of the stability of the plate based on the Hessian criterion.

Analysis of the Hessian stability criterion
Electro-mechanical instability is often considered to occur when the Hessian matrix associated with the second variation of the free energy for the whole system ceases to be positive definite [29]. The rationale of this criterion is that equilibrium corresponds to an extremum of the free energy (and thus its first variation is zero), and that the equilibrium is stable when it corresponds to a minimum of the free energy (and then its second variation is positive).
In different notation and in dimensionless form, the free energy of the whole system, here denoted ψ * , considered in [29] has the form and vanishing of its first variation (for fixed s 1 , s 3 , E 0 ), with s 1 = t 1 /µ, s 3 = t 3 /µ, yields the dimensionless versions of the constitutive relations involving ω * in (19) and (20).  . If, instead, we use E 0 as the independent electric variable, then the corresponding 'energy', denoted ψ, vanishing of the first variation of which yields the constitutive relations in terms of ω in (19) and (20), is given by Note that, on use of (18) in dimensionless form, we have ψ * = ψ − D 0 E 0 , so that ψ is the Legendre transform of ψ * with respect to the conjugate variables E 0 and D 0 related by (30) 4 .
For the free energy ψ * of the whole system to be at a minimum, its second variation must be positive, i.e. the associated Hessian matrix must be positive definite, at a point of equilibrium. The second variations of ψ * and ψ are written compactly as with the subscripts representing partial derivatives. For the equi-biaxial case these become 2 × 2 matrices, given by and we now focus on this case for illustration. It is straightforward to show thatω * D 0 D 0 = −1/ω E 0 E 0 by using the formulas (30) 3,4 . Now the determinants of the Hessians above are given by and on specializing (18) we haveω * (λ, D 0 ) =ω(λ, E 0 ) + D 0 E 0 , from which the following connections, given in [8] in dimensional form, can be obtained: These equations are independent of the specific forms ofω * andω, and so are valid for any choice of (equi-biaxial) energy density function. They have some interesting interpretations, which we now discuss in respect of the neo-Hookean dielectric, for whicĥ and hence Thus here, Note that the maxima of E 0 in Figure 1(b) correspond to det H * = 0, equivalentlŷ ω λλ = 0, which also corresponds to a maximum of s at fixed E 0 . Note that H * is positive definite up to the maxima as E 0 is increased from 0. By contrast, s is monotonic with respect to λ at fixed D 0 andω * λλ > 0, while det H < 0 and H is indefinite, thus defining a saddle point ofω. Note that it would be incorrect to conclude here that the chargecontrolled actuation is unstable, as we now show.
Consider the connectionω * (λ, the first variation of which yields If we now take the second variation then the terms involving δ 2 λ, δ 2 E 0 , δ 2 D 0 cancel and we are left with the quadratic connectioñ and hence, by substituting for δD 0 on the right-hand side of (47), we obtaiñ For stability we require the left-hand side to be positive since this is the second variation of the actual free energy ψ * (so the free energy is minimized), whether we have voltage-control of the deformation (when λ and D 0 are free to vary) or charge-control of the deformation (when λ and E 0 are free to vary).
For fixed E 0 , in a voltage-controlled experiment, this reduces simply toω λλ > 0, and this fails where E 0 is a maximum. For the neo-Hookean dielectric, it reads λ −2 + 5λ −8 − 3E 2 0 > 0, and E 0 = (λ −2 + 5λ −8 )/3 is the plot going through the maxima of each loading curve for different values of the pre-load s, as shown by the dashed curve in Figure 1(b).
For fixed D 0 , in a charge-controlled experiment, the left-hand side is positive ifω * λλ > 0, and for the neo-Hookean dielectric, this reads 1+5λ −6 (1+D 2 0 ) > 0, which holds true for all D 0 . In the case of a perfect dielectric, we have D 0 = λ 4 E 0 , and hence 0 = 4λ 3 δλE 0 +λ 4 δE 0 , so for the right-hand side of (49) to be positive we havẽ which confirms the resultω * λλ > 0 for the neo-Hookean dielectric, and thus, that the second variation of the free energy is always positive.
We can therefore conclude that under charge control, equi-biaxial activation is stable according to the Hessian criterion since we haveω * λλ > 0 for the considered neo-Hookean model. On the other hand, activation under voltage control can become unstable in the Hessian criterion sense, as is well known, since the inequalityω λλ > 0 can fail. The results for the Gent model (not developed here) follow the same pattern.

Incremental stability analysis
To investigate the possibility of geometric instabilities, namely the formation of smallamplitude wrinkles on the faces of the plate, we linearise the governing equations and boundary conditions in the neighbourhood of a large deformation and initial electric field.
We introduce the incremental mechanical displacement u, the incremental nominal stress tensorṪ and the incremental Lagrangian electric field and displacement,Ė L anḋ D L , respectively, all of which are functions of the deformed position x [7]. LetṪ 0 ,Ė L0 anḋ D L0 denote their push-forward forms from the reference to the deformed configuration, as defined byṪ 0 = FṪ ,Ė L0 = F −TĖ L ,Ḋ L0 = FḊ L . These satisfy the governing and the relevant incremental constitutive equations arė where A 0 , A 0 and A 0 are, respectively, fourth-, third-and second-order electroelastic moduli tensors (see [24] for their general expressions), and L is the displacement gradient gradu, u being the incremental displacement, which, by incompressibility, satisfies trL ≡ divu = 0.
From (53) 3 we can introduce the scalar electric potential ϕ such thaṫ We now focus on models of the form for which the relevant components of the moduli tensors reduce to Note that we used the connection p = A 02121 , which is a special case of a general formula given in, for example, equation (9.88) of [7]. The vertical bar between the components of A 0 is used to distinguish the single index (associated with a vector) from the pair of indices associated with a second-order tensor. Now, for brevity, we introduce the notations Then, on elimination ofṗ and use of the incompressibility equation u 1,1 + u 2,2 = 0, the required incremental constitutive equations can be written compactly in the forṁ T 011 =Ṫ 022 + 2(b + c)u 1,1 + 2dϕ ,2 , T 012 = au 2,1 + cu 1,2 − dϕ ,1 ,Ṫ 021 = c(u 1,2 + u 2,1 ) − dϕ ,1 , We now convert the system of equations to a first-order system with six variables based on the Stroh approach. For this purpose we choose the variables u 1 , u 2 , ϕ,Ṫ 021 , T 022 ,Ḋ L02 and consider increments that are sinusoidal in the x 1 direction, i.e. solutions of the form where U 1 , U 2 , Φ, Σ 21 , Σ 22 and ∆ are all functions of kx 2 , k = 2π/L is the wave number and L is the wavelength of the wrinkles. We now arrange the variables so that they all have the same dimensions by defining a Stroh vector η as where U is the 'displacement' vector and S is the 'traction' vector. After a little manipulation, the equations (53) and (58) are cast in the form whereā = a/µ,b = b/µ,c = c/µ andd = d/ √ µε, and a prime denotes differentiation with respect to kx 2 . Similarly to Su et al. [24] we can thus write the equations in Stroh form, i.e. as η = iN η, where N is the Stroh matrix and η is the Stroh vector, defined in (60). Note that the vector η is different from its counterpart in the voltage-controlled case [24], due to the different electric boundary conditions and scalings. In the voltage-controlled case, the incremental electric boundary condition is in terms of the electric potential Φ (which must be zero on the faces), whereas in the charge-controlled case, the incremental electric boundary condition is in terms of the electric displacement ∆. In the present situation the Stroh matrix has the dimensionless form where For the models (23) and (24) for which W (I 1 , I 2 ) depends on only I 1 , i.e. W = W (I 1 ), including the Gent dielectric model (27), in equi-biaxial activation we havē whereW (I 1 ) = W (I 1 )/µ, and henceforth we restrict attention to this specialization. For the Gent dielectric, and we recall that E 0 = λ −4 D 0 . For the neo-Hookean dielectric, the expressions simplify considerably as:W = 1/2 andW = 0.
To investigate the conditions for wrinkling to occur, it is sufficient to calculate the thin-plate and thick-plate limits of the dispersion equation, as the behaviour of a plate with finite thickness lies in between the two [24].
The thin-plate limit is calculated from the Stroh matrix as [22,24] det N 3 = 0, which simplifies here to As in the voltage-controlled case, the thin-plate limit can be separated into symmetric and anti-symmetric modes. Anti-symmetric modes are governed by the equationā −c = 0, as in the voltage-controlled case. For the neo-Hookean and the Gent dielectric models this yields respectively, which is the same as (31) in the absence of pre-stress (s = 0). No symmetric modes are possible as they are governed by the equationb +c + 2d 2 = 0, which has no real solutions in (λ, D 0 ). Likewise, the third factor in (68) yields no solutions.
To calculate the thick-plate limit, we first construct a matrix with the eigenvectors η (j) , j = 1, 2, 3, of N with corresponding eigenvalues with positive imaginary part, stacked as the columns as follows where A and B defined above are 3 × 3 matrices and explicit expressions for the components of η (j) , j = 1, 2, 3, are given as follows for j = 2, 3 and where p 2,3 andW are given by (74) below and (66) 1 , respectively. Then the thick-plate limit is given by Based on the analysis of Stroh (see Ting [26] or Shuvalov [22], for instance), we recall that iBA −1 is Hermitian and the above equation is a single real equation, as distinct from det B = 0, which is a complex equation, although its real and imaginary parts are in proportion. For models with W = W (I 1 ), including that of Gent, equation (72) is a quadratic in where (−) and (+) correspond to p 2 and p 3 , respectively. Note that for the neo-Hookean specialization, sinceW = 1/2,W = 0, we obtain p 2 = 1 and p 3 = λ 3 and the thick-plate limit becomes In the absence of charge (D 0 = 0), this reduces to the classical elastic case and recovers the critical stretch for surface instability under equi-biaxial stretch of Green and Zerna [13], specifically λ = 0.666. + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + ++ + + + + + + + + + + + P P + + + + + + + + We plot the thick-and thin-plate limits, along with the loading curves (31) 1 for the neo-Hookean dielectric model for different values of pre-stress in Figure 3. The loading curves are monotonic, and so the material will not experience the snap-through phenomenon of voltage-controlled actuation [16]. As shown in the previous section, this is connected to to the sign of the second variation of the free energy being always positive.
These theoretical predictions are compared with Finite Element simulations (see Section 5), the results of which are represented by dots in the figure, which also exhibit the stability.
The region between the thick-plate and thin-plate limits is where wrinkling could occur. However, the pre-stretched loading curves do not cross this region, so there is no wrinkling. Charge-controlled dielectric plates are therefore geometrically stable, and will not exhibit wrinkling (provided s > 0). This situation again contrasts with voltagecontrolled plates, which can wrinkle in compression, as here, but also in extension [7,8,23], which is not possible here.

Activation under uni-axial dead load
In order to model the experiments of Keplinger et al. [15], we now consider a plate that is pre-stretched by a uni-axial dead load. A weight is applied in the x 1 -direction and charges on the lateral faces of the dielectric so that an electric field is induced in the x 2 -direction.
In dimensionless form, the loading curves relating the uni-axial stress s, the electric displacement component D 0 and the electric field E 0 to the stretches λ 1 and λ 3 for the neo-Hookean model (27) 1 are given by which lead to expressions for D 0 -λ 1 and E 0 -λ 1 relationships in terms of s (see, for example, [18,14] for details in the voltage-controlled case), namely Plots of D 0 and E 0 versus λ 1 based on (77) for several fixed values of s are shown in Figures 4(a) and 4(b), respectively, as the continuous curves. Notice, in particular, that D 0 is monotonic in λ 1 , while E 0 exhibits maxima, these behaviours being associated with loss of Hessian stability, as we elaborate on below.
It is a simple matter to extend the problem of minimizing the free energy ψ * associated with the whole system from the equi-biaxial to the uni-axial case in order to study material stability. First, we note that for the general biaxial case the second variations of the connection (44), corresponding to (47) and (49) in the equi-biaxial case, arē For the second variations of the free energy of the whole system ψ * to be positive, the 3 × 3 Hessian matrix H * must be positive definite (recall (34) 1 ). According to the equality above, this is equivalent under voltage control (when λ 1 , λ 3 and D 0 are free to vary and E 0 is fixed) toω 11 δλ 2 1 + 2ω 13 δλ 1 δλ 3 +ω 33 δλ 2 3 > 0, for non-zero δλ 1 and/or δλ 3 , i.e. it is equivalent to the leading 2 × 2 minor in H being positive definite.
On the other hand, under charge control (when λ 1 , λ 3 and E 0 are free to vary and D 0 is fixed), the left hand side of the equality (78) tells us that leading 2 × 2 minor of H * should be positive definite for stability, i.e.
for non-zero δλ 1 and/or δλ 3 . The latter inequality always holds for the neo-Hookean dielectric model, sinceω * 11 > 0 and the leading 2 × 2 minor of H * is positive definite, with determinant which factorizes in the equi-biaxial case, with λ 1 = λ 3 = λ, as the first factor coinciding with the corresponding result in the purely equi-biaxial case. Also, using (77) 2 , we find which corresponds toω 11 δλ 2 1 + 2ω 13 δλ 1 δλ 3 +ω 33 δλ 2 for fixed E 0 . This condition means that the leading 2 × 2 minor of H is indefinite, which can hold for fixed E 0 (at least for the neo-Hookean model), and we also have det H < 0.
Plots of E 0 versus λ 1 for the uni-axial case are shown in Figure 4(b) for s = 0, 1, 2, 3, and the connection between E 0 and λ 1 where det H * = 0 is also shown as the dashed curve that passes through the maximum points of E 0 . In Figure 4(a) are shown corresponding plots of D 0 (for s = 0, 1, 2, 3, 4) versus λ 1 , the dashed curve corresponding to where det H * = 0.
In conclusion, for a neo-Hookean dielectric subject to a uni-axial dead load, activation with voltage control can become unstable, but charge-controlled activation is always stable in the sense of the Hessian free energy criterion.
For the study of the geometric stability, we again refer to the limit cases. First, the thin-plate limit, again det N 3 = 0, reduces to The thick-plate limit, is a quadratic in D 2 0 given by Note that these two equations apply for all λ 1 (> 0) and, in particular, they recover the conditions for the equi-biaxial case (69) 1 and (75) when λ 1 = λ 3 = λ. The limit conditions above relate to wrinkles aligned with the direction of the uniaxial load. In Figure 5 we plot the corresponding D 0 -λ 1 curves by solving each condition (85) and (86) together with (76) 2 . The loading curves (77) 1 are also plotted, for different values of uni-axial pre-stress s. X + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + X X X X Figure  . The left-most dashed (blue) curve is the thick-plate limit curve (86) and the other dashed (black) curve is the thin-plate limit (85) curve, equivalent to the hypothetical no-weight curve (s = 0). The shaded region between the thick and thin-plate limits represents values of D 0 and λ 1 for which wrinkling could occur. Because the loading curves for the pre-stressed plate (s > 0) are all monotonic, they will not cross into the wrinkling region, provided the material is pre-stretched, and so wrinkling will not occur in the direction of the uni-axial load. The dots result from Finite Element computations, and follow the theoretical curves closely, although the clamping of the plate creates local, non-homogeneous fields. The main difference with the theoretical predictions is that the simulations eventually breakdown numerically, as indicated by red crosses.
As in the equi-biaxial case, the thin-plate limit is equivalent to the loading curve in the absence of pre-stress (s = 0). The wrinkling zone between the thin-and thick-plate limits is not reached by any of the curves corresponding to a pre-stretch (s > 0), and so the uni-axially pre-stretched plate will not wrinkle in the direction of the load. Note that in the absence of charge (D 0 = 0), we again recover the purely elastic case, where λ 2 = λ −1/2 1 , and the critical stretch for uni-axial surface instability is the Biot value λ 1 = 0.444 [2].
We also investigated wrinkles perpendicular to the direction of the load using the same method. There we looked for wrinkles in the (x 2 , x 3 )-plane, and constructed the Stroh formulation for the variables We then found that the thin-plate condition is identical to the loading curve equation (76) 2 for s = 0, and that the thick-plate limit is On solving this condition together with (76) 2 , no real solutions are found, and so there are no wrinkles perpendicular to the uni-axial load.
In the next section we see that the Hessian and geometric stabilities found from the homogeneous deformation fields can be contradicted by local inhomogeneous effects, as shown in numerical simulations.

Finite Element simulations
To complement the results of the theory, we developed electroelastic Finite Element (FE) models of the equi-biaxial and the uni-axial experiments using the commercial software COMSOL Multiphysics® [4], and coupled the elasticity and electrostatics in two different ways.
In the fully coupled model, COMSOL® uses the second Piola-Kirchhoff stress tensor, denoted P , and implements incompressibility via a volumetric energy function in the form κ(det F − 1) 2 /2, where κ is the initial bulk modulus, taken to be orders of magnitude larger than the shear modulus.
The second way to solve the coupled problem is by considering the effect of the Maxwell stress tensor as a fictitious mechanical boundary condition in the purely elastic problem. Since there are no charges within the volume of a dielectric, it is possible to consider the Maxwell stress as a pressure applied on the external faces of the volume. This adds a boundary traction τ m n to the mechanical problem, where n is the outward normal to the deformed surface of the specimen and τ m is the Maxwell stress tensor We found that both methods lead to the same results, although we noted that imposing the Maxwell stress tensor as a pressure boundary condition seemed to be a slightly more stable method in the uni-axial case.
In the equi-biaxial case, we found no difference between the predictions of the analytical model and those of the FE model, which also displayed stability and could be performed at any level of charge control, see Figure 3.
By contrast, a major difference between the analytical model and the FE model arises in the uni-axial case, because the simulations for the latter eventually break down. We identified the reason for this numerical breakdown to be due to the boundary conditions in the areas close to the clamping playing an initially small but eventually significant role. In the analytical model the strain is homogeneous and the material is free to deform in the transverse x 3 -direction. In the real-world experiments [15] and in the FE numerical model,  the top and bottom parts of the material are clamped and the strain is inhomogeneous in these neighbourhoods, see Figure 6(b). This behaviour is local, however, and the stretch in the direction of the uni-axial tension due to the weight (the x 1 -direction) is almost completely homogeneous, as can be seen in Figure 6(a).
When the stretched plate is electrically activated it expands in area and its thickness reduces. While it is free to expand in the direction of the dead load (the x 1 -direction), the situation in the transverse x 3 -direction is different. There is a central zone where the influence of the clamping is weak, so that the normal stress component P 33 in the x 3 -direction remains close to zero, as in the homogeneous case. On the other hand, the portions of material closer to the clamping areas suffer from the fixed displacement in the x 3 -direction imposed by the clamps. There the application of the uni-axial tension due to the weight increases P 33 , as is clearly visible in Figure 7.
When the plate is progressively activated with an increasing uniform charge distribution on its faces, the stress component P 33 near the clamping zone is progressively reduced until a critical value just below zero is reached: in this configuration, the plate undergoes a lateral compression that makes it buckle in the x 3 -direction [5]. At that point the FE computation breaks down, presumably because the stiffness matrix stops being positive definite and the solver has trouble converging. This phenomenon does not occur in the analytical model and in the equi-biaxial case, as λ 3 is homogeneous then and P 33 is imposed from the boundary condition and is identically equal to zero everywhere.
For larger weights, the levels of the P 33 stress component before activation are higher, making it possible to activate the dielectric plate with a larger value of the electric charge before the instability condition is reached, as can be seen from the dots in Figure 5.
Despite the significant difference in the transverse behaviour between the analytical and numerical model, due to the different boundary condition imposed, there is very good agreement in the results, as can be seen in Figure 5. As long as the FE model stays below the point of negative P 33 , the D 0 -λ curves follow those of the homogeneously deformed analytical model very closely.  Figure 6. The uppermost curve corresponds to the static dead-load condition (D 0 = 0); then an increasing charge activation is performed until the simulation reaches the instability point, where the stress is slightly compressive throughout (lowest curve) and the computation crashes. The colour coding for the levels of P 33 in the simulations goes from about 1 kPa in blue to about 0 kPa in dark orange.

Conclusion
In conclusion, we found that both equi-biaxial and uni-axial modes in the chargecontrol actuation of a dielectric plate are stable, whether the stability analysis is based on a Hessian criterion for the free energy of the whole system, or on the formation of small-amplitude inhomogeneous wrinkles.
By comparing the different Hessian criteria that result from the voltage-and chargecontrol situations, we found that charge-controlled actuation is always stable with respect to the Hessian criterion, in complete contrast to voltage-controlled actuation, which, according to the Hessian criterion, can become unstable.
We also investigated the possibility of small-amplitude wrinkles and found that the wrinkling conditions in the limiting cases of thin and thick plates occur only in compression, whereas it has been shown that wrinkles can exist in extension in the voltagecontrolled case [24]. As a result, charge-controlled actuation, which always occurs in extension, is also geometrically stable, again in contrast to voltage-control actuation.
To account for the difference between a theoretical homogeneous uni-axial deformation and the local inhomogeneous fields created by clamps in practice, we also conducted Finite Element simulations to verify our analytical results. We found complete agreement in the equi-biaxial case and very close agreement in the uni-axial case. So the assumption of homogeneous deformation is well justified for modelling the behaviour of charge-controlled activation of a dielectric plate in equi-biaxial stretch, and in uni-axial stretch when the aspect ratio of the specimen is high.
In the uni-axial case, Finite Element simulations reveal that the fringe effects are localised in a portion of area near the clamping zone and that they do not significantly affect the homogeneous loading curves of the system, although they have a strong effect on the eventual instability of the setup, a possibility that the homogeneous solution cannot capture. We may also argue that the emergence of compressive lateral stresses inside the plate seen in the simulations has a real-world counterpart, and that an equi-biaxial pre-stress leads to larger actuations than a uni-axial pre-stress in practice.
Of course, the plate may become unstable due to other causes than free energy instability or inhomogeneous small-amplitude wrinkles. Other mechanisms include for instance charge localisation [19] or thickness effects [9,30].