Finite element modeling of electro‐viscoelasticity in fiber reinforced electro‐active polymers

In this work, aspects considering material modeling of electro‐mechanical coupling in fiber reinforced electro‐active polymers (EAP) and the corresponding finite element implementation are considered. We propose a constitutive model that takes into account the electro‐viscoelastic behavior of the isotropic matrix and the influence of unidirectional fibers on both the hyperelastic response and the viscous behavior of the whole composite. Two distinct existing models that describe the electro‐mechanical coupling, are demonstrated and implemented, moreover, a numerical link between both models for three‐dimensional continua in terms of tensor calculus, is identified. We propose the extended‐tube model for the elastic response with some of its parameters evolving in response to the electric field, in order to fit electro‐viscoelastic experiments. Regarding the finite element implementation, in addition to the deformation field and the electric potential, two pairs of field variables are introduced on the element level, to enforce quasi‐incompressibility and quasi‐inextensibility. It is shown that using the mixed finite element improves the convergence behavior for the simulation of soft EAP with relatively stiff fibers. Moreover, the choice of the model that expresses the nature of the underlying coupling is shown to noticeably affect the degree of simulated actuation in fiber reinforced actuators.


INTRODUCTION
In the context of soft robotics and biomimetic artificial muscles, electro-active polymer (EAP) are considered as promising candidates among other smart materials, due to their quasi-instantaneous response and their high deformation capability in response to the electric field. Where in some classes of EAP, activation time of less than 1 ms is demonstrated 1,2 and strains up to more than 300% are recorded. 3 Due to those properties, several prototypes of biomimetic actuators and grippers have been based on EAP, [4][5][6][7][8] with fibers embedded either in the EAP itself, or in the passive portion of the structure. For instance, anisotropy is introduced in EAP by aligning fibers in a certain orientation, in order to break the symmetry and to obtain bending deformation of the composite along predefined directions. 7,8 According to Reference 9, the term electrostrictive can be used in general, to indicate materials which demonstrate quadratic proportionality between their mechanical response (stress) and the induced electric field, regardless of the underlying mechanism. Specifically speaking, dielectric elastomer actuators (DEA) are electrostrictive soft polymers, 9 that are coated with conductor electrodes on the surfaces spanning along their relatively small thickness. Electrostriction in DEA is originated in terms of external electrostatic forces that develop through the application of voltage difference between its electrodes. Besides the mechanical behavior in response to the electric field, materials can be characterized based on their polarization behavior. Linear dielectric is a material which is polarized linearly in the electric field. 10 Materials with a stress response that is linear in the square of the electric field (quadratic), and at the same time exhibit polarization which varies linearly in the electric field, are referred to as dielectric materials in the literature. [11][12][13] From a mathematical point of view, the latter mentioned properties imply that the polarization and the stress are related to the electric field and the square of electric field, through a constant electric permittivity, respectively. However, electro-mechanical experiments of some dielectrics have demonstrated that the electric permittivity can vary depending on the prestretch (deformation) of the material. 14 The latter observation implicitly states that in dielectrics, the relation between stress and electric field is not necessarily quadratic, but is at least quadratic or is nonlinear in general, because the electric permittivity can vary during the electro-mechanical coupling process, while the material exhibits continuous shape distortion. Moreover, the polarization can evolve nonlinearly in the electric field, due to the deformation.
The concept of electroelastic deformations is numerically studied during the last century in References 15-19, among others. More recent works have been devoted to express linear and nonlinear electroelasticity. 11,12,[20][21][22][23][24] Moreover, solutions of simple boundary-value problems are provided in References 11,[20][21][22] In this contribution, we consider dielectrically linear materials with deformation-independent electric properties, based on the work of McMeeking and Landis. 12 Furthermore, we base our numerical work of nonlinear electro-mechanical solids with deformation-dependent electric properties, on the developments by Dorfmann, Bustamante, and Ogden. 11,[20][21][22] Dorfmann and Ogden 11 have provided a comprehensive review of electro-elastic deformations, both on the theoretical and the numerical sides. Moreover, using an example of planar actuator with equi-biaxial deformation, they demonstrated the link between the deformation-independent formulation and their proposed nonlinear constitutive model by writing the electric permittivity as a function of the stretch. 11 In this work, we demonstrate the link-point between the two aforementioned models 11,12 in the three-dimensional unrestricted environment using tensor calculus.
The numerical implementation of the electroelastic constitutive equations is demonstrated in Reference 25, where the variational formulation is set and the finite element linearization is performed in the material setting. Zäh and Miehe 26 briefly outlined the finite element formulation in both material and spatial settings. Moreover, they mentioned that implementing the coupled finite element in a spatial setting is preferable from a computational point of view. 26 Locking phenomena or artificial stiffening occurs in finite elements due to the nearly incompressible behavior of rubber-like materials, especially when the structural response is bending dominated. Motivated by the latter aspects, mixed finite elements have been used in References 27 and 28 to simulate electro-viscoelastic experiments and the response of electro-elastic bodies surrounded by free space, respectively.
The electro-viscoelastic coupled response of EAP has been investigated on both, experimental [29][30][31][32] and numerical 30,32-36 levels. A comprehensive characterization of the electro-viscoelastic behavior for VHB 4905 dielectric is conducted in, 30,31 where passive viscoelastic and electro-viscoelastic loading-unloading tests with various loading rates have been performed. Moreover, both passive and coupled single-step relaxation tests are demonstrated. 30,31 Regarding the loading-unloading tests, they have stretched the dielectric with different loading rates, while applying constant voltage difference across the specimen's thickness. It is demonstrated that the electric field influences the viscoelastic response significantly. 30 Furthermore, the authors have proposed and utilized a constitutive model to simulate the experiments and identify the corresponding material parameters. 30 Regarding the dissipative behavior, they set several nonequilibrium branches and adopt the energy function as it is proposed in Reference 37 with strain-type internal variables and a linear evolution law. For modeling the equilibrium response in Reference 30, they use the hyperelastic eight-chain model of Arruda and Boyce 38 with the shear modulus depending quadratically on the electric field. Such ansatz used to describe the dependency of mechanical parameters on electromagnetic fields is demonstrated in Reference 39. However, the consistent finite element linearization of the coupled electro-mechanical Arruda-Boyce model is not outlined in Reference 30. In this work, we extend the hyperelastic extended-tube model 40 to an electro-mechanical formulation, by coupling some of its parameters with the electric field, in order to fit the electro-mechanical experiments presented in References 30,31. Moreover, within a finite element framework, we demonstrate the consistent linearization of the coupled extended tube model. Regarding the time-dependent behavior of the isotropic EAP, we refer to nonlinear viscoelasticity by setting several nonequilibrium branches characterized by the simple Neo-Hookean energy function in terms of the elastic deformation. The evolution law for the viscous deformation is based on the work of Bergström and Boyce, 41 however, with the recent theoretical interpretations, modifications and the corresponding numerical treatment as proposed in  Biological tissues exhibit anisotropic mechanical response due to their naturally oriented internal structure. On the other hand, the response of an isotropic engineering material can be modified, through embedding fibers along predefined directions, thereby leading to a fabricated anisotropic material. As examples of works that have considered the modelling of viscoelastic anisotropic rubber-like materials, we refer to  Regarding the numerics related to the description of viscoelasticity in anisotropic living tissues, either in their passive or active state, we mention. 48,49 All the aforementioned approaches are based on defining structural tensors along predefined directions, then coupling those tensors to the right Cauchy-Green deformation tensor. Thereby, scalar invariants are obtained, 50 which are used as input in the energy function, with mechanical properties assigned as different than those for the isotropic matrix. In this work, we model mechanically transversely isotropic EAP considering the influence of unidirectional fibers on both the hyperelastic and the viscous behavior of the whole electro-mechanical material. Similar to Reference 47, we split the anisotropic energy contribution into elastic and viscous parts. However, based on the recommendations in References 51,52, we do not perform a volumetric-isochoric split of the deformation for the anisotropic contribution. Numerical frameworks for fiber reinforced composites have been proposed in References 45-47, taking into account the directional dependency of both the elastic and the viscous response of composite materials. In order to amend the influence of the time-dependent behavior for fibers, we define a viscous energy function in the logarithmic strain space and use an evolution law for the viscous deformation of the fibers, as it is suggested in Reference 48.
The use of a standard finite element for simulating the behavior of fiber reinforced rubber-like materials might lead to volumetric locking due to the quasi-incompressibility of the material. Moreover, the nearly inextensible response along the stiffer preferred direction can lead to a locking phenomenon or poor convergence behavior. In order to tackle those problems, formulations of mixed finite elements for anisotropic hyperelastic materials are proposed in References 53,54, among others. Dal 53 has proposed an extended Hu-Washizu mixed finite element formulation, referred to it as "Q1P0F0" element and performed a comprehensive study to examine the performance of the newly proposed element. In addition to the displacement field, he defined four additional field variables, which are treated on the finite element level as averagese. 53 Stability of mixed finite elements is assessed by the inf-sup condition, among others. The inf-sup condition can be alternatively referred to as Babuška-Brezzi (BB) condition 55,56 or as the LBB condition, including additionally the name of Ladyzhenskaya. 57 Regarding nonlinear problems, formulations for the LBB condition are generally nontrivial. 58 In a geometrically linear setting, the Hu-Washizu-type Q1P0 mixed finite element does not fulfill the LBB condition. Nevertheless, the aforementioned finite element has demonstrated relatively satisfying robustness for a wide range of practical applications, considering quasi-incompressible materials that exhibit finite deformations. 53,58,59 DEA materials are regarded as strongly soft in general, however, the fibers are relatively stiff. Therefore, in this work, we propose a coupled and mixed electro-mechanical Q1P0F0 finite element, based on the mechanical counterpart. Moreover, we study its behavior compared to an electro-mechanical Q1P0 element. Regarding the setting and the discretization for the coupled finite element, we base our work in a spatial setting and utilize Galerkin's approach, respectively. The paper at hand is organized as follows. Section 2 outlines the basic equations describing the kinematics of the coupled problem and the corresponding balance equations. Section 3 introduces the proposed material model and demonstrates the terms that are needed for the implementation within a fully monolithic finite element. Moreover, in Section 3, special emphasis is devoted to indicate the properties of models for both, dielectrically linear and deformation-dependent nonlinear electro-mechanics, with indicating the link between both. In Section 4, a mixed finite element of the coupled electro-mechanical problem is demonstrated. In Section 5, a set of numerical examples is demonstrated to show the capabilities of the proposed modeling approach and the corresponding numerical tool. Finally, we end the paper with some concluding remarks and appendices.

BASIC EQUATIONS OF ELECTRO-MECHANICS
This section is devoted to provide an overview of kinematics and balance laws that are needed for the description of an electro-mechanical coupled behavior of fiber reinforced bodies, that are capable of undergoing large deformation. Basic principles of continuum mechanics are briefly presented and utilized to portray the starting point toward a specification of the problem. Thereafter, kinematics and balance laws, describing the mechanical response of fiber reinforced materials, are outlined. Finally, the representation of the coupled problem is achieved by demonstrating the kinematics and F I G U R E 1 Mapping from the undeformed to the deformed configuration of a nonlinearly interacting electro-mechanical body balance laws of electric-field sensitive materials, where the link between the mechanical and the electrical response on the kinematical level is illustrated.

Preliminaries
In the context of large deformation, we refer to the undeformed material body  0 and its counterpart in the deformed configuration  t , as it is shown in Figure 1. We define a nonlinear deformation map x = (X, t) at time t ∈  , that projects points from their reference state X ∈  0 to their current configuration x ∈  t . Moreover, we introduce an electric potential (X, t). The deformation gradient is denoted as is a differential operator with respect to the material coordinates X. The boundary  is decomposed in terms of the mechanical boundary conditions as where Dirichlet-type mechanical boundary conditions can be prescribed on  and Neumann-type tractions can be imposed on  t . Furthermore, regarding the electrical counterpart, the surface  is split into which offers the possibility of applying an electric potential on the portion  or prescribing surface charges q on  q , which compose the Dirichlet and the Neumann electrical boundaries, respectively.

Mechanics of deformable continua
The three-dimensional mapping from the reference to the current state within a material body, that is capable of exhibiting large deformation, relies on the following three operators where F is the deformation gradient, cof[F] reads its cofactor, and J is its determinant. Recall that ∇ X [ ⋅ ] reads the gradient with respect to the material coordinates X. Moreover, (X, t) is a deformation map at time t ∈  where x = (X, t). The condition J = det[F] > 0 is mandatory in order to ensure that the deformation map does not permit self-penetration of the material. The defined three operators in Equation (3) denote the mapping of an infinitesimal line, an area element and a volume element from the undeformed to the deformed state, such that in which dX, dA, and dV are an infinitesimal line, area and volume elements at the reference state and dx, da, and dv represent their counterparts in the deformed setting, respectively. The description of fiber reinforced continuum requires the definition of a unit vector f 0 , which denotes the fiber direction in the material configuration. The vector f 0 can be mapped to the spatial configuration through the deformation gradient F, in which with the operator || ⋅ || for a vector magnitude and the fiber direction at the current state f. Furthermore, as strain measures, we refer to the right and the left Cauchy-Green tensors where g and G are the symmetric covariant, Eulerian and Lagrangian metric tensors, respectively. For Cartesian coordinates the tensors g and G are identical to the second-order identity tensor. Considering quasi-static problems, we refer to the balance of linear momentum in terms of three different measures of the stress as with the total first Piola-Kirchhoff stress tensor P, the total Kirchhoff stress tensor and the total Cauchy stress tensor . The operator ∇ X ⋅ [ ⋅ ] is the divergence with respect to the material coordinates X and its counterpart in the spatial setting is written as ∇ x ⋅ [ ⋅ ]. Moreover, the term 0 constitutes the volume specific mechanical body forces in the undeformed configuration and reads its counterpart in the current setting. Note that in what follows, the mechanical body forces are neglected. The stress tensors P, , and are total electro-mechanical stresses, that contain the contribution of electric body forces. 20,22 The relations between the different stress measures read as follows The first relation in Equation (8) applies due to the second identity in Equation (8), where the operator cof[F] links the infinitesimal area elements dA and da as it is shown in Equation (4.2). Moreover, the vectors N and n are unit surface normals at the material and the spatial configuration, respectively. The treatment of the mechanical Neumann boundary conditions is based on the expressions where t reads the spatial traction vector and T denotes a scaled traction vector of the spatial configuration, which are related through tda = TdA. The balance of angular momentum is satisfied in terms of the total stress tensors and through the relations

Electro-statics of deformable continua
This section introduces the electrical balance laws of electro-static problems. Moreover, it denotes the illustration of the kinematical tie between the mechanical and the electrical response within electro-mechanical active bodies. The assumption of electro-mechano-statics applies only in case of the absence of a varying magnetic field. 20 The electric potential can be parameterized in terms of the material and the current coordinates as (X, t) = (x, t), where X ∈  0 and x ∈  t , respectively. Due to the fact that the electric potential is a scalar quantity, we always refer to it as in this manuscript. In the context of electro-mechano-statics, Faraday's law can be expressed in terms of the reference E and the current e electric fields as respectively. The operator ∇ X × [ ⋅ ] denotes the curl of a vector at the undeformed configuration and ∇ x × [ ⋅ ] is its counterpart at the spatial setting. Equations (11.1) and (11.2) allow us to introduce the reference and the current electric fields as negative gradients of the scalar electric potential , where with ∇ X [ ⋅ ] and ∇ x [ ⋅ ] as the gradient operators in the material and the spatial settings, respectively. The relations in Equation (12) imply that The electric response of materials can be characterized by a relation that links the electric field with the electric induction, more commonly known as electric displacement. To this end, we introduce a Cauchy-type electric displacement d and define the generalized equation d v is the electric displacement in vacuum, d m is its counterpart in deformable matter, 0 is the electric permittivity of the vacuum, and p reads the Cauchy-type polarization in matter. The portion d v = 0 e implies that, in vacuum, the relation between the Cauchy-type electric displacement and the electric field takes a linear form in terms of the constant free space permittivity 0 . On the other hand, depending on the choice of the material model, the relation between the electric displacement in matter d m = p and the electric field e, can in general be nonlinear. Further specification of the latter aspect is presented later in this paper. As alternative measures of the Cauchy-type vector d, we define an electric displacement vector D in the reference configuration and a Kirchhoff-type electric displacement d = Jd. Moreover, we introduce a polarization vector P in the material setting and its Kirchhoff counterpart p = Jp. The infinitesimal charges identity D ⋅ dA = d ⋅ da along with Equation (4.2) enables us to express the reference electric displacement D and the reference polarization P as Making use of Equations (12) and (15), the relation between the reference electric displacement, electric field and polarization can be expressed as where D v = J 0 C −1 ⋅ E constitutes the reference electric displacement in vacuum and D m = P denotes its counterpart in deformable matter. The relation between D v and E does not take a linear form as it is shown in Equation (16). However, linearity applies only between the true measures d v and e, as it is demonstrated in Equation (14). Therefore, as a reference to characterize the nature of electric response in materials, we assess the relation between the spatial electric field e and the Cauchy-type electric displacement d. Note that in what follows, we assume that the contribution of electric response in vacuum can be neglected for our problems under consideration, 25,60 where Equation (14) can be reduced to d = d m . Furthermore, Gauss law for electricity can be expressed in terms of the reference D, the Kirchhoff-type d and the Cauchy-type d electric displacements as respectively. 0 f reads the free volume charges per unit reference volume and t f is its spatial counterpart. In this manuscript, we restrict ourselves to cases where free volume charges do not exist ( 0 f = t f = 0). 11 Regarding material surfaces which are in contact with conductors (i.e., electrodes), the electric Neumann-type boundary conditions on the surface  q can be prescribed as follows with Q and q as surface charge densities in the material and the spatial settings, respectively.

MATERIAL MODELING AND CONSTITUTIVE EQUATIONS
In this section, we present a material model which constitutively describes the coupled electro-viscoelastic response of fiber reinforced EAP. The constitutive model expresses the behavior of an electro-active composite material with hyperelastic properties of the isotropic matrix that can evolve in response to an applied electric field. Directional dependency of both the hyperelastic and the viscous response of the material that can arise due to the presence of reinforcing fibers, is taken into account. Regarding the electro-mechanical coupled response, we demonstrate two distinct constitutive models. The first description is based on the assumption that the material follows a quasi-linear dielectric response, with only density-dependent electric permittivity. The second model expresses nonlinear electro-mechanical coupling, where the electric and coupled electro-mechanical properties of the material are deformation-dependent. The material model is in general described and expressions of the constitutive equations, that are necessary for the finite element treatment, are outlined. Thereafter, a detailed specification of the constitutive model is provided. The consistent linearization of the constitutive equations is provided in Appendix A1.

Coupled energy-enthalpy density function
We impose the electric field E as the electrical independent variable of our problem, in addition to the mechanical deformation gradient F. Moreover, we introduce an electro-mechanical coupled enthalpy density function Ψ coup (F, g, E). Assuming that the mechanical properties of EAP can in general be electric field sensitive, we write the mechanical energy density function as Ψ mec ( F, g, f 0 , E ) and express a total energy-enthalpy density functionΨ For more details about the energy-enthalpy based approach, we refer to Miehe et al. 60 In order to treat the quasi-incompressibility of polymeric materials, we decompose the total deformation into a volumetric and a volume-preserving isochoric part with the volumetric deformation gradient F vol , its isochoric counterpart F and the second-order identity tensor 1.
The total deformation gradient is also split into elastic and viscous contributions, F = F e F v . Likewise, the isochoric deformation obeys an elastic-viscous decomposition, In order to describe the anisotropic behavior of fiber reinforced polymers, we introduce the additive decomposition of the total mechanical energy density into isotropic and anisotropic contributions, Ψ mec = Ψ isot + Ψ ani . The isotropic part of the mechanical energy is split into where Ψ vol e (J) is a hyperelastic-volumetric part and Ψ iso e (F, g, E) reads an isochoric contribution of the equilibrium branch.
is the energy density of a nonequilibrium branch i, with a total number of branches n. The anisotropic portion of the mechanical energy density is decomposed into elastic Ψ ani parts as follows The volumetric-isochoric split is employed to the isotropic energy contribution Ψ isot as it is shown in Equation (21), in order to enforce quasi-incompressibility and to treat the pathological volumetric locking of rubber-like materials. However, on the other hand, we do not refer to the volumetric-isochoric decomposition of the anisotropic energy contribution Ψ ani , but rather express it in terms of the total deformation F. This approach is based on the recommendations in References 51,52. In order to impose quasi-incompressibility and quasi-inextensibility, the volumetric part Ψ vol e and the anisotropic contribution Ψ ani e are treated on the finite element level. However, the rest of the energy contributions are processed at Gauss quadrature points.

Constitutive relations and material-tangent moduli
The total Kirchhoff stress is defined as Based on the additive split of the total energy-enthalpy density functionΨ tot ( as it is demonstrated in Equations (19), (21) and (22), the total stress can be additively decomposed as follows The stress contributions vol e and ani e enter the finite element formulation using constraints which are imposed on the finite element level. Therefore, we define the stresses which enter the formulation at the Gauss points in a standard manner aŝ= The Kirchhoff-type electric displacement vector d can be expressed as and it takes the additive form Commonly, the term d is restricted to d = d coup as it is given in most of the previous works. However in this contribution, the term d e iso arises due to the dependency of elastic material parameters on the electric field, as it is discussed in Section 3.1. Full treatment of a monolithic electro-mechanical finite element requires the computation of the following material-tangent moduli C is the fourth-order mechanical tensor, C and C are third-order coupled material tensors and C constitutes the second-order electric tensor. The corresponding Lie derivatives of the coupled problem are demonstrated in Reference 26. Based on the decompositions of the energy-enthalpy function as it is described in Section 3.1, we can split the material tangent-moduli (Equation (28)) into Furthermore, we define the fourth-order tangent moduli for the terms which are treated at the Gauss quadrature points aŝ Regarding the third-order tensors C and C , it is important to note that the computation of the third-order moduli C , as it is shown in Equation (28), is not necessary, due to the fact that the calculation of the third-order tensor C serves for the complete description of the finite element. Therefore, in what follows, we do not refer to the specification of C .

Electro-hyperelastic extended-tube model
In this section, we demonstrate the constitutive treatment of the volumetric Ψ vol e and the isochoric Ψ iso e hyperelastic contributions of the isotropic energy density function Ψ isot . In order to enforce quasi-incompressibility of the material, we choose the following expression for the volumetric part with the bulk modulus of the material . The volumetric Kirchhoff stresses vol e is introduced in terms of the volumetric pressure p as vol The hyperelastic-isochoric part Ψ iso e of the energy density is substituted by the extended-tube model, 40 with material parameters that depend on the electric field. To this end, Ψ iso e is decomposed into cross-linking W e and topological tube constraints L e contributions , with the parameter describing limited chain extensibility and the parameter considering topological constraints. The first invariant of the isochoric Cauchy-Green tensors readsĪ Alternatively, C and b can be defined in terms of the isochoric principal stretches { a , a = 1 to 3}, as where N a and n a are the eigenvectors at the reference and at the current configuration, respectively. Furthermore, G c (E) and G e (E) are cross-linking and topological constraints shear moduli, respectively. In order to describe the dependency of G c (E) and G e (E) on an applied electric field, we first introduce the pure electrical invariant I 6 as in terms of the reference E or the current e electric fields. Thereafter, we define the following phenomenological expressions, where G c0 and G e0 are the base shear moduli. Moreover,Ĝ c andĜ e are material parameters that describe the dependencies of G c and G c on the electric field, respectively. The choice of negative values forĜ c andĜ e indicates that the material exhibits softening in response to an applied electric stimulus. However, if the values ofĜ c andĜ e are positive, the material is assumed to harden due to the application of electric field. The isochoric Kirchhoff stresses for the equilibrium branch can be obtained through the following relations with the fourth-order deviatoric projection tensor P = I 1 − 1∕3 g ⊗ g −1 , and the fictitious Kirchhoff stress tensor e . Note that I 1 denotes a fourth-order identity tensor as it is defined in Appendix A1, in terms of the second-order identity 1.
Based on the split of Ψ iso e into W e and L e as it is shown in Equation (33), we introduce the decomposition e = e W + e L .
The definition of e , as it is shown in Equation (37), leads to the following relations for e W and e L , wherēI The contribution d iso e of the electric displacement is evaluated through the relation d iso e = − e Ψ iso e (C, J, E). Referring to the decomposition of Ψ iso e into W e and L e , as it is demonstrated in Equation (33), we introduce the following relations The volumetric and isochoric fourth-order moduli C vol e and C iso e , the second-order electric tangent-moduli C iso e and the third-order coupled moduli C iso e can be evaluated as it is demonstrated in Appendix A1.

Isotropic nonequilibrium response
For each nonequilibrium branch i, we use the Neo-Hookean energy expression in terms of the first invariant of the elastic-isochoric right Cauchy-Green tensorĪ C with a viscous shear modulus G v,i for each branch. There is a connection between the evolution of the elastic deformation F e i and the corresponding evolution of the inelastic deformation rate tensord v,i . Based on the approach of Bergström and wherėi > 0 reads the effective creep rate. The normalized projection tensor N P,i is defined as with the viscous Kirchhoff stress iso v,i corresponding to the branch i. We adopt the evolution law for the effective creep rate as it is proposed in 42 depending on the Kirchhoff viscous stresses in the sense of v,i and the elastic stretch e,i for the branch i. Moreover, the material parameters are restricted tȯ0/ m > 0, m > 0 and c ≥ 0. For more details about the model in general and its numerical treatment, the reader is referred to References 42,44.

Anisotropic equilibrium response
The complete description of the mechanical response of anisotropic materials is based on the specification of the additional invariants that express energy storage due to reinforcement of soft materials with a single fiber family. For the equilibrium part of the anisotropic energy Ψ ani , we rely on a quadratic function in terms of the fourth invariant I 4 , where with the elastic shear modulus of the fiber G f e . Moreover, Macaulay brackets, where ⟨x⟩ = (|x| + x)∕2, are used to eliminate the contribution of any compressive deformation, as it is assumed that the fibers bear tensile loading only (tension only condition). Subsequently, the anisotropic stress ani e of the equilibrium branch can be expressed as The explicit expression of the anisotropic moduli C ani e is described in Appendix A2.

Anisotropic nonequilibrium response
We utilize the geometrically linear theory to additively decompose the total logarithmic strain f in fiber direction into elastic Consequently, we define the viscous energy contribution of the fiber Ψ ani v in a quadratic form as . (51) f reads the conjugate thermodynamical force with respect to f v , m f controls the nonlinearity of the evolution law, f is the viscosity parameter and̂f constitutes the activation stress. In order to compute the internal variable f v , we use a backward Euler integration scheme to treat the evolution law in Equation (51). We use a discrete time increment Δt = t n+1 − t n , where t n and t n+1 are the previous and current time, respectively. Note that any variable without subscript is meant to be evaluated at time t n+1 . For more details on the computation of f v , the reader is referred to Appendix A2. After the evaluation of the internal variable f v , the stress tensor ani e can be evaluated as with The algorithmic tangent-moduli C ani v of the nonequilibrium part can be evaluated as it is shown in Appendix A2.

Quasi-linear dielectric response
In this section, we consider a model that describes electrically isotropic materials with density-dependent electric permittivity. 12 To this end, we introduce the coupled enthalpy density Ψ coup (C, E) as However, considering that for Cartesian coordinates g = 1, the Kirchhoff-type electric displacement d We can conclude from the first equation in Equation (56) that the relation between the Kirchhoff electric displacement d coup and the electric field at the current configuration e is linear, where the constant electric permittivity links them. However, the second equation in Equation (56) indicates that quasi-linearity applies between e and the Cauchy-type electric displacement d coup due to the density-dependency of the electric permittivitŷ(J). Consequently, the Kirchhoff coupled stress coup takes the form Note that I g −1 reads a fourth-order identity tensor as it is defined in Appendix A1, in terms of the inverse of the metric tensor g −1 . Assuming that the simplification g = 1 applies, the coupled Cauchy stresses can be expressed as Equation (58) with the relation coup = J coup shows that the stress coup changes quadratically with respect to the electric field e. For the sake of completeness, the associated second-order C coup , third-order C coup and fourth-order C coup tensors are computed as it is shown in Appendix A3.

Nonlinearly deformation-dependent electric response
In order to describe the electric response of materials that show nonlinear dependency of their electric properties in response to their exhibited deformation, we refer to 11 and define the following invariants Thereafter, we describe the coupled enthalpy density in terms of the sixth I 6 and seventh I 7 invariants as where c 1 and c 2 are material parameters, that express the polarization and electro-mechanical coupling of the material. Note that c 1 contributes to the control of the polarization only, however, c 2 influences both the polarization and the electro-mechanical coupling. Based on the enthalpy function defined in Equation (60), the relation between the Kirchhoff-type electric displacement d coup and the current electric field e takes a nonlinear form with respect to deformation, where Note that for Cartesian coordinates, where the Eulerian and Lagrangian metric tensors are identical to the second-order identity (g = G = 1), the Cauchy-type electric displacement d coup can be written as whereĉ 1 (J) = J −1 c 1 andĉ 2 (J) = J −1 c 2 are density-dependent material parameters. In order to further understand the meaning of the underlying nonlinearity, we compare the polarization (d coup − e relation) of the nonlinear formulation (Equation (62)), with that of the quasi-linear counterpart (the second relation in Equation (56)). This leads to the conclusion that in Equation (62), the nonlinearly deformation-dependent electric property that describes polarization, takes the form of a second-order tensor̂1 (J, b) = −2ĉ 1 (J) b − 2ĉ 2 (J) b 2 . In other words, the second-order deformation-dependent tensor̂1 (J, b) maps the vector e to the vector d coup . However, for the quasi-linear formulation, the polarization is controlled by the scalar and only density-dependent electric permittivitŷ(J), which appears in Equation (56). In order to assure that the directions of the vector d coup and the vector e are initially consistent, some restrictions on the parameters c 1 and c 2 should be taken into account. The parameter c 1 should be substituted by a negative value with an absolute value larger than the positive value of the parameter c 2 , where (c 1 < 0, c 2 > 0) and − c 1 > c 2 ). The Kirchhoff coupled stress coup reads Consequently, taking into consideration the condition (G = 1), the coupled Cauchy stress can be introduced as Recall that I b denotes a fourth-order tensor in terms of the left Cauchy-Green tensor b, as described in Appendix A1.
Regarding the electro-mechanical coupling ( coup − e relation), we can notice that for the quasi-linear formulation as it is shown in Equation (58), the coupled stress is obtained by scaling the second-order structural tensor [e ⊗ e] in terms of the electric field, by the scalar electric permittivitŷ(J). On the other hand, for the nonlinearly deformation-dependent formulation, Equation (64) demonstrates that the structural tensor [e ⊗ e] is mapped to the stress coup using a fourth-order deformation-dependent tensor̂2 (J, b) = 2ĉ 2 (J) I b . For the implementation of the model in a finite element framework, we introduce the tangent-moduli C coup , C coup and C coup in Appendix A4.

QUASI-INCOMPRESSIBLE AND INEXTENSIBLE ELECTRO-MECHANICAL FINITE ELEMENT
This section is devoted to present a mixed (Q1P0F0) electro-mechanical finite element. Besides considering a deformation field and an electric potential, the proposed coupled scheme relies on the treatment of four extra field variables as averages on the finite element level. First, we outline the main electro-mechanical balance equations in their strong forms and define additional constraints that are needed for the mixed finite element. Then, we obtain the weak forms of the main governing equations using Galerkin's method. Consequently, we demonstrate the linearization of the weak forms in the Eulerian setting, which is needed for a fully monolithic solution scheme. Finally, we present the corresponding finite element discretization.

