On the computational modelling of nonlinear electro-elasticity in heterogeneous bodies at finite deformations

Dielectric elastomer actuators (DEA) have been demonstrated to exhibit a quasi-immediate electro-mechanical actuation response with relatively large deformation capability. The properties of DEA make them suitable to be used in the form of major active components within soft robotics and biomimetic artificial muscles. However, some of the electro-active material properties impose limitations on its applications. Therefore, researchers attempt to modify the structure of the homogeneous DEA material by the incorporation of fillers that possess distinct electro-mechanical properties. This modification of the material’s structure leads to a fabricated inhomogeneous composite. From the point of mathematical material modelling and numerical simulation, we propose a material model and a computational framework using the finite element method, which is capable of emulating nonlinear electro-elastic interactions. We consider a coupled electro-mechanical description with the electric and the electro-mechanical properties of the material assumed to be nonlinearly dependent on the deformation. Furthermore, we demonstrate a coupled ansatz that expresses the electric response as dielectrically quasi-linear with only density-dependent electric permittivity. We couple the electro-mechanical models to the extended tube model, which is a suitable approach for the realistic emulation of the hyperelastic response of rubber-like materials. Thereafter, we demonstrate analytical and numerical solutions of a homogeneous electro-elastic body with the Neo-Hookean material model and the extended tube model to express the hyperelastic response. Finally, we use the finite element method to investigate several heterogeneous configurations consisting of soft DEA matrix filled with spherical stiff inclusions with changing volume fraction and ellipsoidal inclusions with varying aspect ratio.

The manuscript at hand is laid out as follows. Section 2 is devoted to present basic equations that express the balance states of the coupled problem and the corresponding kinematics. In Section 3, a constitutive material model and the linearized terms needed for the corresponding finite element implementation are outlined. In Section 4, an electro-mechanically coupled finite element is described. Section 5 is devoted to present a series of analytical and numerical examples. Finally, the paper is closed with some concluding remarks in Section 6.

General equations of electro-mechanics
In this section, an overview of the general equations that constitute the basis of coupled electro-mechanics is presented. First, principles of continuum mechanics are used to specifically describe inhomogeneous electro-sensitive materials that are subjected to large deformations. Subsequently, the kinematical link between the mechanical and electrical fields and the associated laws of balance are outlined.

Preliminaries
Considering large deformation, an inhomogeneous material compound B is split into two different electro-mechanical bodies with distinct material properties, where B = B I + B M , as it is shown in Fig. 1. We introduce a nonlinear deformation function x = ϕ(X, t) at time t ∈ T that projects points from the material configuration X ∈ B 0 to the spatial setting x ∈ B t . The deformation gradient F = ∇ X ϕ, its cofactor cof[F ] = det[F ]F −T , and the Jacobian J = det[F ] map an infinitesimal line, an area element, and a volume element from the material setting to the spatial configuration, respectively. Furthermore, an electric potential φ(X, t) is defined. The outer boundary ∂B = ∂B M is decomposed in terms of the mechanical and the electrical boundary conditions as: respectively. The portion ∂B ϕ denotes the Dirichlet mechanical boundary and ∂B φ constitutes the electrical counterpart. Furthermore, the Neumann boundaries read ∂B t and ∂B q , where mechanical tractions and surface charges can be prescribed, respectively.

