Bending control and stability 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 electroelasticity, when both the elastic shear modulus and the electric permittivity change with the thickness coordinate. The theory is illustrated for a neo-Hookean electroelastic energy function with the shear modulus and permittivity varying linearly across the thickness. In general the bending angle increases with the potential difference, and this enables the material inhomogeneity to be tuned to control the bending shape. We derive the Hessian criterion that ensures stability of the bent configurations in respect of a general form of electroelastic constitutive law specialized for the considered geometry. This requires that the Hessian remains positive. For the considered model we show that the bent configuration is stable until the voltage reaches the value for which the cross section of the bent configuration forms a complete circle.


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 [8,9], 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 realize shape control by tuning the applied voltage only, without any mechanical input [8,9].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 Selfpropagating High-temperature 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 [8].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 electroelasticity [14,15] to investigate the nonlinear bending behavior and the stability of such plates.
For an FGDE plate the required analysis, summarized 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 [16], Lu et al. [17], Zhao and Wang [18], Zurlo et al. [19], or Su et al. [20,21].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 the positive definiteness of the second variation of the free energy.In particular, we derive the associated Hessian criterion in respect of a general form of free energy for the considered geometry.On this basis we find that the considered bent configurations are stable for the full range of the applied voltage up to the point where the in-plane section of the plate forms a complete circle, the upper limit of the voltage being dependent on the values of the grading parameters.These results, which we summarize in Section 5, have potential for informing the design of high-performance actuators and sensors.

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 the reference configuration, and the faces X 1 = 0, A are coated with flexible electrodes as depicted in Fig. 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 Fig. 1(b), and the resulting deformation is assumed to have a plane strain character in the (X 1 , X 2 ) plane, as in [8,9].In Fig. 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 [22,23], 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 on the inner and outer bent faces, respectively.It follows that the dependence of λ on X 1 is given by 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 [23].Hence λ = ϕr/L.
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 electroelastic 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 ).
The relevant components of the total Cauchy stress tensor are the radial and circumferential components, denoted τ rr and τ θθ , respectively.These satisfy the equilibrium equation d dr (rτ rr ) = τ θθ . (2.6) The stress difference τ θθ − τ rr and Lagrangian electric field E L = λ −1 E are obtained from the formulas (2.7) see, for example, [14,15].
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 [23], where K is a constant to be determined from the boundary conditions, in a similar way to the situation in the purely elastic case [22].
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 strainenergy function and we recall that The boundary conditions (2.9) provide two expressions for K ,

11)
Fig. 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 ϕ.
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 Let φ(X 1 ) denote the electrostatic potential through the thickness.Then, from Maxwell's equation Curl E 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.4).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 Eq.(2.6), it follows that is automatically satisfied [23].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 ∫ r b ra rτ θ θ dr = 0, equivalently which, from (2.8) 2 , yields and for the model (2.10), where the subscript λ signifies the derivative with respect to λ.
Eqs. (2.12) and (2.13), respectively, give D L and K in terms of λ a and λ b .Eq. (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, Eq. (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 (2.25) To take this further we need to specify how µ and ε depend on X 1 , and through the connection (2.4), on λ.We assume a simple linear dependence of µ and ε on X 1 , as in [24,25], here specified by (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 ] . (2.27) Eq. (2.25) can also be evaluated to give (2.28) 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.5), 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 Eqs.(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 , where 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) voltage V = 0.2, Fig. 2  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 [10] and swelling functionally graded hydrogel plates [25].
For a fixed value of the electric parameter α ε = 0.4 and three values of the elastic parameter α µ , results analogous to those in Figs. 2 and 3 are shown in Figs. 4 and 5, respectively (now with V = 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 of V in Fig. 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 when V reaches the approximate value 0.4385.At this point the voltage, however, has not reached its maximum value.Fig. 6 shows plots of V 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 which V reaches a maximum.
As noted above, it is also seen from Figs. 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 values V = 0.2, 0.25 before the complete ring forms.

It is of interest to examine at what point the voltage does in-
deed reach its maximum, and for this purpose V is considered to   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 emphasize 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 then V = √ 1 − λ −4 , which is a monotone function of the stretch [21].Here our bending deformation (2.1) is also a plane strain deformation, but the V − λ 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  α ε .This is illustrated in Fig. 7 for the fixed value α µ = 0.3 and four negative values of α ε with V = 0.3.The bent configurations in this case are generally thinner and longer than those in Figs. 2  and 4.
Fig. 8 illustrates the corresponding radial and circumferential stress distributions, which are quite different from those for the situation with both α µ and α ε positive in that the signs of the stresses are the opposite of those in Figs. 3 and 5.
In Fig. 9, for fixed parameter values α ε = −0.3,α µ = 0.3, the deformation of a plate with increasing voltage is depicted.The plate thins significantly as it bends further, and no maximum of V occurs; the voltage increases monotonically with increasing values of λ a and λ b , which are both larger than 1.
A plot of V versus λ b is shown in Fig. 10.When V reaches the approximate value 0.756 the plate has deformed into a complete circle.The red bullet points identify the values 0.5, 0.6, 0.7, 0.756 of V on the curve corresponding to the bent sectors in Fig. 9.In Fig. 11, for fixed values α ε = −0.3,α µ = −0.3,by contrast, as the voltage increases the plate bends increasingly to the right until the voltage reaches the approximate value 0.5353, at which point a complete circular ring is formed.Note that the shapes in Figs. 2, 4, 7, 9, 11 are not drawn exactly to the correct scale, as the areas should all be the same.
Fig. 12 provides the corresponding plot of V versus λ b , with the red bullet points identifying the values 0.3, 0.4, 0.5, 0.5353 of V on the curve corresponding to the bent sectors in Fig. 11.In this case the voltage has a maximum, but it is beyond the value at which the complete ring is formed.

Stability analysis
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. [26], the second variation of the free energy density of a homogeneously deformed plate can be written in terms of either ω * or ω as ω * λλ (δλ) 2 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 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 ω λλ ω * , as shown in Dorfmann and Ogden [27] and Conroy Broderick et al. [26].
For the quadratic form in (4.1) to be strictly positive, H * must be positive definite, so that For the models with W given by (2.22), we have Thus, for H * to be positive definite, we require ω * λλ > 0, ω λλ > 0, (4.8) The first of these is always satisfied, while the second is µ which can be written in terms of the voltage on use of (2.27).
When α µ = α ε = 0 and the deformation is homogeneous, we have where DL = D L / √ µ a ε a , and we also note that DL = λ 2 ĒL , where 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 load D2 L = λ 4 −1 and Ē2 L = 1−λ −4 so that DL 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 voltage V = ĒL = 1.This is different from the equi-biaxial situation where Ē2 L exhibits a maximum and pull-in instability occurs; see, for example, [16,20,27].
By contrast with (4.11) the inhomogeneity of µ and ε must be accounted for in (4.10).Fig. 13 illustrates the results of calculation of ω λλ in terms of dimensionless ω = ω/µ a with ωλλ plotted against λ for the relevant range of values of λ a and λ b , which is different for each value of the dimensionless voltage V .The (black) right-most curve in each case corresponds to the value of the voltage at which a complete circular ring is formed.Up to this point in each case ω λλ > 0 for each λ ∈ [λ a , λ b ], but its value decreases with the voltage.Thus, each configuration is stable according to the Hessian criterion, and the second variation (4.2) is therefore positive.These results apply for the two considered pairs of values of (α ε , α µ ), but we have found similar results for a range of other combinations of these parameters.In view of the local result there is no need to consider the global counterpart (4.2) for the model in question.

Concluding remarks
In summary, we have investigated the voltage-induced bending response of an FGDE plate, and formulated a method for analyzing its stability based on the 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 α ε .A bent sector is subject to compressive circumferential stresses at both its inner and outer surfaces, with two neutral axes.
On analyzing the stability we found that the bent configurations are always stable according to the Hessian criterion up to the point that a complete circular ring is formed, the voltage at this point depending on the inhomogeneity parameters α ε and α µ .
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 may exhibit a maximum, it cannot be associated with pull-in instability within the stable domain since a complete ring is formed before the maximum is reached.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 of positive definiteness of the Hessian but we have found that such instability does not arise for the considered model.However, this approach to the analysis of instability tells only part of the story since it does not allow for wrinkling types of instability in membranes [19] or inhomogeneous wrinkling deformations [23,27], which can be induced [28] by the negative circumferential stress in the elastomer (Figs.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.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
.26) 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 depicts the resulting bent shape of the plate for a fixed value of the elastic parameter α µ = 0.3 and four values of electric parameter α ε .The plate is bent into a complete circle for α ε = 1.521.Again for α µ = 0.3 and three values of α ε , in Fig. 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 form r = (r − r a )/(r b − r a ).

Fig. 3 .Fig. 4 .
Fig. 3. Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radius r for the dimensionless voltage V = 0.2 and the following values of (α ε , α µ ): red (0.1, 0.3); purple (0.4, 0.3); blue dashed (0.7, 0.3).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5 .Fig. 6 .Fig. 7 .
figure legend, the reader is referred to the web version of this article.)be a function of λ b as the independent variable after elimination of K and D L .Then,

8 .
Plots of (a) the radial stress τ rr and (b) the circumferential stress τ θθ in scaled dimensionless form versus the normalized radius r for the dimensionless voltage V = 0.3 and the following values of (α ε , α µ ): red (−0.3, 0.3); purple (−0.4,0.3); blue dashed (−0.5, 0.3).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)has nontrivial solutions.