Bifurcation of finitely deformed thick-walled electroelastic cylindrical tubes subject to a radial electric field

This paper is concerned with the bifurcation analysis of a pressurized electroelastic circular cylindrical tube with closed ends and compliant electrodes on its curved boundaries. The theory of small incremental electroelastic deformations superimposed on a finitely deformed electroelastic tube is used to determine those underlying configurations for which the superimposed deformations do not maintain the perfect cylindrical shape of the tube. First, prismatic bifurcations are examined and solutions are obtained which show that for a neo-Hookean electroelastic material prismatic modes of bifurcation become possible under inflation. This result contrasts with that for the purely elastic case for which prismatic bifurcation modes were found only for an externally pressurized tube. Second, axisymmetric bifurcations are analyzed, and results for both neo-Hookean and Mooney–Rivlin electroelastic energy functions are obtained. The solutions show that in the presence of a moderate electric field the electroelastic tube becomes more susceptible to bifurcation, i.e., for fixed values of the axial stretch axisymmetric bifurcations become possible at lower values of the circumferential stretches than in the corresponding problems in the absence of an electric field. As the magnitude of the electric field increases, however, the possibility of bifurcation under internal pressure becomes restricted to a limited range of values of the axial stretch and is phased out completely for sufficiently large electric fields. Then, axisymmetric bifurcation is only possible under external pressure.


Introduction
Electroactive elastomeric materials have attracted much interest in the literature in recent years because their properties are suitable for many potential applications in engineering science, such as in the production of actuators, sensors and other devices. One of the geometries of relevance in this context is that of a circular cylindrical tube actuator made of a dielectric elastomer with flexible electrodes on its curved boundaries (see, for example, Pelrine et al. [1]). Aspects of the associated theoretical analysis of this problem have been examined in [2] and in our recent paper [3], where the analysis was based on the tube maintaining its circular cylindrical shape. An important consideration for the effective operation of such an actuator is its ability to maintain its perfect cylindrical shape. In this paper we predict the onset of bifurcation from the cylindrical shape on the basis of the nonlinear theory of electroelasticity for specific material models.
It is well known that a tube made of a rubberlike material inflated by internal pressure may develop a bulge beyond a critical deformation, giving rise to deviations from the perfect cylindrical shape under symmetric loads. Relevant experimental data can be found in [4] and [5], for example, while for a dielectric elastomer tube recent experiments are described in [6]. For thin-and thick-walled cylindrical shells of elastic material under internal and external pressure, extensive bifurcation analyses have been provided in the papers [7,8] by Haughton and Ogden on the basis of the theory of small deformations superimposed on a finite deformation. Corresponding analysis for an electroactive tube, however, is limited. For example, 60 Page 2 of 27 A. Melnikov and R. W. Ogden ZAMP for an electroactive material the bifurcation of a tube subject to just an axial mechanical loading and an axial electric field was examined recently by Su et al. [9], who found that the stability of the tube depends on the electric field, the material constants and the geometrical dimensions of the tube. Prismatic bifurcations have been examined recently by Bortot and Shmuel [10] for a tube with fixed axial extension but no internal pressure. In order to determine when deviations from the perfect cylindrical configuration are possible we use the theory of incremental deformations and electric fields superimposed on an underlying finite deformation and electric field. We note that by including dependence on time the problem of wave propagation in an electroelastic tube can be analyzed, as has been done by Shmuel and deBotton [11], who considered axisymmetric waves for an axially loaded tube with a radial electric field, but no radial mechanical loads. When the wave speed vanishes, corresponding bifurcation results can be obtained.
This paper is structured as follows. In Sect. 2 we summarize the most important ingredients of the general theory of electroelasticity as formulated by Dorfmann and Ogden [12]. The incremental electroelastic constitutive laws are also provided in this section, along with explicit formulas for the electroelastic moduli tensors that are used in subsequent sections.
In Sect. 3, following [3], we consider the underlying problem of the inflation and extension of an electroelastic tube in the presence of a radial electric field as developed in [3]. The electric field is generated by a potential difference between compliant electrodes on the inner and outer curved surfaces of the tube. Expressions are given for the internal pressure in a tube with closed ends and the axial load on its ends in respect of a general form of electroelastic energy function.
Then, in Sect. 4, we present a bifurcation analysis for the electroelastic tube. First, prismatic bifurcations are analyzed, i.e., we seek incremental solutions of the governing equations for which the cross section is independent of the axial location along the tube. For a specific model of a neo-Hookean electroelastic material we find that prismatic modes of bifurcation are possible under inflation in the presence of an electric field, which is in contrast to the result for the purely mechanical situation obtained in [8]. Second, axisymmetric bifurcations are considered. These are the configurations of the tube for which cross sections remain circular, but with the radius depending on the axial location. The results are compared with those in the purely elastic case: depending on the magnitude of the electric field its presence may make the electroelastic tube more susceptible to bifurcation for each of the two energy functions adopted (neo-Hookean and Mooney-Rivlin-based electroelastic models), i.e., bifurcation becomes possible at lower circumferential stretches than is the case for purely elastic materials without electromechanical interactions. But, for sufficiently large fields bifurcation under inflation is not possible and bifurcation requires external pressure.
For each of the prismatic and axisymmetric bifurcations numerical results are provided based on use of the MATLAB code listed in Appendix B of the PhD thesis [13]. Some important aspects of the code, which employs a numerical scheme used in [8], are discussed and explained briefly in [13].

