Instabilities in multilayered soft dielectrics
Katia Bertoldi∗ and Massimiliano Gei† July 15, 2010

Abstract Experimental observations clearly show that the performance of dielectric elastomericbased devices can be considerably improved using composite materials. A critical issue in the development of composite dielectric materials toward applications is the prediction of their failure mechanisms due to the applied electromechanical loads. In this paper we investigate analytically the influence of electromechanical finite deformations on the stability of multilayered soft dielectrics under plane-strain conditions. Four different criteria are considered: i.) loss of positive definiteness of the tangent electroelastic constitutive operator, ii.) existence of diffuse modes of bifurcation (microscopic modes), iii.) loss of strong ellipticity of the homogenized continuum (localized or macroscopic modes), and iv.) electric breakdown. While the formulation is developed for generic isotropic hyperelastic dielectrics, results are presented for the special class of ideal dielectrics incorporating a neo-Hookean elastic response. The effect of material properties and loading conditions is investigated, providing a detailed picture of the different possible failure modes.

Keywords: Electroactive Polymers, Bifurcations, Localization of Deformation, Composite Materials, Finite Elasticity.

1

Introduction

The application of a voltage through electrodes to soft dielectric elastomers deforms them substantially, giving us the opportunity to use this principle to design a new class of actuators. Discovery and development of these materials were first reported in the works of Pelrine and coworkers (Pelrine et al., 1998, 2000). Immediately, they attracted significant interest because of tremendous potential in areas as robotics, aerospace and biomedical engineering. Currently, dielectric elastomers are widely employed to manufacture devices as reliable electrically-driven actuators, manipulators and energy harvesters
School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA; email: bertoldi@seas.harvard.edu; web-page: www.bertoldi.seas.harvard.edu. † Department of Mechanical and Structural Engineering, University of Trento, Via Mesiano 77, I-38123 Trento, Italy; email: massimiliano.gei@unitn.it; web-page: www.ing.unitn.it/∼mgei.
∗

1

(see Bar-Cohen, 2001; Carpi et al., 2008a; Carpi and Smela, 2009, and references cited therein). A significant challenge in the development of devices based on dielectric elastomers is that they often require the application of extremely high voltages as a result of the material low dielectric constant. This represents a clear limitation in their development toward further applications, but both experimental (Zhang et al., 2002; Huang et al., 2004; Carpi et al., 2008b) and analytical (deBotton et al., 2007) investigations showed that composite materials can provide a solution to this critical issue. When stiff and high-permittivity particles are included in a soft elastomeric matrix, the overall dielectric constant of the material increases considerably, while its deformability may be only marginally affected. DeBotton et al. (2007) showed that the use of biphasic laminate dielectric composites can improve the actuation strain up to 50% or, on the other hand, can provide the same actuation with a sensible decrease of the applied voltage. However, to design a new class of optimized devices based on dielectric composites, further investigations are necessary focusing both on their failure under the applied loads and on the effect of inclusions volume fraction, geometry and material properties on their performance. In this work we concentrate on the latter issue, investigating in a systematic way instabilities that develop in finitely deformed multilayered dielectrics. The optimization of the performance of dielectric elastomer actuators is a challenging task due to their multiple failure modes. For single-phase actuators, electromechanical instabilities (unstable thinning of the actuator, local buckling induced by coexistent states, electric breakdown) have been explore both analytically (Zhao and Suo, 2007; Zhao et al., 2007; Moscardo et al., 2008; Liu et al., 2009) and experimentally (Plante and Dubowsky, 2006), providing also design guidelines to prevent potential failures under operating conditions. Moreover, Dorfmann and Ogden (2010) recently investigated surface instabilities for an electroelastic half space. In composite systems instabilities are even more critical and a larger family of potential failure modes must be considered. In this scenario modeling represents a fundamental tool. Motivated by the development of applications, the nonlinear theory of soft dielectrics, first proposed by Toupin (1956), has been recently reviewed and further developed. In particular, we refer to the work of McMeeking and Landis (2005), Dorfmann and Ogden (2005) and Suo et al. (2008), where both the concepts of Maxwell and total stresses and the formulation of constitutive equations for conservative materials have been discussed. In addition, based on this theory a variational formulation has been built and discretized using the finite element method (Vu et al., 2007). In this paper the stability of multilayered hyperelastic dielectric elastomers deforming at large strains is systematically investigated. Four instability criteria for composites are introduced1 , namely i.) loss of positive definiteness of the tangent electroelastic constitutive operator; ii.) existence of diffuse modes of bifurcation (microscopic modes);
It is important to note that throughout the paper we refer generically to i.)–iv.) as ‘instability’ criteria. More precisely, i.)–iii.) detect a possible bifurcation point along the fundamental deformation path, while iv.) corresponds to a failure threshold for the solid.
1

2

iii.) loss of strong ellipticity of single phases and of the homogenized continuum (localized modes or, in the latter case, macroscopic modes); iv.) electric breakdown. The first three criteria follow from the theory of bifurcation and stability for nonlinear elastic solids developed by Hill (1957) and Biot (1965) and applied subsequently to investigate loss of uniqueness of given loading paths in boundary-value problems for both homogeneous solids (Hill and Hutchinson, 1975; Ogden, 1984; Needleman and Ortiz, 1991; Triantafyllidis and Lehner, 1993) and composite materials (Triantafyllidis and Maker, 1985; Geymonat et al., 1993; Bigoni and Gei, 2001; Triantafyllidis et al., 2006; Michel et al., 2007). Electric breakdown is instead specific for dielectric materials that are characterized by a limit maximum value for the intensity of the electric field, beyond which electric discharges may take place. The four instability criteria are then specialized to rank-one layered composites finitely deformed under plane-strain conditions. A detailed analysis of instabilities is reported for a multilayer with two phases made of ideal dielectrics incorporating a neo-Hookean elastic response (Dorfmann and Ogden, 2005). Interestingly, the results clearly show that depending on the heterogeneity contrast between the phases and on the loading conditions different failure modes may occur.

2
2.1

Theory of elastic dielectrics
Basic notation

In this section we summarize the equations governing the nonlinear electrostatic deformation (electrodynamical effects are excluded) of heterogeneous dielectrics following the formulation previously introduced by McMeeking and Landis (2005), Dorfmann and Ogden (2005), deBotton et al. (2007) and Suo et al. (2008). Let us consider an isolated system consisting of a multi-phase electroelastic body and 0 the complemental surrounding space (Fig. 1) and indicate by B0 and B sur = R3 \ B0 the undeformed stress-free configuration of the body and the surrounding space, respectively. 0 We identify with ∂B0 the boundary separating B0 from the surrounding, while ∂Bint denotes the set of all the internal interfaces between heterogeneities in B0 . The application of both mechanical loadings and electric fields deforms quasi statically the body from B0 0 to the current configuration B and interfaces ∂B0 and ∂Bint to ∂B and ∂Bint , respectively. Such deformation is described by the function χ that maps a reference point x0 in B0 to its deformed position x = χ(x0 ) in B. The associated deformation gradient will be denoted by F = ∂χ/∂x0 , while J identifies its determinant, J = det F. If the surrounding space does not consists of vacuum, the deformation χ can be extended to B0 , yielding sur Bsur = R3 \ B = χ(B0 ). sur At this point it is important to note that when electric and/or magnetic interactions are considered, the solution of a boundary-value problem requires the integration of the governing equations over the entire system (i.e. the body and the surrounding space). As a consequence, specific relations for the jumps over ∂B must be satisfied by the relevant 3

surrounding space
1 ¶B int

¶Bv n B ¶B int
2

t

t

++ w++

¶B int ¶Bt n

3

in the current configuration. In the sketch ∂Bint = ∂B 1 ∂B 2 ∂B 3 , while ∂Bv and ∂B t (∂B = int int int ∂B v ∪ ∂B t ) denote the portions of ∂B where displacements and tractions are applied, respectively.

Figure 1: An isolated system consisting of a heterogeneous dielectric body and the surrounding space

quantities entering the governing equations.

2.2

Spatial formulation of the governing equations

Equilibrium of the dielectric body is ensured when the total stress τ , electric displacement D and electric field E satisfy div τ + b = 0, τT = τ, div D = 0, curl E = 0 (in B ∪ Bsur ), (1)

b being the ‘mechanical’ specific body force. Eqns. (1)1,2 are associated with equilibrium of momenta, while eqns. (1)3,4 are the relevant Maxwell’s equations in the quasi– electrostatic limit. Eqn. (1)4 implies that the electric field is conservative, so that it can be written as a function of the electric potential φ(x) E = −gradφ, φ(x) being a continuous function defined over the entire system, since double layers of charge are not considered. The electric field causes the alignment of dipoles inside the body leading to its polarization, which is described by a vector field P, P = D − ǫ0 E. To specify boundary conditions, the subsets ∂B v and ∂B t (∂B v ∪ ∂B t = ∂B, ∂B v ∩ ∂B t = ∅) are introduced, where displacements and surface tractions are prescribed, respectively. Thus the boundary conditions across ∂B are given by (see Fig. 1) [[v]] = 0, [[τ ]]n = t (on ∂B t ), 4 ˜ v=v (on ∂Bv ), (2)

˜ where v denotes the finite displacement field, v and t are prescribed vector fields, n is the current outward normal to ∂B, ω are surface charges2 and [[ · ]] is the jump operator defined on ∂B as [[f ]] = f solid − f sur . Moreover, assuming that no mechanical tractions and surface charges are applied to ∂Bint and that heterogeneities in the body are perfectly bonded, the boundary conditions across ∂Bint separating phases ‘a’ and ‘b’ are given by [[v]] = 0, [[τ ]]n = 0, [[D]] · n = 0, n × [[E]] = 0 (on ∂Bint ), (3)

[[D]] · n = −ω,

n × [[E]] = 0

(on ∂B),

where now the jump operator is defined as [[f ]] = f b − f a and n is the unit normal pointing towards phase ‘a’. In the particular case of the surrounding space consisting of vacuum, the boundary conditions (2)2,4 specialize to τ n = t + τ ∗ n, D · n = −ω + ǫ0 E∗ · n,
∗

(4) denotes quantities

respectively, where τ and D are evaluated in the body, whereas evaluated in the vacuum, 1 τ ∗ = ǫ0 E∗ ⊗ E∗ − E ∗2 I , 2 D∗ = ǫ0 E∗ ,

ǫ0 being the permittivity of vacuum (ǫ0 = 8.85 pF/m). Equation (4)1 clearly shows that the Maxwell stress τ ∗ acting on the boundary separating the dielectric body from vacuum affects the total stress τ in the body.

2.3

Lagrangian formulation of the governing equations
SF = (SF)T , Div D0 = 0, Curl E0 = 0 (in B0 ∪ B 0 ), (5) sur

Integration of eqns. (1) over B and change of variable from x to x0 yield3 Div S + b0 = 0, which represent the lagrangian formulation of the field equations for the electroelastic response of the system. In the equations above S = Jτ F−T is the total first PiolaKirchhoff stress, D0 = JF−1 D is the lagrangian (or nominal) electric displacement, E0 = FT E is the lagrangian (or nominal) electric field4 and Div and Curl denote the
Often charges are accumulated on the surface through electrodes connected to a battery which imposes a constant electric potential. In this case the boundary condition (2)4 must be replaced by a boundary condition setting the appropriate value of φ over the region where the electrodes are applied. 3 Denoting by dV 0 an infinitesimal volume in B 0 , by dA0 an infinitesimal area on ∂B 0 and by n0 its outward unit normal, b0 dV 0 corresponds to the infinitesimal current body force, t0 dA0 to the infinitesimal current contact force and ω 0 dA0 to the infinitesimal current surface charge. Note that dA0 can be expressed in terms of the current quantities using the Nanson’s formula n dA = JF−T n0 dA0 , whereas dV = J dV 0 . 4 We avoid to define a ‘lagrangian polarization vector’ since there is no Maxwell’s equation associated with this field.
2

5

div and curl operators in B 0 . Since the electric field is conservative, in the reference configuration E0 = −Gradφ0 , where φ0 (x0 ) = φ(x) is the lagrangian (nominal) electric potential. The boundary conditions (2) may be rewritten in lagrangian form as [[v0 ]] = 0, [[S]]n0 = t0 (on ∂B 0 ), t ˜ v0 = v0 (on ∂B 0 ), v (6)

