Thermo-Electro-Mechanical Simulation of Electro-Active Composites

In this contribution, a computational thermo-electro-mechanical framework is considered, to simulate coupling between the mechanical, electrical and thermal fields, in nonhomogeneous electro-active materials. A thermo-electro-mechanical material model and a mixed Q1P0 finite element framework are described and used for the simulations. Finite element simulations of the response of heterogeneous structures consisting of a soft matrix and a stiff incluison are considered. The behavior of the composite material is studied for varying initial temperatures, different volume fractions and various aspect ratios of the inclusion. For some of the examples, the response of the structure beyond a limit point of electro-mechanical instability is traced. Regarding the soft matrix of the composite, thermal properties of silicone rubber at normal conditions have been obtained by molecular dynamics (MD) simulations. The material parameters obtained by MD simulations are used within the finite element simulations.


Introduction
Electro-active polymers (EAP) are smart materials that can be favorably used in artificial muscles and soft robotics. That is due to EAP's relatively quasi-instantaneous mechanical actuation in response to an applied electric field [1]. Furthermore, EAP have demonstrated a capability of exhibiting relatively large strains, where in specific types of EAP, strains up to more than 300% are recorded [2]. Many prototypes of soft robotics have been mainly based on EAP, see for instance [3][4][5]. Dielectric elastomer actuator (DEA) is a class of EAP, which demonstrates electrostrictive behavior [6]. Several researchers have investigated the influence of filling DEA's soft matrices with stiff particles that possess relatively large electric permittivity [7][8][9][10]. It has been demonstrated that the intensity of electro-mechanical coupling in particle-filled DEA is enhanced as the volume fraction of embedded particles increases [8,9]. The influence of embedding stiff fillers within soft DEA has been numerically investigated in [10][11][12][13].
Computational modeling of electro-mechanical interactions in EAP relies either on a purely energy-based approach or, alternatively, on the relatively computationally efficient mixed energy-enthalpy-based approach [12]. In the context of the finite element method, the energy-based approach requires setting a vector electric potential (three degrees of freedom) to describe the electrical part of the problem. The mixed energy-enthalpy approach depends on a scalar electric potential (one degree of freedom) to express the electrical sub-problem.
As an example of numerical solutions of the associated problem, a finite element implementation of electro-elasticity at large deformations has been implemented in [23]. Numerical investigations of composites composed of dielectric soft matrix and stiff inclusions have been performed in [12,24]. It has been shown that through driving the simulation using surface charges, the response of EAP beyond the occurrence of electro-mechanical instability can be simulated using the more convenient mixed energy-enthalpy-based approach [12]. A mixed finite element formulation for composites consisting of quasi-incompressible soft matrix and stiff fibers has been proposed in [25,26].
The actuation behavior of EAP is mainly based on coupling between the electrical field and the mechanical field. Nevertheless, temperature changes resulting due to self-induced heating or cooling, and varying ambient conditions can influence EAP's properties and its actuation response. Therefore, the simulation of interactions between thermal, electrical and mechanical fields have been previously considered, through proposing a thermo-electromechanical model [27] and a thermo-electro-mechanical finite element framework [28,29]. The finite element formulation developed in [28] is a standard displacement-based formulation, and it relies on a partially-staggered solution approach.
In this work, finite element modeling of thermo-electro-mechanical interactions in nonhomogeneous electro-active materials is considered. Regarding the constitutive model used in the present work, a material description that takes into account a quasi-linear dielectric response [25] and an entropic thermo-elastic coupling [30] is adopted. A mixed Q1P0 thermo-electro-mechanical framework is implemented and used for the simulations. Regarding the estimation of thermal properties, thermal conductivity, heat capacity and coefficient of thermal expansion of silicone rubber at normal conditions, are estimated using molecular dynamics (MD) simulations, which is an alternative approach to experimental investigations for the prediction of material properties. In this work, the reason of using MD simulations is to avoid conducting time consuming experiments and to demonstrate an efficient coupling of various numerical approaches, where both MD and finite element simulations are performed to carry out the associated study. Numerical simulations of heterogeneous structures with spherical inclusions under varying initial temperature is performed. The effect of varying volume fraction and changing aspect ratio of a stiff inclusion on the overall composite's response are numerically investigated. Both polarization behavior and intensity of electro-mechanical coupling are evaluated. All simulation examples are driven by applying surface charges. For some of the performed analyses, the structural response beyond the point of electro-mechanical instability is traced.
The present paper is organized as follows. Section 2 introduces the constitutive modeling approach and the associated finite element implementation. In Section 3, the method used to determine thermal properties of silicone rubber is described. In Section 4, several numerical examples are demonstrated. Section 5 ends the paper with a brief summary.

