A Thermo-Electro-Viscoelastic Model for Dielectric Elastomers

Abstract

Dielectric elastomers (DEs) are a class of electro-active polymers (EAPs) that can deform under electric stimuli and have great application potential in bionic robots, biomedical devices, energy harvesters, and many other areas due to their outstanding deformation abilities. It has been found that stretching rate, temperature, and electric field have significant effects on the stress-strain relations of DEs, which may result in the failure of DEs in their applications. Thus, this paper aims to develop a thermo-electro-viscoelastic model for DEs at finite deformation and simulate the highly nonlinear stress-strain relations of DEs under various thermo-electro-mechanical loading conditions. To do so, a thermodynamically consistent continuum theoretical framework is developed for thermo-electro-mechanically coupling problems, and then specific constitutive equations are given to describe the thermo-electro-viscoelastic behaviors of DEs. Furthermore, the present model is fitted with the experimental data of VHB4905 to determine a temperature-dependent function of the equilibrium modulus. A comparison of the nonlinear loading-unloading curves between the model prediction and the experimental data of VHB4905 at various thermo-electro-mechanical loading conditions verifies the present model and shows its ability to simulate the thermo-electro-viscoelastic behaviors of DEs. Simultaneously, the results reveal the softening phenomena and the instant pre-stretch induced by temperature and the electric field, respectively. This work is conducive to analyzing the failure of DEs in functionalities and structures from theoretical aspects at various thermo-electro-mechanical conditions.

1. Introduction

Electro-active polymers (EAPs) that deform under electric stimuli are emerging as promising materials for wide applications such as actuators and sensors [1]. Among various types of EAPs, dielectric elastomers (DEs) have attracted more attention in recent years due to their outstanding deformation abilities (about 50–380%), high energy density, fast response speed, light weight, and mechanical compliance [2], and have found their potential applications in bionic robots [3,4], biomedical devices [5,6], energy harvesters [7,8], flexible electronic devices [9,10] and so on. However, DEs may fail in their functionalities and structures when subjected to complicated or extreme environments. For example, a common design of DE film [11,12] or DE tube [13] may wrinkle or buckle because of the deformation induced by the interaction of charges with opposite signs on neighboring surfaces when an electric field is applied along the thickness direction. Besides, temperature fluctuations can influence the rate-dependent stress-strain relation and cause the electric and mechanical instability of DEs by changing their material properties [14,15]. Therefore, studying the thermo-electro-mechanically coupled process of DEs is of great significance in the development of DE materials and devices.

In the past few years, it has been long observed from experiments that the mechanical responses of DEs are highly rate-dependent [16,17]. The viscous dissipation energy, represented by the area between the loading and unloading curves, varies with the stretching rate. Later, a number of experimental tests, including single-step and multi-step relaxation tests and cyclic loading-unloading tests, were performed by Hossain and co-workers [18,19,20] to find that the viscoelastic behavior of the commonly used DE polymer VHB4910 was also sensitive to temperature and electric field. In particular, an increasing temperature can significantly soften the materials. Recently, Mehnert et al. [21] conducted more comprehensive tests, showing that temperature has a more pronounced effect on the viscoelastic behavior of DE polymer VHB4905 than an electric field.

At the same time, some phenomenological or physics-based models have also been established to simulate the viscoelasticity of DEs [22,23,24,25]. However, to the best of our knowledge, most of these models can only consider the thermo-mechanical or electro-mechanical behaviors of DEs. For example, Mehnert et al. [26] proposed an electro-viscoelastic model for DE polymer VHB4905 at room temperature without consideration of temperature fluctuation, while Alkhoury et al. [27] established a thermo-mechanical model to describe the significant viscoelasticity reduction of VHB4910 due to temperature increase, ignoring the effect of an electric field. Very few constitutive models were developed to describe the thermo-electro-viscoelastic behaviors of DEs. Although Mehnert et al. [28] proposed a theoretical model for dielectric elastomers with coupling thermo-electro-mechanical behaviors, for which temperature-dependent scaling functions were introduced to modify the elastic and viscous energy contributions in order to reflect the softening of materials induced by increasing temperature, the scaling functions were inconsistent and lacked a physical basis. Moreover, the used Yeoh-type and Neo-Hookean-type energy functions could not capture the strain-stiffening effect well, and the instant pre-stretch was not simulated in the theoretical model.

Therefore, in this paper, we aim to develop a thermo-electro-viscoelastic model for DEs at finite deformation and verify this model by using available experimental data from the literature. First, a thermodynamically consistent model for thermo-electro-mechanically coupling problems at finite deformation is established by considering the Gauss law and the electric contribution to the energy balance on the foundation of the classical thermo-mechanically coupling theoretical framework. Second, the viscoelasticity of DEs is described by a rheological model consisting of an elastic ground network and several parallel viscous subnetworks, whose elastic deformations are assumed to be incompressible and moduli to be temperature-dependent. The Gent model [29], which takes the strain-stiffening effect into consideration, is employed to represent the equilibrium and non-equilibrium elastic free energies, and Lagrange multipliers are used to impose the elastic deformation incompressibility of both ground networks and subnetworks. Finally, for model verification, we simulate the loading-unloading curves of VHB4905 at various thermo-mechanically and thermo-electric-mechanically coupling conditions and compare theoretical results with experimental data from Mehnert et al. [21] and Liao et al. [20].

The remainder of this paper is organized as follows: In Section 2, a continuum theoretical framework for thermo-electro-mechanically coupling problems is given, and constitutive equations, including the state and evolving equations, are derived. In Section 3, based on the proposed theoretical framework, specific constitutive equations are given and the thermo-electro-viscoelastic behaviors of VHB4905 are modeled. Finally, conclusions are given in Section 4.

2. Theoretical Framework

2.1. Kinematics

Consider a DE body $B_0$ bounded by the surface ${\partial B}_0$, which is defined in a fixed reference configuration and deforms into $B_t$ with the surface ${\partial B}_t$. Then the deformation gradient is given as follows:

$$\mathit{F}={𝛻}_R\mathit{x}$$

(1)

where $\mathit{x}=\mathit{\chi }(\mathit{X},t)$ is a function introduced to map an arbitrary material point $\mathit{X}$ inside $B_0$ into a spatial point $\mathit{x}$ inside $B_t$ and ${𝛻}_R$ is the gradient operator with respect to the coordinates $\mathit{X}$. To describe the viscoelasticity of DE at finite deformation, a one-dimensional analogy rheological model [22,25,30,31,32] is employed, as shown in Figure 1. Here, DE is assumed to comprise of a ground elastic network, with deformation gradient $\mathit{F}$, and $n$ parallel viscous subnetworks, with elastic deformation gradient ${\mathit{F}}_i^e$ and viscous deformation gradient ${\mathit{F}}_i^v$, resulting in the following relation:

$$\mathit{F}={\mathit{F}}_i^e·{\mathit{F}}_i^v, i=1,...,n$$

(2)

The stress due to the deformation of the ground elastic network only depends on the $\mathit{F}$ and cannot be relaxed, while the stress due to the elastic deformation ${\mathit{F}}_i^e$ can be relaxed with time because the viscous deformation gradient ${\mathit{F}}_i^v$ evolves with time.

Accordingly, the right Cauchy–Green deformation tensors, $\mathit{C}={\mathit{F}}^T·\mathit{F}$ and ${\mathit{C}}_i^e={\left({\mathit{F}}_i^e\right)}^T·{\mathit{F}}_i^e$ with the superscript $T$ denoting the transpose operation, are used to measure, respectively, the deformation of the ground network and the elastic deformation of the subnetworks.

2.2. Balance Laws and Entropy Inequality

