rspa.royalsocietypublishing.org

Research Cite this article: Cohen N, Menzel A, deBotton G. 2016 Towards a physics-based multiscale modelling of the electro-mechanical coupling in electro-active polymers. Proc. R. Soc. A 472: 20150462. http://dx.doi.org/10.1098/rspa.2015.0462 Received: 5 July 2015 Accepted: 25 January 2016

Towards a physics-based multiscale modelling of the electro-mechanical coupling in electro-active polymers Noy Cohen1,3 , Andreas Menzel3,4 and Gal deBotton1,2 1 Department of Mechanical Engineering, and 2 Department of

Biomedical Engineering, Ben-Gurion University, Beer-Sheva 84105, Israel 3 Department of Mechanical Engineering, Institute of Mechanics, TU Dortmund, 44227 Dortmund, Germany 4 Division of Solid Mechanics, Lund University, PO Box 118, 22100 Lund, Sweden NC, 0000-0003-2224-640X

Subject Areas: mechanical engineering, mathematical physics Keywords: dielectrics, electro-active polymers, electromechanical coupling, multi-scale analysis Author for correspondence: Noy Cohen e-mail: [email protected]

Owing to the increasing number of industrial applications of electro-active polymers (EAPs), there is a growing need for electromechanical models which accurately capture their behaviour. To this end, we compare the predicted behaviour of EAPs undergoing homogeneous deformations according to three electromechanical models. The first model is a phenomenological continuumbased model composed of the mechanical Gent model and a linear relationship between the electric field and the polarization. The electrical and the mechanical responses according to the second model are based on the physical structure of the polymer chain network. The third model incorporates a neo-Hookean mechanical response and a physically motivated microstructurally based long-chains model for the electrical behaviour. In the microstructural-motivated models, the integration from the microscopic to the macroscopic levels is accomplished by the micro-sphere technique. Four types of homogeneous boundary conditions are considered and the behaviours determined according to the three models are compared. For the microstructurally motivated models, these analyses are performed and compared with the widely used phenomenological model for the first time. Some of the aspects revealed in this investigation, such as the dependence of the intensity of the polarization

2016 The Author(s) Published by the Royal Society. All rights reserved.

Dielectric elastomers (DEs) are materials that deform in response to electrostatic excitation. Owing to their light weight, flexibility and availability, these materials can be used in a wide variety of applications such as artificial muscles [1], energy-harvesting devices [2,3], micropumps [4] and tunable wave guides [5,6], among others. The principle of the actuation is based on the attraction between two oppositely charged electrodes attached to the faces of a thin soft elastomer sheet. Owing to Poisson’s effect, the sheet expands in the transverse direction. Toupin [7], in his pioneering analysis, found that this electromechanical coupling is characterized by a quadratic dependence on the applied electric field. This was later verified experimentally (e.g. [8]). However, DEs have a low energy density in comparison with other actuators such as piezoelectrics and shape memory alloys [1]. Furthermore, their feasibility is limited due to the high electric fields (approx. 100 MV m−1 ) required for a meaningful actuation as a result of the relatively low ratio between the dielectric and the elastic moduli [9,10]. Specifically, common flexible polymers have low dielectric moduli while polymers with high dielectric moduli are usually stiff. Nevertheless, a few recent works suggest that this ratio may be improved. Huang et al. [11] demonstrated experimentally that organic composite EAPs (electro-active polymers) experience more than 8% actuation strain in response to an activation field of 20 MV m−1 . The experimental work of Stoyanov et al. [12] showed that the actuation can be dramatically improved by embedding conducting particles in a soft polymer. In parallel, theoretical works dealing with the enhancement of coupling in composites also hint at the possibility of improved actuation with an appropriate adjustment of their microstructure [5,13–16]. The above findings motivate an in-depth multiscale analysis of the electromechanical coupling in elastic dielectrics which is inherent from their microstructure. In this work, we consider the class of polymer dielectrics. A polymer is a hierarchical structure of polymer chains each of which is a long string of repeating monomers. We start by analysing the behaviour of a single monomer in a chain. Next, the response of a chain is obtained by a first-level integration from the singlemonomer level to the chain level. Finally, the macroscopic behaviour of the polymer is obtained by a higher level summation over all chains. In this work, we use existing constitutive models for the chains and concentrate on the higher level summation from the chain level to the macroscopiccontinuum level. To this end, the physically motivated micro-sphere technique, which enables to extend one-dimensional models to three-dimensional models by appropriate integration over the orientation space [17,18], is exploited. Accordingly, this method lends itself to the characterization of polymer networks since the single chain is often treated as a one-dimensional object which is aligned along the chain’s end-to-end vector [19,20]. The response of a polymer subjected to purely mechanical loadings was extensively investigated at all length scales. A macroscopic-level analysis and models describing the behaviour of soft materials undergoing large deformation, such as polymers, were developed by Ogden [21]. Semi-empirical phenomenological modelling of the response of shape memory polymers has been discussed by Lu & Du [22] and Lu & Huang [23]. The microscopic-level analysis of Kuhn & Grün [24] yielded a Langevin-based constitutive relation and paved the way to various multiscale models such as the three-chain model [25], the tetrahedral model [26,27] and the eight-chain model [28]. Corresponding micro-sphere implementation of the Langevin model was carried out by Miehe et al. [19]. The electric response of dielectrics to electrostatic excitation was examined by Tiersten [29] and Hutter et al. [30], among others, macroscopically as well as through their microstructure. Starting with the examination of a single charge under an electric field, the relations between the

...................................................

1. Introduction

2

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

field on the deformation, highlight the need for an in-depth investigation of the relationships between the structure and the behaviours of the EAPs at the microscopic level and their overall macroscopic response.