The equations of nonlinear electroelasticity
We consider a deformable electrosensitive body which occupies the reference configuration B r with boundary ∂B r , in the absence of mechanical loads and electric fields. Application of an electric field and mechanical loads induces deformation which results in a new configuration B, with boundary ∂B, commonly referred to as the deformed configuration. We label material points in the reference configuration B r by their position vectors X and their images in the current configuration B by their position vectors x. The deformation is described by the vector field χ, which relates the position of a particle in the reference configuration to the position of the same particle in the current configuration: x = χ(X). The deformation gradient tensor, denoted F, is defined by where Grad is the gradient operator defined with respect to X. The quantity defined by J = detF (2) represents the local ratio of the volume in the deformed configuration to that in the reference configuration, and for an incompressible material the constraint must be satisfied at each X.
Along with the deformation gradient we shall use the right and left Cauchy-Green deformation tensors, respectively, defined by

Governing equations and boundary conditions
In the purely static context considered here the reduced forms of Maxwell's equations for a dielectric material are where E denotes the electric field vector and D is the electric displacement vector in the configuration B.
The operators curl and div are defined with respect to x. Exterior to a material body the corresponding fields are denoted E and D , which satisfy the standard relation D = ε 0 E in a nonpolarizable material, where ε 0 is the vacuum permittivity. They also satisfy field equations (5). Associated with (5) are the standard boundary conditions where n is the unit outward normal to ∂B and σ f is the free surface charge on ∂B per unit area.
In the absence of mechanical body forces the mechanical equilibrium equation incorporating electric body forces can be written simply as div τ = 0, where τ is the total Cauchy stress tensor, which depends on the deformation and electric field via a constitutive law, which will be discussed in Sect. 2.3. When there are no intrinsic mechanical couples, which we assume to be the case here, then τ is symmetric. The boundary condition associated with (7) is where ∂B t is the part of the boundary where the traction is prescribed, t a being the mechanical part of the traction and t m = τ m n, which is the additional traction due to the Maxwell stress τ m , calculated from the fields outside the body B. The Maxwell stress is defined by where I is the identity tensor.