Let $\mathit{d}$ and ${\rho }^c$ denote the electrical displacement and the total charge density reckoned in the current configuration. Then, the Gauss law can be written as follows:

$${\int }_{\partial \Omega }\mathit{d}·\mathit{n}da-{\int }_{\Omega }{\rho }^{c}dv=0$$

(3)

where $\mathit{n}$ is the outward unit normal vector of the area element $da$ on the surface $\partial \Omega$ with $\Omega$ denoting an arbitrary domain inside $B_t$, and $dv$ is the volume element of $\Omega$. The above equation can be equivalently written in the reference configuration as follows:

$${\int }_{\partial {\Omega }_0}\mathit{D}·\mathit{N}dA-{\int }_{{\Omega }_0}{\rho }_R^{c}dV=0$$

(4)

where $\mathit{D}$ and ${\rho }_R^c$ are, respectively, the electrical displacement and the total charge density reckoned in the reference configuration, $\mathit{N}$ is the outward unit normal vector of the area element $dA$ on the surface $\partial {\Omega }_0$ with ${\Omega }_0$ denoting an arbitrary domain inside $B_0$, and $dV$ is the volume element of ${\Omega }_0$. The relations $\mathit{D}=j\mathit{d}·{\left({\mathit{F}}^T\right)}^{-1}$ and ${\rho }_R^c=j{\rho }^c$, with $j$ as the determinate of $\mathit{F}$ and the superscript ‘$-1\mathrm{’}$ denoting the inverse of a tensor, can be obtained since $\mathit{n}da=j{\left({\mathit{F}}^T\right)}^{-1}·\mathit{N}dA$ and $dv=jdV$.

Using the divergence theorem, we can respectively write Equations (3) and (4) locally as follows:

$$𝛻·\mathit{d}={\rho }^c, {𝛻}_R·\mathit{D}={\rho }_R^c$$

(5)

where $𝛻$ is the gradient operators with respect to the coordinates $\mathit{x}$, with $𝛻={𝛻}_R·{\left({\mathit{F}}^T\right)}^{-1}$. Furthermore, the electric field $\mathit{e}$, reckoned in the current configuration, and its counterpart $\mathit{E}$, reckoned in the reference configuration, are respectively defined by the following:

$$\mathit{e}=-𝛻\mathsf{\Phi }, \mathit{E}=-{𝛻}_R\mathsf{\Phi }$$

(6)

where $\mathsf{\Phi }$ is the electric potential reckoned in both the current and reference configurations, and $\mathit{E}={\mathit{F}}^T·\mathit{e}$.

Neglecting the inertial effects, the balance laws of force and moment are given in the current configuration as follows:

$$\mathit{\sigma }·𝛻+\mathit{b}=0, \mathit{\sigma }={\mathit{\sigma }}^T$$

(7)

where $\mathit{\sigma }$ is the Cauchy stress tensor and $\mathit{b}$ is the body force per unit volume in $\Omega$, which can be given in the reference configuration as follows:

$$\mathit{P}·{𝛻}_R+{\mathit{b}}_R=0, \mathbf{P}·{\mathbf{F}}^T=\mathit{F}·{\mathit{P}}^T$$

(8)

where $\mathit{P}$ is the first P-K stress tensor and ${\mathit{b}}_R$ is the body force per unit volume in ${\Omega }_0$, with the relations $\mathit{P}=j\mathit{\sigma }·{\left({\mathit{F}}^T\right)}^{-1}$ and ${\mathit{b}}_R=j\mathit{b}$.

Let $\epsilon$ and $q$ denote, respectively, the internal energy density and the heat source per unit volume in $\Omega$, and ${\mathit{j}}^q$ is the heat flux per unit area on $\partial \Omega$. Then, the energy balance law in the current configuration is written as follows:

$$\frac{d}{dt}{\int }_{\Omega }\epsilon dv={\int }_{\partial \Omega }\left(\mathit{\sigma }·\mathit{n}·\dot{\mathit{\chi }}-{\mathit{j}}^q·\mathit{n}-\mathsf{\Phi }\dot{\mathit{d}}·\mathit{n}\right)da+{\int }_{\Omega }\left(\mathit{b}·\dot{\mathit{\chi }}+q+\mathsf{\Phi }{\dot{\rho }}^c\right)dv$$

(9)

whose corresponding form in the reference configuration is as follows:

$$\frac{d}{dt}{\int }_{{\Omega }_0}{\epsilon }_{R}dV={\int }_{\partial {\Omega }_0}\left(\mathit{P}·\mathit{N}·\dot{\mathit{\chi }}-{\mathit{j}}_R^q·\mathit{N}-\mathsf{\Phi }\dot{\mathit{D}}·\mathit{N}\right)dA+{\int }_{{\Omega }_0}\left({\mathit{b}}_R·\dot{\mathit{\chi }}+q_R+\mathsf{\Phi }{\dot{\rho }}_R^c\right)dV$$

(10)

where $\frac{d}{dt}$ or, equivalently, a superposed dot represents the material time derivative ${\epsilon }_R$ and $q_R$ denotes, respectively, the internal energy density and the heat source per unit volume in ${\Omega }_0$, and ${\mathit{j}}_R^q$ is the heat flux per unit area on $\partial {\Omega }_0$. The relations ${\epsilon }_R=j\epsilon$, $q_R=jq$, and ${\mathit{j}}_R^q=j{\mathit{j}}^q·{\left({\mathit{F}}^T\right)}^{-1}$ exist for the above two equations. The first and fourth terms on the right-hand side represent the mechanical work from the surface traction force and the body force, respectively; the second and fifth terms represent the thermal energy, respectively, from the heat flow across the surface and the heat source inside the body; and the third and sixth terms represent the electric work from the charge change on the surface and inside the body, respectively. From Equations (5)–(8) and using the divergence theorem, the local forms of Equations (9) and (10) can be respectively derived as follows:

$$\dot{\epsilon }=\mathit{\sigma }:𝛻\dot{\mathit{\chi }}-𝛻·{\mathit{j}}^q+q+\mathit{e}·\dot{\mathit{d}}$$

(11)

$${\dot{\epsilon }}_R=\mathit{P}:\dot{\mathit{F}}-{𝛻}_R·{\mathit{j}}_R^q+q_R+\mathit{E}·\dot{\mathit{D}}$$

(12)

Let $\eta$ denote the entropy per unit volume in $\Omega$ and $\vartheta$ is the absolute temperature. The entropy inequality is written globally in the current configuration as follows:

$$\frac{d}{dt}{\int }_{\Omega }\eta dv\ge {\int }_{\partial \Omega }\frac{-{\mathit{j}}^q·\mathit{n}}{\vartheta }da+{\int }_{\Omega }\frac{q}{\vartheta }dv$$

(13)

which can also be written in the reference configuration as follows:

$$\frac{d}{dt}{\int }_{{\Omega }_0}{\eta }_{R}dV\ge {\int }_{\partial {\Omega }_0}\frac{-{\mathit{j}}_R^q·\mathit{N}}{\vartheta }dA+{\int }_{{\Omega }_0}\frac{q_R}{\vartheta }dV$$

(14)

where ${\eta }_R$ is the entropy per unit volume in ${\Omega }_0$ with the relation ${\eta }_R=j\eta$. Using the divergence theorem, the local forms of the above entropy inequalities are as follows:

$$\dot{\eta }\ge -𝛻·\frac{{\mathit{j}}^q}{\vartheta }+\frac{q}{\vartheta }$$

(15)

$${\dot{\eta }}_R\ge -{𝛻}_R·\frac{{\mathit{j}}_R^q}{\vartheta }+\frac{q_R}{\vartheta }$$

(16)

Introducing the Helmholtz free energy density $\phi =\epsilon -\vartheta \eta$ and considering the energy balance (11), the inequality (15) becomes the following:

$$\mathit{\sigma }:𝛻\dot{\mathit{\chi }}+\mathit{e}·\dot{\mathit{d}}-\eta \dot{\vartheta }-\dot{\phi }-\frac{1}{\vartheta }{\mathit{j}}^q·𝛻\vartheta \ge 0$$

(17)

Similarly, introducing the Helmholtz free energy density ${\phi }_R={\epsilon }_R-\vartheta {\eta }_R$ with ${\phi }_R=j\phi$ and considering the energy balance (12), the inequality (16) becomes the following:

$$\mathit{P}:\dot{\mathit{F}}+\mathit{E}·\dot{\mathit{D}}-{\eta }_R\dot{\vartheta }-{\dot{\phi }}_R-\frac{1}{\vartheta }{\mathit{j}}_R^q·{𝛻}_R\vartheta \ge 0$$

(18)

which imposes the thermodynamic constraint on DEs. For convenience, we only employ the formulations given in the reference configuration in the following sections.

2.3. Constitutive Equations

In view of thermo-electro-viscoelastic effects of DEs, the Helmholtz free energy density ${\phi }_R$ can be assumed as a function of variables $\left(\mathit{C},{\mathit{C}}_i^e,\mathit{D},\vartheta \right)$, i.e.,

$${\phi }_R={\phi }_R\left(\mathit{C},{\mathit{C}}_i^e,\mathit{D},\vartheta \right)$$

(19)

so that its material time derivative can be written as follows:

$${\dot{\phi }}_R=\frac{\partial {\phi }_R}{\partial \mathit{C}}:\dot{\mathit{C}}+{\sum }_{i=1}^n\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}:{\dot{\mathit{C}}}_i^e+\frac{\partial {\phi }_R}{\partial \mathit{D}}:\dot{\mathit{D}}+\frac{\partial {\phi }_R}{\partial \vartheta }\dot{\vartheta }$$

(20)

The first and second terms on the right-hand side of the above equation can be rewritten as follows:

$$\frac{\partial {\phi }_R}{\partial \mathit{C}}:\dot{\mathit{C}}=2\mathit{F}·\frac{\partial {\phi }_R}{\partial \mathit{C}}:\dot{\mathit{F}}$$

$$\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}:{\dot{\mathit{C}}}_i^e=\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}:\frac{\partial {\mathit{C}}_i^e}{\partial \mathit{F}}:\dot{\mathit{F}}+\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}:\frac{\partial {\mathit{C}}_i^e}{\partial {\mathit{F}}_i^v}:{\dot{\mathit{F}}}_i^v=2{\mathbf{F}}_i^e·\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}·{\left({\mathit{F}}_i^{vT}\right)}^{-1}:\dot{\mathit{F}}-2{\mathit{C}}_i^e·\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}:{\mathit{D}}_i^v$$

(21)

with

$${\mathit{D}}_i^v=\frac{1}{2}\left({\mathit{L}}_i^v+{\mathit{L}}_i^{vT}\right), {\mathbf{L}}_i^v={\dot{\mathit{F}}}_i^v·{\left({\mathit{F}}_i^v\right)}^{-1}$$

(22)

Then, the substitution of Equations (20) and (21) into Equation (18) yields the following:

$$\begin{gathered} \left(\mathit{P}-2\mathit{F}·\frac{\partial {\phi }_R}{\partial \mathit{C}}-2{\sum }_{i=1}^n{\mathit{F}}_i^e·\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}·{\left({\mathit{F}}_i^{vT}\right)}^{-1}\right):\dot{\mathit{F}}-\left(\mathit{E}-\frac{\partial {\phi }_R}{\partial \mathit{D}}\right)·\dot{\mathit{D}}-\left({\eta }_R+\frac{\partial {\phi }_R}{\partial \vartheta }\right)\dot{\vartheta }- \\ \frac{1}{\vartheta }{\mathit{j}}_R^q·{𝛻}_R\vartheta +{\sum }_{i=1}^{n}2{\mathit{C}}_i^e·\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}:{\mathit{D}}_i^v\ge 0 \end{gathered}$$

(23)

In the case that $\mathit{P}$, $\mathit{E}$, and ${\eta }_R$ are independent of $\dot{\mathit{F}}$, $\dot{\mathit{D}}$, and $\dot{\vartheta }$, the first three terms of the above inequality should vanish so that we have the following constitutive relations:

$$\mathit{P}=2\mathit{F}·\frac{\partial {\phi }_R}{\partial \mathit{C}}+2{\sum }_{i=1}^n{\mathit{F}}_i^e·\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}·{\left({\mathit{F}}_i^{vT}\right)}^{-1}$$

(24)

$$\mathit{E}=\frac{\partial {\phi }_R}{\partial \mathit{D}}$$

(25)

$${\eta }_R=-\frac{\partial {\phi }_R}{\partial \vartheta }$$

(26)

Thus, the inequality (23) reduces to the following:

$$-\frac{1}{\vartheta }{\mathit{j}}_R^q·{𝛻}_R\vartheta +{\sum }_{i=1}^n{\mathit{M}}_i^{neq}:{\mathit{D}}_i^v\ge 0$$

(27)

where ${\mathit{M}}_i^{neq}$ is the non-equilibrium Mandel stress tensor [33], defined as follows:

$${\mathit{M}}_i^{neq}=2{\mathit{C}}_i^e·\frac{\partial {\phi }_R}{\partial {\mathit{C}}_i^e}$$

(28)

More specifically, to satisfy the thermodynamic constraint given by inequality (27), the following constitutive equations are deduced:

$${\mathit{j}}_R^q=-\mathit{Y}·{𝛻}_R\vartheta$$

(29)

$${\mathit{D}}_i^v={\mathit{Q}}_i:{\mathit{M}}_i^{neq}$$

(30)

where $\mathit{Y}$ is a second-order tensor and ${\mathit{Q}}_i$ a fourth-order tensor, both of which are positive-definite. Here, Equation (29) is the Fourier heat conduction law and Equation (30) is the rheological viscous flow rule [22,33].

2.4. Heat Conduction

Now, the first Gibbs relation can be obtained by substituting Equations (21), (24)–(26) and (28) into Equation (20), as follows:

$${\dot{\phi }}_R=\mathit{P}:\dot{\mathit{F}}+{\sum }_{i=1}^n{\mathit{M}}_i^{neq}:{\mathit{D}}_i^v+\mathit{E}·\dot{\mathit{D}}-{\eta }_R\dot{\vartheta }$$

(31)

which can also be simultaneously transformed into the second Gibbs relation by using the relation ${\phi }_R={\epsilon }_R-\vartheta {\eta }_R$, i.e.,

$${\dot{\epsilon }}_R=\mathit{P}:\dot{\mathit{F}}+{\sum }_{i=1}^n{\mathit{M}}_i^{neq}:{\mathit{D}}_i^v+\mathit{E}·\dot{\mathit{D}}+\vartheta {\dot{\eta }}_R$$

(32)

Using the energy balance Equation (12), Equation (32) yields the following:

$$\vartheta {\dot{\eta }}_R=-{𝛻}_R·{\mathit{j}}_R^q+q_R+{\sum }_{i=1}^n{\mathit{M}}_i^{neq}:{\mathit{D}}_i^v$$

(33)

which indicates that the total entropy change includes the heat exchange of DE with its environment and its internal dissipation of heat due to viscous effects. From another perspective, it can also be found that the total entropy change includes the reversible entropy exchange (i.e., the right-hand side of the inequality (16)) and the irreversible entropy production (i.e., the left-hand side of the inequality (27)). Then, taking the derivative of Equation (26) with respect to time and substituting the result into Equation (33) gives the following:

$$c\dot{\vartheta }=-{𝛻}_R·{\mathit{j}}_R^q+q_R+{\sum }_{i=1}^n{\mathit{M}}_i^{neq}:{\mathit{D}}_i^v+\vartheta \frac{\partial \mathit{P}}{\partial \vartheta }:\dot{\mathit{F}}-{\sum }_{i=1}^n\vartheta \frac{\partial {\mathit{M}}_i^{neq}}{\partial \vartheta }:{\dot{\mathit{D}}}_i^v+\vartheta \frac{\partial \mathit{E}}{\partial \vartheta }·\dot{\mathit{D}}$$

(34)

where $c$ is the specific heat capacity, defined as $c=-\frac{{\partial }^2{\phi }_R}{\partial {\vartheta }^2}$. This heat conduction equation explicitly indicates that temperature varies with heat exchange, viscous flow, the temperature-stress coupling effect, and the temperature-electricity coupling effect.

Up to now, we have obtained all the governing equations for the thermo-electro-viscoelastic process of DEs, including kinematic Equations (1), (2) and (22), balance Equations (5), (8) and (34), constitutive Equations (24)–(26) and (28)–(30), which can be solved under given boundary and initial conditions. In the next section, we will derive specific constitutive equations by giving specific Helmholtz free energy for DE and simulate viscoelastic phenomena of the DE polymer VHB4905 under various thermo-electro-mechanical loading conditions.

3. Special Cases

VHB4905 is one class of commercially available DEs, whose highly nonlinear stress-strain relations have been extensively investigated at various conditions [20,26,34]. In this section, specific constitutive equations for the coupling thermo-electro-viscoelastic behaviors of VHB4905 are first given and then used to simulate the loading-unloading curves under two different coupling conditions.

3.1. Specific Constitutive Equations

The Helmholtz free energy density is assumed to be additively decomposed, as follows:

$${\phi }_R={\phi }_R^{eq}+{\phi }_R^{neq}+{\phi }_R^E+{\phi }_R^T$$

(35)

where ${\phi }_R^{eq}$ and ${\phi }_R^{neq}$ are, respectively, the equilibrium and non-equilibrium Helmholtz free energy density resulting from the stretching of the ground network and subnetworks, ${\phi }_R^E$ is the Helmholtz free energy density due to the interaction of quasi-static electric charges in VHB4905, and ${\phi }_R^T$ is the Helmholtz free energy density due to temperature fluctuation.

Here, we adopt the Gent model [29] for ${\phi }_R^{eq}$ and ${\phi }_R^{neq}$ in consideration of the strain-stiffening effect that DEs may stiffen sharply when the network chain approaches its extension limit [35], as follows:

$${\phi }_R^{eq}=-\frac{G^{eq}L}{2}\mathit{ln}\left(\frac{L-tr\mathit{C}+3}{L}\right)$$

(36)

$${\phi }_R^{neq}=-{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i}{2}\mathit{ln}\left(\frac{L_i-tr{\mathit{C}}_i^e+3}{L_i}\right)$$

(37)

where the symbol ‘$tr$’ denotes the trace of a tensor, $G^{eq}$ and $G_i^{neq}$ denote, respectively, the equilibrium modulus of the elastic ground network and the nonequilibrium modulus of the ith viscous subnetwork, and $L$ and $L_i$ represent the extension limits of the elastic ground network and ith viscous subnetwork, respectively.

For the sake of simplicity, an isotropic formulation of electrostatics in the current configuration is employed [36]:

$$\frac{{\phi }_R^E}{j}=\frac{1}{2{ϵ}_0{ϵ}_r}\mathit{d}·\mathit{d}$$

(38)

whose corresponding form in the reference configuration can be obtained by using the relation $\mathit{D}=j\mathit{d}·{\left({\mathit{F}}^T\right)}^{-1}$, as follows:

$${\phi }_R^E=\frac{1}{2j{ϵ}_0{ϵ}_r}\left({\mathit{F}}^T·\mathit{F}\right):\left(\mathit{D}\otimes \mathit{D}\right)$$

(39)

where the symbol $\otimes$ denotes the dyadic product defining $\mathit{A}\otimes \mathit{B}$ as $A_m{B}_k$ for any vectors $\mathit{A}$ (with $A_m \mathrm{a}\mathrm{s} \mathrm{i}\mathrm{t}\mathrm{s} \mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{n}\mathrm{t}$) and $\mathit{B}$ (with $B_k \mathrm{a}\mathrm{s} \mathrm{i}\mathrm{t}\mathrm{s} \mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{n}\mathrm{t}$); and ${ϵ}_0$ and ${ϵ}_r$ are, respectively, the vacuum permittivity and the relative permittivity.

For the thermal part of the Helmholtz free energy density, we adopt the following [37]:

$${\phi }_R^T=c_0\left(\vartheta -{\vartheta }_0-\vartheta \mathit{ln}\frac{\vartheta }{{\vartheta }_0}\right)-c_v\frac{{\left(\vartheta -{\vartheta }_0\right)}^2}{2{\vartheta }_0}$$

(40)

which leads to a linear dependency of the specific heat capacity on temperature:

$$c=c_0+c_v\frac{\vartheta }{{\vartheta }_0}$$

(41)

where ${\vartheta }_0$ is a reference temperature, and $c_0$ and $c_v$ are two coefficients of this linear dependency.

Substituting Equations (35)–(38) into (24), (25) and (28), we have

$$\begin{gathered} \mathit{P}=\frac{G^{eq}L}{L-tr\mathit{C}+3}\mathit{F}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}{\mathit{F}}_i^e·{{(\mathit{F}}_i^{vT})}^{-1}+{\left({\mathit{F}}^T\right)}^{-1} \\ ·\left[\frac{1}{j{ϵ}_0{ϵ}_r}\left(\mathit{C}·\mathit{D}\right)\otimes \mathit{D}-\frac{1}{2}\left(\frac{1}{j{ϵ}_0{ϵ}_r}\mathit{C}·\mathit{D}·\mathit{D}\right)\mathit{I}\right] \end{gathered}$$

(42)

$$\mathit{E}=\frac{1}{j{ϵ}_0{ϵ}_r}\mathit{C}·\mathit{D}$$

(43)

$${\mathit{M}}_i^{neq}=\frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}{\mathit{C}}_i^e$$

(44)

with $\mathit{I}$ as the second-order unit tensor.

Next, both the elastic deformations of the ground network and parallel subnetworks are considered to be incompressible [25], for which Lagrange multipliers $\mathsf{\Pi }$ and ${\mathsf{\Pi }}_i$ are introduced to impose these constraint conditions and modify the free energy density as follows:

$${\tilde{\phi }}_R={\phi }_R-\mathsf{\Pi }\left(j-1\right)-{\sum }_{i=1}^n{\mathsf{\Pi }}_i\left(j_i^e-1\right)$$

(45)

where $j_i^e$ is the determinant of the elastic deformation gradient ${\mathit{F}}_i^e$. Replacing ${\phi }_R$ with ${\tilde{\phi }}_R$ in Equations (24) and (28), we can rewrite $\mathit{P}$ and ${\mathit{M}}_i^{neq}$ as follows:

$$\begin{gathered} \mathit{P}=\frac{G^{eq}L}{L-tr\mathit{C}+3}\mathit{F}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}{\mathit{F}}_i^e·{{(\mathit{F}}_i^{vT})}^{-1}+{\left({\mathit{F}}^T\right)}^{-1}·\left[\mathit{E}\otimes \mathit{D}-\frac{1}{2}\left(\mathit{E}·\mathit{D}\right)\mathit{I}\right]- \\ j\mathsf{\Pi }{\left({\mathit{F}}^T\right)}^{-1}-{\sum }_{i=1}^n{j}_i^e{\mathsf{\Pi }}_i{\left({\mathit{F}}^T\right)}^{-1} \end{gathered}$$