Consider the deformation of a hyperelastic dielectric body subjected to electro-mechanical loading from a referential configuration to a current one. In the reference configuration, the body occupies a region V0 ∈ R3 with a boundary ∂V0 . The referential location of a material point is X. In the current configuration, the body occupies a region V ∈ R3 with a boundary ∂V, and we denote the location of a material point by x. The mapping of positions of material points from the reference to the current configurations is x = ϕ(X), and the corresponding deformation gradient is F = ∇X ϕ, where the operation ∇X denotes the gradient with respect to the referential coordinate system. The right and left Cauchy–Green strain tensors are C = FT F and b = F FT . The ratio between the volumes of an infinitesimal element in the current and the reference configurations, J = det F, is strictly positive. In the case of incompressible materials, which are of interest in this work, J = 1. The electric field E acting on the body satisfies the governing equation ∇x × E = 0 in the entire space, where ∇x denotes the gradient with respect to the current coordinate system. Consequently, we can define a scalar field, the electric potential φ, such that E = −∇x φ. The electric displacement field is (2.1) D = ε0 E + P, where ε0 is the permittivity of the vacuum and P is the polarization or the electric dipole-density. We recall that in vacuum P = 0. In the absence of free charges, the electric displacement field is governed by the local equation (2.2) ∇x · D = 0. In [31], the referential counterparts E(0) and D(0) of the electric field and the electric displacement were determined, respectively. Specifically, E(0) = FT E

(2.3)

...................................................

2. Theoretical background

3

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

different electric macroscopic quantities, such as the electric displacement, the electric field and the polarization, and microscopic quantities, such as the free and bound charge densities and the dipoles, were defined and analysed. The study of the response of dielectrics to a coupled electromechanical loading initiated with the pioneering work of Toupin [7], who performed a theoretical analysis at the macroscopic level. Later on, an invariant-based representation of the constitutive behaviour of EAPs was introduced by Dorfmann & Ogden [31]. Subsequently, Ask et al. [32–34], Gei et al. [5] and Jiménez & McMeeking [35] investigated the possible influence of the deformation and its rate on the electromechanical coupling. Thylander et al. [20] made use of a corresponding microsphere technique at the chain level. By employing macroscopic constitutive models for the mechanical and the electrical behaviour of the polymer chains, a few boundary value problems were solved by means of the numerical implementation of the micro-sphere technique and the finite-element method. Initial multiscale physically motivated analyses of the electromechanical response were performed by Cohen & deBotton [36,37]. This work focuses on the implementation of the latter model to EAPs experiencing homogeneous deformations under various types of boundary conditions and comparing their predicted response to the one stemming from the phenomenological models. We begin this work revisiting the common phenomenological models and the recent physically motivated model. We also recall the micro-sphere framework. Next, extending the work of Miehe et al. [19], we include electro-mechanical coupling by employing the micro-sphere technique to compute the macroscopic behaviour of dielectrics with randomly oriented and uniformly distributed dipoles experiencing the macroscopic rotation. In §4, we extend this analysis and determine the response of a polymer according to the physically motivated model under different boundary conditions. These results are compared with the corresponding phenomenological models in order to highlight the differences between the two approaches. The conclusion are gathered in §5.

and D(0) = JF−1 D,

4 (2.4)

(2.5)

In accordance with our assumption that the dielectric solid can be treated as a hyper-elastic material, its constitutive behaviour can be characterized in terms of a scalar electrical enthalpydensity function W. We further assume that W can be decomposed into a mechanical and a coupled contributions, i.e. W(F, E) = W0 (F) + Wc (F, E), where W0 (F) characterizes the material response in the absence of electric excitation and Wc (F, E) accounts for the difference between W with and without electric excitation [36–39]. Accordingly, the polarization is determined via P=−

1 ∂Wc , J ∂E

(2.6)

and the total stress developing in the material can be written as the sum σ = σ m + π + σ v,

(2.7)

1 ∂W(F) T F J ∂F

(2.8)

where σm =

is the mechanical stress due to the deformation of the material, π =E⊗P is the polarization stress stemming from the applied electric field in the dielectric and   σ v = ε0 E ⊗ E − 12 [E · E]I

(2.9)

(2.10)

is the Maxwell stress in vacuum, where I is the second-order identity tensor [38]. We emphasize that this decomposition is purely modelling based as in an experiment the total stress can be measured, but the contributions of the individual components cannot be distinguished. In this work, we consider incompressible materials in which the component of the stress that is energy conjugate to J cannot be determined from the constitutive relation. Accordingly, this workless pressure-like component, namely p I, is to be added to equation (2.7) for the total stress and determined from the equilibrium equations and the boundary condition. Assuming no body forces, the total stress satisfies the local equilibrium equation ∇x · σ = 0.

(2.11)

The electrical boundary conditions are given in terms of either the electric potential or the charge per unit area ρa , such that D · nˆ = −ρa , where nˆ is the outward pointing unit normal. Practically, in EAPs, ρa is the charge on the electrodes. The mechanical boundary conditions are given in terms of the displacement or the mechanical traction t. Owing to the presence of the electric field the total stress in the vacuum outside of the body does not necessarily vanish. Therefore, the mechanical traction at the boundary is [σ − σ v ] · nˆ = t, where the expression for σ v , the Maxwell stress outside the material, is given in equation (2.10) in terms of the electric field in the vacuum [39].

(a) Existing models for the behaviour of dielectrics Within the framework of finite deformation elasticity, the simplest constitutive model is the wellknown neo-Hookean model that requires only one material parameter. The corresponding strain

...................................................

P(0) = JF−1 P.

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

where ∇X × E(0) = 0 and ∇X · D(0) = 0. We note that unlike E(0) and D(0) , the referential polarization is not uniquely defined. In order to ensure that the referential polarization is energy conjugate to the referential electric field such that (1/J)E(0) · P(0) = E · P, we adopt the definition

energy-density function (SEDF) is

5 (2.12)

where μ is the shear modulus and I1 = tr(C) is the first invariant of the right Cauchy–Green strain tensor. This model does not capture the lock-up effect experimentally observed in various polymers which corresponds to a significant stiffening of the material at large strains. Gent [40] proposed a phenomenological constitutive model in which this effect is accounted for. The SEDF for this model,   I1 − 3 μJm ln 1 − , (2.13) W0G (F) = − 2 Jm depends on two parameters, μ and Jm . The latter is the lock-up parameter such that Jm + 3 is the value of I1 at the maximum deformation that can be experienced by the polymer. Specifically, the expression in equation (2.13) becomes unbounded at I1 = Jm + 3, thus capturing the lock-up phenomenon. With regard to the response of the dielectric to electrical excitation, a quadratic dependence of Wc on E is commonly assumed, leading to the relation P = χ E,

(2.14)

where χ is the susceptibility of the material [41]. Consequently, the electric displacement is D = εE,

(2.15)