Lagrangian forms of the electric fields
Lagrangian forms of the electric field vectors are defined by where we recall that J = det F. They satisfy the counterparts of equations (5) in the reference configuration, i.e., CurlE L = 0, DivD L = 0.
Here the operators Curl and Div are defined with respect to X.
In order to obtain a Lagrangian form of the equilibrium equation (7) we introduce the total nominal stress tensor T defined by which is a generalization of the nominal stress tensor in nonlinear elasticity to the electromechanical context. By a standard identity, this yields the Lagrangian form of the electromechanical equilibrium equation (7), i.e., DivT = 0.
The boundary condition associated with (13) can be obtained with the help of the relation which makes use of Nanson's formula nds = JF −T NdS connecting infinitesimal areas ds and dS in the current and reference configurations, n and N being the respective normals to these areas. Then, (8) transforms to where ∂B rt is the pre-image of B t , t A and t M = T M T N (with T M = JF −1 τ m ) being the mechanical traction and the Maxwell traction per unit reference area, respectively.
With the help of Nanson's formula and (10) the boundary conditions (6) convert to Lagrangian form as where N is the unit outward normal to ∂B r and σ F is the free surface charge density per unit area of ∂B r .

Constitutive equations
For the problem discussed in this paper it is convenient to choose D L as the independent electric variable. For mechanically incompressible materials the total stress tensor and the electric field in Lagrangian form are, respectively, where Ω * (F, D L ) is a total energy density function, as defined in [12], and p is a Lagrange multiplier that accommodates the constraint (3). For an incompressible isotropic electroelastic material Ω * depends on F and D L through invariants of the right Cauchy-Green deformation tensor C and D L , i.e., with I 3 = det C = 1 omitted since we are considering incompressible materials. By expressing Ω * in terms of the invariants, then expanding out (17) in terms of the derivatives of Ω * with respect to the invariants, we obtain, on use of (12) with J = 1 and (10) 1 , the total stress τ and electric field E, i.e., where Ω * i stands for ∂Ω * /∂I i , i = 1, 2, 4, 5, 6, and b is the left Cauchy-Green tensor defined in (4) 2 .

Incremental formulation
In this section we summarize the equations governing incremental deformations and electric displacements superimposed on a deformed configuration and an initial electric field. A more detailed discussion of this theory can be found in [14,15].

Incremental equations and boundary conditions.
The increment of each variable is now identified by a superimposed dot. For example,ẋ is the increment in the displacement andḞ = Gradẋ is the corresponding increment in the deformation gradient. The incrementsĖ L ,Ḋ L ,Ṫ must satisfy the incremental governing equations CurlĖ L = 0, DivḊ L = 0, DivṪ = 0.
Outside the material increments in the field vectors are connected byḊ = ε 0Ė and satisfy the equations The incremental forms of the boundary conditions (16) and (15) with J = 1 are whereτ m is the incremental Maxwell stress obtained from (9) and expressed aṡ We now work in terms of the push-forward versions of the increments inĖ L ,Ḋ L andṪ defined by (since J = 1)Ė which are the incremental counterparts of (10) and (12) and also referred to as quantities updated from the reference configuration B r to the deformed configuration B, with updated quantities identified by a zero subscript. A detailed discussion of the concept of updating can be found in [16] for purely mechanical problems, and here we use this approach for updating variables to the deformed configuration for our electromechanical problem.
The governing equations (22) are then updated to and the corresponding boundary conditions are updated to where L =ḞF −1 = gradu, u (=ẋ) being a function of x. The incremental form of (2) is then obtained asJ = Jtr (ḞF −1 ), and for an incompressible material J = 0 and the incremental form of the incompressibility condition (3) becomes Let e 1 , e 2 , e 3 be unit basis vectors in an orthogonal curvilinear coordinate system. Then, in component form, Eq. (29) 3 yields the three scalar equationṡ T 0ji,j +Ṫ 0ji e k · e j,k +Ṫ 0kj e i · e j,k = 0, i = 1, 2, 3, in which summation over repeated indices j and k from 1 to 3 is implied and the notation ,j represents the derivative associated with the jth curvilinear coordinate and is made explicit in Sect. 4 for cylindrical polar coordinates.

