Nonlinear coupled electro-mechanical behavior of a novel anisotropic fiber-reinforced dielectric elastomer

International Journal of Non-Linear Mechanics 119 (2020) 103364

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Nonlinear coupled electro-mechanical behavior of a novel anisotropic fiber-reinforced dielectric elastomer Alireza Ahmadi, Masoud Asgari ∗ Faculty of Mechanical Engineering, K. N. Toosi University of Technology, P. O. Box: 19395-1999, Tehran, Iran

ARTICLE

INFO

Keywords: Fiber-reinforced DE Nonlinear FEM Anisotropic Finite deformation Electroelasticity Blocking force

ABSTRACT Dielectric elastomers are smart materials which produce mechanical actuation when subjected to an electrical field. In this study, we created an efficient framework to model the electromechanical coupling behavior of fiber-reinforced dielectric elastomers using large deformation continuum mechanics, electroelasticity and finite element method. We derived a consistent tangent modulus based on the Cauchy stress definition by incorporating dielectric effect into the tangent modulus. We used the Newton–Raphson method to solve the nonlinear finite element problem at hand with Abaqus and a new user defined material. In the next step, we calibrated the model with experimental tensile tests data on Cotton-reinforced Silicone rubber and check its validity by simulating the exact uniaxial tensile tests and comparing the results to experiments. We further checked the ability of the model to take the dielectric effect into account by comparing simulation results to available experimental data. We showed that the present model can trace the behavior of fiber-reinforced dielectric elastomers with proper accuracy. Also, we presented an example of a fiber-reinforced bending actuator and analyzed the effect of fiber orientation on its actuation and blocking force. We used the Maxwell– Garnett effective medium model to consider the effect of fiber inclusions on the dielectric constant of the material.

1. Introduction Electro-active polymers (EAP) are smart materials capable of producing mechanical output when subjected to external electrical stimuli. EAPs are typically divided into two main groups, Dielectric, and Ionic EAPs. Dielectric elastomers (DE) are the subject of this paper. One of the main drawbacks of DEs is the need for very high voltages, sometimes up to 10 kVs [1]. On the other hand, they do not consume much energy for their actuation [2], they are capable of undergoing large actuation strains, sometimes more than 100% [3], and their mechanical response is quite fast [4,5]. The aforementioned benefits make DEs a good choice for applications like artificial muscles [6,7], sensors [8,9] and other types of actuators like micro pumps [10], pneumatic valves [11], grippers [12,13], tactile actuators [14,15], and underwater vehicles [16]. The most important step in the fabrication of the above-mentioned DE actuators is to investigate whether the material is capable of producing sufficient amount of actuation strain and force. To this end, many researchers investigated the behavior of isotropic DEs in a wide variety of configurations experimentally and analytically. For instance, Guo et al. [17] analyzed a strip-like DE both analytically and experimentally

and revealed that increasing the permittivity and insulation of the elastomer, decreasing the electric insulation of the electrodes and also increasing the intensity of the electric field can lead to electromechanical efficiency of up to 19.8%. When subjected to voltage, DEs experience a decrease in the thickness direction. This decrease causes the electrical field to increase and consequently causes more contraction in the thickness direction. At a specific voltage, this positive feedback causes a tremendous decrease in the thickness direction and makes the DE unstable. This phenomenon is known as the pull-in instability and has been investigated in many previous works [18–20]. One possible way to circumvent this drawback is to use pre-stretched DEs. The effect of pre-stretch on the suppression of the pull-in instability is discussed in some previous works [21,22]. Moreover using pre-stretched elastomers increases the actuation of DEs. As an example, Zhu et al. [23] investigated the effect of pre-stretch on DEs and found that increasing the pre-stretch enhances the actuation and efficiency of DEs. In a more detailed analysis of the pre-stretch effect, He et al. [24] concluded that an optimum amount of pre-stretch leads to a considerable decrease in the amount of drive voltage needed for the actuation. To further enhance the performance of DEs, some researchers investigated fiber-reinforced DEs for specific applications. Lyu and Zhu [25]

∗ Corresponding author. E-mail address: [email protected] (M. Asgari).

https://doi.org/10.1016/j.ijnonlinmec.2019.103364 Received 3 June 2019; Received in revised form 29 September 2019; Accepted 26 November 2019 Available online 29 November 2019 0020-7462/© 2019 Elsevier Ltd. All rights reserved.

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 1. A fiber-reinforced DE patch in the initial and deformed states.

noticed a 13.5% increase in the actuation strain of DE actuator by embedding stiff fibers in the elastomeric matrix. Subramani et al. [26] investigated the electromechanical characteristics of Anisotropic Dielectric Electroactive Polymers with Tunability (ADEPT) and highlighted the great influence of incorporated stiff fibers. Their findings show that ADEPT composites can increase the actuation strain up to almost 35% compared with isotropic DEs. The ratio of the stored mechanical energy to the input electrical energy (known as electromechanical coupling efficiency) was also investigated and an increase of about 30% was discovered by adding fibers to the elastomeric matrix. When DEs are subjected to the electrical field, they expand in the area in both planar directions. In many cases, the desired actuation is in a single direction. In these cases, constraining the actuation in one direction will increase the desired actuation in the other direction due to the incompressibility of the elastomers. This operative and useful idea is the motivation of some studies on development of the fiber-reinforced dielectric elastomers [27–29]. Similar to planar DEs, tubular and spherical DEs are also quite useful actuators and are investigated in many previous works [9,30– 34]. Many of the discussed drawbacks and also applications of the planar DEs are also true for tubular DEs. Fiber stiffening via using fiberreinforced DEs can improve electromechanical performance of tubular actuators as well [25,30,33]. Aside from experimental investigations and analytical modelings, some works are dedicated to the finite element analysis of DEs. The implementation of a finite element model (FEM) for the analysis of the DEs is based on the large deformation continuum mechanics approaches and also on theories of electroelasticity. The concept of electroelasticity is discussed thoroughly in previous papers [35–38]. A large number of previous works aimed at the finite element analysis of the isotropic DEs. For example, Wang et al. [39] used FEM to design and optimize a rotary actuator based on DEs. Li et al. [40] used FEM to design a dielectric elastomer resonator with tunable frequencies and mode shapes. Circular dielectric elastomer actuators were also numerically analyzed by Wissler and Mazza [41]. While isotropic models have been the main focus in finite element analysis of DEs, so far, finite element modeling of fiber-reinforced DEs has been rarely studied by researchers. The work by Sharma and Joglekar [42] on the numerical implementation of the anisotropic dielectric elastomers can be considered one of the pioneers in this field of research. In their research, using the developed numerical framework through a variational approach, the effect of stiff vertical and horizontal fibers on the actuation response and stability of the DEs was studied. The developed framework was further validated by comparing thickness stretch results with an analytical solution. Considering the rapidly growing applications of anisotropic DEs, the present research is aimed at the development of a general finite element framework based on nonlinear continuum mechanics and electroelasticity theories to be capable of modeling fiber-reinforced DEs and implementable in various finite element software. To this end, the paper first proposes a new model for the analysis of anisotropic DEs by

generalizing previous models of fiber-reinforced hyperelastic materials. Then, by using a continuum mechanics approach, the definitions of the Cauchy stress and the consistent tangent modulus including dielectric effect are derived, which lays the groundwork for the numerical implementation of DEs. Moreover, the material model is calibrated by previously published results on uniaxial tensile tests of fiber-reinforced elastomers and the validity of the model is checked by simulation of the tensile tests. We also check the dielectric effect and the validity of the model to implement it by using previously published experimental results on DEs. We then simulate a numerical example of a fiberreinforced DE bending actuator using the proposed framework. Finally, we show that incorporating fibers into an elastomeric bending actuator can lead to almost 15% increase in the actuated displacement and 60% increase in the blocking force. 2. Developing a nonlinear electro-mechanical model for fiberreinforced dielectric elastomers 2.1. Nonlinear continuum mechanics formulation Consider the fiber-reinforced DE patch shown in Fig. 1. The deformation of this continuum is described by 𝑥 = 𝜒(𝑋), where 𝑋 and 𝑥 are a material point position vectors in the reference and deformed states, respectively. 𝜒 is a function which relates the initial and final states of the material. Each fiber has an initial direction vector 𝐴(𝑛) which is changed to 𝑎(𝑛) in the deformed state. The deformation gradient is defined by 𝐹 = 𝑔𝑟𝑎𝑑(𝜒),

(1)

which in the Cartesian coordinate system takes the form 𝐹𝑖𝑗 =

𝜕𝑥𝑖 . 𝜕𝑋𝑗

(2)

The left and the right Cauchy–Green deformation tensors are also defined as 𝐵 = 𝐹𝐹𝑇,𝐶 = 𝐹𝑇𝐹.

(3)

The determinant of the deformation gradient which is a measure of volume change is also defined by 𝐽 = det(𝐹 ).

(4)