where ε = ε0 + χ is the permittivity of the material. It is common practice to assume that χ , and consequently ε, are constants. This assumption is in agreement with the invariant based framework of Dorfmann & Ogden [31]. Furthermore, an experiment carried out by Di Lillo et al. [42] on VHB 4910 polymer film showed that this assumption is fairly accurate. Recent experiments with various types of polymers imply that the permittivity, and therefore the relationship between the polarization and the electric field, is deformation dependent [43–46]. A possible explanation for this dependency of the susceptibility on the deformation is related to the alteration of the inner structure of the polymer [36,37]. In the polymer, the monomers in the chain can move or rotate relative to their neighbours thus providing the chains with a freedom to deform [47]. In order to better understand the response of polymers to an electro-mechanical loading, their microstructure should be accounted for. This can be accomplished in terms of a multiscale analysis consisting of three stages: the first involves the examination of the behaviour of the monomers, the second corresponds to the analysis of the response of the individual chains and the third deals with the polymer behaviour at the continuum level. Treloar [48] and Flory & Rehner [26] carried out such multiscale analyses of polymers subjected to mechanical loading. It was assumed that the directions of the monomers, or links, composing a chain are random, and consequently it was found that the chains are distributed according to a Gaussian distribution. Based on statistical considerations and the laws of thermodynamics, the variation in the entropy of the chain due to its deformation was determined. The overall variation in the entropy of the polymer is computed by summing the entropies of the chains. Remarkably, their result recovered the macroscopic neo-Hookean behaviour. Furthermore, a comparison between the micro- and macro-analyses related the macroscopic shear modulus to the number of chains per unit volume N0 . Specifically, it was found that μ = k T N0 ,

(2.16)

where k and T are the Boltzmann constant and the absolute temperature, respectively. Owing to the assumptions that lead to the use of Gaussian statistics, the lock-up effect was not captured in the above-mentioned analysis. A more rigorous examination of the polymer behaviour by Kuhn & Grün [24] revealed that this phenomenon is a result of the finite extensibility

...................................................

μ [I1 − 3], 2

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

W0nH (F) =

6

where l is the length of a link, nl is the number of links in a chain, r is the distance between the two ends of the chain, and L −1 (•) is the inverse of the Langevin function L (β) ≡ coth(β) −

1 r = . β nl l

(2.18)

Assuming that all chains undergo the macroscopic deformation, i.e. r = F r0 , where r and r0 are the current and referential end-to-end vectors, respectively, the stress associated with a chain is derived from equation (2.17),   √ r0 −1 r L Fˆr0 ⊗ Fˆr0 , = kT n (2.19) σ LC l m r nl l where rˆ0 is a unit vector in the direction of r0 and equation (2.8) is used. Here, r0 = l

√

nl

(2.20)

is the average length of the referential end-to-end vectors [47,49]. The quantity r/nl l ≤ 1 describes the ratio between the end-to-end and the contour lengths of the chain, and from equation (2.20) it follows that



r0 1 r = Fˆr0 · Fˆr0 = Fˆr0 · Fˆr0 √ . (2.21) nl l nl l nl It can be shown that the limit r/nl l → 1 results in L −1 (r/nl l) → ∞, thus capturing the experimentally observed lock-up phenomenon. The lock-up stretch is associated with the chain undergoing the largest extension such that the stretch ratio of its end-to-end vector is λmax =

√

nl .

(2.22)

It is important to note that the first term in the Taylor series expansion of equation (2.19) about r/nl l = 0 reproduces the Gaussian model [24,49]. A few works proposed models that consider specific finite networks of chains. The three-chain model by Wang & Guth [25] examines a network of three chains which are located along the axis of the principal directions of the deformation gradient. Flory & Rehner [26] and Treloar [27] proposed a network of four chains that are linked together at the centre of a regular tetrahedron, and their other ends are located at the vertices of the tetrahedron. The tetrahedron deforms according to the macroscopic deformation while the chains experience different stretches. In the model proposed by Arruda & Boyce [28], eight representative chains in specific directions relative to the principal system of the macroscopic deformation gradient are used to determine the macroscopic behaviour. An anisotropic worm-like chain model in which no inherent alignment between the chosen and the principal coordinate systems is assumed was considered by Kuhl et al. [50]. A multiscale-level analysis of the response of polymers with Gaussian distribution to electromechanical loading was carried out by Cohen & deBotton [36]. In this study, the changes in the magnitudes of the dipolar monomers due to the applied electric field and their rearrangement due to the mechanical deformation were accounted for. This analysis assumed that the polymer chains are composed of a very large number of monomers and that the interaction of the polymer with the electric field is solely through the interactions of the dipolar monomers with the field. Following a model described in Stockmayer [51], Cohen & deBotton [37] considered the class of uniaxial dipoles in which the dipole is aligned with the line segment between the two contact points of a monomer with its neighbours. Thus, taking ξˆ to be the unit vector along this line segment, the first-order estimate of the dipole moment of a monomer with respect to the electric

...................................................

(2.17)

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

of the chains. According to this analysis, the SEDF associated with a polymer chain is      r r L −1 (r/nl l) LC −1 W0 = kTnl L + ln , nl l nl l sinh(L −1 (r/nl l))

field is

7

mu = K[ξˆ ⊗ ξˆ ]E,

(2.23)

K [I − ξˆ ⊗ ξˆ ]E. 2

(2.25)

An expression for the dipole of a chain made out of transversely isotropic dipoles is determined by following the steps which led to the derivation of equation (2.24) for the uniaxial dipoles. The resulting polarization of the polymer is determined by summing the dipoles of the chain in a representative volume element of a volume V R (e.g. [41]) P=

1 mc , VR

(2.26)

i

and the resulting polarization stress is computed via equation (2.9).

(b) The micro-sphere technique Consider a unit sphere whose surface represents the directions of the referential end-to-end vectors. The directional averaging of a quantity • over the unit sphere can be approximated with a discrete summation  m

1 • dA ≈ •(i) w(i) , (2.27) • = 4π A i=1

(i)

where the index i = 1, . . . , m refers to a unit direction vector rˆ0 , •(i) is the value of quantity • in

(i) the direction rˆ0

m w = and w(i) is an appropriate non-negative weight function constrained by Σi=1 i 1 (e.g. [17,19,52]). In general, equation (2.27) can be combined with a more general anisotropic distribution function [53,54]. (i) In the case of polymers, the vectors rˆ0 represent the directions of the end-to-end vectors of the polymer chains or, from a numerical point of view, the integration directions in orientation space. For randomly and isotropically oriented chains, or rather isotropic integration schemes, (i) the vectors rˆ0 satisfy m