Electro-mechanics of deformable continua
The electrical part of the coupled problem is based on two Maxwell equations: the Gauss law for electricity and Faraday's law. In the absence of varying magnetic fields, the associated Maxwell equations can be introduced in terms of the electric displacement and the electric field in the spatial setting, as (2) Fig. 1 Nonlinear mapping from the undeformed to the deformed configuration of an inhomogeneous electro-mechanical body respectively. The operator ∇ x · [•] reads the divergence with respect to the spatial coordinates x and ∇ x × [•] is the curl of a vector at the deformed setting. Moreover, t denotes prescribed free charges per unit current volume [10]. Equation (2.2) enables expressing the spatial electric field e as the negative gradient of a scalar electric potential φ, where Regarding the electrical state at the interface ∂B I , the jump conditions (4) apply [33], where • := [•] out − [•] in and q I denotes the density of free surface charges at the interface ∂B I . Due to the fact that the influence of the surrounding space is neglected in this contribution, the electrical conditions at the outer boundary ∂B M = ∂B can be expressed as (5) with the density of free charges q M at the outer boundary ∂B q t . Regarding the continuity of the scalar electric potential φ on the interface ∂B I , the following condition applies: Moreover, a scalar electric potential φ =φ can be prescribed on the outer Dirichlet boundary ∂B φ . The electrical polarization of materials is expressed mathematically through the generalized relation: with the contribution 0 e as the electric displacement in vacuum and p as the polarization of deformable matter at the current configuration. In the free space, the relation applies, with the constant electric permittivity of the vacuum ( 0 = 8.854 × 10 −12 N/V 2 ), which in turn expresses the electric response in free space as linear. On the other hand, the relation between the polarization and the electric field e can take a nonlinear form. The electric field, electric displacement, and polarization at the reference configuration are linked to their spatial counterparts, through the pull-backs: (8) respectively. As deformation measures, we introduce the right and the left Cauchy-Green tensors as where g is the Eulerian metric tensor and G denotes the Lagrangian metric tensor. Both tensors g and G are covariant and symmetric. The balance of linear momentum is written in terms of the total Cauchy stresses σ , with neglecting the inertial effects as with the volume-specific mechanical body forces ργ at the spatial setting. In the context of coupled electro-mechanics, the total stress tensor σ includes a contribution that is originated due to electric body forces [19,34]. The total first Piola-Kirchhoff stress tensor P can be expressed as a pull-back of the total Cauchy stress tensor σ , where The Dirichlet and Neumann mechanical jump conditions at the interface ∂B I can be expressed as respectively, with the traction t I at the boundary ∂B I t . Moreover, mechanical boundary conditions can be prescribed on the outer boundary ∂B M = ∂B as

Constitutive model for electro-elasticity
This section is devoted to demonstrate a material model that describes the coupled electro-elastic behavior of EAP at large strains. As a first step, we introduce the constitutive equations and the terms needed for the finite element implementation in their general form. Subsequently, we present the material model and the associated algorithmic expressions.

Generalized constitutive description and material-tangent moduli
The coupled material response is expressed by introducing a total energy-enthalpy density functionˆ tot , which is decomposed into a purely hyperelastic energy contribution hypel and an electro-elastic enthalpy density¯ coup as (14) We introduce the multiplicative decomposition of the total deformation gradient F into a volumetric F vol and an isochoric contributionF as follows: where 1 reads the second-order identity tensor. Based on the volumetric-isochoric split of the deformation in Eq. 15, we define an additive split of the hyperelastic energy contribution hypel into where vol denotes a purely volumetric part and iso constitutes a volume-preserving contribution of the energy density.
Referring to the decompositions of the total energy-enthalpy density function as it is shown in Eqs. 14 and 16, the total Cauchy stresses σ can be decomposed and introduced as The electric displacement vector d at the current configuration is given as (18) In order to implement the constitutive model within a fully monolithic finite element framework, we introduce the material-tangent moduli as follows: (19) with the fourth-order tensor Cϕϕ, the third-order coupled material tensors Cϕφ and Cφϕ. Moreover, Cφφ denotes the second-order material tensor. The associated Lie derivatives of the coupled problem are shown in [25]. According to the splits of the total energy-enthalpy density function as demonstrated in Eqs. 14 and 16, the decomposition of the fourth-order material tensor Cϕϕ takes the form: with the volumetric part Cϕϕ vol , the isochoric contribution Cϕϕ iso , and the coupled part Cϕϕ coup .

