Bending control and instability of functionally graded dielectric elastomers

A rectangular plate of dielectric elastomer exhibiting gradients of material properties through its thickness will deform inhomogeneously when a potential difference is applied to compliant electrodes on its major surfaces, because each plane parallel to the major surfaces will expand or contract to a different extent. Here we study the voltage-induced bending response of a functionally graded dielectric plate on the basis of the nonlinear theory of electro-elasticity, when both the elastic shear modulus and the electric permittivity change with the thickness coordinate. The theory is illustrated for a neo-Hookean electro-elastic energy function with the shear modulus and permittivity varying linearly across the thickness. We find that in general the bending angle increases with the potential difference provided the Hessian remains positive, but instability can arise as the potential difference increases and the Hessian vanishes. We derive the Hessian criterion for predicting the instability, and show how the material gradients can be tuned to control the bending shape, and to delay or promote the onset of the instability of the elastomer.


Introduction
Dielectric elastomers are soft active materials capable of undergoing large deformations rapidly in response to an applied potential difference (voltage) across their thickness, and have therefore attracted considerable academic and industrial attention in recent years [1][2][3][4][5]. Shape control, which is eagerly pursued in dielectric elastomers, has considerable potential for applications to, for example, soft robots, energy harvest systems, actuators and sensors [6,7]. Hence, based on the concept of smart-bending actuation [9,14], an intelligent self-controlled switch can be designed to keep the working voltage in a desired range: as the voltage rises or falls, the material bends toward or away from a contact point so as to moderate the voltage automatically.
One clever strategy for designing voltage-responsive bending solids is to use dielectric elastomers with physical properties that vary in the thickness direction. When subject to a voltage through the thickness, a non-uniform deformation is generated in the material, which can result in a global bending. Dielectric-based multilayers composed of pre-stretched dielectric and inactive elastic layers have been proposed to realise shape control by tuning the applied voltage only, without any mechanical input [9,14]. However, discontinuous shear stresses may arise at the interface of the layers, which may result in slide, exfoliation or crack formation during the bending deformation [9][10][11].
To overcome this potential problem, this paper proposes to study functionally graded dielectric elastomers (FGDEs) with material properties varying continuously through the thickness. Experimentally, bulk functionally graded materials (FGMs) can be manufactured by methods such as Powder Metallurgy Technique, Centrifugal Casting and Solid Freeform Technology [12], etc. To produce thin FGMs, techniques such as Physical or Chemical Vapour Deposition (PVD/C-VD), Plasma Spraying and Self-propagating Hightemperature Synthesis (SHS) are more appropriate [13]. Spatial inhomogeneous dielectric properties can be introduced by embedding electroactive particles into a polymer FGM matrix unevenly through a well-established additive manufacturing (3D printing) process [14]. The electric field can be generated by applying a voltage through the flexible electrodes covered on the upper and lower faces of the resulting elastomers.
If the gradients of properties vary in a monotone manner though the plate thickness, we expect the expansion or contraction of each plane in the plate to also vary in the same way; hence, the plate will bend. Here we use the framework of nonlinear electro-elasticity [15,16] to investigate the nonlinear bending behaviour and the stability of such plates.
For an FGDE plate the required analysis, summarised in Section 2, is complicated by the inherent inhomogeneity, but we nonetheless manage to derive analytical formulas governing the voltage-induced bending of a plate made of a neo-Hookean dielectric with linear gradients in its material properties. In Section 3 we obtain numerical results showing the bent shape of the plate and the associated stress distributions. The results vary significantly with the selected values of the grading parameters.
For a plate with uniform properties, the onset of pull-in instability has been examined extensively, see for example the papers by Zhao and Suo [17], Lu et al. [18], Zhao and Wang [19], Zurlo et al. [20], or Su et al. [21,22]. In general, the pull-in instability leads to a dramatic thinning of the plate, thus triggering a giant deformation in the material. Here the corresponding phenomenon would be a sudden increase in the bending angle and, in general, a sudden thinning of the bent plate. However, we find (at least for the geometrical and material parameters chosen as examples) that this type of pull-in instability does not arise, because the maximal configuration of a closed circular ring is reached before the maximum in the voltage-stretch relation.
In Section 4, we examine the stability of the bent configuration on the basis of positive definiteness of the second variation of the free energy, for the considered geometry. In particular we derive the associated Hessian criterion. On this basis we find that the considered bent configurations are either unstable, irrespective of the applied (non-zero) voltage, or stable for a limited range of the applied voltage, depending on the values of the grading parameters. These results, which we summarise in Section 5, have potential for informing the design of high-performance actuators and sensors.
2 Bending response of an FGDE plate 2