Thermo-Electro-Mechanics
This section is devoted to present a constitutive formulation and its associated numerical treatment, which are used to emulate thermo-electro-mechanical interactions at large deformations. Principles of continuum mechanics are demonstrated and employed within a coupled thermo-electro-mechanical environment. An overview of the kinematics and balance laws of the coupled problem is outlined. Subsequently, a constitutive model is presented. Finally, the used finite element formulation is described.

Preliminaries
The nonlinear deformation map x = ϕ(X, t) at time t ∈ T projects points X in the reference configuration B 0 to their counterparts x in the current setting B t . The deformation gradient F = ∇ X ϕ, its cofactor cof[F] = det[F]F −T and its Jacobian J = det[F] project an infinitesimal line element, area element and volume element from the reference setting to the current setting, respectively. The reference density ρ 0 is connected to the current density ρ by ρ 0 = Jρ. Furthermore, a scalar electric potential φ(x, t) and an absolute temperature θ(x, t) are defined. The outer boundary ∂B t is decomposed in terms of the mechanical, electrical and thermal boundary conditions as respectively. ∂B ϕ t , ∂B φ t and ∂B θ t are the surfaces where Dirichlet mechanical, electric and thermal boundary conditions can be prescribed, respectively. The Neumann boundaries, where mechanical tractions t, surface charge densities w and surface heat flow q θ can be prescribed, read ∂B t t , ∂B w t and ∂B q θ t , respectively.

Mechanical Field
The right Cauchy-Green deformation tensor is introduced as where g is the covariant Eulerian metric tensor. The strong form of the balance of linear momentum for quasi-static problems can be defined as with the divergence with respect to the spatial coordinates ∇ x · [ · ], the total Cauchy stresses σ and the mechanical body forces per deformed unit volume ργ mec . The total stresses can be split into mechanical and ponderomotive stresses as σ = σ mec + σ pon . The ponderomotive stresses σ pon express electro-mechanical coupling in the material. For further details, the reader is referred to [20,31,32].

Electrical Field
In the case of electro-statics, where magnetic fields are absent, the electric field at the reference configuration E and the electric field at the current configuration e are defined as where ∇ X [ · ] denotes the differential operator with respect to the reference coordinates X, and ∇ x [ · ] is its counterpart with respect to the current coordinates x. Gauss's law for electricity can be written with its differential form in terms of the current electric where f is the free charges per deformed unit volume. In insulators, there are no free volume charges ( f = 0). The electric displacement in the deformed setting d is connected to its counterpart in the undeformed setting D by The surface charge densities at the undeformed and at the deformed setting are defined as with the unit normal vector N on the surface ∂B 0 and the unit normal vector n on the surface ∂B t .

Thermal Field
The heat flux at the reference configuration Q and its counterpart at the current configuration q are introduced as respectively, where k is a parameter of heat conductivity. The evolution of temperature is expressed by the transient heat equation, which reads where ρc is the heat capacity, ρr is the heat source, ρw ext is the external power, all introduced per unit deformed volume. The term of external power ρw ext can be decomposed into mechanical, electro-mechanical and electrical parts as where l =ḞF −1 is the spatial velocity gradient, d = (l + l T )/2 is its symmetric part andė denotes the rate of change of the current electric field with respect to time. The terms ρw mec ext , ρw coup ext and ρw elec ext express the work done by the mechanical stresses σ mec , the ponderomotive stresses σ pon and the electric displacement d, respectively.

Material Model
The deformation gradient is multiplicatively decomposed into purely volumetric and isochoric (volume-preserving) contributions as F = F volF , where the volumetric contribution is defined by F vol = J 1/3 1 and the isochoric contribution is given byF = J −1/3 F. The isochoric part of the right Cauchy-Green deformation tensor is introduced as C =F T gF. Electro-mechanical and entropic thermo-mechanical couplings of the material are expressed through the energy density function where θ 0 is the reference temperature. The volumetric contribution of the mechanical energy density is introduced as where κ is the bulk modulus. The isochoric part Ψ iso e of the mechanical energy function is specified by the Neo-Hookean model as where G c is the shear modulus andĪC = trC denotes the first invariant of the isochoric right Cauchy-Green tensor. The thermal expansion is expressed in terms of the reference internal energy e 0 (J) = καθ 0 lnJ (14) with the coefficient of thermal expansion α [30]. The thermal function T 0 (θ) is introduced by where ρ 0 c = −θ∂ 2 θθ T 0 (θ) denotes the heat capacity per unit reference volume. The coupled electro-mechanical contribution of the energy function Ψ coup is considered as it is suggested in [12,25], and it is introduced here as where r [−] denotes the relative electric permittivity of the material and 0 = 8.854 × 10 −12 N/V 2 is the electric permittivity of vacuum. The total Cauchy stresses σ and the ponderomotive stresses σ pon can be computed by respectively. The current electric displacement d is introduced as As it is shown in Equation (11), the electro-mechanical contribution Ψ coup (J, C, E) is not scaled with a temperature field. Therefore, the the ponderomotive stresses σ pon and the electric displacement d as defined in Equations (17) and (18), respectively, are temperature independent. Due to that, the external power terms ρw coup ext and ρw elec ext , as defined in Equation (10), are equal to zero.