(46)

$${\mathit{M}}_i^{neq}=\frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}{\mathit{C}}_i^e-j_i^e{\mathsf{\Pi }}_i\mathit{I}$$

(47)

where we have used the remaining unchanged relation (43). Furthermore, the Cauchy stress tensor $\mathit{\sigma }$ can also be obtained by using the relations $\mathit{\sigma }=\frac{1}{j}\mathit{P}·{\mathit{F}}^T$, $\mathit{d}={ϵ}_0{ϵ}_r\mathit{e}$, $\mathit{D}=j{\left(\mathbf{F}\right)}^{-1}·\mathit{d}$, $\mathit{E}={\mathit{F}}^T·\mathit{e}$, and $j=j_i^e=1$ in Equation (46), as follows:

$$\begin{gathered} \mathit{\sigma }=\frac{G^{eq}L}{L-tr\mathit{C}+3}\mathit{F}·{\mathit{F}}^T+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}{\mathit{F}}_i^e·{\mathit{F}}_i^{eT}+{ϵ}_0{ϵ}_r\left[\mathit{e}\otimes \mathit{e}-\frac{1}{2}\left(\mathit{e}·\mathit{e}\right)\mathit{I}\right] \\ -\left(\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i\right)\mathit{I} \end{gathered}$$

(48)

Employing $j=j_i^e=1$ and the multiplicative decomposition of $\mathit{F}$ in Equation (2), the viscous incompressibility, $tr{\mathit{D}}_i^v=\mathbf{0}$ or $j_i^v=1$ with $j_i^v$ as the determinant of ${\mathit{F}}_i^v$, can be deduced, from which the fourth-order tensor ${\mathit{Q}}_i$ in Equation (30) can be taken as follows [25]:

$${\mathit{Q}}_i=\frac{1}{2v_i}\left({\mathit{I}}_4-\frac{1}{3}\mathit{I}\otimes \mathit{I}\right)$$

(49)

where ${\mathit{I}}_4$ is the fourth-order unit tensor and $v_i$ the viscosity of the $i\mathrm{t}\mathrm{h}$ subnetwork. Then, the relaxation time for the $i\mathrm{t}\mathrm{h}$ subnetwork can be defined as ${\tau }_i=\frac{v_i}{G_i^{neq}}$ [31].

3.2. Thermo-Viscoelastic Coupling

In the thermo-viscoelastic experiments of Mehnert et al. [21], VHB4905 samples with the dimensions $130 \mathrm{m}\mathrm{m}\times 10 \mathrm{m}\mathrm{m}\times 0.5 \mathrm{m}\mathrm{m}$ are first heated at different temperatures in the thermal chamber for about 10 min, which has proven to be sufficient for the samples to induce the temperature of the chamber, and then mounted onto the machine Inspekt 5 kN for multi-step relaxation tests and cyclic loading-unloading tests.

During heating, we assume that the deformation of the samples can be neglected due to the relatively small thermal expansion coefficient and the free boundary conditions. Substituting Equations (41) and (29) into Equation (34) and neglecting all the heat sources, the heat conduction equation reduces to the following:

$$\left(c_0+c_v\frac{\vartheta }{{\vartheta }_0}\right)\dot{\vartheta }={𝛻}_R·\left(\mathit{Y}·{𝛻}_R\vartheta \right)$$

(50)

where we consider the isotropic Fourier heat conduction so that $\mathit{Y}=\left({\kappa }_0+{\kappa }_v\frac{\vartheta }{{\vartheta }_0}\right)\mathit{I}$ with ${\kappa }_0$ and ${\kappa }_v$ are two coefficients of the linear temperature-dependent conductivity [38]. The Neumann boundary condition is used to describe the convective heat transfer between the sample and the air in the thermal chamber, given as follows:

$${\mathit{j}}_R·\mathit{n}=h\left({\vartheta }_s-{\vartheta }_t\right)$$

(51)

where $h$ is the convective heat transfer coefficient, ${\vartheta }_s$ and ${\vartheta }_t$ are, respectively, the temperature of the sample surface and the temperature of air in the thermal chamber. Here, the room temperature ($296 \mathrm{K}$) is taken as the reference temperature (i.e., ${\vartheta }_0=296 \mathrm{K}$) and the temperature of air in the thermal chamber is set as ${\vartheta }_t=$ $353 \mathrm{K}$. The initial condition for the temperature of the sample is assumed to be as follows:

$${\left.\vartheta \right|}_{t=0}=296 \mathrm{K}$$

(52)

The material parameters, $c_0,c_v,{\kappa }_0,$ and ${\kappa }_v$, used in the calculation are given in

Table 1. Material parameters.

Parameters	Values
$c_0\left(\left(\mathrm{J}/\left({\mathrm{m}}^3·\mathrm{K}\right)\right)\right)$	$8.86\times {10}^5$
$c_v\left(\mathrm{J}/\left({\mathrm{m}}^3·\mathrm{K}\right)\right)$	$9.18\times {10}^5$
${\kappa }_0\left(\mathrm{W}/\left(\mathrm{m}·\mathrm{K}\right)\right)$	$0.42$
${\kappa }_v\left(\mathrm{W}/\left(\mathrm{m}·\mathrm{K}\right)\right)$	$-0.17$
$h\left(\mathrm{W}/\left({\mathrm{m}}^2·\mathrm{K}\right)\right)$	$20$

, based on Dippel et al. [38] for natural rubber, which resembles VHB4905 in molecular structure. Furthermore, considering that the convective heat transfer coefficient for forced gas convection is about $20~300 \mathrm{W}/\left({\mathrm{m}}^2·\mathrm{K}\right)$, $h=20 \mathrm{W}/\left({\mathrm{m}}^2·\mathrm{K}\right)$ is chosen, as given in Table 1, for the convective heat transfer between the sample and the air in the thermal chamber so that we can roughly estimate the maximum time for reaching the steady state of heat conduction.

The heat conduction inside VHB4905 samples can be simplified as a one-dimensional problem since the dimension in the thickness direction is much smaller than the other two lateral dimensions, which results in negligible convective heat transfer on the lateral surfaces. By solving Equation (50) with the boundary condition (51) and the initial condition (52) through the commercial software COMSOL Multiphysic 5.3, the temperature distribution along the thickness direction during heating is obtained and shown in Figure 2. It can be seen that the temperature along thickness direction is nearly uniform and increases with time until it reaches $353 \mathrm{K}$ (the temperature in the thermal chamber). This is because the great heat conductivity of VHB4905 leads to rapid heat transfer across a small distance between neighboring surfaces (0.5 mm), and forced gas convection between the sample and the air in the thermal chamber results in a steady temperature distribution after about 100 s. The numerical results reveal that 10 min is sufficient for VHB4905 samples to induce the temperature of the thermal chamber.

Next, the corresponding deformation gradient in the mechanical tests can be written as follows:

$$\begin{array}{c}\mathit{F}=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right], {\mathit{F}}_i^e=\left[\begin{array}{ccc}{\lambda }_{i1}^e& 0& 0\\ 0& {\lambda }_{i2}^e& 0\\ 0& 0& {\lambda }_{i3}^e\end{array}\right], {\mathit{F}}_i^v=\left[\begin{array}{ccc}{\lambda }_{i1}^v& 0& 0\\ 0& {\lambda }_{i2}^v& 0\\ 0& 0& {\lambda }_{i3}^v\end{array}\right]\hfill \\ {\mathit{D}}_i^v=\left[\begin{array}{ccc}{\dot{\lambda }}_{i1}^v/{\lambda }_{i1}^v& 0& 0\\ 0& {\dot{\lambda }}_{i2}^v/{\lambda }_{i2}^v& 0\\ 0& 0& {\dot{\lambda }}_{i3}^v/{\lambda }_{i3}^v\end{array}\right]\hfill \end{array}$$