.1 Kinematics and electric field
Consider an FGDE plate of initial length L, thickness A and width H, the latter being assumed to be much longer than its thickness and length. The solid occupies the region 0 ≤ X 1 ≤ A, −L/2 ≤ X 2 ≤ L/2, 0 ≤ X 3 ≤ H in the reference configuration, and the faces X 1 = 0, A are coated with flexible electrodes as depicted in Figure 1(a). The material properties are taken to be inhomogeneous with grading dependent on the thickness coordinate X 1 .
On application of a potential difference V (voltage) across the electrodes the plate is bent into the shape shown in Figure 1(b), and the resulting deformation is assumed to have a plane strain character in the (X 1 , X 2 ) plane, as in [9,14]. In Figure 1(a) we have indicated the dependence on X 1 of the shear modulus µ and permittivity ε, which will be specified in Section 2.2. The material is assumed to be incompressible.
The deformation of the plate is described by the equations as in [23,24], where (r, θ, z) are the cylindrical polar coordinates in the deformed configuration, and r a , r b are the radii of the inner and outer bent surfaces, respectively. The associated deformation gradient, denoted F, has the diagonal form diag(λ −1 , λ, 1) with respect to the cylindrical polar axes, where is the circumferential stretch. It takes the values 3) Figure 1: A block of FGDE with continuously varying shear modulus µ(X 1 ) and permittivity ε(X 1 ) along the thickness, coated with two compliant electrodes. The layer bends when loaded with voltage V through the thickness, with the bending angle ϕ.
where X 1 = X 1 /A is a non-dimensional measure of position through the thickness, µ A and ε A are the shear modulus and permittivity of the solid at the face X 1 = 0, and 0 ≤ β µ < 1, β ε ≤ 1 are two dimensionless parameters characterizing the functionally graded properties of the shear modulus and permittivity of the elastomer, respectively. Subject to a voltage V , the elastomer bends, as depicted in Figure 1(b), according to the deformation [22,23] where (r, θ, z) are the cylindrical coordinates in the current configuration, and r a , r b are the radii of the inner and outer bent surfaces, respectively. The associated deformation gradient is is the circumferential stretch. It takes the values 3 Figure 1: A block of FGDE with continuously varying shear modulus µ(X 1 ) and permittivity ε(X 1 ) through the thickness, coated with two compliant electrodes. The layer bends when loaded with a voltage V across the electrodes, with the bending angle ϕ.
on the inner and outer bent faces, respectively. Note that, by taking X 2 = L/2 and θ = ϕ/2 in (2.1) 2 , we obtain an expression for the bending angle ϕ, namely as given in [24]. Hence λ = ϕr/L. We shall need X 1 in terms of λ later. On use of (2.1) 1 , (2.2) and (2.3) this is written as On application of the voltage V a radial electric field component, denoted E, is generated in the bent configuration of the material, assuming that edge effects can be neglected, and the corresponding electric displacement field component is denoted D. Each of E and D depends on r and is independent of θ and z. Maxwell's equation div D = 0 then reduces to d(rD)/dr = 0, so that rD is constant. The corresponding Lagrangian field, given by D L = F −1 D in general for an incompressible material, reduces to the single component D L = λD, which by (2.2) is therefore a constant.

Constitutive law
For an isotropic electro-elastic material in general we take the energy to be a function of F and D L , but for the considered geometry, deformation and electric field, this reduces to dependence on λ and D L . We denote the energy function by ω * (λ, D L ).

4
The relevant components of the total Cauchy stress tensor are the radial and circumferential components, denoted τ rr and τ θθ , respectively. These satisfy the equilibrium equation The stress difference τ θθ − τ rr and Lagrangian electric field E L = λ −1 E are obtained from the formulas see, for example, [15,16]. Since, from (2.2), λ is proportional to r, while D L is constant, the combination of (2.6) and (2.7) 1 followed by integration, leads to as derived in [24], where K is a constant to be determined from the boundary conditions, in a similar way to the situation in the purely elastic case [23].
Here we assume that the inner and outer surfaces of the bent sector at r a and r b are free of mechanical traction, so that there being no Maxwell stress on these surfaces. For definiteness, we consider an energy function which has the form where W is derived from any isotropic purely elastic strain-energy function and we recall that D = λ −1 D L . The boundary conditions (2.9) provide two expressions for K, and hence an expression for D 2 L , namely where ε a and ε b are the values of ε at r = r a and r = r b , respectively, i.e. for X 1 = 0, A. The form of the function ε(X 1 ) will be exemplified below. It follows that (2.13)