Finite Element Formulation
In order to avoid volumetric locking that can arise due to the material's nearly incompressible behavior, a mixed Q1P0 formulation is adopted and implemented within a thermo-electro-mechanical framework. To this end, the mean dilatation D e and the mean pressure p e are introduced as additional field variables and treated as averages on finite element level [25,33]. The method of weighted residuals as previously used in [24,25,30], is employed to formulate a thermo-electro-mechanical finite element framework. To this end, the main strong forms as given in Equations (3), (5) and (9) are scaled with the test functions δϕ, δφ and δθ, respectively. Thereafter, integration by parts is applied and the Gauss theorem is used. Performing the latter two steps allows to express the associated weak forms as In Equation (19), the surface tractions t, the surface charge density w and the surface heat flow q θ are introduced as t = σ · n on ∂B t t , w = −d · n on ∂B w t and q θ = q · n on ∂B q θ t , respectively. The finite element implementation of the problem is performed through discretizing the weak forms in Equation (19) and the associated incremental forms. Therefore, the undeformed domain B 0 is split into finite elements as B 0 ≈ n e e=1 B e 0 , where n e is the total number of solid elements. To treat the material model within the mixed Q1P0 finite element framework, the total energy density function as given by Equation (11) is additively decomposed into two contributions and is reformulated to whereΨ is considered on each Gauss point. The volumetric part Ψ vol tot is parameterized in terms of the mean dilatation D e and can be written as The average dilatation and average mean pressure are evaluated over each finite element as In Equation (19), the total Cauchy stresses are additively decomposed as σ =σ + σ vol , where the partσ = 2J −1 ∂ gΨ is computed at Gauss points and the volumetric contribution σ vol = p e g −1 is evaluated at each element B e t . The method followed to treat the power term ρw ext within the Q1P0 finite element is detailed in [33]. The isoparametric concept of the finite element method is utilized, where the fields ϕ e , φ e , θ e and the associated test functions are interpolated over the finite element B e t using nodal values and Lagrangian shape functions. Regarding the solution scheme adopted, the thermo-electro-mechanical problem is decomposed into an electro-mechanical and a thermal sub-problem. The electromechanical part of the problem is solved monolithically at a fixed temperature. Thereafter, the thermal sub-problem is solved and attached to the whole system using a staggered transfer of solution fields [28]. The evaluation of both sub-problems is based on an exchange of information after each iteration. The presented material model and finite element formulation are implemented within an in-house finite element program.

Estimation of Thermal Properties
Firstly, polymeric chains, which consist of 418 monomer units, have been constructed by the Moltemplate software [34]. The monomer unit of silicone rubber is shown in Figure 1a. After that, ten polymeric chains have been randomly distributed in a periodic supercell (see Figure 2a) by Packmol software [35]. In the next step, the chains have been crosslinked in a canonical ensemble via an algorithm similar to the one presented in [36]. During this procedure, CH 3 groups participating in the creation of a bond between the polymeric chains are turned into CH 2 groups. The crosslink bridge of silicone rubber is shown in Figure 1b. The degree of crosslinking of the obtained model of silicone rubber is equal to 9.8 %. It is defined as given in [37]. The force field description has been used to describe interactions between the atoms, where the CH 3 group has been modeled as one united atom. The total potential energy of a polymeric system is calculated in this approximation as Non-bonded interactions are modeled only with van der Waals interactions. The cutoff distance is set to 10 Å in all simulations. The Lorentz-Berthelot mixing rules are used to find the missing parameters of the Lennard-Jones potential. The force field parameters have been taken from the OPLS-AA (All Atom) [38] force field, except the parameters of Lennard-Jones potential for the CH 3 group, which has been modeled by the OPLS-UA (United Atom) [38] force field. All simulations are carried out using the LAMMPS [39] software package. After the crosslinking procedure, density of the system has been equilibrated in an isothermal-isobaric ensemble at normal conditions for 200 ps. Before calculation of thermal conductivity by the Green-Kubo method, it has been modeled in a canonical ensemble for 100 ps for equilibration of temperature. By the Green-Kubo formula, the thermal conductivity of an isotropic material can be found as where h is the heat flux calculated as where e i is the total energy of i-th atom, S i is the stress tensor of i-th atom, υ i and υ j are velocities of i-th and j-th atoms. The model of silicone rubber has been simulated in a microcanonical ensemble for 3 ns with a time step of 1 fs. The heat flux has been calculated at each time step and the thermal conductivity has been obtained for each correlation time interval, which is equal to 3 ps.
For calculation of specific heat capacity at constant pressure and coefficient of thermal expansion, the model has been simulated in an isothermal-isobaric ensemble at normal conditions for 100 ps and data points at equilibrium have been taken. The specific heat capacity and the coefficient of thermal expansion have been obtained as derivatives of enthalpy and volume as in [40][41][42] The thermal conductivity obtained by the Green-Kubo method is roughly 0.15 W/m/K (see Figure 3), which is close to experimental value (≈0.17 W/m/K [43]). The specific heat capacity at constant pressure found by the MD simulations is around 1500 J/kg/K. For comparison, in [44], it is approximately 1.2 kJ/kg/K. Calculated coefficient of thermal expansion (2.83 × 10 −4 1/K) is close to measured value (roughly 3 × 10 −4 1/K [45]).