(i)

rˆ0 w(i) = 0

(2.28)

i=1

and m

i=1

(i)

(i)

rˆ0 ⊗ rˆ0 w(i) =

1 I. 3

(2.29)

In view of equations (2.28) and (2.29), the micro-sphere technique naturally lends itself to the calculation of the macroscopic polarization and stress. Specifically, for a polymer with chains composed of nd dipoles and N0 chains per unit referential volume, equation (2.26) may be

...................................................

mt =

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

where K is a material constant. The dipole of a chain composed of nd uniaxial dipoles with an end-to-end vector in the direction rˆ is [36]   16 r K nd I+ [I − 3ˆr ⊗ rˆ] E. (2.24) mc ≈ 3 3 π 2 nd l

If we assume that all chains undergo the macroscopic deformation, then rˆ = F r0 / [F r0 ] · [F r0 ]. In accordance with a second type of dipoles discussed by Stockmayer [51], Cohen & deBotton [37] proposed an expression for a transversely isotropic (TI) dipole, where the dipole is aligned with the projection of the electric field on the plane perpendicular to ξˆ ,

written as P=

8 (2.30)

σ Lm =

N0 LC σ m , J

(2.31)

where we use the notation suggested in equation (2.27). Bažant & Oh [17] demonstrated that a specific choice of 42 directions guarantees sufficient accuracy for the application discussed in their work. We follow this conjecture where the integration directions and the corresponding weight functions are given in table 1 of Bažant & Oh [17]. We note that other integration schemes are available, as demonstrated by Waffenschmidt et al. [55,56] and references therein.

3. Dielectrics with randomly distributed monomers Consider a model of a dielectric composed of n0 monomers per unit referential volume, which are treated as mechanical rods and electric dipoles. The dielectric is subjected to a mechanical deformation, locally represented by F, and an electric field E. We assume that the electric field induced on a monomer by its neighbours is small in comparison with the applied electric field [37]. We examine first a dielectric with uniaxial monomers, the behaviour of which is governed by the quadratic form in equation (2.23). If all of the dipoles experience the macroscopic rotation, i.e. ξˆ = R ξˆ 0 where R = FC−1/2 is a proper rotational tensor, then the polarization according to equation (2.30) is P = n0 mc = n0 KR

42

i=1

n0 K (i) (i) ξˆ 0 ⊗ ξˆ 0 w(i) RT E = E, 3

(3.1)

where equations (2.27) and (2.29) are used. The corresponding polarization stress is π =E⊗P=

n0 K E ⊗ E. 3

(3.2)

In the case of a dielectric composed of n0 TI dipolar monomers per unit referential volume (cf. equation (2.25)), which mechanically act as rigid rods, the same assumptions that led to equation (3.1) lead to P = n0 mc = n0

42  K n0 K (i) (i) E. I − Rξˆ 0 ⊗ ξˆ 0 RT w(i) E = 2 3

(3.3)

i=1

Accordingly, the expression for the polarization stress is given in equation (3.2). We note that the polarization and polarization stress calculated in equations (3.1), (3.3) and (3.2) are identical to the exact expressions obtained by Cohen & deBotton [37].

4. Dielectric elastomers We examine the behaviours of incompressible DEs according to the above-described microstructural model under various homogeneous electromechanical loading conditions and compare the predicted repose to common phenomenological models. To facilitate the comparison, we assume that in the limit of infinitesimal deformations and small electric excitations the models admit the same behaviour. Specifically, we assume that the initial shear modulus and the electric susceptibility are μ = 0.1 MPa and χ = 3 ε0 , as these values are representative values for the moduli of polymers. Whenever the lock-up effect is accounted for, we choose the model parameters such that under purely mechanical biaxial loading the lock-up stretch is λlu = 5.

...................................................

and the macroscopic stress according to equation (2.19) as

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

N0 mc , J

undeformed configuration

9

flexible electrodes DV

Z X

z

x

Figure 1. A sketch of a polymer under an electro-mechanical loading. (Online version in colour.) The precise models and the numerical values assumed for their parameters are as follows: (i) The macroscopic-phenomenological model. The mechanical behaviour is characterized by the Gent model (2.13) with the aforementioned shear modulus and Jm = 47. The electric behaviour is determined according to the linear model (2.15) with the initial permittivity ε = 4ε0 . (ii) The physically motivated microscopic model. The Langevin model (2.19) is used in order to describe the mechanical behaviour, where nl = 25 is chosen to fit the assumed lock-up stretch according to equation (2.22) and N0 = μ/kT. We employ the long-chains model (2.24) with chains that are composed of nd = 100 uniaxial dipoles (equation (2.23)) to characterize the dielectric response of the polymer. The material constant K is chosen such that KN0 nd /3 = 3ε0 to ensure that the referential polarization is identical to the one admitted by the macroscopic model. (iii) The Gaussian model. The neo-Hookean model (2.12) with (2.16) are used, where N0 = μ/kT, in conjunction with the long-chains model (2.24) to characterize the mechanical and the electrical behaviours, respectively. We assume that a chain is composed of nd = 100 uniaxial dipoles (equation (2.23)), where the chosen long-chains model constant is identical to the one determined for the microscopic model. In the following, we examine a thin layer of a polymer whose opposite faces are covered with flexible electrodes with negligible stiffness (figure 1). The electrodes are charged with opposite charges so that the difference in the electric potential induces an electric field across the layer. From a mechanical point of view, we consider four boundary conditions. In the first two cases, different displacements are prescribed at the boundary and consequently the deformation gradient is defined. In the following two cases, we set the traction on the boundaries. We choose a cartesian coordinate system in which the referential electric field is aligned with the y-axis ˆ and calculate the macroscopic polarization according to the microscopic and the Gaussian models via equation (2.30). The polarization stress is computed according to equation (2.9). The mechanical stress is computed via equations (2.13), (2.12) and (2.19) for the macroscopic, the Gaussian and the microscopic models, respectively. The pressure term is determined from the traction free boundaries to which the electrodes are attached and subsequently the total stress is computed via equation (2.7). For convenience, we define the dimensionless normal stress along the xˆ -axis σ¯ = (1/μ)ˆx · σ xˆ and the dimensionless referential electric displacement along