(53)

Here ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are the principle stretches of the deformation gradient $\mathit{F}$; ${\lambda }_{i1}^e$, ${\lambda }_{i2}^e$, and ${\lambda }_{i3}^e$ are the principle stretches of the elastic deformation gradient ${\mathit{F}}_i^e$; and ${\lambda }_{i1}^v$, ${\lambda }_{i2}^v$, and ${\lambda }_{i3}^v$ are the principle stretches of the viscous deformation gradient ${\mathit{F}}_i^v$. In consideration of the equal lateral stretches during the uniaxial tensile test and the incompressible deformation of VHB4905 samples ($j=j_i^e=j_i^v=1$), we have the following:

$${\lambda }_1=\lambda , {\lambda }_2={\lambda }_3={\left(\lambda \right)}^{-\frac{1}{2}}$$

$${\lambda }_1={\lambda }_{i1}^e{\lambda }_{i1}^v=\lambda , {\lambda }_{i2}^e{\lambda }_{i2}^v={\lambda }_{i3}^e{\lambda }_{i3}^v={\left(\lambda \right)}^{-\frac{1}{2}}$$

$${\lambda }_{i1}^e={\lambda }_i^e, {\lambda }_{i2}^e={\lambda }_{i3}^e={\left({\lambda }_i^e\right)}^{-\frac{1}{2}},{\lambda }_{i1}^v={\lambda }_i^v,{\lambda }_{i2}^v={\lambda }_{i3}^v={\left({\lambda }_i^v\right)}^{-\frac{1}{2}}$$

(54)

with $\lambda$, ${\lambda }_i^e$, and ${\lambda }_i^v$, respectively, as the principle stretches of the deformation gradients $\mathit{F}$, ${\mathit{F}}_i^e$, and ${\mathit{F}}_i^v$ along the tensile direction.

Let $P_1$, $P_2$, and $P_3$ denote the components of $\mathit{P}$ along three principal directions, respectively. Substituting Equations (53) and (54) into Equation (46) and then discarding the electric terms, we have the following:

$$P_1=\frac{G^{eq}L\lambda }{L-2{\lambda }^{-1}-{\lambda }^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda \left({\lambda }_i^v\right)}^{-2}}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}-\frac{\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i}{\lambda }$$

$$P_2=P_3=\frac{G^{eq}L{\lambda }^{-\frac{1}{2}}}{L-2{\lambda }^{-1}-{\lambda }^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda }^{-\frac{1}{2}}{\lambda }_i^v}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}-\left(\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i\right){\lambda }^{\frac{1}{2}}$$

(55)

According to the boundary condition, $P_2=P_3=0$ for the uniaxial loading-unloading tests, we can further obtain the following:

$$\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i=\frac{G^{eq}L{\lambda }^{-1}}{L-2{\lambda }^{-1}-{\lambda }^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda }^{-1}{\lambda }_i^v}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}$$

(56)

and

$$P_1=\frac{G^{eq}L(\lambda -{\lambda }^{-2})}{L-2{\lambda }^{-1}-{\lambda }^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i\left[{\lambda \left({\lambda }_i^v\right)}^{-2}-{\lambda }^{-2}{\lambda }_i^v\right]}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}$$

(57)

Here, the incompressibilities of the ground network and the parallel subnetworks have an effect on the stress, and Lagrange multipliers are eliminated according to the boundary conditions.

Let $M_{i1}^{neq}$, $M_{i2}^{neq}$, and $M_{i3}^{neq}$ denote the components of ${\mathit{M}}_i^{neq}$ along three principal directions, respectively. Similarly, substituting Equations (53) and (54) into Equation (47), we have the following:

$$\begin{array}{c}{M}_{i1}^{neq}=\frac{G_i^{neq}{L}_i{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}-{\mathsf{\Pi }}_i\hfill \\ M_{i2}^{neq}=M_{i3}^{neq}=\frac{G_i^{neq}{L}_i{\lambda }^{-1}{\lambda }_i^v}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}-{\mathsf{\Pi }}_i\hfill \end{array}$$

(58)

Then, substitution of Equations (58), (53) and (49) into Equation (30) yields the following:

$$\frac{{\dot{\lambda }}_i^v}{{\lambda }_i^v}=\frac{L_i}{3{\tau }_i\left[L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3\right]}\left[{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}-{\lambda }^{-1}{\lambda }_i^v\right]$$

(59)

which describes the viscous flow in VHB4905 subjected to uniaxial tension. It is worth noting that the elastic incompressibility of the subnetworks has no influence on the viscous flow since only the deviatoric stress, excluding the hydrostatic pressure ${\mathsf{\Pi }}_i$, promotes the viscous flow.

Under the initial conditions, ${\left.P_1\right|}_{t=0}=0$, ${\left.\lambda \right|}_{t=0}={\left.{\lambda }_i^v\right|}_{t=0}=1$, the coupled Equations (57) and (59) can be simultaneously solved to predict the mechanical behaviors under the given material parameters. To determine the material parameters, the least square method is used to fit Equations (57) and (59) with the experimental data of Mehnert et al. [21], who performed multi-step relaxation tests for identifying the elastic responses at different stretches and cyclic loading-unloading tests for acquiring the viscous information at different stretching rates. In each step of a multi-step relaxation test with eight consecutive steps, the sample is stretched quickly with an increased 25% deformation and then held at this deformation state for thirty minutes (the viscous stress is considered to be relaxed almost completely after this period), for which Equations (57) and (59) can be reduced to the equilibrium equation:

$$P_1=\frac{G^{eq}L(\lambda -{\lambda }^{-2})}{L-2{\lambda }^{-1}-{\lambda }^2+3}$$

(60)

Here, we take $L=155$, according to Kollosche et al. [39], and fit Equation (60) with the multi-step relaxation tests [21] to obtain $G^{eq}$ = 15.12 kPa, 13.05 kPa, and 11.84 kPa for the temperatures 296 K, 313 K, and 333 K, respectively. Figure 3 gives a comparison of the equilibrium stress-stretch relations, respectively, from model fitting and experimental data. The experimental data are in good agreement with the fitting curves, indicating that Equation (60) is capable of predicting the equilibrium response of VHB4905. Three discrete points (296, 15.12), (313, 13.05), and (333, 11.84) are further used to fit, as shown in Figure 4, the following exponential function:

$$G^{eq}\left(\vartheta \right)=aexp\left(\frac{b\left(\vartheta -{\vartheta }_g\right)}{\vartheta -{\vartheta }_g-c}\right)$$

(61)

with $a=0.13 \mathrm{k}\mathrm{P}\mathrm{a},b=4.16, \mathrm{a}\mathrm{n}\mathrm{d} c=8.06 \mathrm{K}$. Here, ${\vartheta }_g$ ($233 \mathrm{K}$) is the glass transition temperature of VHB4905. This dependence of the equilibrium modulus on temperature is of the Williams–Landel–Ferry (WLF) type [40], and indicates that increasing temperature would soften VHB4905. This function can be used to calculate $G^{eq}=25.56 \mathrm{k}\mathrm{P}\mathrm{a}$ for 273 K.