where v0 (x0 ) = v(x). The lagrangian form of eqns. (3) can be obtained analogously. We would like to remark that when vacuum surrounds the body, there is no deformation associated with it, so that the definition of lagrangian quantities in B0 has no sur meaning. Therefore, eqns. (4) can be rewritten in lagrangian form as Sn0 = t0 + Jτ ∗ F−T n0 bnd where Fbnd = F|∂B0 . and D0 · n0 = −ω 0 + ǫ0 JF−1 E∗ · n0 , bnd (7)

[[D0 ]] · n0 = −ω 0 ,

n0 × [[E0 ]] = 0 (on ∂B 0 ),

2.4

Constitutive equations

Let us consider a conservative material whose response is described by the free-energy function W = W (F, D0 ). Application of the first and second law of thermodynamics yields (McMeeking and Landis, 2005; Dorfmann and Ogden, 2005; Suo et al., 2008) S= ∂W , ∂F ∂W , ∂D0 τ = 1 ∂W T F , J ∂F ∂W . ∂D0 (8)

whereas the electric field is obtained as E0 = E = F−T (9)

For an incompressible material (i.e. J = 1) eqns. (8) modify as S= ∂W − p F−T , ∂F τ = ∂W T F − p I, ∂F (10)

where p is an arbitrary hydrostatic pressure. Objectivity requires the free energy W (F, D0 ) to be written as a function of the right Cauchy-Green tensor C = FT F and D0 , which are both invariant under arbitrary rigidbody motions. Moreover, for an isotropic hyperelastic dielectric body the free energy W (C, D0 ) can be expressed as a function of the invariants of C (see Dorfmann and Ogden, 2005) I1 = tr C, I2 = 1 (tr C)2 − tr C2 , 2 I3 = det C = J 2 ,

and of three additional invariants depending on D0 , namely I4 = D 0 · D 0 , I5 = D0 · CD0 , 6 I6 = D 0 · C 2 D 0 .

Due both to the lack of available experimental data and to the desire of a simple enough formulation that allows a better understanding of the material response, an uncoupled form for the free energy is often considered W = Welas (I1 , I2 , I3 ) + Wpol (I4 , I5 , I6 ), (11)

where Welas is the strain energy when electric effects are disregarded, whereas Wpol represents the contribution of the polarization of the solid. Several investigations (Eringen and Maugin, 1989; McMeeking and Landis, 2005; Zhao et al., 2007; Zhao and Suo, 2008) showed that the uncoupled free energy (11) well captures the behavior of large classes of soft dielectrics such us ideal dielectrics and electrostrictive materials. In Section 5 results will be presented for an incompressible dielectric (I3 = 1) whose behavior is captured using the free energy µ 1 W nH = (I1 − 3) + I5 , (12) 2 2ǫ where a neo-Hookean form has been adopted for Welas and the dependency of Wpol on I4 and I6 has been neglected. In eqn. (12) µ is the shear modulus in the reference configuration and ǫ = ǫ0 ǫr denotes the dielectric constant of the solid taken to be unaffected by the deformation, so that the free energy (12) predicts a linear proportionality between electric displacement and electric field along any loading path. The adoption of the form (12) is justified by experimental data on typical materials for soft dielectrics (see, e.g., Goulbourne et al., 2005, and Kofod et al., 2003) showing that i.) the simple form of the free energy W nH describes well the electromechanical behavior of these materials up to a maximum tensile stretch of about two, while for highly strained specimens more refined models such those proposed by Ogden and Yeoh must be employed; ii.) the dielectric permittivity ǫ is almost insensitive on the large strain applied to the specimen.

2.5

Incremental boundary-value problem

The mechanics of incremental deformations superimposed upon a given state of finite deformation allows the investigation of instabilities that develop in dielectric bodies subjected to nonlinear electrostatic deformations. Here, in the spirit of Dorfmann and Ogden (2010), we focus on conservative dielectric materials that are isotropic in the reference configuration and satisfy at each stage of deformation eqns. (5)-(6). 2.5.1 Lagrangian formulation

˙ Let us consider a perturbation t0 and ω 0 of the tractions and the surface charges applied ˙ 0 on ∂B that takes the body to a new equilibrium configuration where eqns. (5)-(6) are still satisfied and leaves the body force density b0 unchanged. The incremental problem is governed by ˙ Div S = 0, ˙ Div D0 = 0, ˙ Curl E0 = 0 (in B0 ∪ B 0 ), sur (13)

7

˙ ˙ ˙ where S, D0 and E0 denote the increments of total first Piola-Kirchhoff stress, nominal electric displacement and nominal electric field caused by the perturbation. Moreover, the incremental jump conditions at the external boundary of the body take the form ˙ [[x]] = 0, ˙ ˙ [[S]]n0 = t0 (on ∂B 0 ), t ˙ x=0 (on ∂B 0 ), v (14)

˙ ˙ x = χ(x0 ) denoting the incremental deformation. Similar boundary conditions hold at 0 ∂B int , namely ˙ [[x]] = 0, ˙ [[S]]n0 = 0, ˙ [[D0 ]] · n0 = 0, ˙ n0 × [[E0 ]] = 0 (on ∂B 0 ). int (15)

˙ [[D0 ]] · n0 = −ω 0 , ˙

˙ n0 × [[E0 ]] = 0 (on ∂B 0 ),

When the dielectric body is surrounded by vacuum, the incremental counterpart of eqns. (7) is given by ˙ ˙ ˙ ˙ ˙ Sn0 = t0 + S∗ n0 , D0 · n0 = −ω 0 + D0∗ · n0 , ˙ (16) ˙ ˙ where S and D0 are evaluated in the body, while S∗ = Jτ ∗ F−T and D0∗ = ǫ0 JF−1 E∗ bnd bnd are evaluated in vacuum. Assuming that all incremental quantities are sufficiently small, the constitutive equations (8)1 –(9)1 can be linearized as ˙ ˙ ˙ S = C0 F + B0 D0 ,
∗ ˙ ˙ ˙ E0 = B0T F + A0 D0 ,

0 0 ˙ ˙ ˙0 SiJ = CiJkL FkL + BiJM DM , 0 ˙0 ˙ ˙0 EM = BiJM FiJ + A0 N DN , M

(17)

∗ 0 ˙ ˙ where F = Gradχ and B 0T M iJ = BiJM . It follows from eqns. (8)1 -(9)1 and (17) that the components of the electroelastic moduli tensors C0 , B0 and A0 are given by

0 CiJkL =

∂2W , ∂FiJ ∂FkL

0 BiJM =

∂2W , 0 ∂FiJ ∂DM

A0 N = M

∂2W , 0 0 ∂DM ∂DN

(18)

which imply the symmetries
0 0 CiJkL = CkLiJ ,

A0 N = A0 M . M N

(19)

For incompressible materials the incremental total first Piola-Kirchhoff stress tensor and lagrangian electric field are given by ˙ ˙ ˙ ˙ S = C0 F + pF−T FT F−T − pF−T + B0 D0 , ˙
∗ ˙ ˙ ˙ E0 = B0T F + A0 D0 ,

(20)

where the lagrange multiplier p has been introduced by the incompressibility constraint ˙ ˙ −1 ) = 0. tr(FF

8

2.5.2

Updated lagrangian formulation

Instabilities are often investigated formulating the incremental boundary value problem in an updated lagrangian formulation, where the reference configuration moves and is identified with the current configuration. Push-forward transformations based on linear momentum balance and divergence theorem (see Appendix A) allow the introduction of the incremental updated quantities Σ= 1˙ T SF , J 1 ˙ ˆ D = FD0 . J (21)

ˆ The work-conjugate variables to Σ and D can be easily obtained from the second-order work, namely ˙ ˙ ˙ ˙ Pext = {φ0 ω 0 + t0 · x} dA0.
∂B0

Substitution of eqns. (14)2,4 and application of the divergence theorem yield Pext =
B
0

∪B0 sur

˙ ˙ ˙ ˙ {E0 · D0 + S · F} dV 0 = ˆ ˙ E = F−T E0

B∪Bsur

ˆ ˆ {E · D + Σ · L} dV,

(22)

where

(23)

˙ is the updated incremental lagrangian electric field, u(x) = x, L = grad u and eqns. (15) have been taken into account. Eqn. (22) highlights the analogy between surface charge and surface tractions, displacements and electric potentials. Introduction of the incremental updated quantities into eqns. (13) yields (see Appendix A) div Σ = 0, ˆ div D = 0, ˆ curl E = 0 (in B ∪ Bsur ), (24)

which represent the updated lagrangian formulation of the field equations governing the incremental problem. The incremental boundary conditions relating the jumps in the updated lagrangian formulation to the perturbations in the reference configuration can be written as [[u]] = 0, ˙ [[Σ]]n dA = t0 dA0 (on ∂B t ), ˆ n × [[E]] = 0 u=0 (on ∂B v ), (25)

and

ˆ [[D]] · n dA = −ω 0 dA0 , ˙ [[u]] = 0, [[Σ]]n = 0,

(on ∂B), (on ∂B int ). (26)

ˆ [[D]] · n = 0,

When vacuum surrounds the dielectric body eqns. (25)2,4 reduce to ˙ Σn dA = t0 dA0 + Σ∗ n dA, where ˙ Σ∗ = τ ∗ + τ ∗ tr L|bnd I − (L|bnd )T , 9 ˆ ˙ D∗ = ǫ0 E∗ + tr L|bnd I − L|bnd E∗ , ˆ ˆ D · n dA = −ω 0 dA0 + D∗ · n dA, ˙ (27)

ˆ n × [[E]] = 0

˙ E∗ denoting the increment of the electric field in the vacuum and with L|bnd = (grad u)|∂B . Introduction of eqn. (21) into equations (17) yields ˆ Σ = CL + BD, ˆ ˆ E = BT ∗ L + AD, (28)

where the components of the updated constitutive tensors are given by Ciqkp = 1 0 C FpL FqJ , J iJkL
0 −1 Biqa = BiJM FqJ FM a , −1 −1 Aab = J A0 N FM a FN b . M

(29)

For an incompressible material eqn. (28)1 modifies as ˆ Σ = CL + p LT − p I + BD, ˙ (30)

whereas eqns. (29) are still valid substituting J = 1. Analogously to C0 and A0 the updated constitutive tensors C and A possess major symmetries Ciqkp = Ckpiq , Aab = Aba . Moreover, as pointed out by Dorfmann and Ogden (2010), substitution of eqn. (28)1 into the equality Σ − ΣT = Lτ − τ LT provides Ciqkp + τip δkq = Cqikp + τpq δik , Biqa = Bqia , (31)

for compressible materials, while for incompressible materials eqn. (31)1 becomes Ciqkp + (τip + p δip )δkq = Cqikp + (τpq + p δpq )δik . In Section 5.1 the components of the incremental constitutive tensors defined above will be given for the class of dielectrics whose response is described by the free energy (12).

3

Homogenization of multilayered soft dielectrics

Failure mechanisms occurring in composite materials are often predicted by investigating appropriately homogenized properties of the composite. Thus at this stage, following deBotton (2005) and deBotton et al. (2007), the homogenized response of incompressible rank-one layered soft dielectrics finitely deformed under plane-strain conditions is derived. Here we focus on a biphasic layered structure whose constituents, a and b, have volume fractions ca = h0 a / h0 a + h0 b and cb = 1 − ca , h0 a and h0 b denoting the thicknesses of the two phases in the reference configuration B0 (see Fig. 2). The macroscopic deformation gradient Fav , the macroscopic total first Piola-Kirchhoff stress Sav , the lagrangian electric field E0 av and the lagrangian electric displacement D0 av are defined as Fav = ca Fa +cb Fb , Sav = ca Sa +cb Sb , E0 av = ca E0 a +cb E0 b , D0 av = ca D0 a +cb D0 b , (32)

10

RVE
0 x2

m0 n0
0 x1

h0b b D
0av

a

h0a

Figure 2: Schematic representation of a biphasic multilayered dielectric deformed under plane-strain
conditions.

respectively. In the absence of charges along the interface the electric displacement continuity equation ([[D0 ]] · n0 = 0) can be recast as D0 a · n0 = D0 b · n0 , and can be alternatively expressed as D0 a − D0 b = βm0 , (33)

where n0 is the unit normal to the layers, β is a real parameter and m0 is a unit vector parallel to the layers, such that n0 · m0 = 0. It follows from eqns. (32)4 and (33) that the electric displacement field in the two phases can be written as D0 a = D0 av + cb βm0 , D0 b = D0 av − ca βm0 . (34)

Following a similar approach, deBotton (2005) has shown that interface compatibility ([[F]]m0 = 0) requires Fa = Fav I + α cb m0 ⊗ n0 , Fb = Fav I − α ca m0 ⊗ n0 , (35)

where α is a real parameter. Therefore, the macroscopic free energy of the composite is given by (36) W av = ca W a (Fav , D0 av , α, β) + cb W b (Fav , D0 av , α, β), where α and β are the scalar parameters previously introduced. Their values are obtained imposing continuity of tractions ([[S]]n0 = 0) and tangential components of the electric field (n0 × [[E0 ]] = 0) at the interface. Once α and β are determined, the macroscopic total stress and electric field can be obtained through eqns. (32)2,3 . In Section 5.2 the formulation will be specialized to the case of a layered structure made of two phases whose response is described by the free energy (12). 11

4

Instabilities of multilayered soft dielectrics

Since a limiting factor in the design of composite materials is their failure under the applied loads, following the pioneering work of Hill (1957) the investigation of the stability of composites subjected to purely mechanical loadings has attracted considerable attention (Ogden, 1984, Triantafyllidis and Maker, 1985; Geymonat et al., 1993; Triantafyllidis et al., 2006; Michel et al., 2007). Here, guided by well-established criteria for the ‘pure’ mechanical case, for the first time the onset of instabilities for elastic dielectric composites finitely deformed by the application of both electrical and mechanical loadings is investigated. Three classical bifurcation criteria widely used in the mechanical case are considered and adapted, namely, i.) loss of positive definiteness of the tangent electroelastic constitutive operator, ii.) diffuse-mode (microscopic) instability, iii.) loss of strong ellipticity (macroscopic instability when the analyzed solid is the homogenized continuum). In addition, –iv.)– electric breakdown is investigated, which may represent a critical failure mode when electric fields are considered. Taking as reference the theory developed by Hill (1957) (see also Petryk, 1993) for nonlinear elastic solids, loading processes producing a unique response –the ‘principal equilibrium path’– that can be parameterized in terms of a scalar loading parameter t are considered. For each t ≥ 0, an admissible set of incremental deformation and lagrangian ˙ ˙ electric displacement is defined as a pair (χ, D0 ), with both fields not identically zero, that satisfies the continuity conditions (15) at internal interfaces and such that ˙ χ=0 (on ∂B 0 ), v ˙ [[D0 ]] · n0 = 0 (on ∂B 0 ).