referential electric (0)field and ¯ = (1/√εμ)D(0) · y, ˆ respectively. the y-direction ˆ E¯ (0) = (ε/μ)E(0) · yˆ and D

In the following ¯ (0) are E¯ = (ε/μ)E · yˆ and D ¯ = examples, the current configuration counterparts of E¯ (0) and D √ ˆ respectively, as follows from equations (2.3) and (2.4). (1/ εμ)D · y,

(a) Equibiaxial stretching perpendicular to the electric field In this case, the material is stretched along the axes xˆ and zˆ such that λx = λz = λ. As stated previously, the y-axis ˆ is aligned with the referential electric field and due to the assumed incompressibility λy = 1/λ2 . This setting is common in various experiments with EAPs [43–46].

...................................................

y

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

Y

deformed configuration

(a) 100

(b) 20

10

–

s– –200

D

10

–300 5 –400 –500

5

10

2

15

20

25

0

0

l

5

10 – E

15

20

Figure 2. The dimensionless stress σ¯ (a) and electric displacement D¯ (b) versus the stretch of the transverse plane and the dimensionless electric field according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle symbols).

¯ as functions of λ2 and E, ¯ respectively. The curves with the square Figure 2a,b depicts σ¯ and D symbols correspond to the macroscopic model, the curves with the hollow circle symbols to the Gaussian model, and the curves with the filled circle symbols to the microscopic model. The applied referential electric field is E(0) = 50 MV m−1 . We point out that at λ2 = 1 the dimensionless stress according to the three models is not zero but very small. Figure 2a illustrates the stress increase at the lock-up stretch according to the macroscopic and the microscopic models. As expected, this effect is not observed when the Gaussian model is employed. Furthermore, since the electric field tends to stretch the material in the transverse plane, as long as the prescribed stretch is smaller than the electrically induced stretch, the overall stress is compressive. We emphasize that since the potential is held fixed, as the layer is stretched the current electric field increases and hence also the electromechanically induced stress. This gives rise to different types of loss of stability phenomena which are outside the scope of the current work. The reader is referred to the works, for example, by Dorfmann & Ogden [57], Bertoldi & Gei [58], Rudykh et al. [15] and Shmuel et al. [6]. Only when the prescribed stretches are large enough, does the total stress becomes tensile. In figure 2b, we observe a linear dependence of the electric displacement on the electric field according to the macroscopic model, as follows from equation (2.15) and the assumed constant ¯ the slope permittivity. In accordance with the choice of the dimensionless quantities E¯ and D, of this curve is unity. By contrast, the Gaussian and the microscopic models predict a stronger than linear increase in the electric displacement as we stretch the material. This is a result of the predicted increase in the permittivity due to the stretching of a polymer with chains made up of uniaxial dipoles [36].

(b) Pure shear deformation in the plane of the electric field Once again we analyse a thin layer of a polymer whose opposite faces are covered with flexible electrodes with negligible stiffness. As before the y-axis ˆ is aligned with the electric field, but in this case the deformation of the material along the zˆ -axis is constraint such that λz = 1. The layer is stretched along the xˆ -axis such that λx = λ and the incompressibility condition yields λy = 1/λ. Figure 3a depicts the dimensionless normal stress that develops along the xˆ -axis as a function of the stretch according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle symbols) under the applied referential electric field E(0) = 100 MV m−1 . The inability of the neo-Hookean model to capture the lock-up stretch is again clearly depicted. We also

...................................................

15 –100

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

0

(a) 100

(b) 12

11

80

s

8 –

D

20 0

6

–20

4

–40 2

–60 1

2

3

4 l

5

6

7

2

4

6

8

10

12

–

E

Figure 3. The dimensionless stress σ¯ (a) and electric displacement D¯ (b) versus the axial stretch and the dimensionless electric field according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle symbols).

notice that there is a difference between the lock-up stretches predicted by the macroscopic and the microscopic models. This is because the lock-up stretch according to the Langevin model is determined by the maximum eigenvalue of the deformation gradient as follows from equation (2.22), whereas according to the Gent model it depends on the first invariant of the right Cauchy–Green strain tensor. Treloar [49] presented experimental results demonstrating that polymers lock-up at different values under different types of deformations, and therefore we conclude that in this aspect the Gent model may be a better predictor. A comparison of different models and their calibration according to the data reported by Treloar [49] was carried out by Steinmann et al. [59]. The curves with the square, hollow and filled circle symbols in figure 3b correspond to the macroscopic model, the Gaussian model and the microscopic model, respectively. Here, ¯ on E¯ is illustrated. Owing to the constant permittivity, once again we the dependency of D observe the linear dependency predicted by the macroscopic model. The Gaussian and the microscopic models, which are based on the microscopic long-chains model, predict a change in the permittivity as a result of the mechanical stretch, in agreement with the experimental findings of Choi et al. [43], Wissler & Mazza [44], McKay et al. [45] and Qiang et al. [46]. Since the deformation is dictated by the boundary condition, the Gaussian and the microscopic models predict the same electric behaviour.

(c) Equibiaxial actuation normal to the electric field We once again examine a thin layer of a polymer whose opposite faces are covered with flexible electrodes with negligible stiffness. This time, however, the circumferential boundary of the layer is traction free, thus allowing the layer to expand in the plane transverse to the direction of the electric field in response to the electric excitation. We choose the same system of axes as in §4a and, thanks to the symmetry of the loading, the deformation gradient is diagonal with λx = λz = λ. Owing to the assumed incompressibility λy = 1/λ2 . Figure 4a displays the dimensionless referential electric field E¯ (0) as a function of the induced stretch of the transverse plane according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle symbols). The loss of stability discussed in the works of Dorfmann & Ogden [57], Bertoldi & Gei [58], Rudykh & deBotton [60] and Shmuel et al. [6] is revealed once again. We note that after the peak at E¯ (0) ≈ 0.7, even though the current electric field

...................................................

40 –

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

10

60

(a) 1.0

(b)

12

70

50 40

–

–

E (0)

D (0)

0.6 0.4

30 20

0.2

10 0

5

10

15

20

25

2

l

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

–

E (0)

Figure 4. The dimensionless referential electric field E¯(0) (a) and electric displacement D¯ (0) (b) versus the stretch of the transverse plane and the dimensionless referential electric field according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle symbols).