Hyperelastic material model
The description of the hyperelastic response of the material is achieved by specifying the volumetric part vol and the isochoric part iso of the total energy-enthalpy function. The contribution vol is chosen in a way that allows the simulation of both compressible and quasi-incompressible materials, and is introduced as: where κ is the bulk modulus. The volumetric Cauchy stresses σ vol and the associated moduli Cϕϕ vol can be computed as: (22) with the fourth-order identity tensor I A ij kl = (A ik A jl + A il A jk )/2, where the second-order arbitrary tensor A is defined as A = A ij e i ⊗ e j in terms of the Cartesian basis e i . The extended-tube model [32] is utilized to describe the isochoric part iso of the total energy density function. Taking into regard the limited chain extensibility of the network chains and the topological tube constraints, the isochoric contribution iso is split into a cross-linking contribution W c and a topological tube constraints part L e , where iso C = W c ĪC + L e λ a , with the cross-linking G c and the topological constraints G e shear moduli. Limited chain extensibility is expressed by the material parameter δ and the topological constraints are taken into account through the parameter β. Moreover,ĪC = trC reads the first invariant of the isochoric right Cauchy-Green tensor, withC = J −2/3 C. For the computation of the isochoric stresses σ iso and tangent moduli Cϕϕ iso , we refer to [35,36].

Dielectrically quasi-linear material response
The electric response of the material can expressed as quasi-linear regarding the relation between the current electric field e and the current electric displacement d with density-dependent electric permittivity, by specifying¯ coup as (24) where = m 0 reads the electric permittivity of the material, with m as the associated dielectric constant. Referring to the definition in Eq. 18, the electric displacement d takes the form: with the density-dependent permittivityˆ (J ) = J −1 . The coupled stress σ coup can be expressed as (26) Consequently, using Eq. 19, the associated material-tangent terms are introduced as

Nonlinear electric response
Referring to [21], we describe electro-elasticity with a nonlinear deformation-dependent electrical and electro-mechanical response of the material by introducing the invariants (28) and setting the coupled enthalpy function¯ coup as with the material parameter c 1 that has an influence on the polarization of the material and the parameter c 2 which influences the polarization of the material and describes the corresponding electro-mechanical coupling. Consequently, the electric induction can be introduced as In what follows of this section, we set g = 1 which is the case for Cartesian coordinates. Regarding the electro-mechanical coupling, we introduce the stresses σ coup as For the sake of completeness with regard to the finite element implementation of the model, we refer to Eq. 19 and compute the tangential material tensors as: (32)

Finite element implementation
In this section, we first explain the weak forms of the main governing equations for the coupled electro-mechanical problem under consideration. Thereafter, in the context of a fully monolithic finite element, the linearizations of the weak forms are presented. In a last step, we outline the discretization procedure of the introduced finite element approach.

Weak formulation and its linearization
Following the method of weighted residuals, we multiply the strong forms in Eqs. 2.1 and 10 by the virtual test functions δϕ and δφ, where both test functions vanish at the boundaries ∂B ϕ and ∂B φ , respectively. Subsequently, an integration by parts is applied and the Gauss theorem is used, in order to introduce the weak forms of Eqs. 10 and 2.1 as (33) respectively, where the terms associated with the free volume charges t and the mechanical body forces ργ are not demonstrated. Regarding (33), in the following, we only consider the treatment of the terms within the solid domain B. As a first step toward the consistent linearization of the coupled problem, we define the increments ϕ and φ that correspond to the main state variables ϕ and φ, respectively. Subsequently, we utilize the directional derivative of the weak form expression in Eq. 33.1 along ϕ and its counterpart in Eq. 33.2 along φ, to introduce the linearized forms: with the expressions where L reads the linearization and D denotes the directional derivative of the residual terms. The identities ∇ x δϕ = ∇ X δϕ F −1 and F −1 = −F −1 ∇ x ϕ are utilized to express the incremental form of the gradient term ∇ x δϕ as Consequently, the linearized form of the total Cauchy stresses σ is expressed as with the the fourth-order Cϕϕ and third-order Cϕφ material moduli, as they are defined in Eq. 19. Note that based on the material model presented in Section 3, the symmetries Cϕϕ ij kl = Cϕϕ klij , Cϕϕ ij kl = Cϕϕ jilk , and Cϕφ ij k = Cϕφ jik apply. The relations ∇ x δφ = ∇ X δφ F −1 and F −1 = −F −1 ∇ x ϕ are exploited to describe the incremental expression of the gradient term ∇ x δφ as Thereafter, linearization of the electric displacement d can be evaluated as respectively.