Incremental constitutive equations.
IncrementsḞ andḊ L in the deformation gradient and Lagrange electric displacement induce incrementsṪ andĖ L in the nominal stress and Lagrange electric field, leading to incremental forms of the constitutive laws. From (17) we obtaiṅ where A * , A * , A * , which are, respectively, fourth-, third-and second-order tensors, denote electroelastic moduli associated with the total energy Ω * . Their component forms are written which enjoy the symmetries (37) The vertical bar in the component form of A * separates the first two and the third indices, because the first two indices are associated with a second-order tensor, whereas the third one is associated with a vector. The tensor A * maps a vector into a second-order tensor, whereas its transpose maps a secondorder tensor into a vector. In component form this is expressed as A * αi|β = A * T β|αi . The component forms of the equations in Eq. (35) are theṅ where F −1 αi is defined as (F −1 ) αi . When expressed in terms of the invariants I m , m ∈ {1, 2, 4, 5, 6}, the components in (36) can be expanded as where I is the index set {1, 2, 5, 6}, Ω * n = ∂Ω * /∂I n and Ω * mn = ∂ 2 Ω * /∂I m ∂I n , m, n ∈ {1, 2, 4, 5, 6}. The required expressions for the derivatives of the invariants with respect to F and D L in (39)-(41) are given in "Appendix A.1," and the resulting expressions for the moduli are listed in "Appendix A.2." The updated versions of (35) arė ZAMP Bifurcation of finitely deformed thick-walled electroelastic Page 7 of 27 60 and in component form the connections between the electroelastic moduli tensors (36) and their updated forms are The symmetries (37) carry over to the updated versions of the moduli, while A * 0 gains the symmetry We also record here the useful connection

Extension and inflation of a tube
The geometry of a circular tube is described by cylindrical polar coordinates R, Θ, Z in its reference configuration according to where A and B are the internal and external radii and L is the length of a tube. The deformed geometry of the tube with circular symmetry maintained is described in terms of cylindrical polar coordinates r, θ, z as where a, b and l are the radii and the length in the deformed configuration.
For an incompressible material, on which we focus here, the deformation is given by where λ z is the uniform axial stretch of the tube. The azimuthal stretch is denoted λ and given by λ = r/R, and the deformation gradient is diagonal with respect to the cylindrical polar axes. Then, by the incompressibility condition (3), the radial stretch is λ r = λ −1 λ −1 z , and in terms of the independent stretches λ and λ z the invariants I 1 and I 2 are then 3.1.1. Field variables and boundary conditions. Flexible electrodes are attached to the boundaries R = A and R = B in the reference configuration. A potential difference is then applied between the electrodes along with an axial mechanical load and internal pressure, resulting in the deformed configuration described above and charges on the electrodes with equal magnitudes and opposite signs. Let the charges be denoted Q(a) and Q(b), so that Q(a) + Q(b) = 0. Because of the symmetry it can be assumed that the electric field is purely radial and therefore independent of the coordinates θ and z, i.e., the geometry of the tube is such that end effects can be neglected. Let D r (r) then denote the only component of the electric displacement D. The equation (5) 2 satisfied by D yields the specialization so that rD r is a constant, and hence, by evaluation at the boundaries r = a and r = b, we have By Gauss's theorem there is no field outside the tube so that D = 0. The boundary condition (6) 2 applied on r = a and r = b then yields where σ fa and σ fb are the free surface charge densities on the two boundaries, i.e., For this problem the Lagrange variable D L has just a single component, which we denote by D LR , related by the specialization of (10) 2 to D r by D LR = λ −1 r D r = λλ z D r . From the definitions (19) of the invariants I 4 , I 5 , I 6 we then obtain the specializations For an incompressible isotropic electroelastic material the energy function Ω * is expressed as a function of the invariants I 1 , I 2 , I 4 , I 5 , I 6 and written with these arguments: Ω * (I 1 , I 2 , I 4 , I 5 , I 6 ). Then the only nonzero component of the electric field E, denoted E r , is obtained from the specialization of (21) as and Eq. (5) 1 is satisfied automatically.