At a critical point tcr along the principal equilibrium path a bifurcation may occur, so that two solutions (1 and 2) of the incremental boundary-value problem illustrated in Section 2.5 exist. We indicate by ∆ the differences between the two possible solutions in any field quantity, i.e. ∆(·) = (·)(1) − (·)(2) . Moreover, dead-loading conditions are considered with the applied tractions t0 and the surface charge ω 0 independent of the ˙ deformation, so that t0 = 0 and ω 0 = 0. ˙ ˙ When the incremental boundary conditions are given by eqns. (14) with t0 = 0 and ˙ ∆D0 , ∆E0 } satisfies the incremental ˙ ˙ ˙ ω0 = 0, it is easy to show that the set {∆χ, ∆S, ˙ boundary-value problem. Application of the principle of virtual displacements and of the divergence theorem yields ˙ ˙ ˙ ˙ {∆S · ∆F + ∆E0 · ∆D0 } dV 0 = 0. (37)

B

0

∪B0 sur

˙ ˙ Noting that χ(2) = D0(2) = 0 represents a solution of the considered homogeneous incremental boundary-value problem, it follows from eqn. (37) that a sufficient condition excluding any bifurcation is ˙ ˙ ˙ ˙ {S(t) · F + E0 (t) · D0 } dV 0 > 0 12 (38)

B

0

∪B0 sur

for all admissible sets, where superscript (1) has been omitted for notational simplicity. Since in the reference configuration, identified by t = 0, condition (38) is certainly satisfied, a bifurcation along the loading path occurs as soon as the functional in (38) becomes positive semi-definite, with the equality holding for at least an admissible set ˙ ˙ called primary eigenmode and indicated by (χcr , D0 ), i.e. cr
B0 ∪B0 sur

˙ ˙ ˙ ˙ {Scr (t) · Fcr + E0 (t) · D0 } dV 0 = 0, cr cr

(39)

˙ ˙0 where Scr (t) and Ecr (t) are given by the incremental constitutive equations (17). Although in this study we focus on infinitely extended solids, we note that condition (38) can be easily adapted to investigate dielectric bodies surrounded by vacuum [with boundary conditions specified by eqn. (16)]. We are now in the position to introduce and briefly discuss the four bifurcation criteria investigated in the paper. • Loss of positive definiteness of the tangent electroelastic constitutive operator (PD). Let us start focusing on the pure mechanical problem, so that the critical bifurca˙ ˙ ˙ ˙ tion condition (39) specializes to B0 S(t) · F dV 0 = 0 for at least a χ = χcr = 0, ˙ · F dV 0 > 0 for all other admissible χ. When only dead-load trac˙ ˙ while B0 S(t) 0 ˙ tions are applied to the boundary (i.e. ∂B v = ∅) and F is a homogeneous field, for compressible materials the bifurcation condition can be rewritten as (Hill, 1967) ˙ ˙ ˙ ˙ C0 (t)F · F = 0, for at least a F = Fcr = 0, which corresponds to loss of positive definiteness of the fourth-order tangent constitutive tensor. Therefore, for a homogeneously deformed body with only tractions applied to its boundary, the global bifurcation condition reduces to a local condition that depends only on the principal equilibrium path. Similarly, for a soft dielectric body homogenously deformed applying only deadload tractions/charges to its boundary, introduction of the incremental constitutive equations (17) into the critical bifurcation condition (39) yields ˙ ˙ ˙ ˙ ˙ ˙ C0 (t)F · F + 2B0 (t)D0 · F + A0 (t)D0 · D0 = 0 (40)

˙ ˙ ˙ ˙ for at least a pair (F, D0 ) = (Fcr , D0 ) = 0. cr ˙ and D0 are considered. Therefore, bifurcation is predicted ˙ when homogeneous F when positive definitiveness is lost by the quadratic form in eqn. (40), which can be identified as the ‘tangent electroelastic constitutive operator’. For an incompressible material the incremental total first Piola-Kirchhoff stress tensor is given by (20)1 , so that the bifurcation condition (40) modifies as ˙ ˙ ˙ ˙ ˙ ˙ ˙ C0 (t)F · F + p(t) tr[(FF−1 )2 ] + 2B0 (t)D0 · F + A0 (t)D0 · D0 = 0 ˙ ˙ ˙ ˙ ˙ for at least a pair (F, D0 ) = (Fcr , D0 ) = 0 with tr(Fcr F−1 ) = 0. cr 13 (41)

The critical loading parameter t at which either condition (40) or (41) is satisfied will be denoted by tPD . In Section 4.1.2 the formulation will be specialized to the case of a multilayered consisting of two phases whose response is described by the free energy (12). • Existence of diffuse modes of bifurcations (microscopic modes).

Diffuse buckling modes with a length scale comparable to the size of the heterogeneity may be explored investigating the propagation of small-amplitude perturbations of arbitrary wavelength superimposed on the current state of deformation (Dowaikh and Ogden, 1990; Triantafyllidis and Lehner, 1993; Bertoldi et al., 2008; Gei, 2008). While a real natural frequency corresponds to a propagating wave, a complex natural frequency identifies a perturbation exponentially growing with time. Therefore, the transition between a stable and an unstable configuration identified by the critical loading parameter tMicro is detected when the frequency vanishes and can be investigated using the quasi-static formulation of the incremental problem. For the purely mechanical case the relation between the failure of Hill’s condition (38) and the existence of zero-speed waves has been investigated in detail (see, i.e., Needleman, 1976; Triantafyllidis, 1980), showing that diffuse modes can be alternatively identified at the failure of the critical condition (39) expressing the admissible incremental deformations in Fourier series. For periodic solids of infinite extent, Geymonat et al. (1993) (see also Triantafyllidis and Bardenhagen, 1996, and Triantafyllidis et al., 2006) rigourously showed that in the pure mechanical case the same approach can be used and that it is sufficient to investigate the incremental boundary-value problem of the elementary unit cell of the composite and apply the Floquet-Bloch boundary conditions, which automatically introduce in the governing equations for the unit cell the information regarding the periodicity of the solid. As diffuse modes strongly depend on the microstructure of the unit cell, they are often called microscopic bifurcation modes. The results presented by Triantafyllidis and Maker (1985) were the first obtained for multilayered nonlinear composites applying this method. A relevant result of the study was the close relationship between onset of long-wavelength microscopic modes and localization of deformation of the homogenized continuum. For electroelastic nonlinear solids, the investigation of zero-speed waves has been used by Dorfmann and Ogden (2010) to explore surface instability, revealing that the critical loading parameter is crucially dependent on the magnitude of the applied electric field and on the magnitudes of the electromechanical coupling parameters in the constitutive equations. The microscopic bifurcation analysis for a periodic bilayered nonlinear dielectric solid under plane-strain conditions will be presented in Section 4.1.3.

• Loss of strong ellipticity of single phases and of the homogenized continuum (localized or macroscopic modes) 14

Periodic perturbations do not represent an effective way to describe a critical bifurcation mode corresponding to fields localized along a narrow band with homogeneous deformation and electric field outside the band. In the pure mechanical case localized modes are identified with shear bands (Hill and Hutchinson, 1975) which appear at loss of positive definiteness of the acoustic tensor, also indicated as loss of strong ellipticity5 . In the mechanics of nonlinear composites Geymonat et al. (1993) showed that the loss of strong ellipticity for the homogenized properties coincides with the longwavelength limit of microscopic diffuse modes, providing macroscopic bifurcation modes. The investigation of macroscopic bifurcation modes at critical loading parameter t = tMacro for soft multilayered dielectrics will be presented in Section 4.1.4. • Electric breakdown (EB).

When the electric field in a part of the system reaches a critical level, failure of the whole system may occur due to electric discharges. Similarly to Moscardo et al. (2008), the critical loading parameter tEB corresponds to the attainment of the critical value of electric field, EEB , in one phase.

4.1

Bifurcation analysis for multilayered composite soft dielectrics

In this Section the instability criteria outlined above are specialized to the case of an infinite, incompressible, multilayered body finitely deformed under plane-strain conditions. The solid consists of alternating identical and perfectly bonded incompressible layers of materials a and b with initial thicknesses h0a and h0b in the undeformed, stress-free configuration, which is used as the reference configuration. Note that the formulation presented in the following part of this Section can be generalized to any number of layers and is valid for incompressible phases whose response is described by a free-energy function W (I1 , I2 , I4 , I5 , I6 ) [see Section 2.4]. The results can be easily specialized to a specific material introducing the appropriate incremental moduli given by eqns. (29). 4.1.1 Fundamental pre-bifurcation path

0 Let us introduce a fixed cartesian coordinate system Ox1 x0 x0 with orthonormal basis 2 3 0 vectors ei (i = 1, 2, 3), and with x2 perpendicular to the direction of lamination in the reference configuration (see Fig. 3). A pure homogeneous plane-strain deformation in the x0 –x0 plane is considered, so that the deformation is completely described by 1 2

x1 = λx0 , 1
5

x2 = x0 /λ, x3 = x0 , 2 3

(42)

In the pure mechanical case, a theorem due to van Hove (1947), proved for elasto-plasticity but easily adapted in elasticity, states that failure of strong ellipticity is the critical condition for bifurcation, eqn. (39), for a homogeneous body, homogeneously deformed, subject uniquely to incremental displacements on the boundary.

15

0 where λ is the stretch along direction x1 . Homogeneity and perfect bonding between layers require that all phases ρ share the same longitudinal stretch, namely λρ = λ (ρ = a, b).

RVE x
0 2

a b
0 x1

n =e 2 0 m =e1

0

h0a h0b

D =D2e2

0

Figure 3: Schematic representation of the composite under investigation. The nominal electric displacement vector D0 is taken to be aligned to the x0 direction 2 and to be independent of x0 and x0 . This assumption is compatible, for instance, with 1 3 electrodes placed at the external boundary of the multilayered solid and parallel to the 0 direction of lamination. It follows from eqn. (5)3 that D2 is constant throughout each phase, so that 0 0 0 D1 = 0, D2 = 0, D3 = 0, (43) in the reference configuration, while in the current configuration D1 = D3 = 0,
0 D2 = D2 /λ.

(44)

0ρ 0 Moreover, continuity of the electric displacement vector implies D2 = D2 (ρ = a, b). A consequence of the above assumptions is that the free-energy function W introduced 0 0 ¯ in Section 2.4 depends only on λ and D2 , W = W (λ, D2 ), and that the only non-zero in-plane stress components are related through

S11 −

¯ S22 ∂W = , λ2 ∂λ

τ11 − τ22 = λ

¯ ∂W . ∂λ

(45)

Additionally, the normal component of the total stress is assumed to vanish along direction x0 throughout the body 2 S22 = 0, or τ22 = 0, (46)

so that continuity of tractions at the interface between the two phases requires
ρ S22 = 0 (ρ = a, b).

(47)

Notice that condition (46) is often met in soft dielectric devices where the material can deform and expand under the action of attractive forces between electrodes. Thus, eqns. (45) specialize to ¯ ¯ ∂W ∂W S11 = , τ11 = λ , (48) ∂λ ∂λ 16

while the nominal electric field is given by
0 E2 =

¯ ∂W , 0 ∂D2

0 0 E1 = E3 = 0.

(49)

In each phase the solution of the plane-strain incremental problem must satisfy the updated field equations [see eqns. (24)] Σ11,1 + Σ12,2 = 0, Σ21,1 + Σ22,2 = 0, ˆ ˆ D1,1 + D2,2 = 0, ˆ ˆ E1,2 − E2,1 = 0, (50)

together with incremental jump conditions at all interfaces between layers given by (26). 4.1.2 Loss of positive definiteness of the tangent electroelastic constitutive operator

For the considered principal equilibrium path loss of positive definiteness of the tangent electroelastic constitutive operator can be investigated using directly the constitutive equations (48) and (49). 0 Small perturbations of stretch λ and electric displacement D2 modify the current stress and electric field as ¯ ¯ ∂2W ˙ ∂2W ˙ 0 ˙ D , S11 = λ+ 0 ∂λ2 ∂λ ∂D2 2 ¯ ¯ ∂2W ˙ ∂2W ˙ 0 ˙0 E2 = λ+ D . 0 0 ∂λ ∂D2 ∂(D2 )2 2 (51)

The matrix associated with the linear system (51) is the Hessian of the free energy, which coincides with the tangent electroelastic stiffness at the current state. Eqn. (51) represents the specialization of eqn. (20), with the electroelastic moduli tensors (18) replaced by partial derivatives of the free energy with respect to the independent variables 0 λ and D2 . Under dead-load tractions/charges applied to the whole boundary, uniqueness of the incremental response is ensured by positive definiteness of the Hessian, namely ¯ ∂2W > 0, ∂λ2 ¯ ¯ ∂2W ∂2W > 0 ∂λ2 ∂(D2 )2 ¯ ∂2W 0 ∂λ ∂D2
2

.

(52)

Failure of the former inequality corresponds to a possible mechanical instability, while failure of the latter condition, which involves both stretch and electric displacement, detects an electromechanical instability. Recently Zhao and Suo (2007, 2009) and Zhao et al. (2007) used a similar bifurcation criterion to investigate electromechanical instabilities in homogeneous three-dimensional actuators subjected to loading conditions similar to those described in Section 4.1.1. Their analyses revealed that the critical bifurcation 0 0 point corresponds to a peak in the D2 –E2 curve. We remark that for composites the bifurcation conditions (52) may be critical either for a single phase or for the homogenized continuum.

17

4.1.3

Existence of diffuse-mode bifurcations (microscopic modes)

Here diffuse modes corresponding to inhomogeneous response of the composite with wavelengths (given by 2π/k1 , where k1 is the wavenumber) on the same order of the characteristic length of the heterogeneities are investigated applying the Bloch-Floquet quasi-periodicity conditions. Moreover, the limit case of long-wavelength limit (k1 → 0) is analytically studied, allowing investigation of the relation between microscopic and macroscopic instabilities. To solve the incremental boundary-value problem given by eqns. (50) in each phase ρ, we seek solutions in the form uρ (x1 , x2 ) = υ ρ (x2 ) exp[i k1 x1 ], ˆ Dρ (x1 , x2 ) = ∆ρ (x2 ) exp[i k1 x1 ]. The incompressibility constraint requires
′ i k1 υ1 + υ2 = 0,

p ρ (x1 , x2 ) = q ρ (x2 ) exp[i k1 x1 ], ˙ (53)

(54) (55)

whereas eqn. (50)3 implies

i k1 ∆1 + ∆′2 = 0,

where the dependency of υ and ∆ on x2 and ρ has been omitted for notational simplicity and ( · )′ = d( · )/dx2 . Substitution of eqns. (53)-(55) into eqns. (50)1,2 and (50)4 leads to the following system of ordinary differential equations in each phase ρ
2 ′′ −i k1 q + k1 (C1122 + C1221 − C1111 ) υ1 + C1212 υ1 + i k1 B112 ∆2 + B121 ∆′1 = 0, 2 ′ i q ′ + i k1 C2121 υ2 + (C1221 + C1122 − C2222 ) k1 υ1 + k1 (B121 − B222 ) ∆1 = 0, 2 B121 υ1 + (B112 + B121 − B222 ) k1 υ1 + A11 ∆′1 − i A22 k1 ∆2 = 0.
′′

(56)

Displacement continuity at the interface [eqn. (26)1] reduces to [[υ1 ]] = 0, [[υ2 ]] = 0, (57)

whereas traction continuity [eqn. (26)2 ] implies
′ [[B121 ∆1 + i k1 (C1221 + p) υ2 + C1212 υ1 ]] = 0,

[[B222 ∆2 + i k1 (C1122 − C2222 − p) υ1 − q]] = 0. and electrostatic interface conditions (26)3 and (26)4 are given by [[∆2 ]] = 0,
′ [[A11 ∆1 + B121 υ1 + i B121 k1 υ2 ]] = 0.

(58)

18

In each layer the general solution to eqns. (56) is found to be the sum of six linearly independent partial solutions y(x2 ) = Wa exp[k1 Za x2 ] a− y(x2 ) = Wb exp[k1 Zb x2 ] b y(x2 ) = Wa exp[k1 Za x2 ] a+ (0 < x2 < ha ), (ha < x2 < hb ), (hb < x2 < ha + hb ), (59)

ρ ρ ρ where Zρ = diag zρ , a− , a+ and b are vectors of unknown constants and yρ = [υ1 (υ1 )′ υ2 ∆ρ ∆ρ q ρ ]. Vector zρ and matrix Wρ contain the eigenvalues and eigenvectors of the 1 2 6×6 matrix Vρ with non zero entries

V12 = 1, V25 V31 V45 V54

2 k1 [B121 (B222 − B112 − B121 ) + A11 (C1122 + C1221 − C1111 )], d1 i k1 A11 i k1 (A11 B112 + A22 B121 ), V26 = − , = d1 d1 k2 = −i k1 , V41 = 1 [C1212 (B112 − B222 ) + B121 (C1111 − C1122 + C1212 − C1221 )], d1 i k1 i k1 B121 (B121 B112 + A22 C1212 ), V46 = , =− d1 d1 2 = −i k1 , V62 = i k1 (C1221 + C1122 − C2222 ), V63 = −k1 C2121 ,

V21 =

V65 = ik1 (B121 − B222 ),
2 B121

(60)

with d1 = − A11 C1212 . The quasi-periodic character of the solution along x2 is ensured imposing of the BlochFloquet type relation yρ (x2 + h) = yρ (x2 ) exp[ik2 h] (h = ha + hb ), (61)

with k2 lying in the unit cell of the reciprocal lattice, i.e. 0 ≤ k2 < 2π/h. The wavenumber k2 is frequently termed ‘Bloch parameter’ and sets the shape of modes along the transverse direction. Substituting eqns. (59) into eqns. (57)-(58) and (61) is found that a non-trivial solution to the problem exists when det (K − exp[ik2 h] I) = 0, with K = (Ga )−1 Gb exp[i Zb hb ](Gb )−1 Ga exp[i Za ha ], Gρ = Qρ Wρ , (63) Qρ being a 6×6 matrix with the following non-zero components Q11 = Q23 = Q55 = −Q46 = 1, Q34 = B121 , Q62 = B121 , Q63 = i k1 B121 , Q32 = C1212 , Q33 = i k1 (C1221 + p), Q45 = B222 , (62)

Q41 = i k1 (C1122 + C2222 − p) , Q64 = A11 . 19

It follows from eqn. (62) that the critical loading parameter at the onset of the bifurcation point (tMicro ) corresponds to the first occurrence of an eigenvalue Λ = exp[i k2 h] of unitary magnitude for the matrix K. Since det K = 1 and trK = trK−1 , the eigenvalues of K can be calculated from Λ6 − trK (Λ5 + Λ) + 1 trK2 − (trK)2 (Λ4 + Λ2 ) 2

1 − (trK)3 − 3trK(trK)2 + 2trK3 Λ3 + 1 = 0. 6

(64)

The critical eigenmode is determined observing that the constant vector a− appearing in eqns. (59) corresponds to the critical eigenvector of K, whereas a+ and b are obtained introducing eqns. (59) into eqns. (57)-(58) and (61), a+ = exp[ik2 h] exp[−k1 Za h] a− , b = exp[ik2 h] exp[−k1 Zb h](Gb )−1 Ga a− . The long-wavelength limit. The long wavelength limit k1 → 0 is now considered to establish a relation between microscopic and macroscopic instabilities. Here we extend the formulation of Triantafyllidis and Maker (1985) for the pure mechanical case and ¯ perform a matrix series expansion of K –defined in eqn. (63)– in powers of k1 = k1 h. Collection of the terms of the same order yields ¯2 ¯4 ¯6 trK = 6 + ξ1 k1 + ξ2 k1 + ξ3 k1 , ¯2 ¯4 ¯6 trK2 = 6 + 4 ξ1 k1 + (16 ξ2 + ξ4 )k1 + (64 ξ3 + ξ5 )k1 , ¯2 ¯4 ¯6 trK3 = 6 + 9 ξ1 k1 + (81 ξ2 + 6 ξ4)k1 + (729 ξ3 + ξ6 )k1 , where coefficients ξi (i = 1, . . . , 6) are given in Appendix B. Substituting eqns. (65) into eqn. (64) and recalling that Λ = exp[i k2 h] we obtain
3 2 −ξ1 + 36 ξ1 ξ2 −720 ξ3 + 3 ξ1ξ4 + 12 ξ5 −2 ξ6 + 3(ξ4 + 12 ξ2 −ξ1 )ψ 2 −6 ξ1 ψ 4 −6 ψ 6 = 0, (66)

(65)

with ψ = k2 /k1 . Therefore, a long wavelength bifurcation mode exists when a real solution for eqn. (66) is found. For the purely mechanical case, Geymonat et al. (1993) rigorously showed the equivalence between long wavelength bifurcation modes and localized, macroscopic instability. Here, in Section 5, the solutions of eqn. (66) will be compared with those for localized modes of the homogenized continuum, allowing us to determine a clear relationship between microscopic and macroscopic instabilities. 4.1.4 Loss of strong ellipticity of single phases and of the homogenized continuum (localized or macroscopic modes)

Bifurcation modes consisting of localized fields in a narrow band surrounded by homogeneous electroelastic deformation are now investigated. In general, for composites, this 20

instability may be critical either for each single phase or for the homogenized continuum. As previously recalled, for the latter case, it can be related to long wavelength diffuse modes. The analysis represents an extension of the formulation presented by Hill and Hutchinson (1975) for incompressible materials and can be used for all isotropic, elastic dielectrics that are deformed according to eqns. (42)-(43). The solution for the field equations (50) is sought in the form ˜ u = υ f (n · x), p = qf ′ (n · x), ˙ ˜ ˆ ˜ D = ∆f ′ (n · x), (67)

where f is a continuous and piecewise continuous differentiable function and n is the normal to the band. Incompressibility and eqn. (50)3 require that ˜ υ · n = 0, ˜ ∆ · n = 0, (68)

˜ ˜ so that for the plane-strain problem under consideration the amplitude vectors υ and ∆ are related through ˜ ˜ ∆ = ηυ, (69) η being a real parameter. To obtain the localization condition in a homogeneous phase of the composite, eqns. (28)2 , (30), (67), (68)1 and (69) are substituted in the incremental field equations (50)1,2,4 , providing four equations in q , η, υ1 and υ2 . Elimination of the first three unknowns yields ˜ ˜ ˜ the condition of band localization as Γ6 ν 6 + Γ4 ν 4 + Γ2 ν 2 + Γ0 = 0, (70)

where ν = n2 /n1 and we have assumed n1 = 0 and ν 2 A11 + A22 = 0. Coefficients Γi (i = 0, 2, 4, 6) depend on the incremental moduli as
2 Γ6 = B211 − A11 C1212 ,

Γ4 = −A22 C1212 − A11 (C1111 − 2C1122 − 2C2112 + C2222 ) − 2B211 (B112 + B211 − B222 ), Γ2 = −A11 C2121 − A22 (C1111 − 2C1122 − 2C2112 + C2222 ) + (B112 + B211 − B222 )2 , Γ0 = −A22 C2121 . Therefore, localization of the deformation into a band occurs when a real solution ν∗ of eqn. (70) exists. In this case the amplitude ratio η in eqn. (69) is given by η=
3 ν∗ (B112 + B211 − B222 ) − ν∗ B211 . 2 ν∗ A11 + A22

(71)

For the homogenized continuum, the procedure can be repeated employing the pertinent homogenized constitutive moduli.

21

4.1.5

Electric breakdown

Electric breakdown represents a possible failure mode that occurs when the current electric field reaches a critical value EEB in one of the phases of the composite. Since for the specific geometry and loading conditions considered here the only non-zero component of the electric field is that along direction x2 , electric breakdown takes place at a critical loading parameter tEB when ρ E2 = EEB , either in phase a or b.

5
5.1

Results
Material properties and constitutive laws

In the following sections results are presented for dielectric bodies whose response is captured using the extended neo-Hookean free energy given by eqn. (12). For this specific material the constitutive equations (9)–(10) specialize to 1 S = −p F−T + µF + FD0 ⊗ D0 , ǫ in the lagrangian formulation and 1 τ = −p I + µB + D ⊗ D, ǫ 1 E= D ǫ (73) 1 E0 = CD0 ǫ (72)

in the eulerian formulation6 , a result fully consistent with the general equations (45) and (49); eqn. (73)2 clearly shows the linear proportionality between D and E. The incremental moduli are obtained substituting the free energy (12) into eqns. (18), yielding7 1 0 0 0 (74) CiJkL = µδik δJL + δik DJ DL , ǫ 1 1 0 0 0 A0 N = CM N , BiJM = (FiM DJ + FiS DS δJM ), (75) M ǫ ǫ
We remark that in (72)1 and (73)1 the hydrostatic pressure p is calculated imposing the boundary conditions. Therefore, p can be replaced by an equivalent expression, such as p = p + E · E/2, so that ˜ eqn. (73)1 can be reformulated as 1 1 τ = −˜ I + µB + D ⊗ D − D · D, p ǫ 2ǫ an expression equivalent to that reported by Zhao et al. (2007) in their eqn. (26) and obtained adapting the compressible theory of nonlinear dielectrics to the incompressible case. 7 Eqns. (74)-(75) are obtained using the derivatives of the invariants I1 , I5 with respect to F and D0 ∂I1 = 2F, ∂F ∂I5 = 2D ⊗ D0 , ∂F ∂I5 = 2CD0 , ∂D0 ∂I5 ∂F∂D0
0 0 = 2δJK FiR DR + 2DJ FiK . iJK 6

22

while the use of eqns. (29) provides the updated constitutive tensors 1 Ciqkp = µδik Bpq + δik Dp Dq , ǫ 1 Biqa = (δia Dq + Di δqa ), ǫ 1 Aab = δab . ǫ

For the geometry and loading conditions described in Sections 4.1.1, it follows from eqn. (72) that the only non-zero components of the total first Piola-Kirchhoff stress and nominal electric field are p S11 = µλ − , λ S22 =
0 µ (D2 )2 + − pλ, λ λǫ

S33 = µ − p,

0 E2 =

0 D2 , λ2 ǫ

(76)

whereas the only non-vanishing in-plane entries of the updated incremental constitutive tensors are C1111 = C2121 = µλ2 , 1 B121 = B211 = D2 , ǫ C2222 = C1212 = µ 2 B222 = D2 , ǫ 1 1 2 + D2 , 2 λ ǫ

1 A11 = A22 = . ǫ

Note that the hydrostatic pressure in eqns. (76)1−3 can be determined using eqn. (46), yielding ¯ 1 + D2 p=µ , λ2 so that the longitudinal stress S11 is given by S11 = µ λ − ¯ 1 + D2 λ3 , (77)

0 √ ¯ with D = D2 / µ ǫ. To illustrate the effect of the material parameters µ and ǫ on the electromechanical response of a homogeneous layer, we focus on a sample subjected to a given longitudinal nominal stress S11 . In this case the stretch λ can be obtained from eqn. (77) as a function 0 of the electric displacement D2 and then can be used to calculate the electric field from eqn. (76)4 . In Fig. 4 the responses of two specimens whose response is characterized by two sets of parameters, {µ1 , ǫ1 } and {µ2 , ǫ2 } with µ2 /µ1 = ǫ2 /ǫ1 = 10, are reported. Focusing on 0 the case for {µ1 , ǫ1 }, Fig. 4a shows that the λ − D2 curve can be decomposed into three parts: an initial quadratic response is followed by a linear behavior for moderate values 0 0 √ 0 of D2 , while λ ∼ D2 for D2 / µ1 ǫ1 > 4. Moreover, we observe that the response for the specimen characterized by {µ2 , ǫ2 } is qualitatively similar to that for {µ1 , ǫ1 }. However, for the material characterized by {µ2 , ǫ2 } an electric displacement hundred times larger than that for {µ1 , ǫ1 } must be applied to reach a given stretch λ. In Fig. 4b the evolution of the dimensionless nominal electric field is reported as a function of the applied electric 0 displacement. As expected from eqn. (76)4 , at a given electric displacement D2 a high electric field is reached in the material with lower permittivity. Interestingly, for applied compressive longitudinal tractions (see curve for S11 /µ1 = −1) the nominal electric field

23

3.5 3 2.5 2

1.4

(a)
S11/m1=1 S11/m1=-1
0

1.2

(b)
S11/m1=-1 m1,e1

m1,e1

1 0.8

l
1.5 1 0.5 0 1 2 3 4
0

S11=0

E2 e1/m1 m2,e2

0.6 0.4 0.2

S11/m1=1 S11=0

m2,e2

5

6

7

8

9

10

0

1

2

3

4
0

5

6

7

8

9

10

D2/ m1e1

D2/ m1e1

Figure 4: Electromechanical behavior of a homogeneous layer deformed under plane strain by the
0 application of a transverse electric displacement field D2 and a longitudinal nominal traction S11 . The material response is captured using the extended neo-Hookean free energy (12) and µ2 /µ1 = ǫ2 /ǫ1 = 10.

curve displays a maximum. The peak corresponds to loss of positive definiteness of the tangential electroelastic constitutive operator and occurs at failure of condition (52)2 . However, in a longitudinally compressed layer a buckling instability may occur earlier along the loading path, preventing the structure to reach the peak.

5.2

Homogenized response

Here the general results obtained in Section 3 for rank-one layered dielectrics are specialized to the extended neo-Hookean free energy (12). Initially we consider the case of layers with an arbitrary inclination with respect to the external electric displacement (n0 = e2 , Fig. 2). Then the results are further specialized to the specific geometry and loading conditions considered in this work, characterized by layers perpendicular to the electric displacement field (n0 = e2 , Fig. 3). For the case of a laminated body with layers arbitrarily oriented with respect to the electric displacement field and whose response is described by the free energy (12), parameters α and β introduced in eqn. (36) are given by α= µb − µa Fav n0 · Fav m0 , ca µb + cb µa Fav m0 · Fav m0

µa ǫa − µb ǫb Fav D0 av · Fav m0 µb − µa β= a b + a b D0 av · m0 . (c µ + cb µa ) (ca ǫa + cb ǫb ) Fav m0 · Fav m0 c µ + cb µ a

(78)

The jump in hydrostatic pressure across each interface is obtained by multiplying the traction continuity condition ([[S]]n0 = 0) with the vector (Fav )−T n0 , yielding pa − pb = ǫb − ǫa 0 av 0 2 1 (D · n ) + µa − µb . b ǫa av )−T n0 · (Fav )−T n0 ǫ (F 24 (79)

Results given in eqns. (78) and (79) are consistent with those obtained by deBotton (2005) for neo-Hookean rank-one laminates finitely deformed by the application of only mechanical loadings. For the geometry and loading conditions described in Sect. 4.1.1 (i.e. n0 = e2 and m0 = e1 ), the only non-vanishing components of the homogenized stress and electric field obtained using eqns. (32)2,3 are
av S11 = µav λ −

µav ¯ 1 + (D av )2 , 3 λ ca cb + ǫa ǫb

0 E2 av =

¯ µav D av , ǫHav λ2
0 D2

(80)

where µav = ca µa + cb µb , ǫHav =

−1

,

¯ D av =

µav ǫHav

.

(81)

The analysis of macroscopic modes for the composite requires the macroscopic incremental constitutive tensors Cav , Bav , and Aav . It follows from the introduction of (78) in eqns. (34) and (35) and the use of (32), (20) and (29) that their only in-plane non-zero components are
av C1111 = λ2 µav , av C1221 av C2112

=

¯ ¯ µHav (D av )2 ǫHav µav µav (D av )2 µav av + , C2222 = 2 + , λ2 ǫav λ2 λ λ2 ¯ µHav − µav (D av )2 µav ǫHav − ǫav av av av = + , C2121 = C1221 + C1111 , 2 2 av λ λ ǫ
av C1212 =

(82)

av av B211 = B121 =

¯ D av

µav ǫHav , ǫav λ

av B222 = 2

¯ D av λ

µav , ǫHav

Aav = 11 ca cb + b µa µ

1 , ǫav
−1

Aav = 22

1 ǫHav

,

(83)

where ǫav = ca ǫa + cb ǫb and µHav =

.

(84)

On the other hand, the investigation of loss of positive definiteness of the tangent electroelastic operator requires the specialization of the macroscopic free energy (36) to the specific geometry and loading conditions considered here, leading to ¯ µav 2 (D av )2 + 1 0 ¯ W av = W av (λ, D2 ) = λ + −2 . 2 λ2 (85)

To illustrate the effect of phases volume fraction and contrast in material properties on the overall electromechanical response, we investigate the behavior of a multilayered dielectric deformed by the application of an electric displacement field perpendicular to av the layers, while leaving the body free to expand laterally (i.e. S11 = 0). The analysis is fully characterized by the dimensionless material parameters m = µa /µb and r = ǫa /ǫb . (86)

25

2

2

(a) m=r=10
1.8

c =0.1
1.8

a

(b) m=r=100
c =0.3 c =0.5
a a

c =0
1.6

a

c =0
1.6

a

l
1.4

l
ca=0.7
1.4

c =0.1

a

ca=0.3
1.2 1.2

c =1
1 0 1 2 3 4 5 6 7 8 9 10 1 0 1 2 3 4 5 6
b

a

c =0.5 a c =0.7
7 8 9 10

a

D

b

D

0 ¯ Figure 5: Overall Db –λ response (with D2 / µb ǫb ) for a periodic bilayered dielectric when an electric

displacement field perpendicular to the layers is applied while leaving the body free to expand laterally av (i.e. S11 = 0).

¯ Fig. 5 shows the electromechanical response of the multilayered structure in the D b –λ space with µav 0 ¯ ¯ D b = D2 / µb ǫb = D av . (87) ca µ a b µb + c mr a We note that in both Figs. 5a and 5b for c = 0 the curve corresponds to that for {µ1 , ǫ1 } in Fig. 4a, as in this case the solid is composed uniquely of the soft phase. The effect of phases volume fraction on the material response is investigated for two sets of material parameters, m = r = 10 in Fig. 5a and m = r = 100 in Fig. 5b. For the geometry considered in this study an increase of parameters m and r is found to strongly reduce the longitudinal elongation.

5.3

Nonlinear fundamental electromechanical paths

0 Although a generic state of deformation described by the pair (λ, D2 ) can be reached following multiple loading histories, here we focus on two fundamental nonlinear electromechanical paths: av ˜av • Path A: a uniaxial macroscopic stress S11 = S11 is applied and kept constant while increasing the external electric displacement field (Fig. 6a).

˜ • Path B: a stretch λ = λ is applied and kept constant while increasing the external electric displacement field (Fig. 6b). Here these two nonlinear pre-bifurcation paths are studied assuming a material response defined by the free energy (12), so that stress and electrical displacement fields in each phase of the multilayered body are given by eqns. (76). 26

x

0 2 0 x1

a b

Path A

1) undeformed conf. (l=1)