Numerical Examples
This section is devoted to present several numerical studies of the coupled thermoelectro-mechanical response of heterogeneous structures using the finite element meth-od. An electro-active composite consisting of soft matrix and stiff inclusion is considered. The material parameters are taken as they are shown in Table 1 Mechanical Dirichlet boundary conditions are assigned to mid-planes of the structure, where the motion along each mid-plane's perpendicular direction is constrained. Similar electrical and mechanical boundary conditions have been previously used in [12]. As the geometry and the boundary conditions are symmetric, only one-eighth of the whole composite is analyzed. All considered configurations are discretized using eight-noded Q1P0 finite elements. The discretization of a structure consisting of soft matrix and ellipsoidal stiff inclusion is shown in Figure 4b. In the first example, the effect of varying initial temperature on the response of a composite with spherical inclusion is studied. The second example is meant to examine the influence of varying volume fraction of the inclusion on the composite's response. In the third example, heterogeneous structures with different aspect ratios of the inclusion are studied. The polarization behavior is evaluated through examining the relation between the electric field and the electric displacement. Moreover, the intensity of electro-mechanical coupling is investigated in terms of the relation between the electric field and the deformation. The global response of the composites is evaluated using normalized and averaged electric field E [−], normalized and averaged electric displacement D [−], and averaged deformation gradient F, which are defined by where here r denotes the relative electric permittivity of the soft matrix, G c is the shear modulus of the soft matrix and V is the total volume of the structure.

Varying Initial Temperature
The influence of varying initial temperature θ init on the electro-mechanical coupling and the polarization behavior of a composite, where entropic thermo-elastic coupling takes place, is investigated. The structure considered here has a spherical inclusion with volume fraction f = 1%. The polarization behavior is represented by the relation between the  Figure 5b demonstrates that increasing the initial temperature leads to less deformation of the composite, in response to the same electric field. The temperature distribution within the composite at an applied electric displacement Figure 6a. The distribution of the electric potential φ is depicted in Figure 6b, where it can be noticed that there is a discontinuity in the distribution of φ. Figure 6c demonstrates that the electric fields within the inclusion tend to zero, compared to the fields within the region of soft matrix. That is due to the inclusion's relatively high electric permittivity. The vectors in Figure 6c indicate the direction of the electric field at the reference configuration E z . Figure 6d shows the distribution of the vertical displacement within both, the soft matrix with relatively low shear modulus and the stiff inclusion with relatively high shear modulus.

Varying Volume Fraction of the Inclusion
The influence of different volume fractions f = {1%, 5%, 10% and 15%} of a spherical inclusion with respect to the whole composite, is studied in this section. Figure 7a Figure 10a demonstrate that the influence of r on the D z − E z relation is relatively insignificant. Figure 10b depicts the influence of r on the relation between E z [−] and F zz [−]. It is noticed from Figure 10b that for the same electric field E z [−], the deformation in z− direction increases slightly by increasing r.

Summary
In the present work, finite element analyzes of thermo-electro-mechanical interactions in heterogeneous structures are performed. The considered structures are composed of a soft matrix and a stiff inclusion. The influences of varying initial temperatures, different volume fractions and varying aspect ratios of the stiff inclusion on the overall response of the composite are studied. For the evaluation of the material's response, both the polarization and the electro-mechanical coupling are evaluated. The thermal material parameters of the soft matrix, namely, thermal conductivity, heat capacity and coefficient of thermal expansion are identified based on molecular dynamics (MD) simulations of silicone rubber. In future contributions, it is aimed to implement the presented thermo-electromechanical model within a multi-scale homogenization framework. That allows to analyze multi-physical response of structures used for industrial applications, on multiple scales. Furthermore, it is aimed to validate the computational model using experimental data.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

EAP
Electro-active polymers DEA Dielectric elastomer actuators MD Molecular dynamics