Stress components and mechanical equilibrium.
From (20) the only nonzero (total) Cauchy stress components are given by (49) and (54), the invariants are functions of two independent stretches and I 4 , we now define a reduced energy function in the form with the invariants on the right-hand side given by (49) and (54). Elimination of p then leads to simple expressions for the stress differences, namely Similarly, the expression (55) is simplified to The equilibrium equation (7) reduces to with which we associate boundary conditions on r = a and r = b. Since there is no electric field outside the tube the Maxwell stress (9) vanishes. Thus, the only load on the surfaces r = a and r = b is the internal pressure P on r = a with zero traction on r = b, i.e., Integration of (62) with the use of (60) 1 and boundary conditions (63) yields where ω * λ = ∂ω * /∂λ.
and from (54) 1 with D LR = λλ z D r , it follows from (51), (52), (53) and (48) that ]. Thus, for any given λ z and initial geometry ZAMP Bifurcation of finitely deformed thick-walled electroelastic Page 9 of 27 60 of the tube, (64) provides an expression for the pressure as a function of the internal radius a and the charge on the electrodes. The axial load N on any cross section of the cylinder is given by and a standard manipulation enables this to be re-cast in the form where F denotes the so-called reduced axial load, which removes from N the contribution of the pressure on the ends of a closed-ended tube.

Bifurcation of an electroelastic circular cylinder
In this section, for consistency with the corresponding analysis in the purely elastic situation in [8], we re-order the coordinates r, θ, z as θ, z, r, associated with the stretches λ, λ z , λ r , respectively. We also use the notation λ 1 , λ 2 , λ 3 for these stretches in the same order, and for the associated components of vectors and tensors we adopt the indices 1, 2, 3. The reduced energy function ω * (λ, λ z , I 4 ) is then replaced by ω * (λ 1 , λ 2 , I 4 ), and stress differences (60) are rewritten as where ω * λ , ω * λz denote the derivatives ∂ω * /∂λ, ∂ω * /∂λ z . The unit basis vectors associated with the cylindrical polar coordinates θ, z, r are denoted as e 1 , e 2 , e 3 , and the derivatives (·) ,k in (34) denoted by subscripts with commas become ∂/r∂θ, ∂/∂z, ∂/∂r for k = 1, 2, 3, respectively. For the cylindrical polar coordinates the only nonzero scalar products e i · e j,k in (34) are The incremental displacementẋ = u is written and the matrix of components of L = gradu with respect to the basis vectors e 1 , e 2 , e 3 is where the subscripts θ, z, r without a preceding comma indicate the corresponding partial derivatives. The incremental incompressibility condition (33) specializes to