Path B

~b S11 ~a S11 ~ av ~ a ~b 2) pre-stressed conf. ( S11=caS11+cbS11) ~ 2) pre-stretched conf. (l=l)

electrodes

~ av ~ a ~b 3) electrically actuated conf. ( S11=caS11+cbS11 ) (a)

~ 3) electrically actuated conf. (l=l) (b)

Figure 6: Schematic representation of the two loading paths considered in this study: Path A (a) and
Path B (b).

av ˜av Path A. A constant longitudinal stress S11 = S11 is applied to the multilayered structure, so that ˜av ˜a ˜b S11 = ca S11 + cb S11 , (88) ˜ ˜ where S a and S b are obtained from eqn. (77) as 11 11

¯ 1 + (D b )2 /m r ˜a S11 = µa λ − λ3

,

¯ 1 + (D b )2 ˜b S11 = µb λ − λ3

.

(89)

Substitution of eqns. (89) in eqn. (88) provides an equation that can be solved for the ˜a ˜b unknown stretch λ. The longitudinal stress in each phase, S11 and S11 that in general change along the deformation path, can then be obtained by introducing λ into eqns. (89). Moreover, the nominal electric fields in the two phases follow from eqn. (76)4 and can be written in dimensionless form as
0 ¯ E a = E2 a