increases monotonically, the Gaussian model predicts a decrease in the referential electric field with an increase of the stretch. In an experiment where the referential electric field is controlled (i.e. the electric potential between the electrodes), the macroscopic and the microscopic models predict a jump in the planar stretch. This effect of a transition between two states is known as snap-through [4,61]. The curves with the square, hollow and filled circle symbols in figure 4b correspond to the macroscopic, the Gaussian and the microscopic models, respectively, where the predicted ¯ (0) on E¯ (0) in the direction of the electric field are depicted. Essentially, this dependencies of D plot illustrates the amount of charge per unit referential surface area as a function of the potential difference divided by the initial thickness of the layer. Practically, this setting can be experimentally obtained by connecting the electrodes to an electric source through which the electric potential between the electrodes is controlled. Initially, we observe an increase of the surface charge with an increase of the electric potential. However, beyond the peak at E¯ (0) ≈ 0.7 there is a reversed trend where, at equilibrium, the electric potential drops while the surface charge increases. This occurs in conjunction with the uncontrollable increase in the area of the actuator as shown in figure 4a. From a practical viewpoint, this implies that beyond the peak, in a manner reminiscent of an electrical short-circuit, excessive current flows from the electric source while the electric potential drops. We stress that due to the thinning of the layer the current electric field increases and may result in a failure of the DE due to electric breakdown. We also note that even though the same electric model is used in both the Gaussian and the microscopic models, there is a difference in the relationships between the electric field and the electric displacement. This is due to the different mechanical deformations resulting from the applied electric field according to the two different models.

(d) Uniaxial actuation normal to the electric field We consider a setting reminiscent of the one considered in §4a, but this time the layer is free to expand only along the xˆ -direction. Consequently, the deformation gradient components are λz = 1, λx = λ and λy = 1/λ. Figure 5a shows the dimensionless referential electric field as a function of the stretch according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle

...................................................

60

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

0.8

(b) 80

3.0

70

–

D(0)

–

E (0)

50

2.0 1.5

40 30

1.0

20

0.5

10 1

2

3

4 l

5

6

7

0

0

0.5

1.0

1.5

2.0

2.5

3.0

–

E (0)

Figure 5. The dimensionless referential electric field E¯(0) (a) and electric displacement D¯ (0) (b) versus the stretch and the dimensionless referential electric field according to the macroscopic model (the curve with the square symbols), the Gaussian model (the curve with the hollow circle symbols) and the microscopic model (the curve with the filled circle symbols).

symbols). As mentioned previously, the predicted lock up stretches are λlu = 5 and λlu ≈ 7 according to the microscopic and the macroscopic models, respectively. Unlike the biaxial case described in §4a, this loading does not admit loss of stability according to the macroscopic and the microscopic models. The Gaussian model, that is based on the mechanical neo-Hookean model, reaches a peak at E¯ (0) ≈ 1 and then monotonically decreases. ¯ (0) on E¯ (0) . According to the Gaussian model, represented Figure 5b depicts the dependence of D by the curve with the hollow circle symbols, no significant changes in the surface charge are observed as we initially increase the electric potential difference. However, beyond the peak at E¯ (0) ≈ 1, this model predicts an unstable behaviour as a result of the electrical long-chains model. In order to maintain equilibrium according to the macroscopic and the microscopic models an increase in the charge on the electrodes requires an increase in the potential difference between them. Thus, the Gaussian model is admitting a behaviour that is qualitatively different from the behaviours of the other two models.

5. Concluding remarks We analysed the behaviour of an incompressible polymer undergoing homogeneous deformations according to a physically motivated model under four types of boundary conditions. This was compared with a phenomenological model based on well-known macroscopic constitutive models for the mechanical and the electrical behaviours. The physically motivated microscopic model combines mechanical and electrical models stemming from the microstructure of the polymer. The third model assumes a Gaussian distribution of the chain network and the mechanical and the electrical behaviours are determined in accordance. We comment that the material parameters nd and nl , which denote the number of dipoles and links, respectively, are used as fitting parameters and in this sense the electrical long-chains model and the Langevin model are not consistent. Further investigations in this regard is needed. In the first two representative examples, we apply a referential electric field by setting the potential difference between the electrodes and control the deformation. In the following two examples, we apply a referential electric field and set the traction on the boundaries of the polymer. In order to determine the polarization and the stress according to the microscopic models, we employ the micro-sphere technique. A comparison between the results reveals that the macroscopic and the microscopically motivated models predict different behaviours. Therefore, this work encourages a further and

...................................................

60

2.5

0

13

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

(a) 3.5

the programming. All authors gave final approval for publication.

Competing interests. We have no competing interests. Funding. N.C. would like to thank the financial assistance of the Minerva Foundation. Additionally, partial financial support for this work was provided by the Swedish Research Council (Vetenskapsrdet) under grant no. 2011-5428 and is gratefully acknowledged by A.M. Finally, N.C. and G.dB. wish to acknowledge the support of the Israel Science Foundation founded by the Israel Academy of Sciences and Humanities (grant no. 1246/11).