Prismatic bifurcations
For prismatic bifurcations u, v and w are independent of z. Furthermore, the equation for w decouples from those for u and v and any solution for w does not affect the cross-sectional shape of the tube. It is therefore convenient to set w = 0. The matrix of components of L then specializes to and the incompressibility condition reduces to Equation (73) is satisfied by taking a function φ(θ, r) such that For the updated incrementḊ L0 in the independent variable D L Eq. (29) 2 reduces to which, similarly to (73), is satisfied by defining a function ψ(θ, r) such thaṫ The remaining governing equations are now used to obtain two coupled equations for φ and ψ. First, we note that the equilibrium equation for i = 2 in (34) is satisfied identically, while the equations for i = 1, 3 yieldṪ 011,1 +Ṫ 031,3 + 1 r (Ṫ 031 +Ṫ 013 ) = 0, T 013,1 +Ṫ 033,3 + 1 r (Ṫ 033 −Ṫ 011 ) = 0.
For the considered underlying cylindrical configuration the components ofṪ 0 in the above two equations are given byṪ where the components of the electroelastic moduli tensors A * 0 and A * 0 are obtained by specializing the general expressions given in "Appendix A.2." On substitution of these expressions into (77) and (78) and use of incompressibility condition (73), we obtainṗ ZAMP Bifurcation of finitely deformed thick-walled electroelastic Page 11 of 27 60 respectively, where a prime denotes differentiation with respect to r and the partial derivatives of the componentsḊ L01 andḊ L03 with respect to r and θ are indicated by a prior comma. Second, Eq. (29) 1 reduces to From (42) 2 appropriately specialized we obtaiṅ with the components of A * 0 and A * 0 given in "Appendix A.2." After substitution of (86) into (85) the latter gives By eliminatingṗ r andṗ θ from (83) and (84) by cross-differentiation and then substituting for expressions (74) and (76) in the result and in (87), followed by some rearrangements, we obtain two coupled equations for φ and ψ, namely We now turn to the boundary conditions associated with the latter two equations. First, the boundary condition (32) is specialized toṪ since E = 0, whereṖ is a prescribed constant increment in P . Hence, by (9) and (27), both the Maxwell stress and its increment vanish. In terms of components this giveṡ By using (81), (45) and the boundary conditions (63), the first of these can be written as while on use of (82) and the incompressibility condition L 11 + L 33 = 0 the second can be rewritten as The term inṗ is eliminated by differentiating (93) with respect to θ (which removes the constantṖ ) and using (83). Then by using the expressions for the components L ij ,Ḋ L01 andḊ L03 in terms of the functions φ and ψ this yields the boundary conditions in which we have also used the boundary condition (63) 1 in the form τ 33 = − P on r = a. And from (92), assuming A * 03131 = 0, we obtain similarly The incremental electric field has to satisfy Eq. (29) 1 , which becomes (87) for the problem considered here, and it also has to satisfy the associated boundary condition (30). Since there is no exterior field, the latter reduces toĖ which, from the constitutive component (86) 1 and the expression for the component L 13 with (74) and (76) 2 , can be written in terms of functions φ and ψ as We now have two coupled partial differential equations (88) and (89) for φ and ψ as functions of r and θ and two boundary conditions in each of (94), (95) and (97). No progress is possible theoretically, and in general, the system of equations and boundary conditions needs to be solved numerically. Toward this we consider solutions for φ and ψ of the form φ = rf n (r) sin nθ and ψ = g n (r) sin nθ for consistency and comparison with the corresponding purely elastic situation considered in [8], n being a nonnegative integer. The equations (88) and (89) are now expressed in terms of the functions f n (r) and g n (r) and their derivatives as To obtain (99) we have used the expression which can be obtained from (62) and (45), and its derivative The boundary conditions (94)
For different values of the parameterσ fa and for λ z = 1 the values of λ a = a/A (the azimuthal stretch at the inner boundary) at which prismatic bifurcations become possible for mode number n = 2 are given in Table 1, along with the corresponding value of λ b = b/B (the azimuthal stretch at the outer boundary). The results for the neo-Hookean electroelastic material withσ fa = 0 are almost identical to those reported in [8] for a different strain-energy function. This is consistent with the findings in [8] that under external pressure the values of λ b (equivalently λ a ) are essentially independent of the form of strain-energy function.
For the considered electroelastic material an important difference from the purely elastic case is that the tube may bifurcate into a prismatic configuration under internal pressure (P > 0). In [8] it was reported that prismatic modes are not possible under internal pressure for the material models considered, but only under external pressure (P < 0). This can also be observed here in Table 1 for the caseσ fa = 0: all nondimensional pressures at which prismatic bifurcations are possible are negative. The values of P/μ in Table 1 were calculated using formula (64) and the specialization of (115) in the form A closed-form expression for P was given in [3] for this constitutive model. In [10], prismatic bifurcations in the absence of an internal pressure were analyzed and solutions showing the shapes of different bifurcation modes were obtained for both the Gent and neo-Hookean models.
Note that for any given underlying boundary-value problem we may prescribe either the charge or potential on the boundary, and the potential difference between r = b and r = a, V say, or charge density σ fa , adjusts according to the formula This was given for the material model (115) in [3] (in a slightly different notation). Thus, results for a fixed potential difference can be obtained from Table 1.