¯ ǫa Db √ = 2 , µa λ mr

0 ¯ E b = E2 b

¯ ǫb Db = 2. µb λ

Path B. Devices based on soft dielectrics are usually prestretched before the application of the electric field, since it has been observed experimentally that the prestretch increases 27

considerably the performance of the actuator, mainly increasing the maximum electric ˜ field before breakdown (Pelrine et al., 2000). Here a longitudinal prestretch λ = λ (with ˜ λ > 1) is applied to the multilayer and, since S22 = 0 throughout the solid, the stress component S11 associated with the prestretch in the two phases is given by 1 pre,a ˜ S11 = µa λ − ˜ λ3 , 1 pre,b ˜ S11 = µb λ − ˜ λ3 .

0 Subsequently an increasing electric displacement D2 is applied (note that D2 ≡ D2 as the current stretch does not vary due to incompressibility), so that the longitudinal stresses change as D2 D2 pre,a pre,b a b S11 = S11 − 2 , S11 = S11 − 2 . ˜ ˜ λ3 ǫa λ3 ǫb a b The electric excitation introduces a compressive stress in the bilayer, so that S11 and S11 s0a s0b drop to zero when D2 reaches the intensities D2 and D2 , respectively, where s0a D2 =

˜ µa ǫa (λ4 − 1),