For an incompressible material, the ratio of deformed volume to the initial volume, 𝐽 is unity and hence no volume change occurs. The invariants associated with the right Cauchy–Green deformation tensor are defined as: 𝐼1 = 𝑡𝑟(𝐶), 𝐼2 = 𝑡𝑟(𝑐𝑜𝑓 (𝐶)), 𝐼3 = det(𝐶),

(5)

where 𝑡𝑟(𝐶) and 𝑐𝑜𝑓 (𝐶) are the trace and cofactor of tensor C, respectively. Again, if the material is incompressible, 𝐼3 = 𝐽 2 = 1. For anisotropic materials, other sets of invariants corresponding to the fibers are introduced. Assume two sets of fibers with orientation 2

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

unit vectors defined as 𝐴(1) and 𝐴(2) as shown in Fig. 1. Using the definition of the deformation gradient, in the deformed state, the fibers will have an orientation defined as 𝑎(1) = 𝐹 𝐴(1) and 𝑎(2) = 𝐹 𝐴(2) . Spencer [43] defined the invariants associated with the fibers of an anisotropic material as 𝐼4 = 𝐴(1) .𝐶𝐴(1) , 𝐼5 = 𝐴(1) .𝐶 2 𝐴(1) , 𝐼6 = 𝐴(2) .𝐶𝐴(2) , 𝐼7 = 𝐴(2) .𝐶 2 𝐴(2)

Table 1 Tensile properties of Nitrile rubber composites with longitudinal fibers [46].

(6)

Using 𝑎(1) and 𝑎(2) , the above invariants can be simplified as 𝐼4 = 𝑎(1) .𝑎(1) , 𝐼5 = 𝑎(1) .𝐶𝑎(1) , 𝐼6 = 𝑎(2) .𝑎(2) , 𝐼7 = 𝑎(2) .𝐶𝑎(2) .

(7)

The invariants associated with the electric field are also defined as [36] 𝐼8 = 𝐸.𝐸, 𝐼9 = 𝐸.𝐶 −1 𝐸, 𝐼10 = 𝐸.𝐶 −2 𝐸,

Yield elong. %

9.1 9.2 9.2 9.3 9.3 9.3 9.2 9.3 9.2 9.3

30 40 37 77 27 35 25 32 20 27

agent. silica-resorcinol hexamethylene tetramine.

One of the assumptions for the chosen anisotropic part of the strain energy function is that fibers cannot withstand compression and in that case, should be excluded from the analysis. To find the fibers under compression, one needs to compute 𝐼4 and 𝐼6 . If these values are less than unity for a fiber, then the fiber is under compression. To enforce the incompressibility condition, a penalty method is utilized. To this end, first, the deformation gradient is decomposed multiplicatively into isochoric (unimodular) 𝐹 and volumetric 𝐹𝑣𝑜𝑙. parts by

(9)

1

𝐹 = 𝐹 𝐹𝑣𝑜𝑙. → 𝐹𝑣𝑜𝑙. = 𝐽 3 𝐼, 𝐹 = 𝐽

−1 3

𝐹

(13)

where 𝐼 is the second order unit tensor. Then the SEF is additively decomposed into isochoric and volumetric parts, as well

(10) 𝑊 = 𝑊 (𝐼 𝑚 ) + 𝑊𝑣𝑜𝑙. (𝐽 ) (𝑚 = 1, 2, …),

These terms will be discussed in the following section.

(14)

where 𝐼 𝑚 (𝑚 = 1, 2, …) are invariants calculated using the isochoric part of the deformation gradient. The volumetric part is the penalty function which is responsible for incompressibility enforcement. There are different choices for the penalty function in the literature. In this paper, we use the quadratic form [48]

2.2. The utilized model for fiber-reinforced elastomers For each part of the SEF defined by Eq. (10), there are many choices available in the literature. The most prominent models for the isotropic hyperelastic part are the neo-Hookean model, Mooney– Rivlin model and Ogden model. Based on experimental results [44], the neo-Hookean model is applicable for strains below 20%, whereas the Mooney–Rivlin model is proved to be applicable for strains up to almost 200% [45]. Although the Ogden model has proved to be accurate for strains well beyond 200%, its applicability to fiber-reinforced hyperelastic materials has no use, since the present fibers are not capable of undergoing that amount of strain. Almost all of the previous models used for the analysis of fiber-reinforced DEs utilized the neoHookean model for the isotropic part of the strain energy function. But, as presented in Table 1, there are fibers with yield elongations beyond the applicability of the neo-Hookean model, and as a result, there is a need for a more general and accurate model for the analysis of fiber-reinforced DEs. In this regard, the Mooney–Rivlin model is chosen for the isotropic part of the SEF as 𝑊𝑖𝑠𝑜. = 𝐶10 (𝐼1 − 3) + 𝐶01 (𝐼2 − 3),

Fiber vol. %

– HRHb – HRH – HRH – HRH – HRH

b Hydrated

A hyperelastic material is represented by a strain energy function (SEF) which is a function of the invariants that some of them were defined earlier. Although the chosen SEF should be able to correctly model the behavior of the material, it should also possess a degree of convexity (convexity, quasiconvexity, or polyconvexity) to be properly used by finite element codes for boundary value problems. For an anisotropic hyperelastic material, the SEF is additively decomposed into isotropic and anisotropic parts as 𝑊 = 𝑊𝑖𝑠𝑜. + 𝑊𝑎𝑛𝑖𝑠𝑜. .

BAa

Glass Glass Carbon Carbon Cellulose Cellulose Aramid Aramid Nylon Nylon a Bonding

(8)

where 𝐸 is the nominal (Lagrangian) electric field vector. It should be noted that, 𝐼10 is a higher order invariant and is neglected in many studies. The true (Eulerian) electric field vector 𝑒 is simply the push-forward of the nominal electric field vector and is defined by: 𝐸 = 𝐹 𝑇 𝑒 → 𝑒 = 𝐹 −𝑇 𝐸

Fiber

𝐾 (𝐽 − 1)2 , (15) 2 where 𝐾 is the bulk modulus of the material. Using Eqs. (10), (14) and (15), we define the slightly compressible form of the SEF used to model the anisotropic hyperelastic material by

𝑊𝑣𝑜𝑙. =