Axisymmetric bifurcations
For axisymmetric bifurcations we take v = 0 and u and w to be independent of θ, so that the components of the displacement gradient specialize to and the incompressibility condition then has the form A function φ(z, r) is then introduced such that and Eq. (119) is satisfied. The governing equation (29) 2 reduces to which can be satisfied by defining a function ψ(z, r) such thaṫ For axisymmetric incremental deformations with no dependence on θ and v = 0 the component of equilibrium equation (34) for i = 1 is satisfied identically, while the components for i = 3, 2 specialize, respectively, toṪ 023,2 +Ṫ 033,3 + T 022,2 +Ṫ 032,3 + 1 rṪ wherein the componentṡ are obtained by specializing (42) 1 . Substitution of these into (123) and (124) and use of (119) yieldṡ p z = A * 03232 w rr + (rA * 03232 + A * 03232 )w r /r + A * 02222 w zz + (A * 02233 + A * 03223 )u rz and when these are substituted into (132) we obtain On use of (120) and (122) this is then expressed in terms of φ and ψ as To express (130) and (131) in terms of φ and ψ we eliminate the terms inṗ by cross-differentiation to arrive at p being given by (102). Thus, we have two coupled equations (135) and (136), to which we now append the boundary conditions. Similarly to Sect. 4.1, the electric field boundary condition (30) reduces tȯ which, on application of the specialization of the incremental constitutive law (42) 2 , yields in terms of functions φ and ψ. The traction boundary condition again has form (90) but now specializes tȯ and hence, on using the components (126) and (129), we obtain and, on use of (129) and the incompressibility condition L 11 + L 22 + L 33 = 0, With the help of (120) and the assumption A * 03232 = 0 Eq. (140) can be expressed in terms of φ as rφ zz − rφ rr + φ r = 0 on r = a, b.
(142) For the boundary conditions (141) the terms inṗ are eliminated by differentiating the equations with respect to z and using the expression forṗ z in (131) to obtain, after use of (120) and (122), and some rearrangement with the help of (45) and (62), in terms of functions φ and ψ.
Analogously to (98) we now write φ(z, r) = rf α (r) cos αz and ψ(z, r) = g α (r) cos αz, where α is a constant to be determined by the boundary conditions on the ends of the tube. The equations (135) and (136) can then be expressed, respectively, as and and the corresponding boundary conditions (138), (142) and (143) as and