5
Let φ(X 1 ) denote the electrostatic potential through the thickness. Then, from Maxwell's equation CurlE L = 0 we have E L = −Gradφ, which specializes here to E L = −dφ/dX 1 . The potential difference φ(0) − φ(X 1 ) is the voltage V , which, on use of (2.7) 2 is given by the latter change of variable making use of (2.5). For the model (2.10), this yields and we recall that ε depends on X 1 and hence on λ.
To complete the formulation of the problem, we first note that the resultant force on the lateral end faces at θ = ±ϕ/2 vanishes since, with the help of the boundary conditions (2.9) and equation (2.6), it follows that is automatically satisfied [24]. We then assume that there is no resultant moment on these faces so that the block is bent by the application of a voltage alone. This requires λ λ ∂ω * ∂λ + ω * + K dλ = 0, (2.18) and for the model (2.10), where the subscript λ signifies the derivative with respect to λ. Equations (2.12) and (2.13), respectively, give D L and K in terms of λ a and λ b . Equation (2.15) then determines V also in terms of λ a and λ b , which we write in the form (2.20) Once the form of W is prescribed, equation (2.19), after substituting for D L and K, yields an implicit connection between λ a and λ b , which we write as Then, in principle, V can be determined in terms of λ a (or λ b ). For further specialization we now take W to be the neo-Hookean strain-energy function, which, for the considered plane deformation, has the form where the shear modulus µ is functionally graded through the thickness, i.e. is a function of X 1 . Then the formulas (2.12) and (2.13) specialize to where µ a and µ b are the values of µ(X 1 ) at X 1 = 0, A (r = r a , r b ), and (2.19) becomes To take this further we need to specify how µ and ε depend on X 1 , and through the connection (2.5), on λ. We assume a simple linear dependence of µ and ε on X 1 , as in [25,26], here specified by with µ b = µ a (1 − α µ ) and ε b = ε a (1 + α ε ), where α µ and α ε are two dimensionless constants characterizing the functionally graded properties of the shear modulus and permittivity, respectively, of the material. The expressions for D 2 L and K in (2.23) and (2.24) are then modified accordingly. The parameters α µ and α ε are subject to the restrictions α µ < 1 and α ε > −1 since both µ and ε are positive. Additionally, since the permittivity of all materials is greater than the vacuum permittivity ε 0 , we must have ε a (1 + α ε ) > ε 0 for each α ε > −1. Note that negative values of α µ and α ε reverse the direction of bending, as is illustrated in Section 4.
The integral in (2.15) can now be evaluated to give Equation (2.25) can also be evaluated to give For a given applied voltage V , after D L is substituted into (2.27) and D 2 L and K into (2.28) from (2.23) and (2.24), the resulting equations yield specific forms of the functions f and g in (2.20) and (2.21), which can now in principle be solved simultaneously for λ a and λ b . The dependence of the bending angle ϕ, given by (2.4), on the applied voltage can then be determined. The corresponding radial and circumferential stress components τ rr and τ θθ can also be determined, from (2.8) with the specializations (2.10) and (2.22). Numerical results for specific values of α µ and α ε are provided in the following section.

