3D constitutive modeling of electro-magneto-visco-hyperelastic elastomers: a semi-analytical solution for cylinders under large torsion–extension deformation

The rise of a new class of smart materials known as electro-magnetorheological elastomers (EMREs) requires comprehensive understanding of their electro-magneto-visco-hyperelastic behaviors. The aim of this paper is to develop a generalized three-dimensional (3D) continuum-based framework of the electro-magneto-visco-hyperelastic behaviors of EMREs. The finite strain model is established based on the linear viscoelasticity theory and non-linear electro-magneto-elastic framework. As EMRE devices can be used in a cylindrical shape undergoing shear and normal stresses in many engineering applications like artificial muscles, a boundary-value problem simulating torsion–extension deformations of EMRE cylinders is developed in the finite strain regime and solved semi-analytically. The behaviors of EMRE cylinders under different loading conditions such as purely mechanical loading, purely electric loading as well as full coupling between mechanical, electric and magnetic loading are studied in detail. Influence of different parameters such as electric field, magnetic field, applied strain (-rate) and their coupling on the induced moment and axial force of the EMRE cylinder as well as its relaxation and creep under torsion–extension loading is also examined. It is shown that EMREs have adaptive capability and great potential in applications where the stiffness needs to be controllable. Due to simplicity and accuracy, the model is expected to be used in the future studies dealing with the analysis of EMREs in particular cylinders under torsion–extension developments like 4D printing of artificial EMRE-based cylindrical muscles.


Introduction
Nowadays, smart materials with exceptional properties have gained considerable attention. Shape memory polymers (SMPs) [1], shape memory alloys (SMAs) [2], hydrogels [3], magnetorheological or electrorheological fluid [4], magnetorheological elastomer (MRE) [5], dielectric elastomer [6] and also EMRE [7] are the most common smart materials used for advanced engineering applications. Light-weight, fast response, shape memory, shape changing, stiffening, multi-triggers and etc, are the common fantastic properties of these smart materials [8][9][10][11][12]. EMREs are a class of smart materials that can change their electric, magnetic and mechanical properties in the presence of full coupling between electric, magnetic and mechanical loadings. In other words, adding electric particles such as carbon nanotube and also ferromagnetic particles such as Fe 3 O 4 or carbonyl iron to elastomeric materials makes the matrix multi-trigger smart materials which are known as electro-magneto-elastic materials. Thanks to the properties of these novel materials, they have great potential applications in 4D printing [13], actuators and sensors, energy harvesters, adaptable optics (tunable lenses) and etc [14][15][16], in particular, artificial muscles [17,18]. The emerging properties of EMREs would receive increasing attention in the future in particular due to their potential for additive manufacturing [13,19].
On the numerical and modeling aspects, few studies have been dedicated to model EMREs in the presence of the full coupling between electric, magnetic and mechanical loadings [7,14,15,[20][21][22]. Most researchers investigated dielectrics or MREs, separately and mostly examined their hyperelastic behaviors. Since the base material of EMREs is a rubber-like material, they could endure large deformations. In this way, the general non-linear electro-magneto-elasticity theory provides some constitutive modeling of EMREs in a large deformation regime. In addition, these materials may also be considered as time-dependent materials. However, there are a few studies for modeling of visco-hyperelastic dielectric elastomers and MREs.
Recently, many attempts have been conducted in order to model hyperelastic behaviors of isotropic dielectrics and MREs [14,15,[23][24][25][26][27][28][29] and their anisotropy [30][31][32] for different applications and deformation regimes. In addition, to model time-dependent behaviors of dielectrics or MREs, there are few studies in the open literature that will be reviewed in the following. An electro-viscoelastic constitutive model for a polyurethane-based dielectric at finite strains was initially developed by Ask et al [33]. Using a coupled finite element formulation, they analyzed several boundary value problems numerically. In a complex model, Saxena and his colleagues [34,35] presented a general framework of finite deformation magneto-viscoelasticity. They investigated uni-axial, pure shear and simple shear deformations of MREs decomposing the free energy into elastic, viscoelastic and equilibrium and non-equilibrium magnetic parts. A new model was developed by Vogel et al [36] to investigate the effects of the electric field on the viscoelastic properties of the dielectric elastomers. Considering a geometrically non-linear finite element framework, they also investigated the relaxation and hysteresis behaviors of the dielectrics. To analyze thermo-electro-viscoelasticity of dielectrics, Mehnert et al [37] developed a mathematical framework decomposing the free energy function into four elastic, electric, time-depended and thermal parts. In a particular problem, for a spring-connected dielectric actuator, Wang [38] investigated creep and cyclic behaviors of the actuator via a non-linear three-element viscoelastic model by solving its differential-algebraic system of equations. Bishara and Jabareen [39] implemented a visco-electric-hyperelastic model using a user subroutine UEL in the commercial finite element method (FEM) software of ABAQUS. More recently, Garcia [40] introduced a new framework to model timedepended behaviors of MREs based on the multiplicative decomposition of deformation gradient tensor into elastic and viscous parts. They examined the effects of rate dependencies of mechanical loading and magnetic field on the uniaxial deformation.
The literature review indicates that the time-dependent behaviors of dielectric elastomers have been investigated by some researchers numerically (e.g. [37,41,42],) and experimentally (e.g. [43],). It also reveals the lack of a comprehensive finite-strain constitutive model with capability of simulating electro-magneto-visco-hyperelastic behaviors of EMREs in a three-dimensional (3D) framework. This paper aims at developing a general 3D continuum-based constitutive model of the electro-magneto-visco-hyperelastic behaviors of EMREs in the finite strain regime. For the time-dependent part of the model, the Maxwell-Wiechert model with ten nonequilibrium branches is considered where viscoelastic parameters are calibrated using a relaxation test of a VHB 4910. For the electro-magneto part, a nominal Helmholtz free energy density function adapted from Kumar and Sarangi's work [15,22] and Dorfmann and Ogdens' work [44] is assumed. Regarding the elastic part, Mooney-Rivlin and exp-exp strain energy density functions are considered. Then, inspired by artificial skeletal muscles and their 4D printing, an EMREbased cylinder as a boundary-value problem under torsionextension deformations is established in the finite strain regime and solved semi-analytically. The solution is verified for a pure torsion-extension loading by implementing the commercial FEM software ABAQUS. It is also shown that the governing equations can be non-dimensionalized to eliminate the magnetic parameters. Due to simplicity and accuracy, the model is expected to be an efficient tool in analyzing EMREs in particular cylinders under torsion-extension in future studies and developments like 4D printing of artificial EMRE-based cylindrical muscles.
The paper is organized as follows. In section 2, a general non-linear 3D continuum-based formulation for the electro-magneto-visco-hyperelastic behaviors of the EMRE is developed. In section 2.3, based on the existing experiments in the open literature, the model parameters are calibrated. Next, the non-linear continuum-based formulation of EMREbased cylinders under large torsion-extension is presented and examined semi-analytically in section 2.4. In section 3, the numerical results are presented. Firstly, the verification of the problem is carried out. Next, the effects of the electric field in the absence of the magnetic field on electric and mechanical behaviors of the cylinder are studied. The influence of some parameters like strain rate is investigated and the relaxation and creep behaviors of the cylinder under a fixed electric field are examined. Then, the coupling of the electric, magnetic and mechanical loading are considered. The effects of the magnetic field under a fixed electric field on the mechanical behaviors of the cylinder are investigated via the semi-analytical solution. Finally, in section 4, the summary and conclusion of this work is presented.

Electro-magneto-visco-hyperelastic constitutive modeling
Mechanical response of viscoelastic materials is a combination of elastic and viscous properties [45]. In general, based on the linear viscoelastic theory of time-dependent materials, the relationship between strains and stresses can be expressed as follows [46,47]: (1) in which, σ vh , ε, σ H and s are total viscoelastic stress, total strain tensor, the stress depending on the elastic part and a non-dimensional function which describes the time (or relaxation) phenomenon of the material, respectively. Also, based on the definition of σ H , equation (1) can be used in both small deformation (i.e. linear viscoelastic material) and large deformation (i.e. visco-hyperelastic material). Meanwhile, based on the assumption of the considering dielectric and magnetic effects of elastomers in their strain energy function instead of decomposing stress components, motivated by the linear viscoelastic theory, the relationship between the total Cauchy stress and strain in electro-magneto-viscohyperelastic materials can be expressed as follows. A schematic drawing of the proposed model is also shown in figure 1.
where σ(E, H, t, ε) is the total stress tensor and σ 0 is a stress depending on the strain, electric and magnetic parts. E, H, and t are the electric field, magnetic field and time, respectively. As mentioned before, s is a non-dimensional function which can be commonly defined in Prony series form as: in which s ∞ and s i are non-dimensional material constants corresponding to the equilibrium and instantaneous (viscous) parts, respectively. Furthermore, the values of s ∞ and s i vary between 0 and 1 and also the constraint s ∞ + n ∑ i=1 s i should be satisfied. τ i indicates the corresponding relaxation time or retardation time in ith branch (i = 1, 2, 3, …, n). In section 2.2, the time-dependent part is discussed in detail.

Electric, magnetic and strain-dependent parts
In the current configuration, electrical field variations are E, D, and P which are electric field vector, electric induction or electric displacement vector and polarization density vector, respectively. Similarly, for the magnetic field, the variations H, B, and M are called the magnetic field vector, magnetic induction vector and magnetization density vector, respectively. For condensed materials, these parameters are simply related together as: in which ε 0 and µ 0 denote the electric permittivity and magnetic permeability of free space. The simplified form of equation (4) for linear isotropic media is defined as [48]: where ε r and µ r are the dielectric constant so-called relative electric permittivity and dimagnetic constant so-called relative magnetic permeability. In addition, the macroscopic formulation of Maxwell's equations (i.e. Maxwell's equations in a vacuum) can be written as [49]: div (D) = ρ (Gaus's law of electricity) div (B) = 0 (Gaus's law of magnetism) where ρ, t, and J denote free electric charge density, time and free electric current density. Also, curl and div denote curl and divergence operators with respect to the current position vector x. For a quasi-static state, ∂/∂t vanishes and in a vacuum with free current or electric charge, the simplified forms of Maxwell's equations (6) in the current configuration can be expressed as: Correspondingly, the nominal electric field vector and the magnetic field vector in the reference configuration are denoted by E l and H l , and the relation between them based on the standard kinematic can be expressed as: Adapting equation (7), the quasi-static Maxwell's equations in the reference configuration are expressed as: Curl wherein Curl and Div denote curl and divergence operators in the reference configuration with respect to the reference position vector X.

Constitutive equation of electro-magneto-hyperelastic elastomers.
Following [50,51], the total Cauchy stress tensor σ 0 (or true stress) followed by the Maxwell's concept for an electro-magneto-elastic material with the mechanical Cauchy stress tensor σ Me can be expressed as [52]: (10) in which τ p , τ e , and τ m are known as polarization stress tensor, electrostatic Maxwell stress tensor and magnetic Maxwell stress tensor, and commonly expressed as [48,53]: However, this superposition of elastic, electric and magnetic stresses might not be accurate, especially for large deformations. To overcome this issue, an amended strain energy function in which the total stress can be determined through the nominal Helmholtz free energy density function Ω = Ω ( F, H l , E l ) is introduced. Recently, Kumar and Sarangi [15,22] presented a general form of the amended strain energy function in terms of b and b −1 that are adopted here [15,22]: where b is left Cauchy-green deformation tensor. Unlike the superposition of stresses, it is preferred to superposition nominal Helmholtz free energy density function. Therefore, the general relationship between total Cauchy stress tensor, electric induction vector and magnetic induction vector with strain energy function for an incompressible isotropic electromagneto-elastic material in the current configuration may be written as [48]: wherein p may be interpreted as a hydrostatic pressure (a Lagrange multiplier associated with the incompressibility constrain) and I is the second-order identity tensor. Also, the corresponding total first Piola-Kirchhoff stress tensor T (or nominal stress) based on the electric induction and magnetic induction in the undeformed configuration for incompressible and isotropic materials is expressed as: Since nominal Helmholtz free energy density function is commonly presented in terms of invariants (I 1 : I 9 ), the principal invariants depending on the tensor b and other invariants depending on E l and H l (or quasi invariants) can be defined as [15,22]: wherein, tr,:, and ⊗ represent trace, double-contraction, and outer product, respectively. In this study, for more simplicity, we consider that the electric field and magnetic field are constant and consequently we ignore the direct coupling between electric field and magnetic field (i.e. ignoring I 10 ). In order to calculate the stresses, electric and magnetic induction, based on the previous relations, equations (12)-(15), the explicit form of σ 0 , D and B can be derived by performing some mathematical manipulation as [15,22]: where, Ω i (i = 1 : 9) = ∂Ω/∂I i . More details on the derivation equation (16) can be found in the appendix. In this model, adopting the Mooney-Rivilin model [54] and exp-exp model [55] for the purely mechanical hyperelastic property of the EMRE, and based on the work done by Kumar et al [15,22] and Dorfmann and Ogden [44], a new and complete version of the nominal Helmholtz free energy density function for isotropic electro-magneto-hyperelastic materials is considered as: Now, by using equations (16) and (17) and having specific deformation gradient tensor, T, D, and B tensors can be easily obtained.

Time-dependent part
As discussed already, to model the time-dependent behaviors of EMREs, the linear viscoelastic theory (equation 2) is used where the non-dimensional function s is introduced by Prony series. In fact, by considering the Maxwell-Wiechert rheological model with an equilibrium branch and n number of non-equilibrium branch for the time-dependent part of the model (see, figure 1), the viscoelastic parameters s ∞ and s i can be calculated. Based on the constitutive equations of the Maxwell-Wiechert model, the relaxation modulus of the model is expressed as [56]: (18) in which E 0 is called instantaneous relaxation modulus which Also E ∞ is the equilibrium elastic modulus of the linear viscoelastic model (i.e. linear elastic part replacing in the orange box in figure 1). Thus, the parameter of s i can be expressed as E i /E 0 . Finally, in the following section, the material parameters are calibrated briefly.

Material and model calibration 2.3.1. Electric constants.
In order to determine the electric constants of VHB 4910 acrylic-based dielectric elastomer, Wissler and Mazza [57] presented the deformation-dependent relative electric permittivity ε r = ε r (λ) in equivalent-biaxial deformation loading as shown in figure 2(a). The deformation gradient F and electric field vector E for the experimental condition of the Wissler and Mazzas' work [57] can be expressed as: Based on equations (16), (17) and (19), the electric displacement vector is obtained as: Considering D 3 = ε 0 ε r E 0 , the deformation-dependent relative electric permittivity of the material is derived as: Finally, by fitting equation (21) with experiment data [57], the dielectric parameters are obtained as listed in table 1. Furthermore, the experiment data and the model prediction of the dielectric part are shown in figure 2(a).

Viscoelastic coefficients.
In this division, the timedependent properties of a VHB 4910 acrylic-based dielectric elastomer reported by Hossain et al [58] are considered. Here, in order to determine the viscoelastic parameters of the model, a single-step relaxation test of a standard sample at 20% strain performed by Hossain et al [58] is used. For the calibration purpose, the relaxation modulus is fitted with the experiment data. The relation between stress and relaxation modulus under small deformation for time-dependent materials is written as: Finally, considering appropriate numbers of branches of the Maxwell-Wiechert model (with ten branches), the timedependent parameters are obtained as listed in table 1. In addition, the experiment data and model prediction of the timedependent model are shown in figure 2(b).

Hyperelastic parameters.
In this division, in order to determine hyperelastic parameters of VHB 4910, Hossain et al's experiment [58] is used. They conducted some experiments on the strain-dependent and time-dependent properties of VHB 4910. First of all, the continuum formulation of the large deformation of elastic materials is presented. The deformation gradient for a pure homogeneous deformation can be written as: in which λ x , λ y and λ z are the longitudinal stretch in the x, y and z directions, respectively. By considering the incompressible postulate for VHB 4910 under uniaxial test, they are defined as: Finally, the true stress σ and nominal stress P Me can be calculated from the following non-linear continuum formulations [59]:    1956 where P z is the first Piola-Kirchhoff stress tensor in the z direction, F is the force measured from the experiment, A o is the cross-section of the test sample in the reference configuration and Ω Me is the mechanical strain energy density function related to the hyperelastic behavior of the material. Eventually, by adapting the Mooney-Rivlin model and expexp model for the mechanical strain energy density functions, the corresponding uniaxial true stresses and the uniaxial first Piola-Kirchhoff stresses can be derived presented as: In addition, the experiment data and prediction of the hyperelastic models are shown in figure 2(c). As the state of the art, figure 2 shows the comparison between experimental data and the predictions in the calibration process of the present model. It can be seen that the predictions are in a good agreement with experimental data showing the high accuracy of the calibration. All material parameters can also be found in table 1. It is noted that there is no related experiment on VHB 4910 to calibrate the material parameters of the magnetic part. To settle this issue, equations depending on the coupling between electric, magnetic and mechanical parts are reformulated in non-dimensional forms so that the equations will not depend on the magnetic parameters.

Torsion-extension deformation
As mentioned in section 1, dielectrics, MREs and EMREs are used in wire, rod and tube forms to produce actuating forces and moments. In many engineering applications, the smart devices experience axial-torsional deformations. EMREs also have great potential in 4D printing of smart actuators like artificial muscles. In many circumstances, these actuators like artificial muscles, may undergo different loading regimes like large bending, torsion-extension, etc From a simulation point of view, introducing a simplified model for these special structures is adorable. To this end, a boundary-value problem of the torsion-extension deformation of EMRE cylinders at a finite strain regime is developed here. The deformation mapping for torsion-extension loading in the cylindrical coordinates based on the reference configuration (R, Θ, Z) and current configuration (r, θ, z) for a solid cylinder with outer radius R o can be expressed as [60]: Where γ and τ denote the cylinder axial stretch, and the amount of twist per stretched length unit, respectively. A schematic of the solid cylinder under extension and torsion deformation in magnetic and electric fields is illustrated in figure 3.
The deformation gradient tensor of the torsion-extension loading of a cylinder in the cylindrical coordinate system can be written as: By adapting equation (28), the left Cauchy green tensor and its inverse can be derived as: and the principal invariants read as: The general form of the equilibrium equations in the radial and tangential directions (in the current configuration) is expressed as: where, F r , F θ , a r , a θ and ρ denote external radial body force, external tangential body force, radial acceleration, tangential acceleration and density, respectively. In the absence of external body forces and acceleration, based on the principle of conservation of linear momentum, equation (31) can be rewritten as: Considering the incompressibility condition and integrating equation (32) with respect to R, the following equation can be derived:ˆr Considering the boundary condition at σ rr (R = R rro = 0), one obtains Lagrange multiplier p can be obtained from equations (34) and (16), and by substituting it into equation (16), the stress tensor σ 0 is obtained. Finally, by substituting the stress tensor σ 0 into equation (2), the total Cauchy stress tensor can be readily calculated.

Results and discussion
In order to verify the solution, preliminary results are presented for a purely mechanical loading of the EMRE cylinder through both semi-analytical solution and an FEM in section 3.1. Then in section 3.2, the effects of the electric field in the absence of any magnetic field on different parameters such as axial force, moment, stress components, electric induction, relaxation and creep behaviors of the EMRE are investigated via the semi-analytical solution. Finally, in section 3.3, the effects of the coupling between electric and magnetic fields on the aforesaid parameters are examined.

Verification
It should be mentioned that since there are no appropriate experiments in the literature to calibrate the proposed model, the solution is verified in the absence of any electric and magnetic fields. To this end, Two different mechanical strains (γ = 1.2, θ = π 2 Rad and γ = 1.05, θ = π 4 Rad) are applied to the EMRE-based cylinder in the absence of any electric and magnetic fields and analyzed via commercial FEM software of ABAQUS CAE (6.17) and the semi-analytical solution implemented in Maple (2018) and Mathematica (11.3). It is noted that firstly the complete explicit forms of the stress components were obtained via Mathematica, then using a numerical integration method the total, Cauchy stress components were obtained via Maple. The variation of the moment and the axial force during a time period are shown in figures 4(a) and (b), respectively. As mentioned before, to model the visco-hyperelastic behavior of the EMRE, Maxwell-Wiechert with ten branches is used. Therefore, for modeling viscohyperelastic behavior of the EMRE in ABAQUS, the calibrated ten non-dimensional parameters s i and τ i reported in table 1 were implemented. Also, by choosing Mooney-Rivlin model for the equilibrium (hyperelastic) part, its material parameters were also implemented. In order to determine the optimum required number of meshes, a mesh study was performed and finally 100 000 elements were chosen. To mesh the cylinder, the cross-section of the cylinder was divided into two central and marginal partitions. Linear tetrahedral (i.e. C3D4H: A 4-node linear tetrahedron, hybrid, linear pressure) and quad elements (i.e. C3D8H: An 8-node linear brick, hybrid, constant pressure) were considered for the central and marginal partitions, respectively.
As can be seen in figure 4, the results from FEM and semianalytical solution are in good agreement in particular for the small value of the twist angle. Figure 4(a) shows that at strains of γ = 1.05, θ = π 4 Rad, the maximum difference between FEM and semi-analytical solution for the moment is 2.2% and at strains of γ = 1.2, θ = π 2 Rad, the maximum difference becomes 6.5%. Also, it is concluded from figure 4(b) that under large strains, in particular, a large twist angle, the axial force from the semi-analytical solution has a higher difference with its corresponding FEM. At strains of γ = 1.05, θ = π 4 Rad, the maximum difference between FEM and semianalytical solution for the axial force reads 6.5% and at strains of γ = 1.2, θ = π 2 Rad, the maximum difference reaches to 10.2%. The reason for these differences under a large twist angle seems to be the definition of the deformation gradient tensor used in the semi-analytical approach (see equation 28) that is not exactly similar to that in ABAQUS and may result in some differences at higher strains. It should be noted that at a strain of γ = 1.2, θ = π 2 Rad, the maximum shear strain (ε 23 ) is about 70% that means the applied strain is very large.

Effect of the electric field in the absence of magnetic field
In this section, in the absence of the magnetic field (merely dielectric), the coupling of the electric field and mechanical loading in torsion-extension of VHB 4910-based dielectric at a finite strain regime is examined. In this regard, the effects of electric field on the resultant moment and axial force, stress components, relaxation, creep behaviors and electric induction are investigated. The effects of the strain rate and applied strain at a fixed electric field are examined as well.

Effect of the electric field on the moment, axial force
and stress components. The effects of the uni-axial electric field E = (0, 0, E 0 ) considering on the axial force and external moment are examined as shown in figure 5.
Since the electric field is only applied in the z direction, it has no effects on σ θz . As the resultant moment is calculated based on the shear stress σ θz (35.1)), both σ θz and resultant moment are independent of E 0 , and consequently, they are not changed by varying E 0 , see figure 5(a). It can be conducted from figure 5(b) that increasing the electric field (E 0 ) makes the dielectric stiffer and finally produces higher σ zz , so that the axial force at an electric field of 50 MV m −1 becomes 2.47 higher than that in the absence of the electric field. Therefore, based on equation (35.2), the axial force increases. Furthermore, in order to illustrate the relation between the axial force, the resultant moment and σ zz with the electric field (E 0 ), figure 5(c) is studied. Based on the results presented in this figure, it can be found that increasing the electric field increases the stress level, moment and axial force since the EMRE becomes stiffer. It reveals that EMRE stiffness can be changed and the material can be classified as adaptive materials with controllable stiffness. This phenomenon is unique compared to other soft smart materials and shows the great potential of EMRE for mechanical applications, especially when high forces are required. For instance, it is expected that soft EMRE actuators could grasp objects with higher weight compared to those made of shape memory polymers.

Effect of the strain rate on the stress components.
Since the proposed EMRE model is time-dependent, the strain rate may affect the stress state. To this end, considering the fixed electric field E = (0, 0, 50 MV m −1 ), 20% fixed axial  strain and π/ 2 twist angle, the dependency of σ rr and σ zz to strain rates are examined as shown in figures 6(a) and (b).
It is well known that viscoelastic materials have damping properties. Therefore, the results presented in figure 6 for large strain rates show that the molecular chains do not allow to reorganize in the equilibrium state producing more stress and strength. It is seen that by increasing the axial stretch rate from 0.003 1/s to 0.12 1/s, σ rr and σ zz increase as much as 105.0% and 44.9%, respectively.

Effect of the electric field on the relaxation and creep behavior.
To investigate how ERMEs relieve the stress under a constant strain and how ERMEs strain under constant stress, relaxation and creep analyses are carried out, respectively. In order to examine the effects of the electric field on the relaxation behavior of EMREs, the following strains are applied to the material at different electric fields, and the results of the resultant moment and axial force under the proposed relaxation condition are shown in figures 7(a) and (b), respectively.
As discussed in figure 5, the electric field does not affect the resultant moment of the relaxation test, see figure 8(a). Also, based on equation (39), up to time t 1 , the resultant moment and axial force increase. However, from t 1 to (t 1 + t 2 ) under strains of γ 0 = 1.2, θ 0 = π 2 Rad, both the resultant moment and axial force decrease and start relaxing. As mentioned before, based on the experimental observation, the EMRE  becomes stiffer under high electric fields producing more axial force. Figure 8(a) reveals that the proposed model is also able to replicate this phenomenon.
Next, the effects of the electric field on the creep behaviors of the dielectric are investigated. Moment and axial force for this case are applied to the cylinder as follows: where M 0 = 3000 Nm, N 0 = 20 000 N, t 1 = 20 s, t 2 = 100 s are chosen. The axial stretch and the twist angle are reported as depicted in figures 8(a) and (b), respectively. It should be mentioned that the resultant moment or σ θz are not changed by the varying electric field in the relaxation (see figure 7(a)) and the creep tests. However, since the control parameters are moment and force in the creep test, both axial stretch and twist angle change when the electric field varies, see figure 9. The results show that, due to the hardening effect of the EMRE in the presence of the electric field, increasing the electric field decreases both axial stretch and twist angle, especially as time passes.

Effect of the uni-axial electric field on the electric induc-
tion By applying electric fields, electric inductions or electric displacements are generated where they have a relation with themselves via equations (16.2) and (36.3). Under different uni-axial electric field E = (0, 0, E 0 ), the electric induction in the z direction (D z ) is shown in figure 9.
Considering equation (36.3), it seems that the electric induction D z has a direct relation with the electric field E z . It means that increasing the electric field increases the electric induction. Also, due to the coupling between the electric field (E z ) and mechanical loading (20% axial stretch and π/2 twist angle), the electric induction (D z ) is generated more than purely electric case, see figure 9(b).

Effect of coupling electric and magnetic field
As mentioned before, there are no appropriate experimental data to calibrate the magnetic part of the model. Therefore, governing equations can be non-dimensionalized to eliminate the magnetic parameters in the constitutive relations. Considering uni-axial electric field and magnetic field E = (0, 0, E 0 ) and H = (0, 0, H 0 ) for the torsion-extension deformation of the cylinder (see equation 28), the stress components σ zθ and σ zz can be derived as: For the sake of simplicity and also due to the lack of appropriate experiments, the non-dimensional form of equation (41) can be rewritten by a normalization process as: (42) can be expressed as: The non-dimensional torsional moment and axial force can also be calculated via equation (35).
The effects of coupling between electric field, magnetic field and mechanical loading on the resultant moment, axial force and the relaxation behaviors of the EMRE will be investigated in the following sections.

The effect of coupling between the electric field and
magnetic field on the moment and axial force. In this section, the influence of coupling between electric, magnetic and mechanical loadings on the non-dimensional resultant moment and axial force under fixed strains γ = 1.2, θ = π 2 Rad is studied as shown in figure 10. Henceforth, (-) means that the parameter is expressed in the non-dimensional form.
As it can be seen from figure 10(a), the effect of the purely magnetic loading is more significant than the purely electric loading on the resultant moment. It is also observed that the coupling between electric, magnetic and mechanical loadings has the most significant effect on the resultant moment. It is found that by applying E 0 * = H 0 * = 5 with respect to purely mechanical loading (i.e. E 0 * = H 0 * = 0), the non-dimensional resultant moment is increased about 457.4%. However, for the non-non-dimensional axial force as shown in figure 10(b), the effect of the purely electric loading becomes much more significant about 4600 times than the purely magnetic loading.  It can be observed from figure 11 that the magnetic field has a prominent effect on both the resultant moment and the axial force. It is seen that by increasing the magnetic field H 0 * from 1 to 30, the non-dimensional resultant moment and axial force are increased 8995.5% and 183.1%, respectively. It can be concluded that not only the electric field makes the EMRE stiffer, but also the magnetic field helps make the EMRE much stiffer. Also, to further illustrate the stress state when non-dimensional electric and magnetic fields are varied, 2D contours of the non-dimensional stress components for strains of γ = 1.2, θ = π 2 Rad are depicted in figure 12. The results show that the stress components vary only in the radial direction. It is seen that the stress gradually increases along the radial direction. It is also observed that, since both electric and magnetic fields are applied in the z direction, σ zz is more sensitive compared to other stress components. Furthermore, it can be found that the radial and longitudinal stresses are more sensitive to magnetic and electric fields, respectively.  The results presented in figures 13(a) and (b) show that increasing the magnetic field increases both the resultant moment and the axial force and finally makes the EMRE stiffer. It is seen that by increasing the non-dimensional magnetic field from 1 to 30, the non-dimensional relaxed resultant moment and relaxed axial force increase around 164.6% and 8694.4%, respectively. Figure 13. The relaxation behavior in terms of (a) the moment, and (b) the force at 20% pre-stretch and π/2 twist angle in the different magnetic fields at the constant electric field E 0 * = 1.

Summary and conclusion
The main objective of this paper was to develop a 3D generalized constitutive model to replicate the electro-magnetovisco-hyperelastic behavior of EMREs. Considering a nominal Helmholtz free energy density function depending on the electric, magnetic and strain loading, the total Cauchy stress components were derived. Also, in order to consider timedependency of the EMERs, the total Cauchy stress under large deformations was derived by adopting the linear viscoelastic theory. The visco-hyperelastic and electric parts were calibrated using experimental data. Due to the lack of appropriate experiment data to calibrate the magnetic part of the model, it was successfully non-dimensionalized eliminating the magnetic parameters. As EMREs are commonly used in a cylindrical shape to produce actuating forces and moments, a boundary-value problem of EMRE cylinder under torsionextension at finite deformation range was developed and solved semi-analytically. The solution accuracy was verified by comparing the results with those from a developed FEM for the purely mechanical loading (i.e. visco-hyperelastoc behavior of EMREs). In a series of parametric studies, the cylinder was loaded by applying simultaneously uniaxial stretch and twist, and the resultant moment and axial force as well as the relaxation and creep behaviors were examined under different mechanical, electric and magnetic loads in detail. It was shown that EMREs have adaptive stiffness capability and great potential in mechanical/biomedical applications like adaptive actuators especially when their stiffness needs to be controllable.

Appendix
To establish equations (16.2) and (16.3) (i.e. electric and magnetic induction vectors), the derivatives of the principal invariants with respect to the electric field and magnetic field in the references configurations can be obtained based on the non-linear electro-magneto-mechanics as: Analogous to the above relations, the explicit form of the total Cauchy stress is obtained as: ORCID iDs M Bodaghi  https://orcid.org/0000-0002-0707-944X M Baghani  https://orcid.org/0000-0001-6695-3128