Numerical solution.
In order to proceed further we use the incremental boundary condition so that radial displacements at the ends of the cylinder are not allowed (and the incrementṪ 022 in the axial load also vanishes). Therefore, from the previous relation, (120) 1 and (144) 1 , we obtain where n = 1, 2, 3, . . . is the axisymmetric mode number. We see from (151) that α may be changed either by mode number n or the length of the cylinder L. We fix n = 1 and we perform our analysis for different lengths of the cylinder, recognizing that an increase in the mode number n is captured by a decrease in the value of L. We again use the dimensionless variabler = r/A, and as in the case considered for prismatic bifurcations we use the nondimensional quantityσ fa defined in (111) as an electric parameter that is a measure of the charge on the inner boundary of a tube. We also introduce new nondimensional variableŝ and analogously to (108) we define the variableŝ y 1 (r) =f α (r),ŷ 2 (r) =f α (r),ŷ 3 (r) =f α (r), y 4 (r) =f α (r),ŷ 5 (r) =ĝ α (r),ŷ 6 (r) =ĝ α (r).
Then we again have a first-order system of equations written in the formŷ i = M ijŷj , with the coefficients M ij obtained from and boundary conditionŝ Expressions for the electroelastic moduli in respect of energy function (115) are given in "Appendix B.2," in both dimensional and dimensionless forms. We define initial values for the system (154) in the form where δ ik is the Kronecker delta. Each k (= 1, . . . , 6) in (158) corresponds to the solution y k of the system (154), the general solution of which can be written in the form where c k are constants.
We require the solution (159) to satisfy boundary conditions (155)-(157), and we are interested in solutions (159) for which at least one constant c k is nonzero. Substitution of (159) into (155)-(157) leads to a 6 × 6 determinant of coefficients of c k , vanishing of which represents the bifurcation criterion for this problem. This method was also used for the prismatic bifurcations.
In addition to the neo-Hookean electroelastic model (115) we consider its Mooney-Rivlin generalization in the form where μ 1 ≥ 0 and μ 2 ≤ 0 are material constants satisfying μ 1 − μ 2 = μ. In what follows we set μ 1 = 0.8μ and μ 2 = −0.2μ. The associated dimensional and nondimensional electroelastic moduli and their required derivatives are listed in "Appendix B.2." In Fig. 1, for a relatively thin-walled tube with A/B = 0.85, curves of λ a versus λ z for which the bifurcation criterion is satisfied are plotted (continuous curves) together with the curves corresponding to P = 0 (dashed curves) for several values ofσ fa , separately for two values of L/B. These include results forσ fa = 0, which corresponds to the result obtained in [8] for the neo-Hookean purely elastic material (subject to the factor 2 correction in the values of L/B in figure 3 of [8], which was also noted in [19]). Note that, for each value ofσ fa , P is positive above the associated dashed curve and negative below it. For the smaller values ofσ fa the bifurcation criterion is met under inflation for any fixed value of λ z above a minimum value, but asσ fa increases the range of λ z for which bifurcation is possible decreases and eventually disappears for large enoughσ fa . This can be seen, for example, for the valueσ fa = 1.4 in Fig. 1a, for which the dashed curve is entirely above the bifurcation curve, in which case bifurcation under inflation is not possible. Bifurcation is possible for negative P , which amounts to bifurcation under external pressure. In this and the subsequent figures the uppermost curve corresponds to the purely elastic case (or, equivalently, to the case when there is no potential between the electrodes and thus no electric field). As the electric field is applied and increased, the corresponding curves drop to lower values of λ a .
For the two values of L/B used in Fig. 1, the corresponding results are shown for a thicker-walled tube with A/B = 0.5 in Fig. 2, and with the same values ofσ fa . This shows that the qualitative features are the same as in Fig. 1. The quantitative differences are that, for a given value of λ z , larger values of the azimuthal stretch are required to achieve bifurcation under inflation corresponding to higher pressure,  Fig. 3. Again, the main qualitative features are unchanged, but as λ z increases, larger values of λ a are required for bifurcation with correspondingly larger pressures, and there is no phasing out of bifurcation under inflation for the considered range of values ofσ fa . The nonmonotonicity of the λ a versus λ z curves after the maxima evident in Figs. 1a and 2a is more pronounced in Fig. 3b. However, this does not correspond to any form of snap-through instability since the values of λ z (when there are more than one) associated with a given λ a correspond here to significantly different values of P , whereas snap through between different values of λ z occurs at fixed P . Figure 4 shows representative plots for the Mooney-Rivlin material (160) for A/B = 0.85 analogous to those in Fig. 1. For the purely elastic case considered in [8], the values of μ 1 and μ 2 used were not specified so a direct comparison of our results with those in [8] is not possible. However, our results are qualitatively very similar to those reported in Figure 3 of [8]. For this material there is no phasing out of the bifurcation under internal pressure unless the electric field becomes very large (not shown), but as for the neo-Hookean material increasing electric field strength reduces the values of λ a for which bifurcation can occur for those fixed λ z for which bifurcation is possible under inflation. Moreover, the range of values of λ z for which bifurcation is possible is more limited than for the neo-Hookean model.
The results described herein are for fixed values ofσ fa . Equation (117) gives the corresponding values of the potential difference and can be used to translate the results to those for a fixed potential difference.
In considering the propagation of axisymmetric electroelastic waves in a tube in the absence of pressure Shmuel and deBotton [11] obtained results for vanishing wave speed, which corresponds to bifurcation. Such results can be read off as special cases of the results here corresponding to points in the above figures where the (dashed) zero pressure curves intersect with the corresponding bifurcation curves. Note, however, that in [11] an incremental electric field exterior to the body needed to be accounted for.
The nonzero second derivatives with respect to F are The second derivatives of I 4 , I 5 , I 6 with respect to D L are where C 2 αβ is defined as (C 2 ) αβ . The mixed derivatives of I 1 , I 2 and I 4 with respect to F and D L are equal to zero. For I 5 , I 6 they are

A.2. Electroelastic moduli
We introduce the notationb = I 1 b − b 2 , D (1) = bD and D (−1) = b −1 D and obtain the following expressions for the components of the electroelastic moduli tensors: