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.

1. 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].

During the last century, several contributions have been devoted to the mathematical description of polarization and electro-elastic coupling [14,15,16,17]. Thereafter, more recent works are proposed to derive mathematical expressions describing different forms of electro-elastic interactions [18,19,20,21,22].

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-electro-mechanical 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.

2. 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.

2.1. Preliminaries

The nonlinear deformation map at time projects points in the reference configuration to their counterparts in the current setting . The deformation gradient , its cofactor and its Jacobian project an infinitesimal line element, area element and volume element from the reference setting to the current setting, respectively. The reference density is connected to the current density by . Furthermore, a scalar electric potential and an absolute temperature are defined. The outer boundary is decomposed in terms of the mechanical, electrical and thermal boundary conditions as respectively. , and are the surfaces where Dirichlet mechanical, electric and thermal boundary conditions can be prescribed, respectively. The Neumann boundaries, where mechanical tractions , surface charge densities w and surface heat flow can be prescribed, read , and , respectively.

2.2. Mechanical Field

The right Cauchy-Green deformation tensor is introduced as where 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 , the total Cauchy stresses and the mechanical body forces per deformed unit volume . The total stresses can be split into mechanical and ponderomotive stresses as . The ponderomotive stresses express electro-mechanical coupling in the material. For further details, the reader is referred to [20,31,32].

2.3. Electrical Field

In the case of electro-statics, where magnetic fields are absent, the electric field at the reference configuration and the electric field at the current configuration are defined as where denotes the differential operator with respect to the reference coordinates , and is its counterpart with respect to the current coordinates . Gauss’s law for electricity can be written with its differential form in terms of the current electric displacement as where is the free charges per deformed unit volume. In insulators, there are no free volume charges (. The electric displacement in the deformed setting is connected to its counterpart in the undeformed setting by The surface charge densities at the undeformed and at the deformed setting are defined as with the unit normal vector on the surface and the unit normal vector on the surface .

2.4. Thermal Field

The heat flux at the reference configuration and its counterpart at the current configuration 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 is the heat capacity, is the heat source, is the external power, all introduced per unit deformed volume. The term of external power can be decomposed into mechanical, electro-mechanical and electrical parts as where is the spatial velocity gradient, is its symmetric part and denotes the rate of change of the current electric field with respect to time. The terms , and express the work done by the mechanical stresses , the ponderomotive stresses and the electric displacement , respectively.

2.5. Material Model

The deformation gradient is multiplicatively decomposed into purely volumetric and isochoric (volume-preserving) contributions as , where the volumetric contribution is defined by and the isochoric contribution is given by . The isochoric part of the right Cauchy-Green deformation tensor is introduced as . Electro-mechanical and entropic thermo-mechanical couplings of the material are expressed through the energy density function where is the reference temperature. The volumetric contribution of the mechanical energy density is introduced as where is the bulk modulus. The isochoric part of the mechanical energy function is specified by the Neo-Hookean model as where is the shear modulus and denotes the first invariant of the isochoric right Cauchy-Green tensor. The thermal expansion is expressed in terms of the reference internal energy with the coefficient of thermal expansion [30]. The thermal function is introduced by where denotes the heat capacity per unit reference volume. The coupled electro-mechanical contribution of the energy function is considered as it is suggested in [12,25], and it is introduced here as where denotes the relative electric permittivity of the material and is the electric permittivity of vacuum. The total Cauchy stresses and the ponderomotive stresses can be computed by respectively. The current electric displacement is introduced as

As it is shown in Equation (11), the electro-mechanical contribution is not scaled with a temperature field. Therefore, the the ponderomotive stresses and the electric displacement as defined in Equations (17) and (18), respectively, are temperature independent. Due to that, the external power terms and , as defined in Equation (10), are equal to zero.

2.6. 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 and the mean pressure 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 , the surface charge density w and the surface heat flow are introduced as on , on and on , 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 is split into finite elements as , where 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 is parameterized in terms of the mean dilatation 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 , where the part is computed at Gauss points and the volumetric contribution is evaluated at each element . The method followed to treat the power term within the Q1P0 finite element is detailed in [33]. The isoparametric concept of the finite element method is utilized, where the fields , , and the associated test functions are interpolated over the finite element 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 electro-mechanical 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.

3. 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 groups participating in the creation of a bond between the polymeric chains are turned into CH 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].

Figure 1

(a) Monomer unit of silicone; (b) crosslink bridge between silicone chains.

Figure 2