Strong forms of the balance equations
In addition to the previously defined deformation field (X, t) and the electric potential (X, t), we define four additional field variables that are treated on the element level, where : the dilatation, which replaces the Jacobian J, to describe volume changes, p: pressure-type Lagrange multiplier, used as the conjugate of , f : kinematic variable describing the fiber stretch, which is used instead of I 4 , s f : penalty term, set as the conjugate of f .
In addition to the balance equations that are defined in Equations (7) and (17), we define four additional relations. The six relations in total are introduced in a strong form as Ψ vol e ( ) − p = 0, The total Kirchhoff stress is defined as in terms of the stresŝas it is defined in Equation (25), the variable p and the variable s f .

Weak forms of the balance equations
Following the method of weighted residuals, we scale the main strong forms of the coupled problem with integrable test functions. We multiply Equation (65.1) by the virtual test function and multiply Equation (65.2) by the test function , where both functions vanish at the corresponding boundaries, = 0 on  and = 0 on  . Consequently, we perform an integration by parts and apply the Gauss theorem, which yields the weak forms In what follows, we consider the processing for the integrals of Equation (67) that are originated in the solid domain . The Neumann-type terms, the traction t and the surface charge density q can be imposed using a surface domain. However, for the sake of brevity, in this work we do not demonstrate the treatment of surface finite elements. Furthermore, the additional constraints in Equations (65.3-65.6), can be expressed in their weak forms as where the first two constraints impose quasi-incompressibility and the last two relations are set to enforce quasi-inextensibility. The equations in Equation (68) are imposed in a weak manner on the subdomains  e , where  ≈ ⋃ n e e=1  e with n e as the total number of subdomains. Based on Equations (68.1-68.2), the kinematic variable and its conjugate penalty term p can be introduced within the element  e as The mean dilatation and the mean pressure p are constant within the element  e . In a similar manner, referring to Equations (68.3) and (68.4), the kinematic variable f and the penalty term s f are expressed as . (70) f and s f can be regarded as the average fiber stretch and the average fiber stress, respectively. Moreover, they are constant over the subdomain  e .