Finite element discretization
The finite element method is employed to discretize (33) and (40) with the number of the finite element nodes n en . In analogy to Eq. 41, we discretize the gradients of the increments ϕ e and φ e and the gradients of the test functions δϕ e and δφ e as A discretized form of Eq. 33 can be defined by using the interpolated expressions in Eqs. 41 and 42 as where with the global residual vector R and the operator A that represents the global assembly of R el evaluated at the local nodes, based on the element connectivity. Similarly, we discretize (40) and introduce the global tangent matrix K as where with the element stiffness matrix K el .

Numerical examples
This section is devoted to demonstrate the features of the presented material model within a finite element framework. Examples of both homogeneous and inhomogeneous electro-mechanical bodies are depicted. For all the considered electromechanical simulations, an electric field develops due to the application of a voltage difference between two electrodes that cover the surfaces of the body under consideration. However, the influence of the electrodes is neglected, assuming that their thickness is relatively small compared to the thickness of the body, which is usually the case. For all the finite element analyses, we use four-noded displacement-based tetrahedral finite elements. The first example is meant to demonstrate an analytical solution of the proposed model, where analytical results are compared against solutions that are obtained using the finite element method (FEM). In the second example, we emulate the response of inhomogeneous electro-active bodies that are composed of a soft matrix filled with a stiff inclusion having a spherical shape, where the influence of the inclusion volume fraction on the simulated response is investigated, for both the density-dependent and the nonlinear deformationdependent electro-mechanical models as shown in Sections 3.3 and 3.4, respectively. As a third example, we illustrate the influence of ellipsoidal inclusions with different aspect ratios on the electro-mechanical behavior of the considered structure.

Analytical solution of electro-mechanical Neo-Hookean and extended tube models
We simulate the electro-mechanical response of a homogeneous cubic structure using both analytical and numerical approaches. Regarding the electro-mechanical coupling, we adopt the nonlinear electro-elastic model as shown in Section 3.4. In our analysis, we drive the structure by applying an electric field and estimate the resulting stretch along the direction of the applied electric field. Regarding the analytical solution of an incompressible material with homogeneous state of deformation and one-dimensional electric field, we introduce the right Cauchy-Green deformation tensor C and the electric field E at the reference configuration as with the stretch λ along the direction of the applied electric field E. We define a normalized energy-enthalpy density function For an incompressible Neo-Hookean material subjected to a homogeneous deformation, the hyperelastic energy density function can be written in its normalized form as Considering the extended tube model as described in Eq. 23, we introduce its simplified and normalized counterpart for incompressible materials with homogeneous boundary conditions and write with the ratios R c = G c /G [−], R e = G e /G [−], and the total shear modulus G = G c + G e . Furthermore, the coupled enthalpy density function as defined by Eqs. 28 and 29 can be rewritten in a normalized and reduced form for incompressible materials with homogeneous stress state as where As the coupled electro-mechanical problem under consideration is stress-free where no external mechanical forces are applied, the one-dimensional total first Piola-Kirchhoff stress P vanishes and is introduced as with the hyperelastic P hypel and the coupled P coup stresses. Note that P coup is referred to as the Maxwell stress in some contributions; see for example [9]. We substitute (51) into (52.3) to yield the specific form of the coupled stress Considering the Neo-Hookean material model, the hyperelastic part of the stresses P hypel can be computed through inserting (49) in Eq. 52.2, which yields In order to perform the analytical solution of the coupled Neo-Hookean model, we substitute (53) and (54) into (52.1), which leads to a closed-form expression for the stretch λ in terms of the normalized electric fieldĒ, where Regarding the solution of the extended tube model, we use Eqs. 50 and 52.2 to express the hyper-elastic portion of the stresses as P hypel = P W + P L , Subsequently, we insert (53) and (56)  Regarding the finite element analyses, we specify the bulk modulus parameter κ to emulate quasi-incompressible behavior of the material, such that the initial Poisson ratio yields ν 0 = 0.499 [−]. We discretize a cubic structure with dimensions of L × L × L using (m × m × m)×6 tetrahedral finite elements with m = 7. In order to verify the insensitivity of the finite element formulation toward mesh distortion, we randomly displace three nodes on each of the six inner surfaces of the discretized cubic structure. For further illustration, Fig. 2 shows the distortion at two different inner surfaces of the body. Moreover, the whole discretized structure is demonstrated in Fig. 2c. Figure 3 visualizes a finite element simulation of the considered problem, where a voltage difference between the two surfaces of the cube induces a one-dimensional electric field. As it is demonstrated in Fig. 4, the analytical and the numerical solutions of the considered problem fit together to