(a) Randomly distributed silicone chains in a periodic supercell; (b) the silicone chains after crosslinking procedure; (c) model of silicone rubber.

The force field description has been used to describe interactions between the atoms, where the CH 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 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 is the total energy of i-th atom, is the stress tensor of i-th atom, and 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 × 1/K) is close to measured value (roughly 3 1/K [45]).

Figure 3

Thermal conductivity of silicone rubber at normal conditions as function of correlation time in picosecond.

4. Numerical Examples

This section is devoted to present several numerical studies of the coupled thermo-electro-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. An initial Poisson’s ratio is considered. The simulations are driven by applying electrical Neumann boundary conditions, where surface charges are prescribed on the top and the bottom surfaces of the composite, as demonstrated in Figure 4a. That permits to trace the structural response beyond an electro-mechanical instability point. The surface charge density W is defined in terms of volume-averaged electric displacement , as indicated in Figure 4a. 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.

Table 1

Parameters for the finite element simulations.

soft matrix
G_c=0.0473 MPa,   ϵ_r=3.0[−],  α=2.83×10^−4K^−1, k=0.15N/(Ks),   θ_0=298K,   ρ_0c=1.438N/(mm^2K).
stiff inclusion
G_c=47.3 MPa,     ϵ_r=10^3[−],  α=1.57×10^−5K^−1,k=6.0N/(Ks),   θ_0=298K,    ρ_0c=3.190N/(mm^2K).

Figure 4

(a) Schematic representation of the structure with applied surface charge density W, (b) finite element mesh of the structure showing one-half of the soft matrix and the full inclusion with a volume fraction and an aspect ratio .

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 , normalized and averaged electric displacement , and averaged deformation gradient , which are defined by where here denotes the relative electric permittivity of the soft matrix, is the shear modulus of the soft matrix and V is the total volume of the structure.

4.1. Varying Initial Temperature

The influence of varying initial temperature 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 . The polarization behavior is represented by the relation between the dimensionless electric displacement and the dimensionless electric field , both evaluated along z-direction. Figure 5a shows that as the temperature increases, lower values of electric field result in response to the same applied electric displacement . Electro-mechanical coupling is evaluated using the relation between and the averaged stretch in z-direction . Figure 5b demonstrates that increasing the initial temperature leads to less deformation of the composite, in response to the same electric field.

Figure 5

Plot of (a) averaged and normalized electric displacement and (b) averaged deformation versus averaged and normalized electric field with varying applied initial temperature .

The temperature distribution within the composite at an applied electric displacement is shown in 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 . 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.

Figure 6

Finite element simulation of the heterogeneous structure at and Kelvin, showing the contour map of (a) absolute temperature in Kelvin, (b) electric potential in Volt, (c) electric field in Volt per millimeter and (d) displacement in millimeter.

4.2. Varying Volume Fraction of the Inclusion

The influence of different volume fractions of a spherical inclusion with respect to the whole composite, is studied in this section. Figure 7a,b demonstrate finite element discretizations of a composite with and a composite with , respectively. The results show that increasing the volume fraction f leads to higher polarization intensity in terms of for the same value of the electric field , as demonstrated in Figure 8a. Regarding electro-mechanical coupling, Figure 8b shows that the structure deforms more as the volume fraction of the stiff inclusion f increases under the same electric field .

Figure 7

Finite element mesh of one half of the soft matrix and the whole spherical inclusion with a volume fraction (a) and (b) .

Figure 8

Plot of (a) averaged and normalized electric displacement and (b) averaged deformation versus averaged and normalized electric field with different volume fractions f of a spherical inclusion with Kelvin.

4.3. Varying Aspect Ratio of the Inclusion

The effect of varying aspect ratio of the inclusion on the global response of the composite, is studied. Finite element meshes of composites with and are depicted in Figure 9a,b, respectively.

Figure 9

Finite element mesh of one half of the soft matrix and the whole inclusion with aspect ratio (a) and (b) , where .

Figure 10a shows that as the aspect ratio r increases, lower values of correspond to the same applied . The results in Figure 10a demonstrate that the influence of r on the relation is relatively insignificant. Figure 10b depicts the influence of r on the relation between and . It is noticed from Figure 10b that for the same electric field , the deformation in direction increases slightly by increasing r.

Figure 10

Plot of (a) averaged and normalized electric displacement and (b) averaged deformation versus averaged and normalized electric field with different aspect ratios r of the inclusion, where and Kelvin.

5. 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-electro-mechanical 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.