s0b D2 =

˜ µb ǫb (λ4 − 1).

s0a s0b Note that when D2 > min{D2 , D2 } part of the system is compressed and this may lead to a buckling instability. Although traditionally buckling instabilities have been viewed as an inconvenience, recently they have been exploited to design actuators able to sustain transverse large displacements (Carpi et al., 2008c; Koo et al., 2008).

5.4

Analysis of instabilities

Loss of uniqueness of the incremental response is now investigated for biphasic periodic multilayers whose response is described by the extended neo-Hookean free energy (12), so that the formulation presented in Section 4 specializes as follows: • Loss of positive definiteness of the tangent electroelastic constitutive operator (PD) This criterion is applicable to structures with only tractions/surface charges applied to their boundary, so that it represents a critical condition only for multilayers loaded following type-A loading paths, since along type-B loading paths displacements are prescribed on portion of the boundary. For the effective free energy (85), conditions (52) reduce to µav ¯ 3 + 3(D av )2 + λ4 > 0, λ4 (µav )2 ¯ 3 − (D av )2 + λ4 > 0. λ6

(90)

It is clear that inequality (90)1 is always satisfied, whereas (90)2 fails when λ4 ≤ ¯ (D av )2 − 3, with the equality setting the critical limit represented by the equation ¯ λPD = [(D av )2 − 3]1/4 , 28 (91)

that can be rewritten, using (87), as 1/4 ca µ a cµ +  mr (D b )2 − 3 . ¯ =  av µ 
b b

λPD

(92)

Introduction of the critical condition (92) into (80)1 yields ¯ av EPD = ¯ Db 3µav ¯ (D b )2 − ca µ a cb µ b + mr ,

(93)

while substituting eqn. (91) into eqn. (80)1 we obtain 4  ¯ S11,PD = − 3 = −4  λPD  cb µ b + −3/4 ca µ a mr (D b )2 − 3 ¯ , 

µav

(94)

¯ ˜av where S11 = S11 /µav . Eqn. (94) shows that the tangent moduli of the homogenized continuum looses positive definiteness only when longitudinal homogenized compressive stresses are applied. Moreover, the results reveal that the critical condition ¯ ¯ (91) corresponds to a peak in the D b –E av curve. Depending on the way the electric excitation is applied to the system, the peak may or may not correspond to loss of stability. When the charge on the electrodes is controlled along the loading path, it is possible to go beyond the peak maintaining the homogeneous quasi-static response. On the other hand, when the voltage is controlled, and as a consequence ¯ the value of the electric field, the system cannot sustain an electric field E av higher than that at the peak. Finally, we remark that in all our calculations the failure of eqns. (90) has been monitored both in each single phase and in the homogenized continuum. For the considered geometry and loading conditions, we have found that loss of positive definiteness of the electroelastic tangent operator occurs earlier along the loading path for the homogenized continuum than for each single phase. • Existence of diffuse bifurcation modes: microscopic modes (Micro)

Diffuse-mode bifurcations are investigated numerically as described in Section 4.1.3. The only simplification introduced in the analysis by the choice of the extended neoρ Hookean free energy is that the eigenvalues zi of matrices Vρ (ρ = a, b) satisfy the characteristic equation
2 k1 − (z ρ )2 2 2 k1 λ4 − (z ρ )2 = 0,

(95)

providing six real values for each phase. 29

In the long-wavelength limit (k1 → 0) the critical stretch is obtained solving eqn. (66), yielding 1/4 ǫHav µHav ¯ λ∞ = 1 − av + (D av )2 1 − av , (96) Micro µ ǫ ¯ which can be easily reformulated in terms of D b using eqn. (87). Note that when 0 D = 0, eqn. (96) reduces to the well-known expression reported by Triantafyllidis and Maker (1985), namely λ∞ Micro|D0 =0 µHav = 1 − av µ
1/4

.

(97)

• Loss of strong ellipticity of single phases and of the homogenized continuum: macroscopic modes (Macro) While band localization within each layer for the considered material model is excluded [eqn. (70) does not provide any real solution ν], for the homogenized continuum substitution of the incremental moduli (83)-(84) in the same equation leads to the same critical loading parameter reported in eqn. (96), so that λMacro = λ∞ . Therefore, a remarkable result of this analysis is that along elecMicro tromechanical loading paths localization is found to correspond to diffuse modes of infinite wavelengths. Similar results were found by Triantafyllidis and Maker (1985) and Geymonat et al. (1993) for the purely mechanical case. • Electric breakdown (EB)

For the geometry and loading conditions described in Sect. 4.1.1, eqn. (73)2 and results of Fig. 4 show that the electric field is significantly higher in the phase with lower permittivity. Therefore, taking b as the phase characterized by shear modulus and permittivity lower than those of phase a, electric breakdown may occur first in phase b. Indicating by EEB the critical value of the electric field, introduction of eqn. (44)3 into eqn. (73)2 leads to the critical condition λEB = ¯ µb D b , ǫb EEB (98)

¯ ¯ ¯ in the D b –λ plane. Alternatively, in the D b –E av space it can be expressed as
b ¯ av = ǫ E µb 2 EEB 1 av ¯b. D µ a a c µ cb µ b + mr

(99)

Here we investigate electrical breakdown using a set of parameters typical for silicones (EEB = 350 MV/m, ǫb = 3, µb = 100 kPa, see Kornbluh and Pelrine, 2008). r

30

For the considered multilayered structures the critical curves corresponding to the four instability criteria investigated in this paper are reported in Figs. 7–11 with continuous lines; black lines correspond to loss of positive definiteness of the tangent electroelastic constitutive operator (PD), blue lines correspond to microscopic modes, while red and brown lines denote macroscopic mode and electric breakdown limits, respectively. Moreover, in the figures loading paths of both type A and type B are reported to visualize the critical point along them. ¯ In Figs. 7, 9 and 10 the critical conditions are sketched in the D b –λ space, while Figs. b ¯ av ¯ 8 and 11 show the results in the D –E diagram where, using (80)2 ,
0 ¯ E av = E2 av

¯ D av ǫHav = 2 . µav λ

(100)

The background color in Figs. 7 and 10 identifies the dimensionless longitudinal stress av ¯ S11 = S11 /µav . In Figs. 8, 9 and 11 the gray shaded areas denote the region of allowable states for type-A loading paths, while for type-B loading paths the material is stable both in the gray and in the green areas. The figures clearly highlight the important effect on the failure modes of both the volume fraction of phases and the contrast in material properties. All these effects are investigated separately below. Effect of phases volume fraction. To explore the influence of phases volume fraction on the failure modes of the bilayer, for a contrast in material properties set by m = r = 10 we investigate four different values of ca , namely ca = 0.05, 0.2, 0.4, 0.6, and the results are reported in Figs. 7-8. Focusing on the four instability criteria, we observe that • for a very low value of ca (Fig. 7a, where ca = 0.05), microscopic instability represents the limit of the region of allowable states at low values of the applied ¯ ¯ ¯ electric displacement field D b (D b < 1 in the inset in Fig. 7a), while for higher D b macroscopic bifurcation and loss of positive definiteness of the tangent electroelastic constitutive operator (PD) are critical; • for moderate values of the volume fraction ca (Fig. 7b, where ca = 0.2), only bandlocalization (Macro) and PD criteria set the limit for instability over the considered ¯ range of D b ; • for higher values of ca (Fig. 7c,d), macroscopic instability is again dominant for low ¯ ¯ ¯ D b (approximately D b < 5), while for D b > 5, microscopic instability represents the critical bifurcation; ¯ • electric breakdown is never critical for the considered range of D b and ca . To illustrate the transition between macroscopic and microscopic instabilities occurring for high ca we focus on the case of Fig. 7c (m = r = 10 and ca = 0.6). Fig. 9a shows 31

m=10, r=10
2

(a) ca=0.05
1.8

PD Macro

0 -20

S11

2

(b) c =0.2
1.8

a

0

S11

Macro l=1.7
~

l=1.7
1.6

~

-10

-40

S11=-0.7 S11=0
0.88

EB

-60 -80

1.6

-20

S11=0 S11=-0.5

EB
-30

l

1.4

l 1.4
-100

l
-120 1.2

-40

1.2

l=1.1
1

~

PD l=1.1
~

0.82

Micro
0.76

-140

-50

Macro
0.7

D
1

b

-160
1.5

1 -60 0.8 0
2

0.8 0 2 4

0

0.5

-180 6 8 10

-70 2 4 6 8 10

D
2

b

D
0

b

(c) ca=0.4
1.8

S11

Micro
-5
1.8

(d) c =0.6
l=1.7
1.6
~

a

S11
0

l=1.7
1.6

~

Micro

-2 -4

-10

EB S11=1 S11=0 PD l=1.1
~

-6 -8 -10

l

1.4

S11=0.5

S11=-0.3 Macro EB

-15

l

1.4

1.2

-20

1.2

l=1.1
1

~

PD
-25
1

-12

Macro

-14 -16

-30 0.8 0 2 4
0.8

-18
0 2 4 6 8 10

D

b

6

8

10

D

b

Figure 7: Bifurcation of a bilayered dielectric composite subjected to a transverse electric displacement field: critical conditions in the longitudinal stretch (λ)–dimensionless electric displace0 ¯ ment (Db = D2 / µb ǫb ) space for m = r = 10. Four phase volume fractions ca are considered, a c = 0.05, 0.2, 0.4, 0.6. Solid lines correspond to Micro (blue), Macro (red), PD (black) and electric ¯ ˜av breakdown (EB, brown) instability criteria. Background color corresponds to S11 = S11 /µav . Shortdashed (green) curves denote type-A loading paths, while long-dashed (magenta) curves denote type-B loading paths. Critical states are indicated with circular markers.