Figure 5 gives a comparison of the experimental data with the model fitting and the model prediction at different stretching rates and a fixed temperature of 273 K. We have found that $n=3$ results in well-converged results as well as physically realistic parameter values. With $G^{eq}=25.56 \mathrm{k}\mathrm{P}\mathrm{a}$ and $L=L_i=155$ ($i=1, 2, 3$) determined, the non-equilibrium modulus $G_i^{neq} (38.44 \mathrm{k}\mathrm{P}\mathrm{a}, 94.23 \mathrm{k}\mathrm{P}\mathrm{a}, 94.96 \mathrm{k}\mathrm{P}\mathrm{a})$ and the relaxation time ${\tau }_i (413.32 \mathrm{s}, 5.43 \mathrm{s}, 1.65 \mathrm{s})$ are obtained by fitting Equations (57) and (59) with the experimental data at stretching rates $\left|\dot{\lambda }\right|=0.03/\mathrm{s}$ and $\left|\dot{\lambda }\right|=0.05/\mathrm{s}$ in Figure 5. Then, the model prediction of $\left|\dot{\lambda }\right|=0.1/\mathrm{s}$ by using all determined parameters agrees well with the experimental data, and each loading-unloading cycle shows typical viscoelastic behavior, which validates the employed viscoelastic model.

Figure 6 also presents a comparison of the experimental data with the model fitting and the model prediction at different stretching rates and temperatures. Here, we assume that the relaxation times and the extension limits are temperature-independent, so that only the dependence of the non-equilibrium moduli on temperature is obtained via the model fitting. The experimental data at the stretching rates $\left|\dot{\lambda }\right|=0.025/\mathrm{s}$ and $\left|\dot{\lambda }\right|=0.05/\mathrm{s}$ are used to determine the non-equilibrium moduli, and $\left|\dot{\lambda }\right|=0.1/\mathrm{s}$ is used to validate the model fitting. It can be seen that the results from the model fitting and the model prediction are close to the experimental data, and the viscoelastic effect at the higher temperature is more significant, indicating that the strategy is reasonable. We assume that the dependence of the non-equilibrium moduli on temperature is also a WLF type like Equation (61) and fit this type with the obtained four discrete points, yielding the following fitting functions:

$$\begin{array}{c}{G}_1^{neq}\left(\vartheta \right)=5.54\times 10^{-2}\mathit{exp}\left(\frac{4.74\left(\vartheta -233\right)}{\vartheta -233-11.04}\right)\hfill \\ G_2^{neq}\left(\vartheta \right)=6.79\times 10^{-3}\mathit{exp}\left(\frac{6.87\left(\vartheta -233\right)}{\vartheta -233-11.17}\right)\hfill \\ G_3^{neq}\left(\vartheta \right)=7.89\times 10^{-3}\mathit{exp}\left(\frac{7.21\left(\vartheta -233\right)}{\vartheta -233-9.34}\right)\hfill \end{array}$$

(62)

The fitting curves of the non-equilibrium modulus versus temperature and the corresponding fitted points are depicted in Figure 7, where the fitting curves match roughly with these points. With these assumptions and relations, all parameters at 353 K (see

Table 2. Material parameters at different temperatures.

	273 K	296 K	313 K	333 K	353 K
$G^{eq}\left(\mathrm{k}\mathrm{P}\mathrm{a}\right)$	23.44	15.12	13.05	11.82	11.07
$G_1^{neq}\left(\mathrm{k}\mathrm{P}\mathrm{a}\right)$	38.44	16.83	13.72	11.55	10.20
$G_2^{neq}\left(\mathrm{k}\mathrm{P}\mathrm{a}\right)$	94.23	27.39	23.90	13.11	13.29
$G_3^{neq}\left(\mathrm{k}\mathrm{P}\mathrm{a}\right)$	94.96	41.04	20.57	25.41	19.49
${\tau }_1\left(\mathrm{s}\right)$	413.32	413.32	413.32	413.32	413.32
${\tau }_2\left(\mathrm{s}\right)$	5.43	5.43	5.43	5.43	5.43
${\tau }_3\left(\mathrm{s}\right)$	1.65	1.65	1.65	1.65	1.65
$L,L_1,L_2,L_3$	155	155	155	155	155

) are obtained and substituted into Equations (57) and (59) for predicting the loading-unloading curves at different stretching rates.

Figure 8 shows a comparison of the loading-unloading curves from the model prediction and the experimental data at $\left|\dot{\lambda }\right|=0.025, 0.05,\frac{0.1}{\mathrm{s}},$ and 353 K. It can be seen that the predicted loading-unloading curves from the present model match reasonably well with the experimental results, and the peak stress of each curve increases with the stretching rate. The comparison indicates that the present model and its fitting parameters are able to describe the viscoelastic effect of VHB4905 at different stretching rates and temperatures.

3.3. Thermo-Electro-Viscoelastic Coupling

In the thermo-electro-viscoelastic experiments of VHB4905 by Mehnert et al. [21], material samples with the dimensions $70 \mathrm{m}\mathrm{m}\times 100 \mathrm{m}\mathrm{m}\times 0.5 \mathrm{m}\mathrm{m}$, coated with carbon conductive grease, acting as compliant electrodes, on its two largest surfaces, are first heated for fifteen minutes to reach a target temperature in the closed thermal chamber. Then, an electric field is applied on the two largest surfaces of the heated samples. Finally, the loading-unloading tests on the samples with different electric fields are conducted.

Let $E$ denote the electric field along thickness direction reckoned in the reference configuration. Substitution of Equation (53) into Equation (46) yields the following:

$$P_1=\frac{G^{eq}L{\lambda }_1}{L-{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda }_1{\left({\lambda }_{i1}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-\frac{\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i}{{\lambda }_1}-{ϵ}_0{ϵ}_r\frac{E^2}{2{\lambda }_1{\lambda }_3^2}$$

$$P_2=\frac{G^{eq}L{\lambda }_2}{L-{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda }_2{\left({\lambda }_{i2}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-\frac{\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i}{{\lambda }_2}-{ϵ}_0{ϵ}_r\frac{E^2}{2{\lambda }_2{\lambda }_3^2}$$

$$P_3=\frac{G^{eq}L{\lambda }_3}{L-{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda }_3{\left({\lambda }_{i3}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-\frac{\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i}{{\lambda }_3}+{ϵ}_0{ϵ}_r\frac{E^2}{2{\lambda }_3^2}$$

(63)

From the boundary condition $P_3=0$, we can obtain the following:

$$\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i=\frac{G^{eq}L{\lambda }_3^2}{L-{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-{{ϵ}_0ϵ}_r\frac{E^2}{2{\lambda }_3^2}$$

(64)

Further considering the incompressibility constraints ${\lambda }_1{\lambda }_2{\lambda }_3={\lambda }_{i1}^v{\lambda }_{i2}^v{\lambda }_{i3}^v=1,$ and the boundary conditions $P_2=0$, we obtain the following:

$$\begin{gathered} P_1=\frac{G^{eq}L\left({\lambda }_1-{\lambda }_1^{-3}{\lambda }_2^{-2}\right)}{L-{\lambda }_1^2-{\lambda }_2^2-{\lambda }_1^{-2}{\lambda }_2^{-2}+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i\left[{\lambda }_1{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_1^{-3}{\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2\right]}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_1^{-2}{\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2+3} \\ -{ϵ}_0{ϵ}_r{{\lambda }_2^2{\lambda }_{1}E}^2 \end{gathered}$$

$$P_2=\frac{G^{eq}L\left({\lambda }_2-{\lambda }_1^{-2}{\lambda }_2^{-3}\right)}{L-{\lambda }_1^2-{\lambda }_2^2-{\lambda }_1^{-2}{\lambda }_2^{-2}+3}+{\sum }_{i=1}^n\frac{G_i^{neq}{L}_i\left[{\lambda }_2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_1^{-2}{\lambda }_2^{-3}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2\right]}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_1^{-2}{\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2+3}-{ϵ}_0{ϵ}_r{\lambda }_2{\lambda }_1^2{E}^2=0$$

(65)

Next, substituting Equation (53) into Equation (47), we have the following:

$$M_{i1}^{neq}=\frac{G_i^{neq}{L}_i{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-{\mathsf{\Pi }}_i$$

$$M_{i2}^{neq}=\frac{G_i^{neq}{L}_i{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-{\mathsf{\Pi }}_i$$

$$M_{i3}^{neq}=\frac{G_i^{neq}{L}_i{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}}{L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_3^2{\left({\lambda }_{i3}^v\right)}^{-2}+3}-{\mathsf{\Pi }}_i$$

(66)

Substituting Equations (49), (53) and (66) into Equation (30) and using the incompressibility constraints ${\lambda }_1{\lambda }_2{\lambda }_3={\lambda }_{i1}^v{\lambda }_{i2}^v{\lambda }_{i3}^v=1$, we can obtain the viscous flow rule:

$$\begin{gathered} \frac{{\dot{\lambda }}_{i1}^v}{{\lambda }_{i1}^v}=\frac{L_i}{3{\tau }_i\left[L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_1^{-2}{\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2+3\right]}\left[{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2} \\ -\frac{{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}+{{\lambda }_1^{-2}\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2}{2}\right] \end{gathered}$$

$$\begin{gathered} \frac{{\dot{\lambda }}_{i2}^v}{{\lambda }_{i2}^v}=\frac{L_i}{3{\tau }_i\left[L_i-{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}-{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2}-{\lambda }_1^{-2}{\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2+3\right]}\left[{\lambda }_2^2{\left({\lambda }_{i2}^v\right)}^{-2} \\ -\frac{{\lambda }_1^2{\left({\lambda }_{i1}^v\right)}^{-2}+{{\lambda }_1^{-2}\lambda }_2^{-2}{\left({\lambda }_{i1}^v{\lambda }_{i2}^v\right)}^2}{2}\right] \end{gathered}$$

(67)

Here, it is worth noting that the hydrostatic pressure ${\mathsf{\Pi }}_i$ is eliminated by the mathematical operation since only the deviatoric stress is considered to promote the viscous flow in Equation (49).

Furthermore, an electric field prior to the mechanical test leads to an extension of the sample in the transverse direction, so we have ${\lambda }_3={\lambda }_{pre}$ and ${\lambda }_1={\lambda }_2=\frac{1}{\sqrt{{\lambda }_{pre}}}$ with ${\lambda }_{pre}$ being the free stretch in the thickness direction after the application of an electric field. The electric field in the sample is considered to reach a steady and uniform state quickly once the electric potential difference between the two largest surfaces is applied. Thus, the electric field can be calculated as $E=\frac{∆\varphi }{l}$, where $∆\varphi$ and $l$ are respectively the electric potential difference and the thickness of the VHB4905 sample. Because the relaxation times of the 1st, 2nd, and 3rd subnetworks are, respectively, 413.32 s, 5.43 s, and 1.65 s, we assume that the 1st subnetwork is completely elastic and the 2nd and 3rd subnetworks are completely viscous prior to the mechanical test. Therefore, using $P_1=0$ in the case of free expansion, Equation (65) can be reduced to the following:

$$\frac{\left(G^{eq}+G_1^{neq}\right)L\left(1/{\lambda }_{pre}-{\lambda }_{pre}^2\right)}{L-2/{\lambda }_{pre}-{\lambda }_{pre}^2+3}=0$$

(68)

which can be solved to obtain ${\lambda }_{pre}$. When the sample starts to be stretched, the initial conditions are given by the following:

$${\lambda }_1={\lambda }_2={\lambda }_{21}^v={\lambda }_{22}^v={\lambda }_{31}^v={\lambda }_{32}^v=\frac{1}{\sqrt{{\lambda }_{pre}}}, {\lambda }_{11}^v={\lambda }_{12}^v=1$$

(69)

A finite difference approach is employed to obtain the numerical results of the evolving stress $P_1$ with the stretch ${\lambda }_1$ by solving the coupled Equations (65) and (67) under the initial conditions (69).

Figure 9 gives a comparison of the evolving stress $P_1$ with the stretch ${\lambda }_1$ at $\left|{\dot{\lambda }}_1\right|=0.2/\mathrm{s}$ from the model prediction and the experimental data at different temperatures and electric fields. The equilibrium and non-equilibrium moduli for $\vartheta =296 \mathrm{K}$ are listed in Table 2 and for $\vartheta =333 \mathrm{K}$ are calculated via the relations (61) and (62) as $G^{eq}=12.36 \mathrm{k}\mathrm{P}\mathrm{a}$, $G_1^{neq}=12.24 \mathrm{k}\mathrm{P}\mathrm{a}$, $G_2^{neq}=17.38 \mathrm{k}\mathrm{P}\mathrm{a}$, and $G_3^{neq}=20.43 \mathrm{k}\mathrm{P}\mathrm{a}$. Here, we assume that the material parameters would not change with the electric field in consideration of the relatively small effect of the electric field on the viscous response in the experiments. It can be seen from Figure 9a that the loading-unloading curve with an applied electric potential difference of 6 kV (or an electric field $12\times {10}^6 \mathrm{V}/\mathrm{m}$) does not start exactly at the stretch of 1 but slightly above 1 (the pre-stretch of ${\lambda }_1$ is $1.0383$), which captures well the electric field-induced deformation prior to the mechanical deformation in the experiment. The curve without the electric potential is slightly below the experimental data because we have neglected the slight influence of the compliant electrodes on the viscoelastic response investigated by Mehnert and Steinmann [34]. Similarly, in Figure 9b, the pre-stretch of ${\lambda }_1$ for 323 K is $1.0540,$ and the distinction of the loading-unloading curves with and without an applied electric potential difference is relatively significant. Furthermore, it can be observed from Figure 9a,b that temperature fluctuation has a significantly more pronounced effect on the viscoelastic response than electric field, and the influence of electric field at high temperatures becomes clearly visible. The predicted loading-unloading curves match well with the experimental data, indicating that the present model is able to describe the thermo-electro-viscoelastic response of VHB4905 and has great potential in engineering applications of DEs.

4. Conclusions

In this paper, a thermodynamically consistent model at finite deformation for the thermo-electro-viscoelastic coupling behaviors of DEs is proposed and verified by the comparison of the nonlinear loading-unloading curves from the model prediction and the experimental data of VHB4905. The major novelty of the present work lies in the following aspects: First, we have proposed a thermo-electro-viscoelastic model for DEs at finite deformation, which can describe the nonlinear response of DEs at various stretching rates, temperatures, and electric fields. Especially, the elastic deformation incompressibility of both the ground network and a few parallel subnetworks is considered by introducing Lagrange multipliers, for which the distinct effects of the incompressibility on the force equilibrium and the viscous flow are clearly presented. The incompressibility conditions directly change the magnitude of the stress but have no influence on the viscous flow because only the deviatoric stress, excluding the hydrostatic pressure, promotes the viscous flow. Second, the WLF-like dependence of the moduli on temperature is found by fitting the model with the experimental data, contributing to the calculation of the moduli at different temperatures for theoretical modeling. Third, we also simulate the instant pre-stretch due to the interaction of quasi-static charges after the application of an electric field, which has not been modeled in the previous work.

The numerical results reveal that increasing temperature can soften DEs significantly, and the influence of the electric field on the mechanical response of DEs with respect to temperature is slight. Thus, this work provides a guide on modeling the thermo-electro-mechanical coupling behaviors of DEs and can help analyze the failure of DEs in functionalities and structures. Simultaneously, this theoretical model can be applied to different materials, such as piezoelectric composites, which will be our future study subject.