References 1. Bar-Cohen Y. 2001 EAP history, current status, and infrastructure. In Electroactive polymer (EAP) actuators as artificial muscles, ch. 1 (ed. Y Bar-Cohen), pp. 3–44. Bellingham, WA: SPIE Press. 2. McKay T, O’Brien B, Calius E, Anderson L. 2010 An integrated, self-priming dielectric elastomer generator. Appl. Phys. Lett. 97, 062911. (doi:10.1063/1.3478468) 3. Springhetti R, Bortot E, deBotton G, Gei M. 2014 Optimal energy-harvesting cycles for loaddriven dielectric generators in plane strain. IMA J. Appl. Math. 79, 929–946. (doi:10.1093/ imamat/hxu025) 4. Rudykh S, Bhattacharya K, deBotton G. 2012 Snap-through actuation of thick-wall electroactive balloons. Int. J. Nonlinear Mech. 47, 206–209. (doi:10.1016/j.ijnonlinmec.2011. 05.006) 5. Gei M, Springhetti R, Bortot E. 2013 Performance of soft dielectric laminated composites. Smart Mater. Struct. 22, 104014. (doi:10.1088/0964-1726/22/10/104014) 6. Shmuel G, Gei M, deBotton G. 2012 The Rayleigh–Lamb wave propagation in dielectric elastomer layers subjected to large deformations. Int. J. Nonlinear Mech. 47, 307–316. (doi:10.1016/j.ijnonlinmec.2011.06.013) 7. Toupin RA. 1956 The elastic dielectric. J. Ration. Mech. Anal. 5, 849–915. (doi:10.1512/iumj. 1956.5.55033) 8. Kofod G, Sommer-Larsen P, Kornbluh R, Pelrine R. 2003 Actuation response of polyacrylate dielectric elastomers. J. Intell. Mater. Syst. Struct. 14, 787–793. (doi:10.1177/104538903039260) 9. Pelrine R, Kornbluh R, Joseph J, Heydt R, Pei Q, Chiba S. 2000 High-field deformation of elastomeric dielectrics for actuators. Mater. Sci. Eng. C 11, 89–100. (doi:10.1016/S09284931(00)00128-4) 10. Pelrine R, Kornbluh R, Pei Q, Joseph J. 2000 High-speed electrically actuated elastomers with strain greater than 100%. Science 287, 836–839. (doi:10.1126/science.287.5454.836) 11. Huang C, Zhang QM, deBotton G, Bhattacharya K. 2004 All-organic dielectric-percolative three-component composite materials with high electromechanical response. Appl. Phys. Lett. 84, 4391–4393. (doi:10.1063/1.1757632) 12. Stoyanov H, Kollosche M, McCarthy DN, Kofod G. 2010 Molecular composites with enhanced energy density for electroactive polymers. J. Mater. Chem. 20, 7558–7564. (doi:10.1039/ c0jm00519c) 13. Tian L, Tevet-Deree L, deBotton G, Bhattacharya K. 2012 Dielectric elastomer composites. J. Mech. Phys. Solids 60, 181–198. (doi:10.1016/j.jmps.2011.08.005) 14. Galipeau E, Ponte Castañeda P. 2012 The effect of particle shape and distribution on the macroscopic behavior of magnetoelastic composites. Int. J. Solids Struct. 49, 1–17. (doi:10.1016/j.ijsolstr.2011.08.014) 15. Rudykh S, Lewinstein A, Uner G, deBotton G. 2013 Analysis of microstructural induced enhancement of electromechanical coupling in soft dielectrics. Appl. Phys. Lett. 102, 151905. (doi:10.1063/1.4801775) 16. Lopez-Pamies O. 2014 Elastic dielectric composites: theory and application to particle-filled ideal dielectrics. J. Mech. Phys. Solids 64, 61–82. (doi:10.1016/j.jmps.2013.10.016) 17. Bažant P, Oh BH. 1986 Efficient numerical integration on the surface of a sphere. ZAMM-J. Appl. Math. Mech./Z. Angew. Math. Mech. 66, 37–49. (doi:10.1002/zamm.19860660108)

...................................................

Authors’ contributions. All authors carried out the research and the analysis jointly and N.C. was responsible for

14

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

a more rigorous investigation into the relation between the two scales aimed at deepening our understanding of the micro–macro relations and the mechanisms which control the actuation in EAPs. Moreover, analysis of this type may open the path to the design and manufacturing of polymers with microstructures that enable to improve the electromechanical coupling.

15 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

18. Carol I, Jirasek M, Bažant ZP. 2004 A framework for microplane models at large strain, with application to hyperelasticity. Int. J. Solids Struct. 41, 511–557. (doi:10.1016/S00207683(03)00416-5) 19. Miehe C, Göktepe S, Lulei F. 2004 A micro-macro approach to rubber-like materials–part i: the non-affine micro-sphere model of rubber elasticity. J. Mech. Phys. Solids 52, 2617–2660. (doi:10.1016/j.jmps.2004.03.011) 20. Thylander S, Menzel A, Ristinmaa M. 2012 An electromechanically coupled micro-sphere framework: application to the finite element analysis of electrostrictive polymers. Smart Mater. Struct. 21, 094008. (doi:10.1088/0964-1726/21/9/094008) 21. Ogden RW. 1997 Non-linear elastic deformations. New York, NY: Dover Publications. 22. Lu H, Du S. 2014 A phenomenological thermodynamic model for the chemo-responsive shape memory effect in polymers based on Flory–Huggins solution theory. Polym. Chem. 5, 1155–1162. (doi:10.1039/C3PY01256E) 23. Lu H, Huang WM. 2013 On the origin of the Vogel–Fulcher–Tammann law in the thermoresponsive shape memory effect of amorphous polymers. Smart Mater. Struct. 22, 105021. (doi:10.1088/0964-1726/22/10/105021) 24. Kuhn W, Grün F. 1942 Beziehungen zwischen elastischen Konstanten und Dehnungsdoppelbrechung hochelastischer Stoffe. Kolloid-Zeitschrift 101, 248–271. (doi:10. 1007/BF01793684) 25. Wang MC, Guth E. 1952 Statistical theory of networks of non-gaussian flexible chains. J. Chem. Phys. 20, 1144–1157. (doi:10.1063/1.1700682) 26. Flory PJ, Rehner J. 1943 Statistical mechanics of cross-linked polymer networks I. rubberlike elasticity. J. Chem. Phys. 11, 512–520. (doi:10.1063/1.1723791) 27. Treloar LRG. 1946 The elasticity of a network of long-chain molecules.-III. Trans. Faraday Soc. 42, 83–94. (doi:10.1039/tf9464200083) 28. Arruda EM, Boyce MC. 1993 A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. J. Mech. Phys. Solids 41, 389–412. (doi:10.1016/00225096(93)90013-6) 29. Tiersten HF. 1990 A development of the equations of electromagnetism in material continua, vol. 39. Springer Tracts in Natural Philosophy. New York, NY: Springer. 30. Hutter K, van de Ven AAF, Ursescu A. 2006 Electromagnetic field matter interactions in thermoelasic solids and viscous fluids, 2nd edn, vol. 710. Lecture Notes in Physics 710. Berlin, Germany: Springer. 31. Dorfmann A, Ogden RW. 2005 Nonlinear electroelasticity. Acta Mech. 174, 167–183. (doi:10.1007/s00707-004-0202-2) 32. Ask A, Menzel A, Ristinmaa M. 2012 Phenomenological modeling of viscous electrostrictive polymers. Int. J. Nonlinear Mech. 47, 156–165. (doi:10.1016/j.ijnonlinmec.2011.03.020) 33. Ask A, Menzel A, Ristinmaa M. 2012 Electrostriction in electro-viscoelastic polymers. Mech. Mater. 50, 9–21. (doi:10.1016/j.mechmat.2012.01.009) 34. Ask A, Menzel A, Ristinmaa M. 2015 Modelling of viscoelastic dielectric elastomers with deformation dependent electric properties. Procedia IUTAM 12, 134–144. (doi:10.1016/ j.piutam.2014.12.015) 35. Jiménez S, McMeeking RM. 2013 Deformation dependent dielectric permittivity and its effect on actuator performance and stability. Int. J. Nonlinear Mech. 57, 183–191. (doi:10.1016/j.ijnonlinmec.2013.08.001) 36. Cohen N, deBotton G. 2015 The electromechanical response of polymer networks with longchain molecules. Math. Mech. Solids 20, 721–728. (doi:10.1177/1081286514550574) 37. Cohen N, deBotton G. 2014 Multiscale analysis of the electromechanical coupling in dielectric elastomers. Eur. J. Mech. A/Solids 48, 48–59. (doi:10.1016/j.euromechsol.2014.02.022) 38. McMeeking RM, Landis CM. 2005 Electrostatic forces and stored energy for deformable dielectric materials. J. Appl. Mech. Trans. ASME 72, 581–590. (doi:10.1115/1.1940661) 39. McMeeking RM, Landis CM, Jimenez SMA. 2007 A principle of virtual work for combined electrostatic and mechanical loading of materials. Int. J. Nonlinear Mech. 42, 831–838. (doi:10.1016/j.ijnonlinmec.2007.03.008) 40. Gent AN. 1996 A new constitutive relation for rubber. Rubber Chem. Technol. 69, 59–61. (doi:10.5254/1.3538357) 41. Blythe T, Bloor D. 2008 Electrical properties of polymers, 2nd edn. Cambridge, UK: Cambridge University Press.