Composites with spherical inclusions
We study the electro-mechanical response of inhomogeneous bodies consisting of a soft matrix and a stiff spherical inclusion with different volume fractions f of the inclusion. In order to express the global electric response of the considered with the total reference volume V of the composite. Moreover, in order to evaluate the actuation or the intensity of the electro-mechanical coupling, we study the relation betweenẼ [−] and the averaged deformation gradientF of the whole structure using , the material parameters are chosen such that they denote base electric and coupled properties of the material which are equivalent to the assumed properties using the parameters = 0 with respect to matrix and = 1000 0 with regard to the inclusion, for the only density-dependent counterpart (Section 3.3). Therefore, we utilize the parameters c 1 = − 0 and c 2 = 0.5 0 to express the coupled response in the matrix and we set the parameters as c 1 = −1000 0 and c 2 = 500 0 to describe the electro-mechanical interaction within the stiff inclusion. The chosen parameters indicate that the inclusion is stiffer and has higher electric permittivity than the matrix. The electro-mechanical analysis is driven by applying an electric potential difference. Moreover, the mechanical boundary conditions are assigned such that the boundaries of the full composite are unconstrained. As the applied boundary conditions are symmetric, we simulate only one-eighth of the whole structure. Note that similar mechanical boundary conditions are adopted in [9]. The structures are discretized using tetrahedral finite elements with geometry adaptive sizing of the elements. The minimum and the maximum sizes of the used elements are the same for all composites with different volume fractions of the spherical inclusion. Finite element meshes of the composite with f = 5% and f = 20% are illustrated in Fig. 5, where the whole inclusion and one-half of the matrix are shown. Regarding the analysis results, Fig. 6 depicts the deformed state of the inhomogeneous structures with varying volume fractions of the spherical inclusion. It can be observed that the distribution of the electric voltage is disturbed due to the existence of the high electric permittive inclusion. Moreover, the electric fields tend to zero within the inclusion compared to the fields within the region of the soft matrix. See the distribution of the  electric potential φ and the electric field E in Fig. 6. In [9], the same coupled model as in Section 3.3 is utilized along with a hyperelastic Neo-Hookean material model to perform stability analyses of similar configurations as the ones presented in this section. It is discussed in [9] that an analysis of the structural configuration under consideration using the coupled ansatz as in Section 3.3, with the simulation being driven by an electric fieldẼ z , allows for performing the analysis up to the point where electro-mechanical instability occurs. Therefore, the authors have suggested that in order to trace the response after the occurrence of instability, the electro-mechanical response of the structure should be controlled by applying an electric displacementD z [9]. Moreover, they [9] have demonstrated that for the considered configurations with the assumption that the material exhibits a dielectrical quasi-linear electro-mechanical behavior as described in Section 3.3, the instability point is reached for levels of the normalized electric fieldẼ z < 1. In this contribution, we use the dielectrical quasi-linear model as shown in Section 3.3 and the nonlinear electro-elastic description as expressed in Section 3.4, both coupled to the extended tube model to perform our study. Moreover, as stability analysis is not our topic, we drive the numerical simulations by applying an electric fieldẼ z within a rangeẼ z ≤ 0.45, before all the considered structures become electro-mechanically unstable. The influence of the inclusion size with respect to the soft matrix on the total actuation intensity and the global electric response of the composite is illustrated in Fig. 7. We can observe from theẼ z -F z relation as it is shown in Fig. 7a that the actuation capability gets improved as the relative size of the inclusion is increased in terms of f . Moreover, the   Fig. 7b indicate that as f increases, larger electric displacementsD z result in response to the same appliedẼ z . Figure 8 illustrates the results of the simulations for composites assuming that the electroelastic coupling is based on the description as in Section 3.4 for a relatively large range of electric fields (i.e.,Ẽ z > 1). Regarding the simulation of the composites assuming nonlinear dependency of the electric properties on deformation, electro-mechanical instability does not occur at least at low levels of the electric field (i.e.,Ẽ z < 1), which is not the case for the earlier considered dielectrical quasi-linear materials. From a pure engineering point of view, instability at relatively low levels of electric loading does not take place in electrically nonlinear materials because the material adapts its electric and coupled properties in response to the deformation, which leads to overpassing the instability point. Figure 8a depicts that at low levels of the electric fieldẼ z , the actuation intensity is in general enhanced as f is set larger. However, at higher levels of the electric field, the actuation behavior is diminished as the relative size of the inclusion increases. On the other hand, the results of the electric response as shown in Fig. 8b imply that as f increases, larger electric displacementsD z result with regard to the same appliedẼ z . In order to compare the results of the dielectrical quasi-linear material and its electrically  nonlinear counterpart, we extract the results of Fig. 8, in which the maximum applied electric field isẼ z = 0.45 [−] and demonstrate them in Fig. 9. We can notice from Fig. 9a that a general trend of actuation improvement is observed, which is caused by increasing the volume fraction of the inclusion. However, comparing Figs. 7a to 9a leads to the conclusion that dielectrically quasi-linear materials are more sensitive toward varying volume fractions of the spherical inclusion, compared to their counterparts with nonlinearly deformation-dependent electro-mechanical properties.