both the evolution of λ (on the left vertical axis) and of the critical wavenumbers k1,cr ¯ and k2,cr (on the right vertical axis) as a function of D b . At low values of the applied ¯ b < 5) macroscopic instability prevails. The electric displacement field (approximately D results clearly show that at this stage the wavenumber k1,cr of the diffuse mode tend to zero, consistently with the long wavelength limit analysis. Moreover, the macroscopic instability analysis detects the formation of a band of localized fields orthogonal to the layers (n = e1 ). Thus, since η = 0 [eqn. (71)], only the deformation localizes into the band, while the electric field is continuous. ¯ Interestingly, for larger values of the applied electric displacement field (D b > 5) the microscopic instability analysis detects diffuse modes characterized by a critical Bloch

32

m=10, r=10
2 1.8 1.6 1.4 1.2 2

PD

Macro

PD
1.8 1.6

EB

l=1.7 l=1.1 S11=-0.7 S11=0
~

~

1.4 1.2

Macro S11=-0.5 S11=0

Eav

1 0.8 0.6 0.4 0.2 0 0

E

av

1 0.8

l=1.7
0.6 0.4

~

(a) c =0.05
2 4 6 8 10

a

0.2 0 0 2

l=1.1
4 6

~

(b) ca=0.2
8 10

D
2 1.8 1.6 1.4 1.2

b

D
2 1.8

b

PD

EB

EB
1.6 1.4

PD

Macro Micro S11=0.5 l=1.7
~

1.2

Eav

1 0.8 0.6 0.4 0.2 0 0 2

Eav
S11=-0.3

Macro l=1.1
~

1

Micro S11=0

0.8 0.6 0.4
~

S11=1 l=1.7

l=1.1
4 6

~

(c) c =0.4
8 10

a

0.2 0 0 2 4 6

(d) ca=0.6
8 10

D

b

D

b

Figure 8: Bifurcation of a bilayered dielectric composite subjected to a transverse electric displace0 ¯ ment field: critical conditions in the dimensionless average electric field (E av = E2 av ǫHav /µav )– 0 ¯ dimensionless electric displacement (Db = D2 / µb ǫb ) space for m = r = 10. Four phase volume a a fractions c are considered, c = 0.05, 0.2, 0.4, 0.6. Solid lines correspond to Micro (blue), Macro (red), PD (black) and electric breakdown (EB, brown) instability criteria. Short-dashed (green) curves denote type-A loading paths, while long-dashed (magenta) curves denote type-B loading paths. Critical states are indicated with circular markers while the gray (gray-green) shaded area highlights the region of stable type-A (type-B) homogeneous paths. In (c) and (d) the two regions coincide.

parameter k2,cr different from zero. Therefore, phases a and b at bifurcation deform ¯ out-of-phase, as shown in Fig. 9b and 9c, where the bifurcation modes at D b = 5 and ¯ 10 are reported. Differently, when microscopic instability is critical at low values of D b b ¯ (D < 1, see Fig. 7a), the analysis predicts critical diffuse modes with a vanishing Bloch cr parameter (k2 = 0), so that both phases a and b at bifurcation deform in-phase. This result is consistent with that obtained in the pure mechanical case by Triantafyllidis and Maker (1985), where only in-phase bifurcation modes were found for neo-Hookean 33

2

7

(a) c =0.6, m=10, r=10
1.8

a

(b)
Db=5, lcr=1.18, k1,cr=3.25, k2,cr=3.55
2h

Micro

k2,cr
Macro

6

1.6

5
h

l

4 1.4

k1,cr

k1,k2
3

0 0 h 2h 3h 4h 5h

1.2
2h

(c)
Db=10, lcr=1.83, k1,cr=3.60, k2,cr=5.63

2 1
h

(b)
0.8 0 5

(c)

1

10

0 15

0 0 h 2h 3h 4h 5h

D

b

Figure 9: Bifurcation modes of a bilayered dielectric composite subjected to a transverse electric

¯ ¯ displacement field. (a): critical dimensionless wavenumbers k1,cr = k1,cr h and k2,cr = k2,cr h for ca = 0.6 ¯ b = 5 and 10. and m = r = 10. (b), (c): graphical representations of modes for D

materials. Turning now the attention on the loading paths introduced in Section 5.3, we observe that for type-A loading paths the equilibrium is stable when specimens are subjected ¯av to a tensile longitudinal stress (S11 > 0) for ca < 0.5. On the other hand, for type-B ˜ loading paths we note that depending on the applied prestretch λ and on ca , the critical mode may correspond either to a microscopic or a macroscopic instability. Finally, it is ˜ interesting to remark that an increase of the prestretch λ not only leads to a change of the bifurcation mode from a localized band to a diffuse mode, but also to an expansion of the region where the structure is stable. Effect of contrast in material properties. Although composite dielectrics characterized by a wide range of contrast in material properties have been fabricated (Zhang et al., 2002; Huang et al., 2004), here the effect of contrast in material properties on instabilities is explored focusing on m = r = 10 and m = r = 100. For the same volume fractions as for m = r = 10, results pertaining to the case m = r = 100 are reported in Figs. 10-11, revealing that • loss of positive definiteness of the tangent electroelastic constitutive operator (PD) does not represent a limiting condition for the region of allowable states in the ¯ considered range of D b and ca ; • electric breakdown represents a critical condition at large values of ca ;

¯ • independently of the volume fraction of phase a, band localization is critical for D b b ¯ approximately lower than 3, while for larger values of D microscopic instabilities become dominant. 34

m=100, r=100
2

(a) ca=0.05
1.8

Micro

0

S11

2

(b) ca=0.2
1.8

Micro l=1.8
~

S11
0 -2

l=1.8 Macro

~

-10

Macro
-20

1.6

S11=1

1.6

EB S11=1

-4 -6 -8

l
1.4

l
S11=0 EB
1.2
~

-30

1.4

PD

-40

S11=0.5
1.2

PD
-10 -12

l=1.1
1
-50

l=1.1
1 -14

~

0

5

D
2

b

10

15

0

5

D S11
2

b

10

15

(c) c =0.4
l=1.8
1.8
~

a

Micro

1 1.8 0

(d) ca=0.6
l=1.8
~

Micro

1.5 1

S11

EB
1.6 -1 1.6

0.5 0

l
1.4

S11=1 EB S11=0.5 Macro

-2 -3 -4

l
1.4

S11=1

-0.5 -1

Macro
1.2
~

1.2

S11=0.5 l=1.1

-1.5 -2

l=1.1
1 -5 -6 0 5 10 15 0 1

~

-2.5 -3 5 10 15

Db

Db

Figure 10: Bifurcation of a bilayered dielectric composite subjected to a transverse electric displacement field: critical conditions in the longitudinal stretch (λ)–dimensionless electric displace0 ¯ ment (Db = D2 / µb ǫb ) space for m = r = 100. Four phase volume fractions ca are considered, ca = 0.05, 0.2, 0.4, 0.6. Solid lines correspond to Micro (blue), Macro (red), PD (black) and electric ¯ ˜av breakdown (EB, brown) instability criteria. Background color corresponds to S11 = S11 /µav . Shortdashed (green) curves denote type-A loading paths, while long-dashed (magenta) curves denote type-B loading paths. Critical states are indicated with circular markers.

When tractions are controlled (type-A loading paths), an increase of m and r has a dramatic effect on the stability of the multilayer. Focusing on ca = 0.05 and 0.2, Figs. ¯ 7a, b show that for m = r = 10 the solid loaded by a longitudinal stress S11 ≥ 0 can deform homogeneously up to a large longitudinal stretch λ, whereas for m = r = 100 (Figs. 10a, b) the critical point is reached at a much smaller longitudinal stretch. For type-B loading paths the situation is completely different, since an increase of parameters m and r allows the application of a slightly higher electric displacement prior to bifurcation, specially for weakly prestretched composites with a low volume fractions ca .

35

m=100, r=100
1 1

Macro
0.9 0.8 0.7 0.6

Micro S11=0 l=1.8 S11=1
~

0.9

Macro
0.8

EB
0.7 0.6

Eav

0.5 0.4 0.3 0.2 0.1 0 1 0.9 0.8 0.7 0.6 0 5
~

Eav
l=1.1

Micro
0.5 0.4 0.3 0.2

l=1.8 S11=0.5 S11=1

~

l=1.1

~

(a) ca=0.05 D
b

0.1 0 0 5

(b) ca=0.2 Db
10 15

10

15

1 0.9

EB

Macro
0.8 0.7 0.6

EB Macro

Eav

0.5 0.4 0.3 0.2 0.1 0 0
~

Eav 0.5
Micro S11=1 S11=0.5
0.4 0.3 0.2 0.1 0 0 5
~ ~

Micro

l=1.8 l=1.1

~

l=1.8

S11=1 S11=0.5

(c) ca=0.4
10 15

l=1.1

(d) ca=0.6 Db
10 15

5

D

b

0 ¯ ment field: critical conditions in the dimensionless average electric field (E av = E2 av ǫHav /µav )– b 0 b ǫb ) space for m = r = 100. Four phase volume ¯ dimensionless electric displacement (D = D2 / µ fractions ca are considered, ca = 0.05, 0.2, 0.4, 0.6. Solid lines correspond to Micro (blue), Macro (red), PD (black) and electric breakdown (EB, brown) instability criteria. Short-dashed (green) curves denote type-A loading paths, while long-dashed (magenta) curves denote type-B loading paths. Critical states are indicated with circular markers, while the gray shaded area highlights the region of stable homogeneous paths.

Figure 11: Bifurcation of a bilayered dielectric composite subjected to a transverse electric displace-

We conclude noting that an increase of m and r reduces the occurrence of bandlocalization and PD instabilities, promoting microscopic bifurcations and electric breakdown failure modes.

36

6

Conclusions

Soft composite dielectrics show great potential in the design of smart devices based on electrosensitive polymers, since their use allows the application of lower voltages as a result of the increased overall dielectric constant. A critical issue related to their development is the prediction of both global and local instabilities that may occur at the macro or at the micro-scale, respectively. Instabilities can be investigated analytically extending to soft dielectrics the well-known theory of bifurcation and stability of nonlinear elastic solids. In this paper a detailed analysis of instabilities of layered soft dielectric composites under plane deformations has been developed. Four instability criteria have been formulated, namely: i.) loss of positive definiteness of the tangent electroelastic constitutive operator of the homogenized continuum (PD); ii.) existence of diffuse modes of bifurcation (microscopic instability); iii.) existence of localized solutions of the homogenized solid (microscopic instability, corresponding to the loss of ellipticity of the homogenized continuum); iv.) electric breakdown. The critical conditions associated with these four criteria have been obtained for a periodic bilayered composite obeying an extended form of the neo-Hookean free-energy function and subjected to a transverse electric displacement field. The results provide for the first time an analysis of failure mechanisms which may occur in a heterogeneous dielectric under external electrical excitation. The effect of phases volume fraction and contrast in material properties on instabilities has been investigated, revealing that for low volume fractions of the stiffer phase –a relevant case in the application–, an increase of the contrast in material properties of the two phases promotes a shift from PD to microscopic instabilities. Moreover, while for constituents not strongly dissimilar band-localization is critical, for high-contrast phases microscopic instabilities and electric breakdown are the dominant failure modes for composite soft dielectrics.

Acknowledgements
K.B. acknowledges Harvard University for financial support and thanks Robert McMeeking for his generous hospitality at UCSB. M.G. gratefully acknowledges the financial support of PRIN grant no. 2007YZ3B24 “Multi-scale Problems with Complex Interactions in Structural Engineering”, financed by Italian Ministry of University and Research.

References
[1] Y. Bar-Cohen (Ed). Electroactive Polymer (EAP) Actuators as Artificial Muscles. SPIE Press, Bellingham, Wa., 2001.

37