Numerical results
For numerical purposes we work in terms of the dimensionless quantities: and equations (2.20) and (2.21) take the dimensionless forms These are the dimensionless forms of (2.27) and (2.28) after substitution for D L and K, Thus,f andḡ depend only on λ a , λ b and the parameters α µ and α ε , i.e. they are independent of µ a and ε a . For an initially rectangular FGDE plate with aspect ratio L/A = 10, subject to a (dimensionless) voltageV = 0.2, Figure 2 depicts the resulting bent shape of the plate for a fixed value of the elastic parameter α µ = 0.3 and three values of electric parameter  α ε . Again for α µ = 0.3 and three values of α ε , in Figure 3 the distributions of the radial and circumferential stress components through the plate are shown in dimensionless forms τ rr = τ rr /µ a ,τ θθ = τ θθ /µ a versus the radial coordinate in the formr = (r − r a )/(r b − r a ).  Figures 10 and 9, respectively. In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε .
As is clear from Figure 10, as the magnitude of the voltage increases the bent sector forms a complete circular sector. At this point the voltage, however, does not reach its maximum value that is possible for the two cases illustrated in Figure 18. The crosses marked in Figure 18 are not accessible. It was also seen from Figures 3 and 9 that the As expected intuitively, we find that increasing α ε makes the plate more susceptible to bending. The stress components vary continuously through the thickness, with the circumferential stress two orders of magnitude larger than that of the radial stress. Interestingly, at least for examples considered here, the bent sector has two neutral surfaces (where τ θθ = 0). Additionally we note that while on both its inner and outer faces the circumferential stress is compressive, these faces are nonetheless in extension as the circumferential stretches are greater than 1 there. Similar phenomena were observed for the bending of purely elastic layered plates [11] and swelling functionally graded hydrogel plates [26].
For a fixed value of the electric parameter α ε = 0.4 and three values of the elastic parameter α µ , results analogous to those in Figures 2 and 3 are shown in Figures 4 and 5, respectively (now withV = 0.25). In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε .
As is clear from the larger value ofV in Figure 4, as the magnitude of the voltage increases the circular sector becomes more and more bent. As the voltage increases further the sector will eventually form a complete circular ring, which, for the values α ε = 0.4, α µ = 0.4, occurs whenV reaches the approximate value 0.4385. At this point the voltage, however, has not reached its maximum value. Figure 6 shows plots ofV versus λ b for different values of α ε and α µ , with the points at which a complete circular ring is formed identified by circles, while the crosses mark the (non-accessible) points at whichV reaches a maximum.
As noted above, it is also seen from Figures 3 and 5 that the inner and outer faces of the bent sector are subject to compressive stress (τ θθ < 0) during the bending deformation for the lower valuesV = 0.2, 0.25 before the complete ring forms.
It is of interest to examine at what point the voltage does indeed reach its maximum, For a fixed value of α ε = 0.4 and three values of α µ , results analogous to those in Figures 2 and 3 are shown in Figures 10 and 9, respectively. In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε .
As is clear from Figure 10, as the magnitude of the voltage increases the bent sector forms a complete circular sector. At this point the voltage, however, does not reach its maximum value that is possible for the two cases illustrated in Figure 18. The crosses marked in Figure 18 are not accessible. It was also seen from   Figures 10 and 9, respectively. In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε .
As is clear from Figure 10, as the magnitude of the voltage increases the bent sector forms a complete circular sector. At this point the voltage, however, does not reach its maximum value that is possible for the two cases illustrated in Figure 18. The crosses marked in Figure 18 are not accessible. It was also seen from Figures 3 and 9 that the inner and outer faces of the bent sector are compressed (τ θθ < 0) during the bending deformation for the lower valueV = 0.3 before the complete sector forms. 9 Figure 3: Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radiusr for the dimensionless voltageV = 0.2 and the following values of (α ε , α µ ): For a fixed value of α ε = 0.4 and three values of α µ , results analogous to those in Figures 2 and 3 are shown in Figures 4 and 5, respectively. In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε . Figure 5: Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radiusr for the dimensionless voltageV = 0.2 and the following values of (α ε , α µ ): red (0.1, 0.3); purple (0.4, 0.3); blue dashed (0.7, 0.3).
For a fixed value of α ε = 0.4 and three values of α µ , results analogous to those in Figures 2 and 3 are shown in Figures 10 and 9, respectively. In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε .
As is clear from Figure 10, as the magnitude of the voltage increases the bent sector forms a complete circular sector. At this point the voltage, however, does not reach its maximum value that is possible for the two cases illustrated in Figure 18. The crosses marked in Figure 18 are not accessible. It was also seen from Figures 3 and 9 that the 9 -0.06 -0.8 Figure 5: Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radiusr for the dimensionless voltageV = 0.2 and the following values of (α ε , α µ ): red (0.1, 0.3); purple (0.4, 0.3); blue dashed (0.7, 0.3).
For a fixed value of α ε = 0.4 and three values of α µ , results analogous to those in Figures 2 and 3 are shown in Figures 10 and 9, respectively. In general, whether bending is enhanced or retarded by a change in α µ depends on the value of α ε .
As is clear from Figure 10, as the magnitude of the voltage increases the bent sector forms a complete circular sector. At this point the voltage, however, does not reach its maximum value that is possible for the two cases illustrated in Figure 18. The crosses marked in Figure 18 are not accessible. It was also seen from Figures 3 and 9 that the inner and outer faces of the bent sector are compressed (τ θθ < 0) during the bending deformation for the lower valueV = 0.3 before the complete sector forms. As is clear from Figure 6, as the magnitude of the voltage increases the bent sector forms a complete circular sector. At this point the voltage, however, does not reach its maximum value that is possible for the two cases illustrated in Figure 14. The crosses marked in Figure 14 are not accessible. It was also seen from Figures 2 and ?? that the inner and outer faces of the bent sector are compressed (τ θθ < 0) during the bending deformation for the lower valueV = 0.3 before the complete sector forms. and for this purpose V is considered to be a function of λ b as the independent variable after elimination of K and D L . Then,

4)
-0.6 Figure 5: Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radiusr for the dimensionless voltageV = 0 .3 and the following values of (α ε , α µ ): red (0 .4, 0.2); purple (0 .4, 0.4); blue dashed (0.4, 0.6). Since, however, the maximum of V does not occur it is necessary to examine when instability might arise, which, for example, might occur for given values of the parameters α µ and α ε as the voltage increases. This is now examined based on an energy approach. has nontrivial solutions. However, because the maximum of V does not occur before the plate is bent into a full circle, it is necessary to examine when instability might arise, which, for example, might occur for given values of the parameters α µ and α ε as the voltage increases. This is examined in the next section, based on an energy approach.
As the maximum of V does not occur in the relevant range of deformations we emphasise that, in contrast to what happens when a homogeneous dielectric plate, equi-biaxially deformed under an applied voltage, the phenomenon of pull-in instability does not appear for a functionally graded plate under plane strain. We recall that in the homogeneous plane strain expansion of a neo-Hookean dielectric plate, there is also no pull-in instability, because thenV = √ 1 − λ −4 , which is a monotone function of the stretch [22]. Here our bending deformation (2.1) is also a plane strain deformation, but theV − λ b relationship is not monotone.
Because of the dependence of µ and ε on X 1 there is a lack of symmetry with respect to X 1 = 0 and X 1 = A. Thus, bearing in mind the limitations on the values of α µ and α ε mentioned in Section 2.2, it is possible to consider negative values of α µ and/or α ε , which leads to a bending response opposite to that shown in Figures 2 and 4. This is illustrated in Figure 7 for the fixed value α µ = 0.3 and three negative values of α ε withV = 0.3.

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, D L ) = ω(λ, E L ) + D L E L . From the analysis of Conroy Broderick et al. [28], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L . For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart of this, which requires integration over X 1 ∈ [0, A] and X 2 ∈ [−L/2, L/2], is, in terms of ω,  Figure 8 illustrates the corresponding radial and circumferential stress distributions, which are quite different from those for the situation with both α µ and α ε positive. In particular, τ rr is negative while there is only one neutral surface τ θθ = 0. Correspondingly, the circumferential stressτ θθ is monotonic, but it has both positive and negative values and the integral in Eq. (2.17) vanishes. The situation is very similar for a fixed negative value of α ε with different values of α µ or for a fixed pair of values of α ε and α µ with different voltages.

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, D L ) = ω(λ, E L ) + D L E L . From the analysis of Conroy Broderick et al. [28], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L . For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart of this, which requires integration over X 1 ∈ [0, A] and X 2 ∈ [−L/2, L/2], is, in terms of ω,

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, D L ) = ω(λ, E L ) + D L E L . From the analysis of Conroy Broderick et al. [28], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as In Figure 9, for fixed parameter values α ε = −0.2, α µ = 0.3, the deformation of a plate with increasing voltage is depicted. Initially the plate bends to the right and then thins significantly as it bends further, but then as the voltage increases more the bending is reversed and further thinning occurs. No maximum ofV occurs; the voltage increases monotonically with increasing values of λ a and λ b , which are both larger than 1. 12

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, D L ) = ω(λ, E L ) + D L E L .
From the analysis of Conroy Broderick et al. [28], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L . For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart of this, which requires integration over X 1 ∈ [0, A] and X 2 ∈ [−L/2, L/2], is, in terms of ω, Figure 11: Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radiusr for the dimensionless voltageV = 0.3 and the following values of (α ε , α µ ): red (−0.

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, From the analysis of Conroy Broderick et al. [28], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L . For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart of this, which requires integration over X 1 ∈ [0, A] and 12 Figure 8: Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radiusr for the dimensionless voltageV = 0.3 and the following values of (α ε , α µ ): red (−0. is reversed and further thinning occurs. No maximum ofV occurs; the voltage increases monotonically with increasing values of λ a and λ b , which are both larger than 1.

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, In Figure 10, for fixed parameter values α ε = −0.6, α µ = 0.3, by contrast, the plate bends increasingly to the right without significant thinning until the voltage reaches the approximate value 0.639, at which point a complete circular ring is formed, while the voltage increases monotonically with increasing values of λ a and λ b , which are both larger than 1.

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, From the analysis of Conroy Broderick et al. [25], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, D L ) = ω(λ, E L ) + D L E L . From the analysis of Conroy Broderick et al. [25], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L . For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart of this, which requires integration over X 1 ∈ [0, A] and X 2 ∈ [−L/2, L/2], is, in terms of ω, We focus initially on the local form of the second variation (4.1). The Hessian matrix associated with the left-hand side of (4.1), denoted H * , is given by which can also be written ω λλ ω * D L D L , as shown in Dorfmann and Ogden [24] and Conroy Broderick et al. [25]. monotonically with increasing values of λ a and λ b , which are both larger than 1. Note that the shapes in Figures 2,4,7,9,11 are not drawn exactly to the correct scale, as the areas should all be the same. Figure 12 provides the corresponding plot ofV versus λ b , with the red bullet points identifying the values 0.1, 0.3, 0.45, 0.6, 0.639 ofV on the curve.

Onset of instability
Consider the energy functions ω * (λ, D L ) given by (2.10) with the specialization (2.20) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, D L ) = ω(λ, E L ) + D L E L .
From the analysis of Conroy Broderick et al. [25], the second variation of the free energy of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L , the corresponding subscripts representing partial derivatives. For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart of this, which requires integration over