Composites with ellipsoidal inclusions
We investigate the influence of the aspect ratio r of ellipsoidal inclusions on the electro-mechanical behavior of inhomogeneous bodies, utilizing the electro-mechanical model with deformation-dependent electric properties; see Section 3.4. We consider a fixed volume fraction f = 10% and varying aspect ratios r= {1.00, 1.25, 1.50, 1.75, and 2.00 [−]}. Using a geometry adaptive meshing scheme, all composites with different aspect ratios of the ellipsoidal inclusion are discretized using tetrahedral finite elements. Moreover, we apply the same minimum and maximum sizes of the finite elements for all the considered examples. Figure  Regarding the results of the analyses, Fig. 11a shows that at low electric field levels, the actuation is slightly improved as the aspect ratio r increases. However, the actuation intensity at higher levels of the electric field is reduced as r is set larger. Nevertheless, relatively, no significant influence of changing the aspect ratio on the actuation behavior is found. Furthermore, we can notice from Fig. 11b that higher electric displacements are obtained for the same electric field as the aspect ratio of the inclusion r is larger. The finite element simulations for composites with r= {1.00, 1.25, 1.50, and 1.75 [−]} are shown in Fig. 12.

Conclusion
In the presented work, we have proposed a constitutive material description and the associated finite element implementation, which are together adequate for the simulation of the coupled electro-mechanical response of inhomogeneous EAP at large deformations. The modelling framework includes the hyperelastic extended tube model coupled to two distinct electromechanical mathematical descriptions. The capabilities of our suggested framework are demonstrated through analytical and numerical examples, with emphasis on heterogeneous bodies. In future contributions, we aim to employ the presented model to study material and structural electro-mechanical instabilities in EAP.