Consistent linearization
First, with regard to the main state variables and of the coupled electro-mechanical problem, we introduce their increments as Δ and Δ , respectively. Thereafter, we refer to the directional derivatives of the weak forms in Equation (67 The operator L and the operator D, denote the linearization and the directional derivatives of the residual terms, respectively. We exploit the relations ∇ x = ∇ X F −1 and ΔF −1 = −F −1 ∇ x Δ in order to express the incremental form of the nonlinear term ∇ x that appears in Equation (67.1) In what follows, we exploit the symmetry of the stresses , the symmetries of the fourth order tangent-moduli C ijkl = C klij and C ijkl = C jilk , and the symmetry of the third-order moduli C ijk = C jik . In analogy to the definition of the total Kirchhoff stresses, as it is shown in Equation (66), the linearized expression of the Kirchhoff stresses can be written as with the material tangent-moduli C and C , as it is given in Equation (28). Moreover, the increments Δp and Δs f are specified in Appendix B1. Regarding the nonlinear gradient term ∇ x that appears in Equation (67.2), we make use of the identities ∇ x = ∇ X F −1 and ΔF −1 = −F −1 ∇ x Δ , to evaluate its incremental form as Furthermore, the linearization of the Kirchhoff-type electric displacement can be expressed as where both C and C are defined in Equation (28). Based on the element averages p and s f as introduced in Equations (69.2) and (70.2), respectively, we can redefine the the Kirchhoff stresses in Equation (66) as witĥas it is defined in Equation (25)

Finite element discretization
Finally, we discretize Equation (67) and the associated incremental forms, using the finite element method. To this end, the domain under consideration  e is decomposed into solid elements, such that  ≈ ⋃ n e e=1  e , where n e is the total number of finite elements. Exploiting the isoparametric concept of the finite element method, we interpolate the deformation field e and the electric potential e and their test functions over the element  e , in terms of nodal values and continuous shape functions as with n en as the number of nodes within the element  e . Similar to Equation (78), the gradients for the test functions and the gradients for the increments of the variables e and e can be discretized as Furthermore, in order to simplify the discretization of ΔG , we introduce the averaged terms, The interpolated representations in Equations (78) and (79) can be incorporated in Equation (67), in order to express the global residual vector in a discrete form. For more details on the discretized residual vector, the reader is referred to Appendix B2. Furthermore, a discretized representation of the global tangent matrix can be obtained as it is described in Appendix B2. In this work, the finite element is specified as an eight-node brick element.

NUMERICAL EXAMPLES
In this section, we demonstrate the capabilities of our electro-mechanical constitutive model together with the proposed finite element through a set of numerical examples. As a norm for all the presented electro-mechanical simulations, we assume that the electric field develops through the thickness of dielectric elastomers upon the application of voltage difference between electrodes which are attached to the surfaces of the elastomer. Nevertheless, for the sake of simplicity, we do not model the electrodes, assuming that their thickness is relatively negligible. In the first numerical example, we simulate a planar actuator, using the two models that describe the electric response and are demonstrated in Sections 3.7 and 3.8. Next, we utilize an electro-mechanical patch test of a fiber reinforced sample to study the convergence behavior of the proposed electro-mechanical Q1P0F0 and compare it to its Q1P0 counterpart. As a third example, we simulate a set of electro-viscoelastic experiments in order to identify the material parameters of the proposed electro-viscoelastic extended-tube model. The fifth example is devoted to demonstrate the influence of fiber reinforcement on the actuation behavior of viscoelastic dielectric elastomer. Finally, we demonstrate the robustness of our numerical framework through simulations for viscoelastic and fiber reinforced EAP-based actuators, which exhibit relatively complex deformation states.

Linear versus nonlinear electric response
The difference between the earlier presented formulation for quasi-linear electric response, see Section 3.7, and its nonlinearly deformation-dependent counterpart, as it is shown in Section 3.8, is demonstrated through using both models to simulate the coupled response of a thin dielectric film. We consider a thin plate with dimensions of 25 × 25 × 0.1 mm, which is discretized by 30 × 30 × 2 elements. A voltage difference of 16, 000 Volts is applied stepwise within 100 time steps using a quadratic sinusoidal function. We emulate the response of three different cases: 1. Quasi-linear formulation, see Section 3.7, with electric permittivity = 4.7 0 . Homogeneous Dirichlet boundary conditions are assumed to apply on the deforming plate. Note that the electro-mechanical parameters , c 1 , c 2 are chosen such that the base electric properties for both formulations are identical. It is assumed that the plate is an elastic Neo-Hookean solid, with parameters = 24 MPa, G c0 = 0.0473 and G e0 = = 0, for the extended-tube model, reducing it to the Neo-Hookean model. Moreover, it is assumed that mechanical material parameters do not vary in response to the electric field, wherê G c =Ĝ e = 0.
The plate with deformation-dependent electric properties shows a slope variation of its polarization curve (d z − e z relation), as its shown in Figure 2(A). The deformation-dependent electric properties allow for a planar elongation up to 69%, for an applied voltage difference of 16, 000 Volts. However, the deforming plate with quasi-constant electric permittivity (quasi-linear formulation), demonstrates a linear polarization curve, as it is illustrated by Figure 2 Figure 2(A) shows that, when the deformation of the plate is totally constrained, the polarization curve of the plate with deformation-dependent electric properties renders the same behavior as the one with quasi-constant electric permittivity (quasi-linear formulation), however, without showing material instability.

Electro-mechanical dual clamped patch
The performance of the mechanical quasi-inextensible and quasi-incompressible finite elements in References 53,54 is tested using a two-dimensional dual clamped test. In this contribution, we use a three-dimensional dual clamped patch to demonstrate the convergence behavior of the proposed electro-mechanical quasi-inextensible and quasi-incompressible element. Similar to the examples in References 53,54, we consider a reinforced sample with fibers aligned to an angle of 60 • by setting f 0 = [0.5, √ 3∕2, 0]. The sample has dimensions of 10 × 10 × 0.5 mm and it is fixed at its both ends, as it is shown in Figure 3(A). Both, total mechanical traction of t 0 = 0.05 N/mm 2 and total electric potential difference Δ = 10,000 Volts, as it is illustrated in Figure 3(A), are monolithically applied. We evaluate the displacement at the middle of the specimen's end, as it is shown in Figure 3(A). The performance for both, the electro-mechanical Q1P0F0 and Q1P0 elements, is assessed through varying the number of elements m = {4, 8, 16, 32, 64 or 128} along the length and the width of the sample while keeping the number of elements along the thickness constant (n = 2), see Figure 3(B). For a given number of elements m, if the numerical solution diverges within 500 steps, the test is considered as divergent. For the patch test, we consider the material as Neo-Hookean with electric field independent parameters = 23.64 MPa and G c0 = 0.0473 MPa, where the initial Poisson's ratio equals to 0.499. The fiber stiffness is varied, such that G f e = {5, 10, 10 2 or 10 4 } MPa. Moreover, the coupled parameters read c 1 = −20 0 and c 2 = 10 0 . Regarding the fibers, we consider a set of simulations with applying tension only condition, and another set assuming that fibers can bear both tension and compression. Figure 4 demonstrates that for the test with tension only condition  Figure 5, shows that the behavior of the Q1P0 finite element is more robust than in the case when the tension only condition is imposed, as divergence can be seen only for G f e = 10 4 at m = 64. We conclude from the results in Figures 4 and 5 that the electro-mechanical Q1P0F0 element is more robust than its Q1P0 counterpart, when the stiffness of the fiber is relatively high. However, for the cases with the tension only condition imposed, the ill-conditioning of the Q1P0 element is more pronounced than in cases where the tension only condition is not assumed.  Figure 6. Note that as it is shown in Table 1 In order to identify the viscous and the electro-mechanical parameters, we use both passive and electro-mechanical loading-unloading experiments, which are presented in Reference 30. The the specimen has dimensions of 100 × 70 × 0.5 mm. Moreover, the specimen is stretched at a strain ratė= 0.1 s −1 . The loading-unloading experiment is conducted for several times, with applying various constant potential difference Δ = {0, 2, 3, 4, 5, and 6} kV across the thickness of the loaded specimen. A finite element model of the full specimen, with both mechanical and electrical loadings, is illustrated in Figure 7. However, for the sake of efficiency, we make use of the symmetric boundary conditions and simulate one-eighth of the specimen size only.
For the fitting process, we use two nonequilibrium branches and identify the corresponding viscous parameters using the passive viscoelastic expirement, where Δ = 0 kV. Thereafter, we identify the coupled electro-mechanical parameters, using the loading-unloading expirement with Δ = 6 kV. Furthermore, we simulate loading-unloading experiments with Δ = {0, 2, 3, 4 and 5} kV. Consequently, we compare the experimental data and the simulation results for the maximum resulted force. The experimental data and the simulation results are illustrated in Figure 8.

Fiber reinforced planar actuator
A planar actuator with the dimensions 25 × 25 × 0.2 mm is discretized by 24 × 24 × 4 elements, as it is illustrated in Figure 11. The mechanical parameters for the isotropic matrix are taken as it is shown in Table 1. However, the coupled parameters read c 1 = −6 0 , c 2 = 3 0 andĜ c =Ĝ e = 0 N/(Vmm) 2 . Furthermore, three different sets of parameters for the fibers are given in Table 2. Note that the three sets of parameters that are assumed to describe the anisotropic response of the material as it is shown in Table 2 are chosen in a way that allows us to investigate the capability of the model for a range of parameters, without necessarily referring to specific materials that are commonly used as fiber reinforcement in all the assumed cases. First, we investigate the influence of anisotropic viscosity on the actuator performance, using the set 2 of parameters, as it is shown in Table 2. We apply a linearly evolving cyclic electric stimulus as it is domenstrated in Figure 9(A) with a time step Δt = 0.01 s. Furthermore, we simulate the electro-mechanical response of the actuator for two different cases by assuming that fibers do not affect the viscous behavior of the structure (elastic fibers), and, another case, where anisotropic viscosity is incorporated (viscoelastic fibers). Figures 9(B) and 9(C) demonstrate that in general, instead of deforming equally along the two planar directions, the fiber reinforced planar actuator exhibits higher stretch in the perpendicular direction (y) than in the direction where fibers are aligned (x). Figure 9(B) depicts that the incorporation of anisotropic viscous effects reduces the stretch in the x-direction, where the fibers are oriented. Furthermore, the actuator exhibits higher elongation in y-direction as it is shown in Figure 9(C). The influence of fiber reinforcement on the behavior of the planar actuator is studied in terms of three sets of parameters for the fibers as it shown in Table 2. The sets of parameters are labeled by 1 to 3, where the set with number one describes a relatively soft class of fibers, and the set labeled by the number three belongs to the most stiff class of fibers. Figures 10(A) and 10(B) depict that as the stiffness of the fibers increase, the actuator shows higher stretch along the perpendicular y-direction, at the expense of the elongation along the x-direction, where the fibers are oriented. Figure 10(A) shows that when the stiffness of the fibers is extremely high compared to the stiffness of the isotropic matrix (parameters set 3), the planar actuator does not deform in the x-direction, however, it shows the highest elongation in y-direction, compared to the cases with less stiff fibers (parameters sets 1 and 2). For more illustration, see Figure 11.

Fiber reinforced actuators
As a prototype for this example, we use a bilayered composite with a fiber reinforced layer that is sandwiched between two electrodes and a passive isotropic layer as it is demonstrated in Figure 12. The plate has dimensions of 50 × 25 × 1.0 mm and is discretized with 40 × 20 × 4 elements. The parameters are taken as given in Table 3, unless otherwise mentioned. We simulate the coupled response of the structure with the setup as it is shown in Figure 12 for various fiber angles = {0 • , 20 • , 45 • and 90 • }. Note that the angle is measured within the plane x − y, as illustrated in Figure 12(A). A cyclic electric potential with a quadratic sinusoidal function, as illustrated in Figure 17(A), is applied to the bottom face of the structure, while the potential at the middle of the composite's thickness is held as zero. We use a time step Δt = 0.0025 s. When the fiber reinforced layer is electrically activated, a force is exerted. Depending on the orientation of the fiber layout, the force causes the composite to deform in a certain direction. Figures 13,14,15, and 16 demonstrate the We have sought to investigate the influence of the different models which are used to describe the electro-mechanical coupled response on the actuation mechanism of bending-twisting structures. To this end, we perform simulations of the composite with = 20 • , using both the quasi-linear formulation of the electric response (Section 3.7) and the for quasi-linear and nonlinearly deformation dependent formulations and (C) with and without considering electric field dependency of the hyperelastic parameters nonlinearly deformation dependent model (Section 3.8). Thereafter, the degree of actuation in terms of the underlying nature of coupling is assessed through studying the absolute magnitude of the displacement. Figure 17(B) demonstrates that the deformation intensity is noticeably dependent on the choice of the model, which is used to describe the pattern of coupling ( coup − e relation). In case of assuming deformation-dependent electric properties of the DEA, the structure is shown to deform with lower intensity than in case of using the quasi-linear model, where the electric permittivity is only density-dependent, however, with no highly nonlinear dependency on the deformation. Furthermore, we perform simulations with and without considering the dependency of the mechanical parameters on the electric field. This comparison can be performed through simulating the response of the structure with two different conditions, first with setting the parametersĜ c andĜ e as it is shown in Table 3, and second with substituting both parameters as zero. Figure 17(C) depicts that the deformation of the structure is slightly higher when the mechanical properties of the material are assumed to evolve in response to the electric field, compared to the case where the underlying dependency is neglected. The latter effect is due to the observed influence of the electric field on the strain-force relation, as it is shown in Figure 8, where the force corresponding to a certain strain intensity reduces as the applied electric potential is set higher.

CONCLUSION
In the presented work, we have proposed a constitutive model and a finite element framework suitable for the simulation of electro-viscoelastic fiber reinforced bodies. The electro-mechanical coupling in isotropic EAP matrices is taken into regard with depicting the influence of viscoelastic anisotropy on the coupled response of the composite material. An electro-mechanical extended tube model is proposed and the corresponding material linearization is demonstrated. The extended tube model with its material parameters set as electric field dependent is used to fit some electro-mechanical experiments of isotropic VHB material, as it is shown in Figures 6 and 8. Nevertheless, we look forward for further experimental investigations on the coupled behavior of materials that show typical S-shaped force-strain response. Based on those experiments, the capabilities of the coupled extended tube model can be more widely studied. Two distinct models, 11,12 that describe the coupled response of electro-active materials are demonstrated and processed for the use within finite element framework. Furthermore, a link between both descriptions is illustrated, for a three-dimensional deformable continuum. We have shown that for dielectrically quasi-linear material, 12 the polarization and the coupling stress constitute of the electric field e and the structural electric field tensor [e ⊗ e], scaled by the scalar density-dependent electric permittivitŷ(J) of the material, see Equations (56) and (58), respectively. On the other hand, the nonlinear deformation-dependent formulation 11 regards the polarization and the coupling stress as mappings of the electric field e and the tensor [e ⊗ e], using deformation-dependent second-order and fourth-order tensors, respectively. The latter relations are depicted in Equations (62) and (64). Moreover, the polarization second-order tensor̂1 (J, b) and the coupling fourth order-tensor̂2 (J, b), are defined in Section 3.8. The implications of both models are depicted through a numerical example in Section 5.1. In Section 4, we have proposed a mixed and coupled finite element to tackle limitations that arise due to nearly incompressibility and inextensibility of the material. Our element can be regarded as the extension of the mechanical Q1P0F0 that is proposed in Reference 53. Moreover, the capability of the finite element is briefly demonstrated in Section 5.2. Finally, in Section 5.4 and Section 5.5, we have studied the influence of mechanical anisotropy on the actuation behavior of soft actuators, in a quantitative and a qualitative sense.
The fourth-order identity tensor I A is defined in terms of an arbitrary second-order tensor A as In terms of the arbitrary second-order tensors A, B and the Cartesian basis e i , the upper and the lower dyadic products read The isochoric Eulerian moduli C iso e can be computed as ( P T ∶ e ⊗ g −1 + g −1 ⊗ e ∶ P ) , where C e reads the fictitious mechanical tangent-moduli of the equilibrium branch. Referring to the decomposition as in Equation (33), we define the additive split C e = C Moreover, the derivative 2Ī ] . (A9) In case of equal principal stretches a → b , the second term of Equation (A7) is computed using l'H•spital's rule The electric tangent-moduli C iso e = 2 ee Ψ iso e (C, J, E) can be evaluated as The third-order electro-mechanical moduli C iso e can be evaluated with the help of the fictitious third-order tensor C e as respectively. Note thatb = F ⋅ F T .

B.1: Finite element linearization
The terms ΔG and ΔG are expressed as respectively. The incremental forms of the averaged pressure Δp and the averaged fiber stretch Δs f are constant within the subdomain  e and they read Δp = Ψ vol e ( ) Δ , where Δ = 1 V e ∫  e Jg −1 ∶ g∇ x Δ dV, The material-tangent moduli C can be rewritten in the following form withĈ as it is expressed in Equation (30). Moreover, the fourth-order tensor V is defined as (B4)

B.2: Finite element discretization
A discretized form of the global residual vector R reads where R = ∫  e ∇ x N I ⋅ dV,