Onset of instability
Consider the energy function ω * (λ, D L ) given by (2.10) with the specialization (2.22) for the considered plane strain situation. We denote by ω(λ, E L ) the corresponding function with the independent variable E L , related to ω * by ω * (λ, From the analysis of Conroy Broderick et al. [28], the second variation of the free energy density of a homogeneously deformed plate can be written in terms of either ω * or ω as where δλ, δD L and δE L , respectively, are variations in λ, D L and E L , the corresponding subscripts representing partial derivatives. For the inhomogeneous deformation considered here, this is also the local form of the second variation, and its global counterpart, which requires integration over X 1 ∈ [0, A] and X 2 ∈ [−L/2, L/2], is, in terms of ω, We focus initially on the local form of the second variation (4.1). The Hessian matrix associated with the left-hand side of (4.1), denoted H * , is given by which can also be written ω λλ ω * D L D L , as shown in Dorfmann and Ogden [27] and Conroy Broderick et al. [28].
For the quadratic form in (4.1) to be strictly positive, H * must be positive definite, so that ω * λλ > 0, det H * ≡ ω λλ ω * D L D L > 0. (4.5) For the models with W given by (2.22), we have Thus, for H * to be positive definite, we require ω * λλ > 0, ω λλ > 0, (4.8) and, in view of (4.7) 2 , the global form of the second variation (4.2) is positive if ω λλ > 0 for each point of the domain of the integral.
On use of the connection E L = ε −1 λ −2 D L , the expressions (2.26) and the formula (2.27), we find that ω λλ is given by and similarly, the expression for ω * λλ is given by (4.9) When α µ = α ε = 0 and the deformation is homogeneous, we have whereD L = D L / √ µ a ε a , and we also note thatD L = λ 2Ē L , whereĒ L = E L ε a /µ a . Then, it is clear that ω * λλ > 0, while det H * can become negative as D L increases from zero, in which case stability according to the Hessian criterion is lost. Since there is no applied loadD 2 L = λ 4 − 1 andĒ 2 L = 1 − λ −4 so thatD L increases indefinitely with λ whileĒ L increases up to an upper limit 1, i.e. a homogeneous plate in plane strain can only support a maximum non-dimensional voltageV =Ē L = 1. This is different from the equi-biaxial situation whereĒ 2 L exhibits a maximum and pull-in instability occurs; see, for example, [17,21,27].
Turning back to the inhomogeneous problem. In considering the inhomogeneity, for stability the integral expression in (4.2) must be positive for all possible choices of δλ and δE L with at least one of them non-zero [29]. Since the second term in the integrand is non-negative and δE L can be chosen to be zero, stability requires the integral to be positive for all non-zero δλ. If, however, ω λλ < 0 for part of the range of λ then by choosing δλ to be only non-zero for those λ where ω λλ < 0, the second variation would be negative, and we would have instability for such values ofV . Thus, stability is guaranteed when ω λλ > 0 for every value of λ ∈ [λ a , λ b ].
For the examples in Figures 2 and 4, with positive values of α µ and α ε , we found (not shown here) that ω λλ is negative for every value of λ ∈ [λ a , λ b ] whenV exceeds 0.1, and for smaller values ofV it is only positive for a very small range of values of λ near λ a . Thus, such configurations are unstable according to both the local Hessian criterion and its global counterpart. We also found this instability for other positive values of α µ and α ε .
These interesting and quite surprising results therefore led us to consider if the situation was different for negative values of α µ and α ε . It turns out, depending on the actual values of α µ and α ε , that the related configurations are stable for a range of values ofV .
In Figures 13(a) and (b), respectively, ω λλ and ω * λλ are illustrated for a negative electric parameter α ε = −0.2 and positive elastic parameter α µ = 0.3, with several different values ofV (refer to the bent shape in Figure 7(b) forV = 0.3, for example). The range of values of λ a and λ b changes as the deformation increases with the voltage. Both ω λλ and ω * λλ remain positive asV increases except that for a very small range of values ofV near 0.85 ω λλ is negative with maximum magnitude about 0.1. In this case, the plate is unstable since, as has been shown in the discussion related to Eq. (4.11), stability requires that ω λλ be positive throughout the plate. With reference to the changing curvature shown in Figure 9 we mention that asV reaches a value between 0.903 and 0.904 the plate becomes straight and the curvature reverses for larger values ofV , but then the configuration becomes unstable (both ω λλ and ω * λλ become negative). situation whereĒ 2 L exhibits a maximum and snap-through and pull-in instability occur; see, for example, [14,18,24].
Turning back to the inhomogeneous problem, we focus on the local Hessian. For the examples in Figures 2 and 4, with positive values of α µ and α ε it has been found that ω λλ is negative as soon as a voltage is applied for each value of λ ∈ [λ a , λ b ], while ω * λλ is positive for the lower values of λ ∈ [λ a , λ b ] but negative for the upper range of values of λ. Thus, such configurations are unstable according to the Hessian criterion. This has also been found for other positive values of α µ and α ε .
We therefore investigate negative values of α µ and α ε , and it turns out that the related configurations, depending on the actual values of α µ and α ε are stable for a range of values ofV . In Figures 13(a) and (b), respectively, ω λλ and ω * λλ are illustrated for α ε = −0.2 with α µ = 0.3 with several different values ofV (refer to the bent shape in Figure 7(b) for V = 0.3). The range of values of λ a and λ b changes as the deformation increases with the voltage. Both ω λλ and ω * λλ remain positive asV increases except for a very small range of values ofV near 0.85 when the negative value has maximum magnitude about 0.1. With reference to the changing curvature shown in Figure 9 we mention that asV reaches a value between 0.903 and 0.904 the plate becomes straight and the curvature reverses for larger values ofV , but then the configuration becomes unstable (both ω λλ and ω * λλ become negative). ω λλ /µ a ω * λλ /µ a λ λ Figure 13: Plots of (a) ω λλ /µ a and (b) ω * λλ /µ a versus λ for the following values ofV reading from left to right: 0.1, 0.3, 0.5, 0.65, 0.75, 0.8, 0.85, 0.88, 0.895, 0.9, with α ε = −0.2 and α µ = 0.3. Note that range of values of λ ∈ [λ a , λ b ] is different in each case.
In Figures 14(a) and (b), respectively, ω λλ and ω * λλ are illustrated for α ε = −0.6 with α µ = 0.3 with several different values ofV (refer to the bent shape in Figure 11(b) for V = 0.3). ForV = 0.1, 0.2, 0.25, ω λλ is positive and hence, by (4.7), the integral in (4.2) is positive, so the second variation of the free energy function is positive and the configuration is stable. This is also the case for values ofV up to nearly 0.3, but ω λλ becomes negative for a small range of values of λ whenV = 0.3 and for a larger range with increasingV , as illustrated by the blue and purple curves in Figure 14 In Figures 14(a) and (b), respectively, ω λλ and ω * λλ are illustrated for α ε = −0.6 with α µ = 0.3 with several different values ofV (refer to the bent shape in Figure 11(b) for V = 0.3). ForV = 0.1, 0.2, 0.25, ω λλ is positive and hence, by (4.7), the integral in (4.2) is positive, so the second variation of the free energy function is positive and the configuration is stable. This is also the case for values ofV up to nearly 0.3, but ω λλ becomes negative for a small range of values of λ whenV = 0.3 and for a larger range with increasingV , as illustrated by the blue and purple curves in Figure 14(a). On the other hand, ω * λλ is positive for all λ andV .
We emphasise that when ω λλ > 0 for λ ∈ [λ a , λ b ] the configuration in question is stable, and because then ω * λλ > ω λλ there is no need to consider the sign of ω * λλ in assessing We emphasize that when ω λλ > 0 for λ ∈ [λ a , λ b ] the configuration in question is stable, and since then ω λλ < ω * λλ there is no need to consider ω * λλ in assessing stability. More generally, in considering the inhomogeneity, for stability the integral expression in (4.2) must be positive for all possible choices of δλ and δE L with at least one of them non-zero. Since the second term in the integrand is non-negative and δE L can be chosen to be zero, stability requires the integral λ b λ a λω λλ (δλ) 2 dλ (4.12) to be positive for all non-zero δλ. If, however, ω λλ < 0 for part of the range of λ then by choosing δλ to be only non-zero for those λ where ω λλ < 0 the second variation would be negative, and we have instability for such values ofV . Thus, there is only stability when ω λλ > 0 for every value of λ ∈ [λ a , λ b ].
I'm still not 100% convinced by this, but I haven't yet found any contradiction in the variational literature.

Concluding remarks
In summary, we have investigated the voltage-induced bending response of an FGDE plate, and formulated a method for analyzing the onset of its instability based on the loss of positive definiteness of the Hessian. The 'bendability' of the plate can be tuned by carefully selecting the grading properties, as illustrated here for the linear grading, by choosing appropriate values of the inhomogeneity parameters α µ and α ε . Depending on the signs of these parameters, a bent sector may be either (a) subject to compressive circumferential stresses at both its inner and outer surfaces, with two neutral axes, when both parameters are positive, or (b) compressive circumferential stress on one surface and tensile stress on

Concluding remarks
In summary, we have investigated the voltage-induced bending response of an FGDE plate, and formulated a method for analyzing the onset of its instability based on the loss of positive definiteness of the Hessian.
We found that the 'bendability' of the plate can be tuned by carefully selecting the grading properties, as illustrated here for linear gradients, by choosing appropriate values of the inhomogeneity gradients α µ and α ε . Depending on the signs of these parameters, a bent sector may be either (a) subject to compressive circumferential stresses at both its inner and outer surfaces, with two neutral axes, when both parameters are positive, or (b) subject to compressive circumferential stress on one surface and tensile stress on the other, with a single neutral axis, when one of the parameters is negative. On analysing the stability we found that configurations associated with (a) are always unstable, while configurations based on (b) are stable for certain combinations of α µ and α ε or, for other combinations, stable provided the voltage is not too large.
In contrast with the situation for a homogeneous plate subject to equi-biaxial deformations, where pull-in instability occurs, in the present plane strain problem there in no such pull-in instability. In particular, even though the voltage-stretch curve exhibits a maximum, it cannot be associated with pull-in instability within the stable domain.
The analysis herein is based on a particular (neo-Hookean based) constitutive model and a simple form of the material grading properties. It can be expected that the results will be somewhat different, at least quantitatively, for different model choices, but the framework presented here can accommodate more general models.
In this paper we have focused only on the possibility of instability arising from the loss 18 of positive definiteness of the Hessian. This approach to analysis of stability tells only part of the story of stability since it does not allow for wrinkling types of instability in membranes [20] or inhomogeneous wrinkling deformations [24,27], which can be induced [30] by the negative circumferential stress in the elastomer (Figures 3, 5, 8).
Applications where gradient properties could be beneficial are in the design of sensors and actuators and complicated motions of soft robots, for example.