𝑊 = 𝑊 𝑖𝑠𝑜. + 𝑊 𝑎𝑛𝑖𝑠𝑜. + 𝑊 (𝐽 ) = 𝐶10 (𝐼 1 − 3) + 𝐶01 (𝐼 2 − 3) } ∑ 𝑘1 { 𝐾 + exp[𝑘2 (𝐼 𝑓 − 1)2 ] − 1 + (𝐽 − 1)2 2𝑘 2 2 𝑓 =4,6

(16)

Now, consider a uniform dilation as 𝐹 = 𝜆𝐼. The tensors and invariants are then calculated as 𝑇

𝐽 = 𝜆3 , 𝐹 = 𝐼, C = F F = 𝐼, 𝐼 1 = 𝐼 2 = 3, 𝐼 4 = 𝐼 6 = 1.

(17)

Putting these invariants into Eq. (16) yields 𝐾 3 (𝜆 − 1)2 , (18) 2 which is isotropic, and as a result, the anisotropy of the material is lost. Note that it was expected for all isochoric parts of the SEF including the term 𝐶10 (𝐼 1 − 3) + 𝐶01 (𝐼 2 − 3) to vanish due to the uniform dilation. In order to overcome this drawback, the anisotropic part of the model is modified by Nolan et al. [49] as 𝑊 =

(11)

where 𝐶10 and 𝐶01 are material constants and are found by calibration. Also, 𝐼1 and 𝐼2 are the first and second invariants of the right Cauchy–Green deformation tensor and were defined by Eq. (5). Following Holzapfel et al. [47], the anisotropic part of the strain energy function is chosen to be: ∑ 𝑘1 { } 𝑊𝑎𝑛𝑖𝑠𝑜. = exp[𝑘2 (𝐼𝑓 − 1)2 ] − 1 , (12) 2𝑘 2 𝑓 =4,6

𝑊 = 𝑊 𝑖𝑠𝑜. + 𝑊𝑎𝑛𝑖𝑠𝑜. + 𝑊 (𝐽 ) = 𝐶10 (𝐼 1 − 3) + 𝐶01 (𝐼 2 − 3) ∑ 𝑘1 { } 𝐾 + exp[𝑘2 (𝐼𝑓 − 1)2 ] − 1 + (𝐽 − 1)2 2𝑘2 2 𝑓 =4,6

where 𝑘1 > 0 is a stress-like parameter and 𝑘2 > 0 is a dimensionless parameter to be found by calibration via available published results. 𝐼𝑓 (f = 4, 6) are invariants related to the fibers which were defined by Eqs. (6) and (7).

(19)

In this modified version, the anisotropic part of the strain energy function is calculated using the total right Cauchy–Green deformation tensor instead of its isochoric counterpart. 3

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 2. Fiber-reinforced elastomer test specimen configuration (dimensions in mm).

Fig. 3. Results of tensile tests for Cotton-reinforced Silicone rubber specimens. Source: Data from Peel [50].

as an iterative operator, the final result (if converged) is correct if the definition of the Cauchy stress is correct. In order to derive the definition of stress, namely the Cauchy stress, the first and second derivatives of the strain energy function with respect to strain should be calculated. The first derivative is used in the definition of the second Piola–Kirchhoff stress tensor 𝑆𝑖𝑗 as 𝑆𝑖𝑗 = Fig. 4. Fiber-reinforced DE strip and its element.

𝜕𝑊 𝐶=𝐼+2𝐸 ∗ 𝜕𝑊 ←←←←←←←←←←←←←←←←←←→ ← 𝑆𝑖𝑗 = 2 𝜕𝐸𝑖𝑗∗ 𝜕𝐶𝑖𝑗

= 2( In order to incorporate the dielectric effect into the model, the strain energy function defined by Eq. (19) is augmented by the term 𝑊𝑒𝑙𝑒𝑐. which is defined by [51,52]

𝜕𝑊 𝑖𝑠𝑜. 𝜕𝑊𝑎𝑛𝑖𝑠𝑜. 𝜕𝑊𝑣𝑜𝑙. 𝜕𝑊𝑒𝑙𝑒𝑐. + + + ), 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗

(21)

where 𝐸𝑖𝑗∗ is the Lagrangian strain tensor. Using the definitions for each part of the strain energy function, one can arrive at the following results

−1 𝜀𝐽 𝐶 −1 ∶ [𝐸 ⊗ 𝐸], (20) 2 where 𝜀 is the permittivity of the material. It is worth mentioning that by assuming that the anisotropy in mechanical properties is due to the fiber reinforcement of the elastomer, the electric characteristics due to the electrodes could be assumed independent of mechanical anisotropy, as these are independent parts. Implementing the abovementioned models into a finite element package requires the definitions of stress and tangent modulus. The solution of a nonlinear finite element problem is mostly calculated by the Newton–Raphson method. It is proved that the Newton–Raphson method possesses a quadratic convergence rate provided that the tangent modulus is calculated suitably. If this is not the case, the solution may not converge during the normal amount of iterations allowed by FEM software. It is noteworthy that since the tangent modulus serves

−2 −4 𝜕𝑊 𝑖𝑠𝑜. 𝜕𝐼 𝜕𝐼 𝜕 𝜕 = 𝐶10 1 + 𝐶01 2 = 𝐶10 ( (𝐽 3 𝐼1 )) + 𝐶01 ( (𝐽 3 𝐼2 )), 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗

𝑊𝑒𝑙𝑒𝑐. =

(22) Noting that 𝜕𝐼1 𝜕𝐼2 𝜕𝐽 1 = 𝐽 𝐶𝑖𝑗−1 , = 𝛿𝑖𝑗 , = 𝐼1 𝛿𝑖𝑗 − 𝐶𝑖𝑗 , 𝜕𝐶𝑖𝑗 2 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗

(23)

where 𝛿𝑖𝑗 is the Kronecker delta. By combining Eqs. (22) and (23), we arrive at −2 −1 −4 −2 𝜕𝑊 𝑖𝑠𝑜. = 𝐶10 𝐽 3 [ 𝐼1 𝐶𝑖𝑗−1 + 𝛿𝑖𝑗 ] + 𝐶01 𝐽 3 [ 𝐼2 𝐶𝑖𝑗−1 + 𝐼1 𝛿𝑖𝑗 − 𝐶𝑖𝑗 ] (24) 𝜕𝐶𝑖𝑗 3 3

Similarly for 𝑊𝑎𝑛𝑖𝑠𝑜. , ∑ 𝜕𝑊𝑎𝑛𝑖𝑠𝑜. 𝜕𝐼𝑓 𝜕𝑊𝑎𝑛𝑖𝑠𝑜. = ( )( ) 𝜕𝐶𝑖𝑗 𝜕𝐼𝑓 𝜕𝐶𝑖𝑗 𝑓 =4,6 4

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 5. DE elements with fibers having complementary angles, (a) fibers at 𝜃 < 45◦ , (b) fibers at 𝜃 ′ = 90◦ − 𝜃.

=

∑

) (𝑓 ) 𝑘1 (𝐼𝑓 − 1) exp[𝑘2 (𝐼𝑓 − 1)2 ]𝐴(𝑓 𝑖 𝐴𝑗 ,

(25)

According to Eq. (31), in order to define the tangent modulus, 𝜕𝑆 we need to calculate 𝜕𝐶𝑖𝑗 . This is done in the following section by 𝑘𝑙 separating different parts of the second Piola–Kirchhoff stress tensor similar to Eq. (21)

𝑓 =4,6

where, 𝐴(𝑓 ) (𝑓 = 4, 6) are the orientation vectors in the reference configuration for the family of fibers which are associated with the invariant 𝐼𝑓 . Also, note that, 𝜕𝐼𝑓

𝑖𝑠𝑜.

𝑆𝑖𝑗 = 𝑆 𝑖𝑗 + 𝑆𝑖𝑗𝑎𝑛𝑖𝑠𝑜. + 𝑆𝑖𝑗𝑣𝑜𝑙. + 𝑆𝑖𝑗𝑒𝑙𝑒𝑐.

𝜕 ) 𝜕𝐶𝑘𝑙 (𝑓 ) (𝐴(𝑓 ) 𝐶 𝐴(𝑓 ) ) = 𝐴(𝑓 𝐴𝑙 𝑘 𝜕𝐶 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 𝑘 𝑘𝑙 𝑙 𝑖𝑗 ) 1 ) ) (𝑓 ) (26) = 𝐴(𝑓 [ (𝛿𝑖𝑘 𝛿𝑗𝑙 + 𝛿𝑖𝑙 𝛿𝑗𝑘 )]𝐴(𝑓 = 𝐴(𝑓 𝑖 𝐴𝑗 . 𝑘 𝑙 2 The volumetric and electric parts are also calculated as follows =

𝜕𝑊𝑣𝑜𝑙. 𝜕𝑊𝑣𝑜𝑙. 𝜕𝐽 𝐾 =( )( ) = 𝐽 (𝐽 − 1)𝐶𝑖𝑗−1 , 𝜕𝐶𝑖𝑗 𝜕𝐽 𝜕𝐶𝑖𝑗 2

Taking a derivative of Eq. (32) with respect to the 𝐶𝑘𝑙 yields 𝑖𝑠𝑜.

𝜕𝑆𝑖𝑗 𝜕𝐶𝑘𝑙

=

𝜕𝑆 𝑖𝑗

+

𝜕𝐶𝑘𝑙

𝜕𝑆𝑖𝑗𝑎𝑛𝑖𝑠𝑜. 𝜕𝐶𝑘𝑙

+

𝜕𝑆𝑖𝑗𝑣𝑜𝑙. 𝜕𝐶𝑘𝑙

+

𝜕𝑆𝑖𝑗𝑒𝑙𝑒𝑐. 𝜕𝐶𝑘𝑙

.

−2 { 1 1 −1 1 −1 −1 = 2𝐶10 𝐽 3 𝐼 𝐶 −1 𝐶 −1 − (𝐶𝑘𝑙 𝛿𝑖𝑗 + 𝐶𝑖𝑗−1 𝛿𝑘𝑙 ) + 𝐼1 (𝐶𝑖𝑘 𝐶𝑙𝑗 𝜕𝐶𝑘𝑙 9 1 𝑖𝑗 𝑘𝑙 3 6 } −4 { 4 −1 +𝐶𝑖𝑙−1 𝐶𝑗𝑘 ) + 2𝐶01 𝐽 3 𝐼 𝐶 −1 𝐶 −1 9 2 𝑖𝑗 𝑘𝑙 2 2 −1 −1 − 𝐼1 (𝐶𝑘𝑙 𝛿𝑖𝑗 + 𝐶𝑖𝑗−1 𝛿𝑘𝑙 ) + (𝐶𝑖𝑗−1 𝐶𝑘𝑙 + 𝐶𝑘𝑙 𝐶𝑖𝑗 ) 3 3 } 1 1 −1 −1 −1 −1 + 𝐼2 (𝐶𝑖𝑘 𝐶𝑙𝑗 + 𝐶𝑖𝑙 𝐶𝑗𝑘 ) − (𝛿𝑖𝑘 𝛿𝑗𝑙 + 𝛿𝑖𝑙 𝛿𝑗𝑘 ) 3 2

𝜕𝑆𝑖𝑗𝑖𝑠𝑜.

(34)

𝜎𝑚𝑛 =

It is more convenient to use the definition of the 𝑊𝑎𝑛𝑖𝑠𝑜. in the calculation of the

(29)

𝑓 =4,6

𝜕𝑆𝑖𝑗𝑎𝑛𝑖𝑠𝑜.

+𝜀[𝑒𝑚 𝑒𝑛 − 21 𝛿𝑚𝑛 𝑒𝑙 𝑒𝑙 ]

𝜕𝐶𝑘𝑙

where 𝐵 is the left Cauchy–Green deformation tensor and 𝑎(𝑓 ) (𝑓 = 4, 6) are the fiber orientation vectors in the deformed configuration for the family of fibers associated with the invariant 𝐼𝑓 . Note that the abovementioned procedure for the electric part resulted in the derivation of the well-known Maxwell stress tensor which is the last term in Eq. (29). The next step to implement the model by a finite element code is to define the tangent modulus. Using the consistent form of the tangent modulus will ensure having less problems with convergence issues and is also numerically more economical. The tangent modulus 𝙲 is defined as a tensor which relates the Jaumann rate of the Kirchhoff stress and the rate of deformation tensor 𝐷 as follows ∇(𝐽 𝑎𝑢𝑚𝑎𝑛𝑛)

𝜏

= 𝜏̇ + 𝜏.𝜔 − 𝜔.𝜏 = 𝐽 𝙲 ∶ 𝐷,

𝑎𝑛𝑖𝑠𝑜. 𝜕𝑆𝑖𝑗

𝜕𝐶𝑘𝑙

as

=2

∑ 𝜕 𝜕𝑊𝑎𝑛𝑖𝑠𝑜. 𝜕𝐼𝑓 𝜕 𝜕𝑊𝑎𝑛𝑖𝑠𝑜. ( )=2 ( ) 𝜕𝐶𝑘𝑙 𝜕𝐶𝑖𝑗 𝜕𝐶𝑘𝑙 𝜕𝐼𝑓 𝜕𝐶𝑖𝑗 𝑓 =4,6

=2

2 ∑ 𝜕 2 𝑊𝑎𝑛𝑖𝑠𝑜. 𝜕𝐼𝑓 𝜕𝐼𝑓 𝜕𝑊𝑎𝑛𝑖𝑠𝑜. 𝜕𝐼𝑓 + ) ( 2 𝜕𝐶𝑖𝑗 𝜕𝐶𝑘𝑙 𝜕𝐼𝑓 𝜕𝐶𝑖𝑗 𝜕𝐶𝑘𝑙 𝜕𝐼𝑓 𝑓 =4,6

(35)

The last term in Eq. (35) is proved to be zero, 𝜕𝐼𝑓2 𝜕𝐶𝑖𝑗 𝜕𝐶𝑘𝑙

𝜕𝐶𝑝𝑟 (𝑓 ) 𝜕 𝜕𝐼𝑓 𝜕 ( )= (𝐴(𝑓 ) 𝐴 ) 𝜕𝐶𝑖𝑗 𝜕𝐶𝑘𝑙 𝜕𝐶𝑖𝑗 𝑟 𝜕𝐶𝑘𝑙 𝑝 𝜕 = (𝐴(𝑓 ) 𝐴(𝑓 ) ) = 0 (𝑓 = 4, 6). 𝜕𝐶𝑖𝑗 𝑘 𝑙 =

(36)

This is because 𝐴(𝑓 ) (𝑓 = 4, 6) are orientation vectors in the reference configuration and are constant. Simplifying Eq. (35) leads to

(30)

where 𝜏 and 𝜔 are Kirchhoff stress tensor and spin tensor respectively. This definition is similar to the definition used by Abaqus in its user-defined material models. Following the procedure provided by Nguyen and Waas [53], the following expression can be derived for the definition of 𝙲 according to Eq. (30) 𝙲𝑚𝑛𝑟𝑠 =

(33)

After long but straightforward calculations, each part of the Eq. (33) is calculated as follows

(27)

−1 𝜕𝐶𝑝𝑞 𝜕𝑊𝑒𝑙𝑒𝑐. 1 𝜕𝐽 −1 = − 𝜀𝐸𝑝 𝐸𝑞 [ 𝐶𝑝𝑞 + 𝐽 ] 𝜕𝐶𝑖𝑗 2 𝜕𝐶𝑖𝑗 𝜕𝐶𝑖𝑗 1 −1 −1 −1 −1 −1 − (𝐶𝑝𝑖 𝐶𝑞𝑗 + 𝐶𝑝𝑗 𝐶𝑞𝑖 )]. (28) = − 𝜀𝐽 𝐸𝑝 𝐸𝑞 [𝐶𝑖𝑗−1 𝐶𝑝𝑞 4 Combining Eqs. (24) (25), (27), and (28) along with Eq. (21), the definition of the second Piola–Kirchhoff stress tensor can be found. By pushing this tensor forward, we calculate the Cauchy stress tensor as

1 2 1 𝐹 𝐹 𝑆 = {𝐶 (− 𝐼 𝛿 + 𝐵 𝑚𝑛 ) 𝐽 𝑚𝑖 𝑛𝑗 𝑖𝑗 𝐽 10 3 1 𝑚𝑛 +𝐶01 (− 32 𝐼 2 𝛿𝑚𝑛 + 𝐼 1 𝐵 𝑚𝑛 − 𝐵 𝑚𝑘 𝐵 𝑛𝑘 ) ∑ ) (𝑓 ) + 𝑘1 (𝐼𝑓 − 1) exp[𝑘2 (𝐼𝑓 − 1)2 ]𝑎(𝑓 𝑚 𝑎𝑛 } + 𝐾(𝐽 − 1)𝛿𝑚𝑛

(32)

𝜕𝑆𝑖𝑗𝑎𝑛𝑖𝑠𝑜. 𝜕𝐶𝑘𝑙

=2

∑

) (𝑓 ) (𝑓 ) (𝑓 ) 𝑘1 (1+2𝑘2 (𝐼𝑓 −1)2 ) exp[𝑘2 (𝐼𝑓 −1)2 ]𝐴(𝑓 𝑖 𝐴𝑗 𝐴𝑘 𝐴𝑙 . (37)

𝑓 =4,6

Similarly for the volumetric and electric parts, we arrive at the following equations 𝜕𝑆𝑖𝑗𝑣𝑜𝑙.

𝜕𝑆𝑖𝑗 2 1 𝐹 𝐹 𝐹 𝐹 + (𝜎𝑚𝑟 𝛿𝑛𝑠 + 𝜎𝑛𝑠 𝛿𝑚𝑟 + 𝜎𝑚𝑠 𝛿𝑛𝑟 + 𝜎𝑛𝑟 𝛿𝑚𝑠 ). (31) 𝐽 𝑚𝑖 𝑛𝑗 𝑟𝑘 𝑠𝑙 𝜕𝐶𝑘𝑙 2

𝜕𝐶𝑘𝑙 5

=

1 𝐾𝐽 2

{ } −1 −1 −1 −1 − (𝐽 − 1)[𝐶𝑖𝑘 𝐶𝑙𝑗 + 𝐶𝑖𝑙−1 𝐶𝑗𝑘 ] , (2𝐽 − 1)𝐶𝑖𝑗−1 𝐶𝑘𝑙

(38)

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 6. Experimental and FEM results of tensile tests for a specimen with fibers at 0◦ .

Fig. 7. Experimental and FEM results of tensile tests for a specimen with fibers at ±15◦ .

𝜕𝑆𝑖𝑗𝑒𝑙𝑒𝑐. 𝜕𝐶𝑘𝑙

where 𝑒 is the Eulerian electric field vector as defined by Eq. (9). It should be noted that the CTM defined by Eq. (40) possesses both 𝑚 ↔ 𝑛 and 𝑟 ↔ 𝑠 minor symmetries as expected.

{ 1 −1 −1 −1 −1 −1 −1 = − 𝜀𝐽 𝐸𝑝 𝐸𝑞 𝐶𝑘𝑙 [𝐶𝑖𝑗−1 𝐶𝑝𝑞 − (𝐶𝑝𝑖 𝐶𝑗𝑞 + 𝐶𝑝𝑗 𝐶𝑖𝑞 )] 4

−1 [𝐶 −1 𝐶 −1 + 𝐶 −1 𝐶 −1 ] −𝐶𝑝𝑞 𝑖𝑘 𝑗𝑙 𝑖𝑙 𝑗𝑘 −1 𝐶 −1 + 𝐶 −1 𝐶 −1 ] + 𝐶 −1 [𝐶 −1 𝐶 −1 + 𝐶 −1 𝐶 −1 ] −𝐶𝑖𝑗−1 [𝐶𝑝𝑘 𝑗𝑞 𝑞𝑙 𝑝𝑙 𝑞𝑘 𝑝𝑘 𝑖𝑙 𝑝𝑙 𝑖𝑘 −1 [𝐶 −1 𝐶 −1 +𝐶𝑝𝑖 𝑗𝑘 𝑞𝑙

.

(39)

3. Model calibration using experimental data

−1 ] + 𝐶𝑗𝑙−1 𝐶𝑞𝑘

In order to use the proposed model for the fiber reinforced elastomers, the model coefficients should be calibrated with proper experimental data. Experiments often use the stress (force)-strain (stretch) curves to characterize the behavior of materials being tested. Often, the stress directly measured in experiments is the nominal stress which is the force divided by unit cross-sectional area in the reference configuration. The calibration process for hyperelastic material models was thoroughly discussed by Ogden et al. [54] for different kinds of experiments. An example of such calibration for uniaxial tension of isotropic silicone rubber strip can be found in Amabili et al. [55]. For simple tension of an incompressible material, the nominal stress is calculated as 𝑑 𝑊̂ 𝑡= , (41) 𝑑𝜆

} −1 𝐶 −1 + 𝐶 −1 𝐶 −1 ] + 𝐶 −1 [𝐶 −1 𝐶 −1 + 𝐶 −1 𝐶 −1 ] + 𝐶𝑖𝑞−1 [𝐶𝑝𝑘 𝑝𝑗 𝑗𝑙 𝑝𝑙 𝑗𝑘 𝑖𝑘 𝑞𝑙 𝑖𝑙 𝑞𝑘 Finally, using Eq. (31) along with Eqs. (34), (37), (38) and (39) yields the definition of the consistent tangent modulus (CTM) for an anisotropic dielectric elastomer 4 4 (𝐶 𝐼 + 4𝐶01 𝐼 2 )𝛿𝑚𝑛 𝛿𝑟𝑠 − (𝐶 + 2𝐶01 𝐼 1 )(𝐵 𝑚𝑛 𝛿𝑟𝑠 + 𝐵 𝑟𝑠 𝛿𝑚𝑛 ) 9𝐽 10 1 3𝐽 10 2 8 + (𝐶10 𝐼 1 + 2𝐶01 𝐼 2 )(𝛿𝑚𝑟 𝛿𝑛𝑠 + 𝛿𝑚𝑠 𝛿𝑛𝑟 ) + 𝐶 (𝐵 𝐵 𝛿 + 𝐵 𝑟𝑘 𝐵 𝑠𝑘 𝛿𝑚𝑛 ) 3𝐽 3𝐽 01 𝑚𝑘 𝑛𝑘 𝑟𝑠 ∑ 4 4 + 𝐶01 (𝐵 𝑚𝑛 𝐵 𝑟𝑠 − 𝐵 𝑚𝑟 𝐵 𝑛𝑠 ) + 𝑘 (1 + 2𝑘2 (𝐼𝑓 − 1)2 ) 𝐽 𝐽 𝑓 =4,6 1

𝙲𝑚𝑛𝑟𝑠 =

) (𝑓 ) (𝑓 ) (𝑓 ) × exp[𝑘2 (𝐼𝑓 − 1)2 ]𝑎(𝑓 𝑚 𝑎𝑛 𝑎𝑟 𝑎𝑠

,

+𝐾[(2𝐽 − 1)𝛿𝑚𝑛 𝛿𝑟𝑠 − (𝐽 − 1)(𝛿𝑚𝑟 𝛿𝑛𝑠 + 𝛿𝑚𝑠 𝛿𝑛𝑟 )] 1 − 𝜀𝑒𝑘 𝑒𝑘 (𝛿𝑚𝑛 𝛿𝑟𝑠 − (𝛿𝑚𝑟 𝛿𝑛𝑠 + 𝛿𝑚𝑠 𝛿𝑛𝑟 )) 2 +𝜀(𝛿𝑟𝑠 𝑒𝑚 𝑒𝑛 + 𝛿𝑚𝑛 𝑒𝑟 𝑒𝑠 − 𝛿𝑚𝑠 𝑒𝑟 𝑒𝑛 − 𝛿𝑚𝑟 𝑒𝑠 𝑒𝑛 − 𝛿𝑛𝑠 𝑒𝑚 𝑒𝑟 − 𝛿𝑟𝑛 𝑒𝑠 𝑒𝑚 ) 1 + (𝜎𝑚𝑟 𝛿𝑛𝑠 + 𝜎𝑛𝑠 𝛿𝑚𝑟 + 𝜎𝑚𝑠 𝛿𝑛𝑟 + 𝜎𝑛𝑟 𝛿𝑚𝑠 ) 2

−1

where 𝑊̂ = 𝑊 (𝜆, 𝜆 2 ) is the strain energy function of the material and −1 𝜆1 = 𝜆, 𝜆2 = 𝜆3 = 𝜆 2 are the principal stretches. In the present study, we calibrate the material model using experimental stress–stretch data on simple uniaxial tension of Cotton-reinforced Silicone rubber specimens with fibers at different directions. As stated before, the role of the strain energy function is to exhibit the real behavior of the elastomer

(40) 6

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 8. Experimental and FEM results of tensile tests for a specimen with fibers at ±30◦ .

Fig. 9. Experimental and FEM results of tensile tests for a specimen with fibers at ±45◦ .

as precise as possible. To this end, this research calibrates the model for each fiber angle separately instead of using multi-objective curvefitting. During multi-objective curve-fitting, a set of coefficients is found that gives a relatively small error for each fiber angle while calibrating the model for each fiber angle separately will give the most identical curve to the experimental curve. In other words, 𝑘1 and 𝑘2 in Eq. (12) are assumed to be a function of fiber angle. The calibration process starts by deriving an equation for nominal stress as a function of axial stretch. During a uniaxial tension, the deformation gradient can be written as ⎡𝜆 ⎢ ⎢0 𝐹 =⎢ ⎢ ⎢0 ⎣

0 1 √ 𝜆 0

0 ⎤ ⎥ 0 ⎥ ⎥ 1 ⎥ √ ⎥ 𝜆⎦

𝐼2 = 2𝜆 +

0 1 𝜆 0

(42)

0⎤ ⎥ 0⎥. 1 ⎥⎥ 𝜆⎦

1 sin2 (𝜃𝑓 )] (47) 𝜆2 In order to calibrate the model, first the Mooney–Rivlin coefficients are found by fitting Eq. (47) with the uniaxial stress–stretch curve of an elastomeric matrix and then by using fiber-reinforced specimens, the coefficients 𝑘1 and 𝑘2 are found. The fitting is done by minimizing the following objective function × [2𝜆 cos2 (𝜃𝑓 ) −

(43)

𝜒2 =

Then, using Eqs. (5) and (6), the invariants 𝐼1 , 𝐼2 , 𝐼4 and 𝐼6 will be easily calculated as 𝐼 1 = 𝜆2 +

2 𝜆

(45)

2 0 0⎤⎧ ⎧cos(𝜃 )⎫ ⎡𝜆 ⎫ ⎥ ⎪cos(𝜃𝑓 )⎪ 𝑓 ⎪ ⎢ 1 ⎪ 0 0 ⎢ ⎥ 𝐼4 = 𝐼6 = 𝐴𝑓 .𝐶𝐴𝑓 = ⎨ sin(𝜃𝑓 ) ⎬ . ⎨ sin(𝜃𝑓 ) ⎬ 𝜆 1 ⎥⎥ ⎪ 0 ⎪ ⎪ 0 ⎪ ⎢⎢ 0 ⎩ ⎭ ⎣0 ⎩ ⎭ 𝜆⎦ 1 = 𝜆2 cos2 (𝜃𝑓 ) + sin2 (𝜃𝑓 ) (46) 𝜆 where 𝜃𝑓 (𝑓 = 4, 6) is the angle of fibers corresponding to the invariant 𝐼𝑓 with respect to the tension direction. Using Eqs. (44)–(46) along with Eqs. (10)–(12) and (41), the definition of the nominal stress during uniaxial tension can be derived as ∑ 1 1 𝑡 = 2𝐶10 (𝜆 − ) + 2𝐶01 (1 − )+ 𝑘1 (𝐼𝑓 − 1) exp{𝑘2 (𝐼𝑓 − 1)2 } 𝜆2 𝜆3 𝑓 =4,6

The corresponding right Cauchy–Green deformation tensor will be calculated as ⎡ 𝜆2 ⎢ 𝐶 = 𝐹𝑇𝐹 = ⎢ 0 ⎢ ⎢0 ⎣

1 𝜆2

𝑛 ∑ (𝑡 − 𝑡exp . )2𝑖

(48)

𝑖=1

where 𝑡exp . is the experimental nominal stress and 𝑛 is the number of data points. To preserve the convexity and stability of the strain energy

(44) 7

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

characterize the behavior of the material containing fibers with angles greater than 45◦ . In this section, we show that for dielectric elastomers, this drawback can be overcome by noting the fact that unlike uniaxial test specimens which undergo lateral compression during axial tension, a DE element stretches in all planar directions. Consider a strip of fiber-reinforced DE with fiber angles ±𝜃 as shown in Fig. 4 and take an element out of it which contains two separate fibers. Now consider two different configurations of the DE element having complementary fiber angles. When a voltage is applied to these elements, they will exhibit horizontal and vertical stretches which are named (𝜆1 , 𝜆2 ) and (𝜆′1 , 𝜆′2 ) as shown in Fig. 5. If the configuration (a) is rotated 90◦ clockwise or counterclockwise, the configuration (b) is achieved and the stretches are equal in such a way that 𝜆′1 = 𝜆2 and 𝜆′2 = 𝜆1 . Hence, it can be concluded that these two configurations are actually the same. As a result, the material properties of a DE element containing fibers at ±𝜃 are exactly the same as of a DE element containing fibers at ±(90 − 𝜃). Again it should be emphasized that the reasoning used here is allowed only because a DE element subjected to voltage is under positive stretch in all planar directions. 4. Validation of the implemented model The validation process for the derived formulations consists of two steps. First, the isotropic and anisotropic parts of the strain energy function are validated by simulating tensile tests which were used earlier to calibrate the model and we show that the finite element results have good agreement with experimental data, and that the model works correctly. In the second step, the validity of the model to incorporate the dielectric effect is checked by simulating an isotropic DE actuator and comparing the results to experimental data. In this way, the contributions of different terms to the strain energy function are checked and since no coupling effects between the dielectric and anisotropic parts are assumed, the whole formulation can be considered validated. The validation proves the correctness of the derived formulations in Section 2.2 and their implementation in the finite element software is justified. Note that due to lack of proper experimental data on fiber-reinforced DEs, we used a part-by-part basis for the validation of the formulations. It would be better and more accurate to use experimental data on the fiber-reinforced material activated electrically. In that way the coupling effects between the dielectric and anisotropic parts in the SEF would also be taken into account.

Fig. 10. Schematic of experimental setup [56], (a) Pre-stretched state, (b) Actuated state. Table 2 Calibration results of 𝑘1 and 𝑘2 for each fiber angle. Fiber angle ◦

0 ◦ ±15 ◦ ±30 ◦ ±45

𝑘1 (MPa) 78.93 12.08 3.07 1.88

𝑘2 ≈ 0.1

function, 𝑘1 and 𝑘2 should both be positive [47]. Similarly, 𝐶01 ≥ 0 and 𝐶10 + 𝐶01 ≥ 0 should also be satisfied during the curve-fitting procedure [45]. The experimental data for uniaxial tensile tests of Cotton-reinforced Silicone rubber specimens were taken from Peel [50] and used for model calibration. Fig. 2 shows the test specimen configuration for the fiber-reinforced elastomer. Note that the fiber volume fraction is 51.8%. The fiber angles with respect to the tension direction are 0◦ , ±15◦ , ±30◦ and ±45◦ . The results of tensile tests are shown in Fig. 3. As it is clear in this figure specimens with small fiber angles are relatively stiffer compared with other specimens. It should be noted that the results which are used for material calibration and depicted in Fig. 3 are limited to stretches before the onset of fiber failure which was characterized by a softening in the stress–stretch curve of the specimens. In the next step, the model coefficients are calibrated for the elastomeric matrix i.e. Silicone rubber and the Mooney–Rivlin parameters are found to be 𝐶10 = 0.12780 MPa and 𝐶01 = 0.00139 MPa. Then, by using these parameters, the 𝑘1 and 𝑘2 coefficients are calibrated for each fiber angle and the results are shown in Table 2. We found that the 𝑘2 coefficient is almost 0.1 for all fiber angles and can be considered constant. For other fiber angles, one can simply find the 𝑘1 coefficient by interpolation. It should be emphasized that the fiber orientations used for model calibration in this section were limited to the angles less than 45◦ . This is because for fiber angles greater than 45◦ , one can simply use Eq. (46) and show that 𝐼4 and 𝐼6 are less than unity at some stretches and the fibers are compressed during axial tension of the specimen due to lateral compression. Hence, uniaxial tests are not sufficient to

4.1. Model validation without dielectric effect Using the specimen configuration in Fig. 2 and material properties calculated in Section 3, the tensile tests are simulated and the results are compared to experimental data. We use 12 000, 8-node linear brick elements with reduced integration and enhanced hourglass control options. Also, the bulk modulus in Eqs. (29) and (40) was set to 𝐾 = 100 × 𝜇 = 100 × 2(𝐶10 + 𝐶01 ), where 𝜇 is the initial shear modulus of the elastomer. Note that a higher amount of bulk modulus leads to less compressibility and more accurate results, but may cause convergence problems. For more information about choosing the bulk modulus, the readers are referred to [57]. The results of the simulation for different specimens are shown in Figs. 6–9. Note that the colors on specimens depict the amplitude of the Cauchy stress in tensile direction at the final stretch. As it is shown, the model is able to follow the real behavior of the Cotton-reinforced Silicone rubber with reasonable accuracy. On the other hand, we mention that no convergence issues were noticed during the simulation and the consistency of the derived tangent modulus can be concluded. It should be underscored that here, simulation of tensile tests for verification of the calibration procedure was done due to lack of proper experimental data on fiber-reinforced elastomers, and should be avoided if suitable data are available. 8

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 11. Pre-stretched state, (a) Experiment [56] (b) FEM.

used 1200 elements to model the DE strip. Figs. 11 and 12 show the DE strip before and after actuation. Note that the colors on the simulated strip depict the amplitude of the Cauchy stress in the vertical direction. Fig. 13 shows the simulation results of the DE strip actuation. As it is clear in this figure, the numerical results have good accuracy when compared to experimental data. For stretches close to 3.0 (200% strain), the numerical results have a very small deviation from experiments. This is because as stated in Section 2.2, the Mooney–Rivlin model cannot accurately capture the real behavior of hyperelastic materials for strains around and more than 200%. Note that the initial vertical stretch was previously set to 𝜆𝑝1 = 2.36 using a dead load. Note also that anisotropic mechanical properties are assumed to be related to the fiber reinforcement in the elastomer core. So, the electric characteristics due to the electrodes could be assumed independent of mechanical anisotropy as these are independent parts. Based on this point, the coupling between elastic anisotropy and electric properties has not been taken to account.

4.2. Model validation with dielectric effect In order to show the ability of the model to take the dielectric effect into account, a DE actuation was simulated and the results were compared to the experimental data obtained by Wang et al. [56]. A piece of 3M™ VHB 4905 strip with initial dimensions 𝐿1 = 5, 𝐿2 = 60, 𝑡 = 0.5 mm was first pre-stretched two times in the horizontal direction so that its length was changed to 𝜆𝑝2 𝐿2 = 120 mm. Then a dead load 𝑃1 was added in the vertical direction to produce an initial stretch of 𝜆𝑝1 = 2.36. In the next step, a 5 kV voltage was applied to the electrodes on the sides of the DE strip and the vertical actuation was monitored. Fig. 10 shows the schematic of the experimental setup of the DE strip in the pre-stretched and actuated states. Note that as derived earlier, the model works with the electrical field and not the voltage. The nominal electrical field 𝐸 is related to the applied voltage 𝑉 by: 𝐸=

𝑉 , 𝑑0

(49) 5. Effect of anisotropy on fiber-reinforced DE actuation

where 𝑑0 is the distance between electrodes in the reference configuration. We emphasize that even when DEs are pre-stretched, Eq. (49) holds true and the reference configuration should be utilized to define the distance between electrodes. This is because as shown by Eq. (9), the true electrical field is automatically calculated by the deformation gradient and the effect of pre-stretch on the electrical field is taken into account. The material parameters for the VHB 4905 strip are found to be 𝐶10 = 0.016 MPa and 𝐶01 = 0.0073 MPa [58] and the permittivity of F the material is set to 𝜀 = 𝜀𝑟 𝜀0 = (4.7) × (8.85 × 10−12 m ) [59], where 𝜀𝑟 and 𝜀0 are relative and vacuum permittivities, respectively. Similar to previous simulations, the bulk modulus was also set to 𝐾 = 100𝜇. We

In this section, we analyze the effects of fiber angles on the actuation of a DE bending actuator. To this end, a piece of Cotton-reinforced Silicone rubber strip with dimensions 100 × 10 × 0.25 mm was electroded on both sides and attached to an un-electroded Silicone rubber strip with the same dimensions, as depicted in Fig. 14. Then, the electrodes were subjected to a 4 kV voltage. After applying the voltage, the electroded strip tends to expand in area while the un-electroded strip does not. In this way, both strips bend to produce the desired expansion for the electroded strip. This is how a bending DE actuator works. The material properties for Cotton-reinforced Silicone rubber were calculated in Section 3. To account for the effect of fibers on the 9

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 12. Actuated state, (a) Experiment [56] (b) FEM.

Fig. 13. Comparison of numerical and experimental results for DE vertical actuation.

was earlier said to be 51.8%, the effective relative permittivity of the Cotton-reinforced Silicone rubber composite is calculated to be 3.23. Note that dielectric constants vary with humidity, temperature, frequency [63] and also deformation [64]. In this regard, the readers are encouraged to use more accurate data to calculate the effective permittivity of the materials. The simulation results for the actuation of the simple and fiberreinforced DE bending actuators are shown in Fig. 15 with contours showing the magnitude of the von Mises stress. As shown in this figure, and as predicted, small fiber angles yield smaller actuations in the vertical direction. This is because the bending deformation is a result

dielectric constant of the material, we use the approximate Maxwell– Garnett model [60]. In this model, the effective dielectric constant of the composite 𝜀𝑒𝑓 𝑓 is calculated based on the dielectric constants of the matrix 𝜀𝑚 , inclusions 𝜀𝑖 , and the volume fraction of the inclusions 𝛿𝑖 as follows: 𝜀𝑒𝑓 𝑓 − 𝜀𝑚 𝜀 − 𝜀𝑚 ( ) = 𝛿𝑖 ( 𝑖 ) (50) 𝜀𝑒𝑓 𝑓 + 2𝜀𝑚 𝜀𝑖 + 2𝜀𝑚 Using data from the literature, the relative permittivity of Silicone rubber and Cotton are found to be 3.3 [61] and 3.18 [62], respectively. Using these values, and also the fibers volume fraction which 10

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 17. Deformation of simple and fiber-reinforced constrained actuators, (a) Base state, (b) Fibers at 0◦ , (c) Fibers at ±45◦ , (d) No fibers, (e) Fibers at 90◦ .

of the axial expansion in the upper strip. Small fiber angles are stiff enough in the axial direction, so that they allow small axial stretch and consequently smaller bending deformations. To grasp the effect of adding Cotton fibers to the Silicone rubber strip, we depict distribution of the vertical deformation on the DE strip in Fig. 16. We find that adding fibers at 90◦ caused a 15.1% increase in the vertical actuation compared with fiberless DE. This is because in this configuration, the lateral extension is more constrained compared

Fig. 14. Schematic of a fiber-reinforced DE bending actuator.

Fig. 15. Bending actuation of simple and fiber-reinforced DEs, (a) Base state, (b) Fibers at 0◦ , (c) Fibers at ±45◦ , (d) No fibers, (e) Fibers at 90◦ .

Fig. 16. Distribution of the vertical deformation for simple and fiber-reinforced DE strips.

11

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364

Fig. 18. Blocking force for the simple and fiber-reinforced DE actuators.

with simple DE bending actuator and this leads to higher stretches in the horizontal direction and consequently, the actuator bends more.

6. Conclusion

The blocking force is an important characteristic of every actuator, especially for bending actuators that have strong potential as biomimetic grippers. In simple words, the blocking force is the amount of force an actuator can exert when it is constrained in the direction of actuation. For the bending actuator under investigation, Fig. 17(a) shows the constraint that produces the blocking force. Fig. 17(b–e) depict the deformation of different actuators when they are constrained in the vertical direction at the tip and are subjected to a 4 kV voltage. Without a constraint at the tip, the strips would bend downward as was previously depicted in Fig. 15, but when the constraint is applied, the strips deform to a different curved configuration. This is because unlike the lower strip, the upper strip in the DE which is subjected to a voltage is expanding in the area. This causes the actuator to bend so that the upper strip has the opportunity to expand. When the vertical displacement at the tip of the strip is constrained, it does not change the behavior of the DE actuator. The upper strip still wants to expand and this explains the curved configurations seen in Fig. 17. Simulations show that actuators with small and moderate fiber angles (i.e. 0◦ and ±45◦ ) experience a very small deformation. For other actuators (i.e. fiberless DE and the actuator with fibers at 90◦ ), a significant bending deformation occurs. These behaviors are similar to the unconstrained deformations which were discussed previously.

This research presented a validated platform for finite element analysis of simple and fiber-reinforced dielectric elastomers. In order to implement the behavior of a hyperelastic material in a finite element code, the definition of Cauchy stress and consistent tangent modulus is desired. While the correctness of the Cauchy stress definition will ensure the correctness of results, an accurate tangent modulus will ensure least convergence issues during implementation. By adding the dielectric effect in the tangent modulus, we derived a consistent tangent modulus for the model at hand. On the other hand, a new idea for fiber-specific model calibration was also presented and its applicability was tested by comparing the simulation results to experimental data on uniaxial tests. If proper experimental data on fiber-reinforced elastomers are available, the readers are encouraged to use those data to verify the calibration procedure. This would be more accurate than the method used in this paper. We also checked the ability of the derived formulations to correctly trace the behavior of a DE actuator with pre-stretch and compared it to experimental results. Finally, an example of a fiber-reinforced DE bending actuator was simulated. The effect of fiber inclusions on the dielectric constant of the material was taken into account by using the Maxwell–Garnett effective medium model. The simulation results show that the DE with fibers at 90◦ shows an almost 15.1% increase in the actuation distance and 60% increase in blocking force compared with fiberless DE.

To compare the blocking force of the actuators, the applied voltage is increased from zero to 4 kV and the blocking force is reported as a function of voltage in Fig. 18. As shown in this figure, actuators with fibers at 0◦ and ±45◦ have smaller blocking forces. To find the reason for this trend, attention should be paid to the change in the geometry of the strips. When voltage is increased, a curved beam is formed and the simulations show that the radii of curvature in these curved beams are directly proportional to the applied voltage. In a previous research [65], it was demonstrated that the stiffness of a curved beam is proportional to the radius of curvature and its powers. Hence, it can be concluded that those beams that are bent more exert higher blocking forces because their stiffness is higher.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Ethical Statement: Article was conducted according to ethical standards References [1] K. Wang, G. Ouyang, X. Chen, H. Jakobsen, Engineering electroactive dielectric elastomers for miniature electromechanical transducers, Polymer Reviews 57 (2016) 369–396. [2] J. Biggs, K. Danielmeier, J. Hitzbleck, J. Krause, T. Kridl, S. Nowak, E. Orselli, X. Quan, D. Schapeler, W. Sutherland, J. Wagner, Electroactive polymers: developments of and perspectives for dielectric elastomers, Angew. Chem. Int. Ed. Engl. 52 (2013) 9409–9421. [3] X. Zhao, Z. Suo, Theory of dielectric elastomers capable of giant deformation of actuation, Phys. Rev. Lett. 104 (2010) 178302. [4] R. Pelrine, R. Kornbluh, Q. Pei, J. Joseph, High-speed electrically actuated elastomers with strain greater than 100%, Science 287 (2000) 836–839. [5] K. Jia, T. Lu, T.J. Wang, Response time and dynamic range for a dielectric elastomer actuator, Sensors Actuators A 239 (2016) 8–17.

Although the fibers at small orientation angles increase the stiffness of the structure, we emphasize that the stiffness is a function of both material properties and the geometry. The large bending deformations occurred in the fiberless DE and the one with 90◦ fibers increased their bending stiffness and as a result, their final stiffness is even higher than the actuators with fibers at 0◦ and ±45◦ . The simulation results show that the DE with 90◦ fibers shows an almost 60% increase in the blocking force compared with fiberless DE. 12

A. Ahmadi and M. Asgari

International Journal of Non-Linear Mechanics 119 (2020) 103364 [36] A. Dorfmann, R. Ogden, Nonlinear electroelasticity, Acta Mech. 174 (2005) 167–183. [37] C. Trimarco, Material electromagnetic fields and material forces, Arch. Appl. Mech. 77 (2007) 177–184. [38] C. Trimarco, On the Lagrangian electrostatics of elastic solids, Acta Mech. 204 (2009) 193. [39] N. Wang, H. Guo, B. Chen, C. Cui, X. Zhang, Design of a rotary dielectric elastomer actuator using a topology optimization method based on pairs of curves, Smart Mater. Struct. 27 (2018). [40] Y. Li, I. Oh, J. Chen, Y. Hu, A new design of dielectric elastomer membrane resonator with tunable resonant frequencies and mode shapes, Smart Mater. Struct. 27 (2018). [41] M. Wissler, E. Mazza, Modeling of a pre-strained circular actuator made of dielectric elastomers, Sensors Actuators A 120 (2005) 184–192. [42] A.K. Sharma, M. Joglekar, A numerical framework for modeling anisotropic dielectric elastomers, Comput. Methods Appl. Mech. Engrg. 344 (2019) 402–420. [43] A.J.M. Spencer, Constitutive theory for strongly anisotropic solids, in: Continuum Theory of the Mechanics of Fibre-Reinforced Composites, Springer, 1984, pp. 1–32. [44] C.W. Macosko, R.G. Larson, Rheology: principles, measurements, and applications, 1994. [45] N. Kumar, V.V. Rao, Hyperelastic Mooney-Rivlin model: determination and physical interpretation of material constants, Parameters 2 (2016) 01. [46] J. O’connor, Short-fiber-reinforced elastomer composites, Rubber Chem. Technol. 50 (1977) 945–958. [47] G.A. Holzapfel, T.C. Gasser, R.W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elast. Phys. Sci. Solids 61 (2000) 1–48. [48] J. Simo, R. Taylor, Penalty function formulations for incompressible nonlinear elastostatics, Comput. Methods Appl. Mech. Engrg. 35 (1982) 107–118. [49] D. Nolan, A. Gower, M. Destrade, R. Ogden, J. McGarry, A robust anisotropic hyperelastic formulation for the modelling of soft tissue, J. Mech. Behav. Biomed. Mater. 39 (2014) 48–60. [50] L.D. Peel, Fabrication and mechanics of fiber-reinforced elastomers (Ph.D. thesis), Department of mechanical engineering, Brigham Young University, Utah, United States, 1998. [51] A. Büschel, S. Klinkel, W. Wagner, Dielectric elastomers–numerical modeling of nonlinear visco-electroelasticity, Internat. J. Numer. Methods Engrg. 93 (2013) 834–856. [52] A. Ask, A. Menzel, M. Ristinmaa, On the modelling of electro-viscoelastic response of electrostrictive polyurethane elastomers, in: IOP Conference Series: Materials Science and Engineering, IOP Publishing, 2010, p. 012101. [53] N. Nguyen, A.M. Waas, Nonlinear, finite deformation, finite element analysis, Z. Angew. Math. Phys. 67 (2016) 35. [54] R. Ogden, G. Saccomandi, I. Sgura, Fitting hyperelastic models to experimental data, Comput. Mech. 34 (2004) 484–502. [55] M. Amabili, P. Balasubramanian, I.D. Breslavsky, G. Ferrari, R. Garziera, K. Riabova, Experimental and numerical study on vibrations and static deflection of a thin hyperelastic plate, J. Sound Vib. 385 (2016) 81–92. [56] Y. Wang, B. Chen, Y. Bai, H. Wang, J. Zhou, Actuating dielectric elastomers in pure shear deformation by elastomeric conductors, Appl. Phys. Lett. 104 (2014) 064101. [57] D. Dassault Systèmes, Abaqus analysis user’s guide, in, Technical Report Abaqus 6.14 Documentation, Simulia Corp, 2016. [58] S. Son, N. Goulbourne, Dynamic response of tubular dielectric elastomer transducers, Int. J. Solids Struct. 47 (2010) 2672–2679. [59] T. Vu-Cong, C. Jean-Mistral, A. Sylvestre, Impact of the nature of the compliant electrodes on the dielectric constant of acrylic and silicone electroactive polymers, Smart Mater. Struct. 21 (2012) 105036. [60] O. Levy, D. Stroud, Maxwell Garnett theory for mixtures of anisotropic inclusions: Application to conducting polymers, Phys. Rev. B 56 (1997) 8035. [61] C. Löwe, X. Zhang, G. Kovacs, Dielectric elastomers in actuator technology, Adv. Eng. Mater. 7 (2005) 361–367. [62] F.S.C. Mustata, A. Mustata, Dielectric behaviour of some woven fabrics on the basis of natural cellulosic fibers, Adv. Mater. Sci. Eng. 2014 (2014). [63] R. Patel, N. Gupta, Effect of humidity on the complex permittivity of epoxy-based nanodielectrics with metal oxide fillers, Int. Trans. Electr. Energy Syst. 23 (2013) 846–852. [64] S.M. Jiménez, R.M. McMeeking, A constitutive law for dielectric elastomers subject to high levels of stretch during combined electrostatic and mechanical loading: Elastomer stiffening and deformation dependent dielectric permittivity, Int. J. Non-Linear Mech. 87 (2016) 125–136. [65] R. Palaninathan, P. Chandrasekharan, Curved beam element stiffness matrix formulation, Comput. Struct. 21 (1985) 663–669.

[6] L. Liu, C. Zhang, M. Luo, X. Chen, D. Li, H. Chen, A biologically inspired artificial muscle based on fiber-reinforced and electropneumatic dielectric elastomers, Smart Mater. Struct. 26 (2017). [7] H. Wang, S. Qu, Constitutive models of artificial muscles: a review, J. Zhejiang Univ. Sci. A 17 (2016) 22–36. [8] A. Khan, F.R. Khan, H.S. Kim, Electro-active paper as a flexible mechanical sensor, Actuator Energy Harvest. Transducer: Rev. Sensors (Basel) 18 (2018). [9] J. Chavanne, Y. Civet, Y. Perriard, Modelling of a dielectric electroactive polymer tubular shape sensor for pressure measurements, 2016. [10] F.A. Mohd Ghazali, C.K. Mah, A. AbuZaiter, P.S. Chee, M.S. Mohamed Ali, Soft dielectric elastomer actuator micropump, Sensors Actuators A 263 (2017) 276–284. [11] M. Hill, G. Rizzello, S. Seelecke, Development and experimental characterization of a pneumatic valve actuated by a dielectric elastomer membrane, Smart Mater. Struct. 26 (2017). [12] K.-R. Heng, A.S. Ahmed, M. Shrestha, G.-K. Lau, Strong dielectric-elastomer grippers with tension arch flexures, in: Electroactive Polymer Actuators and Devices (EAPAD) 2017, International Society for Optics and Photonics, 2017, p. 101631Z. [13] S. Shian, K. Bertoldi, D.R. Clarke, Dielectric elastomer based grippers for soft robotics, Adv. Mater. 27 (2015) 6814–6819. [14] S. Mun, S. Yun, S. Nam, S.K. Park, S. Park, B.J. Park, J.M. Lim, K.U. Kyung, Electro-active polymer based soft tactile interface for wearable devices, IEEE Trans. Haptics 11 (2018) 15–21. [15] H. Phung, C.T. Nguyen, T.D. Nguyen, C. Lee, U. Kim, D. Lee, J.-d. Nam, H. Moon, J.C. Koo, H.R. Choi, Tactile display with rigid coupling based on soft actuator, Meccanica 50 (2015) 2825–2837. [16] F. Berlinger, M. Duduta, H. Gloria, D. Clarke, R. Nagpal, R. Wood, A modular dielectric elastomer actuator to drive miniature autonomous underwater vehicles, 2018. [17] S. Guo, X. Jia, K. Hashimoto, T. Li, Elastic characteristics of a strip-like electro-active elastomer actuator, J. Intell. Mater. Syst. Struct. 29 (2018) 2570–2580. [18] A.K. Sharma, M.M. Joglekar, Effect of anisotropy on the dynamic electromechanical instability of a dielectric elastomer actuator, Smart Mater. Struct. 28 (2019). [19] B. Li, H. Chen, J. Qiang, S. Hu, Z. Zhu, Y. Wang, Effect of mechanical pre-stretch on the stabilization of dielectric elastomer actuation, J. Phys. D: Appl. Phys. 44 (2011). [20] Y. Su, W. Chen, M. Destrade, Tuning the pull-in instability of soft dielectric elastomers through loading protocols, Int. J. Non-Linear Mech. (2019). [21] J.S. Plante, S. Dubowsky, Large-scale failure modes of dielectric elastomer actuators, Int. J. Solids Struct. 43 (2006) 7727–7751. [22] G. Kofod, The static actuation of dielectric elastomer actuators: how does pre-stretch improve actuation?, J. Phys. D: Appl. Phys. 41 (2008) 215405. [23] Y. Zhu, H. Zhou, H. Wang, Electromechanical characteristic analysis of a dielectric electroactive polymer (DEAP) actuator, Smart Mater. Struct. 24 (2015). [24] L. He, J. Lou, J. Du, J. Wang, Finite bending of a dielectric elastomer actuator and pre-stretch effects, Int. J. Mech. Sci. 122 (2017) 120–128. [25] L.Z. Lyu, S.J. Zhu, The electromechanical behavior of dielectric elastomer actuator stiffened by fiber, Key Eng. Mater. 765 (2018) 12–15. [26] K.B. Subramani, E. Cakmak, R.J. Spontak, T.K. Ghosh, Enhanced electroactive response of unidirectional elastomeric composites with high-dielectric-constant fibers, Adv. Mater. 26 (2014) 2949–2953. [27] T. Lu, J. Huang, C. Jordi, G. Kovacs, R. Huang, D.R. Clarke, Z. Suo, Dielectric elastomer actuators under equal-biaxial forces, uniaxial forces, and uniaxial constraint of stiff fibers, Soft Matter 8 (2012). [28] J. Huang, T. Lu, J. Zhu, D.R. Clarke, Z. Suo, Large, uni-directional actuation in dielectric elastomers achieved by fiber stiffening, Appl. Phys. Lett. 100 (2012) 211901. [29] K.B. Subramani, R.J. Spontak, T.K. Ghosh, Influence of fiber characteristics on directed electroactuation of anisotropic dielectric electroactive polymers with tunability, Compos. Sci. Technol. 154 (2018) 187–193. [30] L. He, J. Lou, J. Du, H. Wu, Voltage-driven nonuniform axisymmetric torsion of a tubular dielectric elastomer actuator reinforced with one family of inextensible fibers, Eur. J. Mech. A Solids 71 (2018) 386–393. [31] D. McCoul, Q. Pei, Tubular dielectric elastomer actuator for active fluidic control, Smart Mater. Struct. 24 (2015). [32] S. Son, N.C. Goulbourne, Finite deformations of tubular dielectric elastomer sensors, J. Intell. Mater. Syst. Struct. 20 (2009) 2187–2199. [33] L. He, J. Lou, J. Du, Voltage-driven torsion of electroactive thick tubes reinforced with helical fibers, Acta Mech. 229 (2018) 2117–2131. [34] E.M. Mockensturm, N. Goulbourne, Dynamic response of dielectric elastomers, Int. J. Non-Linear Mech. 41 (2006) 388–395. [35] D. Vu, P. Steinmann, Nonlinear electro-and magneto-elastostatics: material and spatial settings, Int. J. Solids Struct. 44 (2007) 7891–7905.

13