[2] K. Bertoldi, M.C. Boyce, S. Deschanel, S.M. Prange and T. Mullin. Mechanics of deformation-triggered pattern transformations and superelastic behavior in periodic elastomeric structures. J. Mech. Phys. Solids 56, 2642-2668, 2008. [3] D. Bigoni and M. Gei. Bifurcation of a coated, elastic cylinder. Int. J. Solids Structures 38, 5117-5148, 2001. [4] M.A. Biot. Mechanics of Incremental Deformations. Wiley, New York, 1965. [5] F. Carpi and E. Smela (Eds). Biomedical Applications of Elecroactive Polymer Actuators. John Wiley & Sons Ltd., Chichester, UK, 2009. [6] F. Carpi, D. De Rossi, R. Kornbluh, R. Pelrine and P. Sommer-Larsen (Eds). Dielectric Elastomers as Electromechanical Transducers. Elsevier, Oxford, UK, 2008a. [7] F. Carpi, G. Gallone, F. Galantini and D. De Rossi. Silicone-poly(hexylthiophene) blends as elastomers with enhanced electromechanical transduction properties. Adv. Funct. Mat. 18, 235-241, 2008b. [8] F. Carpi, G. Fantoni, G. Frediani and D. De Rossi. Buckling actuators with integrated displacement sensor. In Dielectric Elastomers as Electromechanical Transducers, F. Carpi, D. De Rossi, R. Kornbluh, R. Pelrine and P. Sommer-Larsen (Eds). Elsevier, Oxford, UK, 2008c, pp. 132–140. [9] G. deBotton. Transversely isotropic sequentially laminated composites in finite elasticity. J. Mech. Phys. Solids 53, 1334-1361, 2005. [10] G. deBotton, L. Tevet-Deree and E.A. Socolsky. Electroactive heterogeneous polymers: analysis and applications to laminated composites. Mech. Adv. Mat. Struct. 14, 13–22, 2007. [11] A. Dorfmann and R.W. Ogden. Nonlinear electroelasticity. Acta Mech. 174, 167–183, 2005. [12] A. Dorfmann and R.W. Ogden. Nonlinear electroelastostatics: incremental equations and stability. Int. J. Eng. Sci. 48, 1–14, 2010. [13] M.A. Dowaikh and R.W. Ogden. Interfacial waves and deformations in pre-stressed elastic media. Proc. R. Soc. Lond. A 433, 313-328, 1991. [14] A.C. Eringen and G.A. Maugin. Electrodynamics of Continua. Springer, New York, 1989. [15] M. Gei. Elastic waves guided by a material interface. Eur. J. Mechanics-A/Solids 27, 328–345, 2008. [16] G. Geymonat, S. M¨ller and N. Triantafyllidis. Homogenization of nonlinearly elastic mau terials, microscopic bifurcation and macroscopic loss of rank-one convexity. Arch. Rat. Mech. Anal. 122, 231–290, 1993. [17] N.C. Goulbourne, E. Mockensturm and M. Fracker. A nonlinear model for dielectric elastomer membranes. J. Appl. Mech. 72, 899–906, 2005.

38

[18] R. Hill. On uniqueness and stability in the theory of finite elastic strain. J. Mech. Phys. Solids 5, 229–241, 1957. [19] R. Hill. Eigenmodal deformations in elastic/plastic continua. J. Mech. Phys. Solids 15, 371-386, 1967. [20] R. Hill and J.W. Hutchinson. Bifurcation phenomena in the plane tension test. J. Mech. Phys. Solids 23, 239–264, 1975. [21] C. Huang, Q.M. Zhang, G. deBotton and K. Bhattacharya. All-organic dielectricpercolative three-component composite materials with high electromechanical response. Appl. Phys. Lett. 84, 4391–4393, 2004. [22] G. Kofod, P. Sommer-Larsen, R. Kornbluh and R. Pelrine. Actuation response of polyacrylate dielectric elastomers. J. Intell. Mater. Syst. Struct. 14, 787–793, 2003. [23] J.C. Koo, H.R. Choi, K. Jung, J. Nam, Y. Lee and S. Lee. A new braille display system design using a polymer-based soft actuator tactile display. In Dielectric Elastomers as Electromechanical Transducers, F. Carpi, D. De Rossi, R. Kornbluh, R. Pelrine and P. Sommer-Larsen (Eds). Elsevier, Oxford, UK, 2008, pp. 239–248. [24] R. Kornbluh and R. Pelrine. High performance acrylic and silicone elastomers. In Dielectric Elastomers as Electromechanical Transducers, F. Carpi, D. De Rossi, R. Kornbluh, R. Pelrine and P. Sommer-Larsen (Eds). Elsevier, Oxford, UK, 2008, pp. 33–42. [25] Y. Liu, L. Liu, K. Yu, S. Sun and J. Leng. An investigation on electromechanical stability of dielectric elastomers undergoing large deformation. Smart Mater. Struct. 18, 095040, 2009. [26] R.M. McMeeking and C.M. Landis. Electrostatic forces and stored energy for deformable dielectric materials. J. Appl. Mech. 72, 581–590, 2005. [27] J.C. Michel, O. Lopez-Pamies, P.P. Castaneda and N. Triantafyllidis. Microscopic and macroscopic instabilities in finitely strained porous elastomers. J. Mech. Phys. Solids 55, 900–938, 2007. [28] M. Moscardo, X. Zhao, Z. Suo and Y. Lapusta. On designing dielectric elastomer actuators. J. Appl. Phys. 104, 093503, 2008. [29] A. Needleman. Necking of pressurized spherical membranes. J. Mech. Phys. Solids 24, 339–359, 1976. [30] A. Needleman and M. Ortiz. Effect of boundaries and interfaces on shear-band localization. Int. J. Solids Structures 28, 859–877, 1991. [31] R.W. Ogden. Non-Linear Elastic Deformations. Ellis Horwood, Chichester, 1984. [32] R. Pelrine, R. Kornbluh and J. Joseph. Electrostriction of polymer dielectrics with compliant electrodes as a means of actuation. Sens. Act. A 64, 77-85, 1998. [33] R. Pelrine, R. Kornbluh, Q. Pei and J. Joseph. High-speed electrically actuated elastomers with strain greater than 100%. Science 287, 836-839, 2000.

39

[34] H. Petryk. Theory of bifurcation and instability in time-independent plasticity. In Bifurcation and Stability of Dissipative systems, Q.S. Nguyen (Ed). Springer-Verlag, Wien, 1993, pp. 95–152. [35] J.S. Plante and S. Dubowsky. Large-scale failure modes of dielectric elastomer actuators. Int. J. Solids Structures 43, 7727–7751, 2006. [36] Z. Suo, X. Zhao and W.H. Greene. A nonlinear field theory of deformable dielectrics. J. Mech. Phys. Solids 56, 467–486, 2008. [37] R.A. Toupin. The elastic dielectric. Arch. Rat. Mech. Anal. 5, 849–915, 1956. [38] N. Triantafyllidis. Bifurcation phenomena in pure bending. J. Mech. Phys. Solids 28, 221–245, 1980. [39] N. Triantafyllidis and S. Bardenhagen. The influence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models. J. Mech. Phys. Solids 44, 1891–1928, 1996. [40] N. Triantafyllidis and F.K. Lehner. Interfacial instability of density-stratified two-layer systems under initial stress. J. Mech. Phys. Solids 41, 117–142, 1993. [41] N. Triantafyllidis and B.N. Maker. On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites. J. Appl. Mech. 52, 794–800, 1985. [42] N. Triantafyllidis, M.D. Nestorovic and M.W. Schraad. Failure surfaces for finitely strained two-phase periodic solids under general in-plane loading. J. Appl. Mech. 73, 505–516, 2006. [43] L. van Hove. Sur l’extension de la condition de Legendre du calcul des variations aux integrales multiples a plusieurs fonctions inconnues. Proc. Sect. Sci. K. Akad. van Wetenschappen, Amsterdam, vol. 50, 1947, pp. 18-23. [44] D.K. Vu, P. Steinmann and G. Possart. Numerical modelling of non-linear electroelasticity. Int. J. Numer. Meth. Engrg. 70, 685–704, 2007. [45] Q.M. Zhang, H. Li, M. Poh, F. Xia, Z.-Y. Cheng, H. Xu and C. Huang. An all-organic composite actuator material with a high dielectric constant. Nature 419, 284-289, 2002. [46] X. Zhao and Z. Suo. Method to analyze electromechanical stability of dielectric elastomers. Appl. Phys. Lett. 91, 061921, 2007. [47] X. Zhao, W. Hong and Z. Suo. Electromechanical hysteresis and coexistent states in dielectric elastomers. Phys. Rev. B 76, 134113, 2007. [48] X. Zhao and Z. Suo. Electrostriction in elastic dielectrics undergoing large deformation. J. Appl. Phys. 104, 123530, 2008. [49] X. Zhao and Z. Suo. Electromechanical instability in semicrystalline polymers. Appl. Phys. Lett. 95, 031904, 2009.

40

Appendix A - Push-forward operations defining quantities in the updated lagrangian formulation. The integration of (14)2,4 over ∂B and the change of variables from reference to current configurations, observing that F is continuous across the interface ∂B, provide the following push-forward transformations: ˙ [[S]]n0 dA0 =
∂B
0

˙ 1 [[S]] FT n dA = J ∂B

[[Σ]]n dA,
∂B

(A.1)

and
∂B
0

˙ [[D0 ]] · n0 dA0 =

1 ˙ [[D0 ]] · FT n dA = J ∂B

ˆ [[D]] · n dA,
∂B

(A.2)

which yield the identities (21). The integrals in (A.1) and (A.2) represent the total incremental force and charge in the system, respectively, that vanish as incremental body forces and volume charge are null. Therefore, application of the divergence theorem to (A.1) and (A.2) yields the incremental field equations (24)1,2 in the updated lagrangian formulation. The updated lagrangian formulation of eqn. (24)3 is obtained considering an arbitrary integration path Γ0 in the reference configuration and an infinitesimal fibre dl0 tangent to Γ0 . In the current configuration Γ = χ(Γ0 ) and dl = F dl0 , so that ˙ E0 · dl0 =
Γ

Γ0

˙ E0 · F−1 dl =
Γ

ˆ E · dl.

(A.3)

(A.3) provides (in a way alternative to eqn. (22)) the relationship (23). Since the electric field is conservative, the integral (A.3) evaluated along a close path Γ0 is null. Therefore application c of Stokes’ theorem to (A.3) provides eqn. (24)3 . Appendix B - Coefficients of the expansion (65)

(ca )2 ˜ a 2 (cb )2 ˜ b 2 ˜b ˜a tr(Z ) + tr(Z ) + ca cb tr AZ A−1 Z , 2 2 (ca )4 ˜ a 4 (cb )4 ˜ b 4 ca cb ˜a ˜b ξ2 = tr(Z ) + tr(Z ) + tr (Z )2 A−1 (Z )2 24 24 4 a (cb )3 b (ca )3 c c ˜b ˜a ˜b ˜a + tr A(Z )3 A−1 Z + tr AZ A−1 (Z )3 , 6 6 a )6 b )6 b ca )3 (c (c (c ˜a ˜b ˜b ˜a tr(Z )6 + tr(Z )6 + tr A(Z )3 A−1 (Z )3 ξ3 = 720 720 36 ca (cb )5 cb (ca )5 ˜b ˜a ˜b ˜a + tr A(Z )5 A−1 Z + tr AZ A−1 (Z )5 120 120 a )4 (cb )2 (ca )2 (cb )4 (c ˜b ˜a ˜b ˜a + tr A(Z )2 A−1 (Z )4 + tr A(Z )4 A−1 (Z )2 , 48 48 ξ1 = ˜ ˜ ξ4 = (ca cb )2 tr AZ A−1 Z
b a 2

(B.1)

˜b ˜a − tr A(Z )2 A−1 (Z )2 )

,

(B.2)

41

ξ5 = (ca )2 cb

1 1 ˜b ˜a ˜b ˜a ˜b ˜a ˜b ˜a tr A(Z )3 A−1 Z AZ A−1 Z + tr A(Z )2 A−1 Z A(Z )2 A−1 Z 3 4 7 2 1 ˜b ˜b ˜a ˜a ˜b ˜a − tr A(Z )4 A−1 (Z )2 + (ca )4 cb tr AZ A−1 (Z )3 AZ A−1 Z 12 3 1 7 ˜b ˜a ˜b ˜a ˜b ˜a + tr AZ A−1 (Z )2 AZ A−1 (Z )2 − tr A(Z )2 A−1 (Z )4 4 12 3 ˜b ˜b ˜a ˜b ˜a ˜a tr AZ A−1 Z A(Z )2 A−1 (Z )2 − tr A(Z )3 A−1 (Z )3 , + ca cb
4 4 b b a b a b b a b

(B.3)

ξ6 = (ca )2 cb

˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜ 5 tr A(Z )2 A−1 Z AZ A−1 Z A(Z )2 A−1 + 3 tr A(Z )2 A−1 Z A(Z )2 A−1 Z
a b

a

˜ ˜ ˜ − 8 tr A(Z )2 A−1 (Z )2 A(Z )2 A−1 + (ca )4 cb
2

˜a ˜b ˜b ˜a ˜b ˜a ˜b ˜a ˜b ˜a 5 tr (Z )2 AZ A−1 Z AZ AZ A−1 Z + 3 tr AZ A−1 (Z )2 AZ A−1 (Z )2
a

˜ ˜ −8 tr A(Z )2 A−1 (Z )4

b

˜b ˜a ˜b + (ca cb )3 tr (AZ A−1 Z )3 − 13 tr A Z ,

3

˜a A−1 Z

3

˜b ˜a ˜b ˜a + 12 tr AZ A−1 Z A(Z )2 A−1 (Z )2 ˜ρ with Z = Zρ /k1 and A = (Ga )−1 Gb .

(B.4)

42