16 ...................................................

rspa.royalsocietypublishing.org Proc. R. Soc. A 472: 20150462

42. Di Lillo L, Schmidt A, Carnelli DA, Ermanni P, Kovacs G, Mazza E, Bergamini A. 2012 Measurement of insulating and dielectric properties of acrylic elastomer membranes at high electric fields. J. Appl. Phys. 111, 024904. (doi:10.1063/1.3676201) 43. Choi HR et al. 2005 Effects of prestrain on behavior of dielectric elastomer actuator. In Smart structures and materials, pp. 283–291. Bellingham, WA: SPIE. (doi:10.1117/12.599363) 44. Wissler M, Mazza E. 2007 Electromechanical coupling in dielectric elastomer actuators. Sensors Actuators A Phys. 138, 384–393. (doi:10.1016/j.sna.2007.05.029) 45. McKay TG, Calius E, Anderson IA. 2009 The dielectric constant of 3m vhb: a parameter in dispute. In SPIE smart structures and materials + nondestructive evaluation and health monitoring, p. 72870P. International Society for Optics and Photonics. 46. Qiang J, Chen H, Li B. 2012 Experimental study on the dielectric properties of polyacrylate dielectric elastomer. Smart Mater. Struct. 21, 025006. (doi:10.1088/0964-1726/21/2/025006) 47. Flory PJ. 1953 Principles of polymer chemistry. Ithaca, NY: Cornell University Press. 48. Treloar LRG. 1943 The elasticity of a network of long-chain molecules. I. Trans. Faraday Soc. 39, 36–41. (doi:10.1039/tf9433900036) 49. Treloar LRG. 1975 The physics of rubber elasticity. Oxford, UK: Clarendon Press. 50. Kuhl E, Menzel A, Garikipati K. 2006 On the convexity of transversely isotropic chain network models. Philos. Mag. 86, 3241–3258. (doi:10.1080/14786430500080296) 51. Stockmayer WH. 1967 Dielectric dispersion in solutions of flexible polymers. Pure Appl. Chem. 15, 2816. (doi:10.1351/pac196715030539) 52. Menzel A, Waffenschmidt T. 2009 A microsphere-based remodelling formulation for anisotropic biological tissues. Phil. Trans. R. Soc. A 367, 3499–3523. (doi:10.1098/rsta.2009.0103) 53. Alastrué V, Martínez MA, Doblaré M, Menzel A. 2009 Anisotropic micro-sphere-based finite elasticity applied to blood vessel modelling. J. Mech. Phys. Solids 57, 178–203. (doi:10.1016/ j.jmps.2008.09.005) 54. Alastrué V, Martínez MA, Menzel A, Doblaré M. 2009 On the use of non-linear transformations for the evaluation of anisotropic rotationally symmetric directional integrals. application to the stress analysis in fibred soft tissues. Int. J. Numer. Methods Eng. 79, 474–504. (doi:10.1002/nme.2577) 55. Waffenschmidt T, Menzel A, Kuhl E. 2012 Anisotropic density growth of bone-a computational micro-sphere approach. Int. J. Solids Struct. 49, 1928–1946. (doi:10.1016/ j.ijsolstr.2012.03.035) 56. Ostwald R, Bartel T, Menzel A. 2014 A Gibbs-energy-barrier-based computational microsphere model for the simulation of martensitic phase-transformations. Int. J. Numer. Methods Eng. 97, 851–877. (doi:10.1002/nme.4601) 57. Dorfmann A, Ogden RW. 2010 Nonlinear electroelastostatics: incremental equations and stability. Int. J. Eng. Sci. 48, 1–14. (doi:10.1016/j.ijengsci.2008.06.005) 58. Bertoldi K, Gei M. 2011 Instabilities in multilayered soft dielectrics. J. Mech. Phys. Solids 59, 18–42. (doi:10.1016/j.jmps.2010.10.001) 59. Steinmann P, Hossain M, Possart G. 2012 Hyperelastic models for rubber-like materials: consistent tangent operators and suitability for treloar’s data. Arch. Appl. Mech. 82, 1183–1217. (doi:10.1007/s00419-012-0610-z) 60. Rudykh S, deBotton G. 2011 Stability of anisotropic electroactive polymers with application to layered media. Z. Angew. Math. Phys. 62, 1131–1142. (doi:10.1007/s00033-011-0136-1) 61. Goulbourne N, Mockensturm E, Frecker M. 2005 A nonlinear model for dielectric elastomer membranes. J. Appl. Mech. 72, 899–906. (doi:10.1115/1.2047597)