Viscoelastic Effects on the Response of Electroelastic Materials

Abstract

Electroelastic materials, as for example, 3M VHB 4910, are attracting attention as actuators or generators in some developments and applications. This is due to their capacity of being deformed when submitted to an electric field. Some models of their actuation are available, but recently, viscoelastic models have been proposed to give an account of the dissipative behaviour of these materials. Their response to an external mechanical or electrical force field implies a relaxation process towards a new state of thermodynamic equilibrium, which can be described by a relaxation time. However, it is well known that viscoelastic and dielectric materials, as for example, polymers, exhibit a distribution of relaxation times instead of a single relaxation time. In the present approach, a continuous distribution of relaxation times is proposed via the introduction of fractional derivatives of the stress and strain, which gives a better account of the material behaviour. The application of fractional derivatives is described and a comparison with former results is made. Then, a double generalisation is carried out: the first one is referred to the viscoelastic or dielectric models and is addressed to obtain a nonsymmetric spectrum of relaxation times, and the second one is the adoption of the more realistic Mooney–Rivlin equation for the stress–strain relationship of the elastomeric material. A modified Mooney–Rivlin model for the free energy density of a hyperelastic material, VHB 4910 has been used based on experimental results of previous authors. This last proposal ensures the appearance of the bifurcation phenomena which is analysed for equibiaxial dead loads; time-dependent bifurcation phenomena are predicted by the extended Mooney–Rivlin equations.

1. Introduction

Recently, viscoelastic models have been proposed to give an account of the dissipative behaviour of the electroelastic materials [1,2,3,4,5,6,7,8,9,10,11,12]. This is justified because the range of temperatures at which these materials are being currently used is close to the glass transition temperature. Indeed, above the glass transition temperature, the materials in a rubber-like state still exhibit some rest of the viscoelastic behaviour. As a consequence, in these materials, the response to an external mechanical or electrical force field implies relaxation processes towards a new state of thermodynamic equilibrium. The characteristic time required for these relaxation processes is called the relaxation time. The literature on the subject concerning mechanical and dielectric behaviour is extensive (see, for example, [13,14,15,16,17,18]). The more typical mechanical experiments are called creep or stress relaxation according to the input type (stressing or stretching). The same occurs in the case of the dielectric relaxation or retardation experiments.

In order to represent the most elementary relaxation or retardation processes, different types of rheological or dielectric models have been proposed. These models consist of a short number of springs and dashpots for the mechanical case or capacitors and resistances, for the dielectric counterpart. These models are appealing in the sense that they allow easy visualisation of the corresponding mechanical or electrical relaxation processes. However, they are only a crude representation of the actual behaviour of the system under study. In fact, most of the proposed models are related to systems governed by a single relaxation time. This problem may be solved [13,16] by adding, in series or parallel, depending on the problem, a finite set of these viscoelastic or dielectric models. This is equivalent to assuming a discrete distribution of relaxation times. However, instead of a discrete distribution of relaxation times, it seems more realistic to consider a continuous distribution or spectra of these times. In fact, it is well known that, in practice, viscoelastic and dielectric materials, as is the case of polymers, exhibit a distribution of relaxation times instead of a discrete set of relaxation times [13,15].

The molecular origin for this distribution may be due to the different microscopic moieties involved in the process, or else to a distribution of the shape and size of the defects associated with the neighbourhood surrounding the molecular moieties. The macroscopic result is, in any case, a broad relaxation spectrum of these relaxation times. For this reason, in the present approach, a continuous distribution is considered.

On the other hand, in the theoretical background, a substantial modification of the classical model equations for the purely electroelastic materials is necessary in order to consider the viscoelastic and dielectric relaxation effects both of which are intrinsically dissipative.. In this sense, an approach based on classical irreversible thermodynamics has been proposed by Suo et al. [4]. Indeed, these authors assume that the response of an elastomer to mechanical and electrical force fields is delayed in time due to the relaxation mechanisms associated with the motion of the structural moieties. The basic idea is to introduce into the free energy density convenient additional parameters, taking the form of internal variables, as in the classical irreversible thermodynamics (CIT), describing the different degrees of freedom associated with the dissipative relaxation processes. In the same way, the analysis of the electromechanical instability requires taking into account, together with the kinematic variables, the new internal ones just introduced. The basic lines of this approach have been followed by other authors [6,8]. The theory has been tested by assuming a rheological model governed by a single relaxation time in the framework of a neo-Hookean model. In particular, the pertinent viscoelastic data for the VHB 4910 are obtained from earlier papers by Wissler and Mazza [19] and Planté and Dubowsky [20]. In the approach used by Suo et al. [4], only mechanical dissipation has been taken into account. This is justified on the basis that the mechanical relaxation times are several orders of magnitude larger than the dielectric counterpart retardation times and, for this reason, the viscoelastic relaxation process takes place at several orders of magnitude larger than the dielectric ones. Consequently, the highest dielectric relaxation in the time domain, which is the space charge relaxation, has been completely relaxed. In the present approach, the viscoelastic data are obtained, as explained later, from the paper by Wang et al. [7], and incompressible materials are assumed. In any case, it should be mentioned that the dielectric relaxation processes can be addressed in the same way as the counterpart viscoelastic ones.

Recently, Ghosh and Lopez-Pamies [9] have proposed a two-potential framework in order to model the viscoelastic–dielectric elastomers. In this case, both the viscoelastic and dielectric relaxations are accounted for. The functional form of these two potentials allows electromechanical coupling and creep and relaxation processes after, respectively, applying and removing the mechanical and electrical inputs. The first potential is used for free energy and the second for dissipation. In this second case, the reduced dissipation inequality, for the case of isothermal processes, imposes additional constraints on the dissipation potential. The constitutive model used in [9] had been previously proposed by Lopez-Pamies [21], and it is expressed in terms of only the first invariant of the right Cauchy–Green deformation tensor. As in [4], a single relaxation time has been assumed for the viscoelastic as well as for the dielectric relaxation.

On this background, attention is paid to enlarge the viscoelastic and dielectric effects in electroelastic materials to consider the case of a distribution of relaxation times. For this purpose, the simplified approach used in Ref. [4] is followed in the present paper. Our formal strategy is to introduce the classical linear differential equations governing the relaxational behaviour of nonlinear viscous terms instead of the classical Newtonian ones. These nonlinear terms take the formal structure of fractional derivatives [22] which physically correspond to constant phase elements (CPEs), as defined in Appendix A in the structure of the equations governing the relaxational (retardational) processes. It should be noted that some progress in this direction has been made by Xiao et al. [23]. The constant phase elements induce the appearance of distributions (spectra) of relaxation (retardation) times. Moreover, since these distributions do not exhibit in practice symmetrical shape, the dynamic equations should be modified accordingly. Specifically, these modifications result in the introduction of two constant phase elements, respectively, for the low- and high-frequency sides of the spectrum.

Accordingly, the plan of the paper is as follows: after the introductory remarks in this section, in Section 2, the background of the paper is developed. In Section 3, the formalism of the fractional derivatives is introduced in the viscoelastic or dielectric models in order to give an account of a distribution of relaxation times. Moreover, the formal treatment in terms of fractional derivatives is developed. In Section 4, a double generalisation is made. The first one refers to viscoelastic or dielectric models and is addressed in order to obtain a nonsymmetric spectrum of relaxation times. The second is the adoption of the more realistic Mooney–Rivlin equation, instead of the classical neo-Hookean model, for the stress–strain relationship of the elastomeric material. Consequently, the theory developed in [4] is modified accordingly for the cases of compliant and floating electrodes. It should be noted that the adoption of the Mooney–Rivlin model equation ensures the appearance of the bifurcation phenomena. In Section 5, a practical example including a discussion on the dynamic mechanical measurements is analysed. In particular, the acrylic elastomer VHB 4910 from 3M is considered as usual. Section 6 is devoted to the numerical strategy for the calculations. In particular, the case where the system is subjected to an external electric force field but not to mechanical loadings is studied. A comparison between the results obtained by considering neo-Hookean and Mooney–Rivlin materials governed by a single relaxation time and the present results for the distribution of relaxation times is also made. Section 7 is devoted to analysing the effect of viscoelastic relaxation on the bifurcation processes appearing in the material under study. Section 8 and Section 9 are, respectively, concerned with a discussion of the results and the conclusions of the present research.

2. Theoretical Background

Elastomers show a nonlinear hyperelastic behaviour superimposed with dissipative viscoelastic and dielectric processes. The presence of these dissipative processes is due to the fact that the temperatures at which these materials are in use remain, in many cases, close to the glass transition temperature at which the viscoelastic and dielectric effects are still present. For example, let us consider the case of the VHB 4910, one of the materials that have merited preferential attention by part of the researchers due to its applicability. According to our experimental data [24], the glass transition temperature of the VHB 4910 measured by differential scanning calorimetry (DSC) at 20 °C/min is close to –37 °C. Although the tensile stress–strain analysis is carried out at 25 °C, it is still possible to detect experimentally some viscoelastic effects in the material behaviour. In any case, dynamic mechanical data should be the most convenient ones to analyse the present problem, as discussed later.

It is well known that the glass transition is detected in dynamic mechanical experiments by a drastic diminution in the viscoelastic modulus of pure amorphous polymers of about three decades in the value of the elastic modulus. However, in the case of the VHB 4910, only a diminution of about half a decade is observed. This fact can be due to the chemical structure of the material, or to the residual character of the viscoelastic effects. Moreover, an additional experimental problem appears due to the fact that at the temperatures above the glass transition the material is too soft. This fact frequently precludes obtaining reliable and reproducible dynamic results. In these conditions, the dynamic mechanical data should be derived from alternative techniques. In practice, many of the dynamic results are obtained indirectly from uniaxial extension experiments giving the extensional (Young) modulus as a function of the time or, alternatively, from the shear modulus. These facts are probably in the origin of the dispersion of the dynamic viscoelastic data at hand, as will be discussed later.

On the other hand, the classical nonlinear elastic models are time independent and, consequently, do not give an account of the viscoelastic behaviour which is frequency (time) dependent. As mentioned above, Suo et al. [4] have recently proposed an approach to give an account of the dissipative viscoelastic behaviour of the dielectric elastomers based on the CIT. On this basis, they have constructed a theory that, as a particular case, includes the annihilation of the classical hessian determinant to determine the stability or instability conditions. The aforementioned theory requires new ingredients to give an account of the irreversible processes suffered by the material during the relaxation processes. Then, a set of internal variables are included as a part of the Helmholtz free energy density. Subsequently, a differential equation to represent the kinetics of the actual irreversible processes appears. This differential equation needs to be solved in conjunction with the constitutive stress–strain relations derived from the free energy density. Suo et al. rightly state that ‘there is considerable flexibility in choosing kinetic models to fulfil the thermodynamic inequality’, which, in fact, is derived as a consequence of the second law of thermodynamics. In any case, the choice should be experimentally consistent. On this basis, an order one kinetic has been chosen to illustrate the general approach. A first-order kinetic, which is represented by a first-order differential equation, is equivalent to considering a single relaxation time associated with the relaxation process. This relaxation time is related to the dashpot appearing in the classical solid standard model (Zener model). It should be stressed that two types of standard solid models are possible. One is a spring in parallel with a Maxwell element, that is, a spring with a dashpot in series (Figure 1a), and it is used in this paper as the basis of the more complex models. The second one is a spring in series with a Kelvin–Voigt element, which is a spring in parallel with a dashpot (Figure 1b), and it has been used by other authors [2,7]. The first model is related to the case of a strain input (stress relaxation experiment), whereas the second one corresponds to the case in which a stress input is used (creep experiment). The parameters of both models are obviously related [15], and similar models have been proposed for the case of dielectric relaxation processes [18].

In this respect, the simplest dielectric model corresponds to the results of the classical Debye theory of the dielectric relaxation based on the rotational Brownian motion [25], which considers a single relaxation time. In fact, the Debye equivalent electric circuit includes two capacitors representing the passive elements for the storage of the electrical energy and a resistor to give an account of the dissipative process. The mechanical equivalent circuit in the electromechanical analogy is, of course, the above-mentioned Zener model consisting of two springs and a dashpot. However, the single relaxation time model initially proposed by Debye requires modification in order to fit the experimental results of the materials under study because of the actual presence of a distribution of relaxation times. A first alternative model was proposed by Cole and Cole [26] after studying the dielectric properties of some alcohols, showing non-Debye behaviour. This model implies the substitution of the resistor in the electric circuit by a CPE, which is an empirical impedance function (see Appendix A) and has a parabolic form in the time domain as can be found in the literature [16,18]. It is noticeable that the introduction of a CPE induces the appearance of a distribution of relaxation times. However, in the case of polymeric materials, the experimental data and the corresponding relaxation time spectra show a nonsymmetric shape. As a consequence, Havriliak and Negami (HN) [27] have modified the Cole–Cole model (CC) [26] to represent accurately the actual dielectric relaxation data of most of the polymers. The HN model predicts in fact a nonsymmetric relaxation time spectrum in agreement with the nonsymmetric shape of the experimental data and has been extensively used in practice. As an alternative to the HN equation, a biparabolic model has been proposed by Huet [28,29], which includes two CPEs. Each one of these CPE accounts for, respectively, the low- and high-frequency sides of the spectra. It should be noted that the parameters of the HN equation are closely related to those of the biparabolic equation [18] (p. 205). The equivalent mechanical model is easy to obtain, taking into account that the permittivity is, in the electromechanical analogy, mechanical compliance and assuming that the corresponding mechanical counterpart is an electric modulus. It should be noted that in mechanical terms the adoption of a CPE is equivalent to assume, instead of a simple Newtonian behaviour, a non-Newtonian one for the mechanical dissipative element. On this basis, a modification of the kinetic model used in [4] is pertinent.

3. Formal Analysis

The starting point to the formal analysis of the viscoelastic behaviour is the standard solid model proposed in [4] (Figure 1a), which is formed by a spring in parallel with a Maxwell element, that is, a spring in series with a dashpot. However, due to the fact that the viscoelastic primary data are obtained from uniaxial tensile modulus, the parameters of the standard solid model refer to the relaxed and unrelaxed tensile moduli and the corresponding tensile viscosity (respectively, E_R, E_U, and η_E).

The differential equation governing this model is given by

$$\left({\mathrm{E}}_{\mathrm{R}}+{\mathsf{\tau }}_{\mathrm{E}}{\mathrm{E}}_{\mathrm{U}}\frac{\mathrm{d}}{\mathrm{dt}}\right)\mathsf{ϵ}=\left(1+{\mathsf{\tau }}_{\mathrm{E}}\frac{\mathrm{d}}{\mathrm{dt}}\right)\mathsf{\sigma }$$

(1)

where $\mathsf{\sigma }$ and $\mathsf{ϵ}$ are, respectively, the stress and the strain, ${\mathsf{\tau }}_{\mathrm{E}}={\mathsf{\eta }}_{\mathrm{E}}/\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)$ is the characteristic relaxation time, ${\mathsf{\eta }}_{\mathrm{E}}$ (${\mathsf{\eta }}_{\mathrm{E}}>0$) is the viscosity coefficient associated with the dashpot, and ${\mathrm{E}}_{\mathrm{R}},{ \mathrm{E}}_{\mathrm{U}}$ are, respectively, the relaxed and unrelaxed tensile modulus. In forming this differential equation (see Ref. [15]), the following relation between the stress in the Maxwell element in parallel and the strain in the dashpot is used:

$${\mathsf{\sigma }}_2={\mathsf{\eta }}_{\mathrm{E}}\dot{\mathsf{ϵ}}$$

(2)

where a superscript dot indicates time derivative. Equation (2) corresponds to the classical Newtonian behaviour in the dashpot.

The stress response of the model given by Equation (1) in the time domain to a strain input given by $\mathsf{ϵ}={\mathsf{ϵ}}_0\mathrm{H}\left(\mathrm{t}\right)$, where $\mathrm{H}\left(\mathrm{t}\right)$ is the unit Heaviside function (step function), is given by

$$\mathsf{\sigma }\left(\mathrm{t}\right)={\mathsf{ϵ}}_0\left[{\mathrm{E}}_{\mathrm{R}}+\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)\exp \left(-\mathrm{t}/{\mathsf{\tau }}_{\mathrm{E}}\right)\right]$$

(3)

This seems to be the type of response depicted in Figure 2b of [4].

Then, the corresponding relaxation modulus in the time domain is given by

$$\mathrm{E}\left(\mathrm{t}\right)={\mathrm{E}}_{\mathrm{R}}+\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)\exp \left(-\mathrm{t}/{\mathsf{\tau }}_{\mathrm{E}}\right)$$

(4)

and the dynamic modulus in the frequency domain is obtained from the Laplace transform of both sides of (4) and is defined by

$${\mathrm{E}}^{*}=\mathrm{sE}\left(\mathrm{s}\right)={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right){\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}}{1+{\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}}={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)}{1+{\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{-1}}$$

(5)

where $\mathrm{s}=\mathrm{j}\mathsf{\omega }$, $\mathrm{E}\left(\mathrm{s}\right)=\mathsf{\sigma }\left(\mathrm{s}\right)/{\mathsf{ϵ}}_0$, and ${\mathrm{E}}^{*}=\mathsf{\sigma }\left(\mathrm{s}\right)/\mathsf{ϵ}\left(\mathrm{s}\right)$ is the defined dynamic modulus.

This is the classical single relaxation time dynamic response in the frequency domain corresponding to the standard solid model.

However, in order to give an account of the observed distribution of relaxation times experimentally observed, there are several possibilities, but in the present case, a slightly more sophisticated model is considered, as the one shown in Figure 2, where a mechanical CPE instead of a dashpot is used. The empirical mechanical impedance of this CPE may be expressed in terms of the parameters shown in Figure 2 as

$${\mathrm{E}}_{\mathrm{CPE}}^{*}=\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right){\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{\mathsf{\alpha }}$$

(6)

where $\mathsf{\alpha }\in \left(0,1\right)$.

Taking into account Equation (6), the total mechanical impedance of the scheme, represented in Figure 2, may be found as the addition of the impedance of the upper branch of the scheme in that figure and the inverse of the addition of the inverses of the impedance of the two elements of the lower branch of the same figure.

$${\mathrm{E}}^{*}={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right){\mathrm{E}}_{\mathrm{CPE}}^{*}}{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)+{\mathrm{E}}_{\mathrm{CPE}}^{*}}={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right){\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{\mathsf{\alpha }}}{1+{\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{\mathsf{\alpha }}}={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)}{1+{\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{-\mathsf{\alpha }}}$$

(7)

Since $\mathrm{s}$ represents in the time domain the first-order derivative with respect to the time, ${\mathrm{s}}^{\mathsf{\alpha }}$ can be interpreted in the present context as a fractional derivative of order $\mathsf{\alpha }$. The introduction of the constant $\mathsf{\alpha }$ is enough to create a distribution of relaxation times corresponding to the viscoelastic model given by Equation (7).

It can be shown that this result should be obtained by assuming a new relation between the stress and the rate of strain in the dashpot as follows:

$${\mathsf{\sigma }}_2={\mathsf{\eta }}_{\mathrm{E}}\frac{{\mathrm{d}}^{\mathsf{\alpha }}{\mathsf{\epsilon }}_2}{{\mathrm{dt}}^{\mathsf{\alpha }}}$$

(8)

where ${\mathsf{\sigma }}_2$ and ${\mathsf{\epsilon }}_2$ are, respectively, the stress in the bottom branch of the scheme shown in Figure 2, and the strain in the CPE of this Figure, and ${\mathsf{\eta }}_{\mathrm{E}}$ is a pseudo-viscosity.

Equation (8) may be considered as corresponding to a ‘fractional dashpot’ where $\mathsf{\alpha }=0, 1$ refer, respectively, to the ideal plasticity and the Newtonian liquid. In our case, being $0<\mathsf{\alpha }<1$ the obtained result in the time domain corresponds to the situation of nonexponential stress relaxation, as observed in the actual experiments.

Of course, more complex models different from that given by Equation (8) should be assumed as it will be shown later. Then, taking into account Equation (8), the differential Equation (1) is modified, and thereafter, the pertinent algebra can be reformulated according to (see Appendix B)

$$\left({\mathrm{E}}_{\mathrm{R}}+{\mathrm{E}}_{\mathrm{U}}{\mathsf{\tau }}_{\mathrm{E}}{}^{\mathsf{\alpha }}\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}\right)\mathsf{ϵ}=\left(1+{\mathsf{\tau }}_{\mathrm{E}}{}^{\mathsf{\alpha }}\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}\right)\mathsf{\sigma }$$

(9)

where ${\mathsf{\tau }}_{\mathrm{E}}{}^{\mathsf{\alpha }}=\frac{{\mathsf{\eta }}_{\mathrm{E}}}{{\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}}$, being ${\mathsf{\eta }}_{\mathrm{E}}$ the pseudo-viscosity, and $\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}\mathrm{f}\left(\mathrm{t}\right)$ denotes the fractional derivative of order $0<\mathsf{\alpha }\le 1$ of a function $\mathrm{f}\left(\mathrm{t}\right)$.

It is a well-known fact that there are several kinds of fractional derivatives [22]. In this work, we have chosen the Caputo derivative because its properties make it possible to consider initial conditions in terms of the value of the function $\mathrm{f}\left(\mathrm{t}\right)$ as is required in our context. The intuitive meaning of Caputo derivative lies in the concept of fractional integral, termed Liouville fractional integral, which is defined by

$$({\mathrm{J}}_0^{\mathsf{\alpha }}\mathrm{f})\left(\mathrm{t}\right)=\frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}_0^{\mathrm{t}}{\left(\mathrm{t}-\mathrm{s}\right)}^{\mathsf{\alpha }-1}\mathrm{f}\left(\mathrm{s}\right)\mathrm{ds}$$

(10)

where $\mathsf{\Gamma }\left(\mathsf{\alpha }\right)$ is the Euler’s Gamma function. Then, the Caputo derivative is defined by

$$\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}\mathrm{f}\left(\mathrm{t}\right)=({\mathrm{J}}_0^{\mathrm{n}-\mathsf{\alpha }}{\mathrm{f}}^{\left(\mathrm{n}\right)})\left(\mathrm{t}\right)$$

(11)

where $\mathrm{n}=-\left[-\mathsf{\alpha }\right],$ (being $\left[-\mathsf{\alpha }\right],$ the greatest integer less than or equal to $(-\mathsf{\alpha }$), and ${\mathrm{f}}^{\left(\mathrm{n}\right)}\left(\mathrm{t}\right)$ denotes de $\mathrm{n}$-th ordinary classical derivative of $\mathrm{f}\left(\mathrm{t}\right)$. Intuitively, in order to obtain the Caputo fractional derivative of $\mathrm{f}\left(\mathrm{t}\right)$ of order $\mathsf{\alpha }$, first, it is necessary to differentiate $\mathrm{n}$ times the function $\mathrm{f}\left(\mathrm{t}\right)$ and then to integrate the exceeding part. This latter term, consisting of calculating an integral (thus involving the evaluation of the function $\mathrm{f}\left(\mathrm{t}\right)$ over its integration domain), allows the interpretation of the Caputo derivative as a ‘memory’ operator since it takes into account the physical ‘history’ described by the function $\mathrm{f}\left(\mathrm{t}\right)$. In numerous physical contexts (see, for example, Refs. [30,31]), the Caputo derivative has demonstrated to successfully capture the ‘past history’ of the corresponding phenomena under study described via the function $\mathrm{f}\left(\mathrm{t}\right)$. Furthermore, it is easy to check, by means of the application of integration by parts the formula to $({\mathrm{J}}_0^{\mathsf{\alpha }}\mathrm{f})\left(\mathrm{t}\right)$, that the Caputo derivative becomes the ordinary derivative when $\mathsf{\alpha }$ is a positive integer, thus retaining the physical meaning of the classical derivative.

In our case, we consider that $0<\mathsf{\alpha }\le 1$; thus, the Caputo derivative of this specific order is given by

$$\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}\mathrm{f}\left(\mathrm{t}\right)=\frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)} {{\displaystyle \int }}_0^{\mathrm{t}}\frac{{\mathrm{f}}^{\prime }\left(\mathrm{s}\right)}{{\left(\mathrm{t}-\mathrm{s}\right)}^{1-\mathsf{\alpha }}}\mathrm{ds}, 0<\mathsf{\alpha }\le 1, \mathrm{t} >0.$$

(12)

Differential models in viscoelasticity including fractional derivatives were suggested by several authors [32,33,34,35,36,37] and used in connection with the viscoelastic or dielectric time-temperature superposition [38].

It is noteworthy that the dielectric counterpart of the mechanical model (Appendix A) consists of a capacitor with capacitance ${\mathrm{C}}_1$ in parallel with a series coupling of a capacitor with capacitance ${\mathrm{C}}_2$ and a CPE, as represented in Figure 3.

4. Generalisation for Non-Symmetric Relaxation Spectrum

As mentioned above, Equations (7) and (A5) predict a symmetric distribution of relaxation times [18], which is unlikely for polymeric materials. For this reason, a biparabolic model with two constant phase elements is proposed to represent the dynamic shear viscoelastic data [18]. The corresponding mechanical model is shown in Figure 4.

In a similar way to how is conducted in Equation (6), one can write

$${\mathrm{E}}_{1\mathrm{CPE}}^{*}=\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right){\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{\mathsf{\alpha }}, {\mathrm{E}}_{2\mathrm{CPE}}^{*}=\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right){\mathsf{\delta }}^{-1}{\left({\mathrm{s}\mathsf{\tau }}_{\mathrm{E}}\right)}^{\mathsf{\beta }}$$

(13)

for, respectively, the low- and high-frequency side of the spectrum. This explains why only two continuous phase elements at least are necessary.

Thus, after the pertinent calculations similar to those of Appendix B, the dynamic tensile modulus can be expressed as

$${\mathrm{E}}^{*}=\mathrm{sE}\left(\mathrm{s}\right)={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)}{1+{\left({\mathrm{s}\mathsf{\tau }}_1\right)}^{-\mathsf{\alpha }}+\mathsf{\gamma }{\left({\mathrm{s}\mathsf{\tau }}_1\right)}^{-\mathsf{\beta }}}={\mathrm{E}}_{\mathrm{R}}+\frac{\left({\mathrm{E}}_{\mathrm{U}}-{\mathrm{E}}_{\mathrm{R}}\right)}{1+{\left({\mathrm{s}\mathsf{\tau }}_1\right)}^{-\mathsf{\alpha }}+{\left({\mathrm{s}\mathsf{\tau }}_2\right)}^{-\mathsf{\beta }}}$$

(14)

where $0<\mathsf{\beta }<\mathsf{\alpha }<1$, ${\mathsf{\tau }}_1=\mathsf{\tau },{ \mathsf{\tau }}_2={\mathsf{\tau }\mathsf{\gamma }}^{-1/\mathsf{\beta }}$, and $\mathsf{\gamma }>0$ is a parameter associated with the nonsymmetric shape of the spectrum [39].

At this point, it is important to note that the viscoelastic parameters appearing in Equation (14) refer to the dynamic tensile modulus, whereas the elastic parameters that appear in the free energy density correspond to the shear modulus. In general, the relationship between the tensile and the shear viscoelastic moduli is given by

$${\mathrm{E}}^{*}=2{\mathsf{\mu }}^{*}\left(1+{\mathsf{\nu }}^{*}\right)$$

(15)

where ${\mathsf{\nu }}^{*}$ is the viscoelastic Poisson ratio. For incompressible materials, the Poisson ratio is commonly taken as ${\mathsf{\nu }}^{*}=0.5$ due to the incompressibility assumption. For this reason, in those that follow and in order to relate the dynamic viscoelastic tensile modulus to the dynamic shear modulus, the following equation is taken:

$${\mathrm{E}}^{*}=3{\mathsf{\mu }}^{*}$$

(16)

In these conditions, in those that follow, a similar expression to that of Equation (14) is taken for the shear modulus, that is,

$${\mathsf{\mu }}^{*}={\mathsf{\mu }}_{\mathrm{R}}+\frac{\left({\mathsf{\mu }}_{\mathrm{U}}-{\mathsf{\mu }}_{\mathrm{R}}\right)}{1+{\left({\mathrm{s}\mathsf{\tau }}_1\right)}^{-\mathsf{\alpha }}+{\left({\mathrm{s}\mathsf{\tau }}_2\right)}^{-\mathsf{\beta }}}$$

(17)

The use of the shear modulus $\mathsf{\mu }$ is justified because from the Hooke’s law and the relation between the stretching $\mathsf{\lambda }$ and the deformation $\mathsf{ϵ}$ for relatively small deformations, the classical result $\mathsf{\mu }={\mathrm{Nk}}_{\mathrm{B}}\mathrm{T}$ is obtained, where N is the number of chains of the network, ${\mathrm{k}}_{\mathrm{B}}$ is the Boltzmann constant and T the absolute temperature. As will be seen later (see Equations (31) and (32)), the shear modulus of the material is closely related to the parameters ${\mathsf{\mu }}_{\mathrm{Ui}}$ and ${\mathsf{\mu }}_{\mathrm{Ri}}$ appearing in Equation (18a,b) for the free energy density.

A similar model to the one presented in Equation (17) should be assumed for the dielectric permittivity of the elastomer by taking into account that the permittivity is a compliance function. However, as mentioned above, the dielectric relaxation process takes place very quickly in comparison with the viscoelastic one and, in these conditions, the dielectric relaxation process is in equilibrium, that is, the material is fully dielectrically relaxed when the viscoelastic process starts. On the other hand, it is assumed that possible conductive relaxation processes in the elastomer, which currently are proportional to ${\mathsf{\omega }}^{-\mathrm{s}}, (0<\mathrm{s}<1)$, are significantly active at frequencies considerably lower than the viscoelastic relaxation processes, and, for this reason, no overlapping between these two processes is expected. In Equation (17), $\mathsf{\alpha },\mathsf{\beta }$ and ${\mathsf{\tau }}_{1,2}$ are, respectively, the fractional-order of the $\mathsf{\xi }$ derivative, and the two relaxation times.

In order to obtain a better fit of the experimental results, a modified Mooney–Rivlin model has been chosen for the free energy density instead of the neo-Hookean one. Accordingly, the free energy density for the case in which viscoelastic effects are present may be written as

$$\begin{array}{c}\mathrm{W}=\frac{{\mathsf{\mu }}_{\mathrm{R}1}}{2}\left({\mathsf{\lambda }}_1^2+{\mathsf{\lambda }}_2^2+{\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}-3\right)+\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{2}\left({\mathsf{\lambda }}_1^{-2}+{\mathsf{\lambda }}_2^{-2}+{\mathsf{\lambda }}_1^2{\mathsf{\lambda }}_2^2-3\right)\\ +\frac{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}{2}\left({\mathsf{\lambda }}_1^2{\mathsf{\xi }}_{11}^{-2}{\mathsf{\xi }}_{12}^{-2}+{\mathsf{\lambda }}_2^2{\mathsf{\xi }}_{21}^{-2}{\mathsf{\xi }}_{22}^{-2}+{\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}{\mathsf{\xi }}_{11}^2{\mathsf{\xi }}_{12}^2{\mathsf{\xi }}_{21}^2{\mathsf{\xi }}_{22}^2-3\right)\\ +\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{2}\left({\mathsf{\lambda }}_1^{-2}{\mathsf{\xi }}_{11}^2{\mathsf{\xi }}_{12}^2+{\mathsf{\lambda }}_2^{-2}{\mathsf{\xi }}_{21}^2{\mathsf{\xi }}_{22}^2+{\mathsf{\lambda }}_1^2{\mathsf{\lambda }}_2^2{\mathsf{\xi }}_{11}^{-2}{\mathsf{\xi }}_{12}^{-2}{\mathsf{\xi }}_{21}^{-2}{\mathsf{\xi }}_{22}^{-2}-3\right)+\frac{{\tilde{\mathrm{D}}}^2}{2\mathsf{\epsilon }}{\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}\end{array}$$

(18a)

In the present approach, it is assumed that the permittivity $\mathsf{\epsilon }$ does not depend on the strain. As it has been discussed in another study [40], the polarity of the acrylic VHB 4910 elastomer corresponds to a type C polymer (in the classical Stockmayer´s terminology [41]), in which the dipoles causing the main dielectric activity are in the lateral chains and, according to it, they are separated from the backbone by flexible segments, and therefore, its permittivity should be scarcely affected by stretching.

Equation (18a) is a generalisation of Equation (17) of the Ref. [4] useful for our purposes. In this equation, ${\mathsf{\mu }}_{\mathrm{Ri}}$ and ${\mathsf{\mu }}_{\mathrm{Ui}}$ with $\mathrm{i}=1,2$ are the relaxed and unrelaxed shear moduli for the two branches of Figure 4. The first two terms on the right-hand side represent the elastic energy of the spring at the top of the scheme given in Figure 4. The third and fourth terms in this Figure correspond to the bottom branch of the before mentioned scheme. Moreover, ${\mathsf{\xi }}_{\mathrm{ij}}, \mathrm{i},\mathrm{j}=1,2$ are the stretches due to the two dashpots. Two subindexes are used to describe the stretching in these two dashpots: the first one refers to each one of the two in-plane stretches, whereas the second refers to the two CPE existing in the model. The last term represents the electrostatic energy of the material, with $\tilde{\mathrm{D}}$ being the nominal dielectric displacement. The relationship between the nominal dielectric displacement and the nominal electric field is given by $\tilde{\mathrm{D}}={\mathsf{\epsilon }\mathsf{\lambda }}_1^2{\mathsf{\lambda }}_2^2\mathrm{Ẽ}$, where $\mathsf{\epsilon }$ is the dielectric permittivity, and $\mathrm{Ẽ}$ is the nominal electric field. In the forthcoming calculations, it is assumed that the permittivity is not dependent on the stretches. Of course, many other empirical or semiempirical models should be checked, but for the present purposes, the chosen approach suffices.

It is noted that Equation (18a) is valid for the case of compliant electrodes (Figure 5a), that is, for the case in which the electrodes are deposited on the two faces of the sample, by ion sputtering, for example. However, in some circumstances, for example, nematic liquid crystal elastomers [42], the constrained geometry is avoided in order to eliminate anchoring effects in the boundaries of the sample. To accomplish this purpose, a gap between the sample and the electrodes is adjusted using adequate spacers in such a way that the distance between the electrodes is maintained systematically larger than the thickness of the sample (Figure 5b). The gap should be large enough to prevent pull-in phenomena.

In that case, the Maxwell stress tensor between the sample slab and the electrodes should be taken into account, and the free energy density, taking into account the viscoelasticity, is now given by

$$\begin{array}{c}\mathrm{W}=\frac{{\mathsf{\mu }}_{\mathrm{R}1}}{2}\left({\mathsf{\lambda }}_1^2+{\mathsf{\lambda }}_2^2+{\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}-3\right)+\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{2}\left({\mathsf{\lambda }}_1^{-2}+{\mathsf{\lambda }}_2^{-2}+{\mathsf{\lambda }}_1^2{\mathsf{\lambda }}_2^2-3\right)\\ +\frac{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}{2}\left({\mathsf{\lambda }}_1^2{\mathsf{\xi }}_{11}^{-2}{\mathsf{\xi }}_{12}^{-2}+{\mathsf{\lambda }}_2^2{\mathsf{\xi }}_{21}^{-2}{\mathsf{\xi }}_{22}^{-2}+{\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}{\mathsf{\xi }}_{11}^2{\mathsf{\xi }}_{12}^2{\mathsf{\xi }}_{21}^2{\mathsf{\xi }}_{22}^2-3\right)\\ +\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{2}\left({\mathsf{\lambda }}_1^{-2}{\mathsf{\xi }}_{11}^2{\mathsf{\xi }}_{12}^2+{\mathsf{\lambda }}_2^{-2}{\mathsf{\xi }}_{21}^2{\mathsf{\xi }}_{22}^2+{\mathsf{\lambda }}_1^2{\mathsf{\lambda }}_2^2{\mathsf{\xi }}_{11}^{-2}{\mathsf{\xi }}_{12}^{-2}{\mathsf{\xi }}_{21}^{-2}{\mathsf{\xi }}_{22}^{-2}-3\right)+\left({\mathsf{\epsilon }}^{-1}-{\mathsf{\epsilon }}_0^{-1}\right)\frac{{\tilde{\mathrm{D}}}^2}{2}{\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}\end{array}$$

(18b)

where ${\mathsf{\epsilon }}_0$ is the permittivity of the fluid where the sample is floating, which usually is the air.

Let us assume that geometry under consideration refers to a squared plate of electroviscoelastic material subjected to two pairs of forces or deformations in the principal directions. An electric field is applied perpendicular to the plate.

For equibiaxial stretches and in order to simplify calculations, it is assumed that

$${\mathsf{\lambda }}_1={\mathsf{\lambda }}_2=\mathsf{\lambda },{ \mathsf{\xi }}_{11}={\mathsf{\xi }}_{21}={\mathsf{\xi }}_1,{ \mathsf{\xi }}_{12}={\mathsf{\xi }}_{22}={\mathsf{\xi }}_2$$

(19)

Then, the following equations for the stresses and the electric field together with two kinetic equations are found (see Appendix C).

$$\begin{array}{l}{\mathsf{\sigma }}_1=\left({\mathsf{\mu }}_{\mathrm{R}1}+{\mathsf{\lambda }}_2^2{\mathsf{\mu }}_{\mathrm{R}2}\right)\left({\mathsf{\lambda }}_1-{\mathsf{\lambda }}_1^{-3}{\mathsf{\lambda }}_2^{-2}\right)\\ \hspace{1em}\hspace{1em}+\left[\left({\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}\right)+{\mathsf{\lambda }}_2^2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}\left({\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}\right)\right]\left({\mathsf{\lambda }}_1{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}-{\mathsf{\lambda }}_1^{-3}{\mathsf{\lambda }}_2^{-2}{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)-{\mathsf{\epsilon }}^{-1}{\tilde{\mathrm{D}}}^2{\mathsf{\lambda }}_1^{-3}{\mathsf{\lambda }}_2^{-2}\\ {\mathsf{\sigma }}_2=\left({\mathsf{\mu }}_{\mathrm{R}1}+{\mathsf{\lambda }}_1^2{\mathsf{\mu }}_{\mathrm{R}2}\right)\left({\mathsf{\lambda }}_2-{\mathsf{\lambda }}_2^{-3}{\mathsf{\lambda }}_1^{-2}\right)\\ \hspace{1em}\hspace{1em}+\left[\left({\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}\right)+{\mathsf{\lambda }}_1^2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}\left({\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}\right)\right]\left({\mathsf{\lambda }}_2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}-{\mathsf{\lambda }}_2^{-3}{\mathsf{\lambda }}_1^{-2}{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)-{\mathsf{\epsilon }}^{-1}{\tilde{\mathrm{D}}}^2{\mathsf{\lambda }}_2^{-3}{\mathsf{\lambda }}_1^{-2}\\ \mathrm{Ẽ}={\mathsf{\epsilon }}^{-1}{\tilde{\mathrm{D}}\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}\end{array}$$

(20a)

for Equation (17), and

$$\begin{array}{l}{\mathsf{\sigma }}_1=\left({\mathsf{\mu }}_{\mathrm{R}1}+{\mathsf{\lambda }}_2^2{\mathsf{\mu }}_{\mathrm{R}2}\right)\left({\mathsf{\lambda }}_1-{\mathsf{\lambda }}_1^{-3}{\mathsf{\lambda }}_2^{-2}\right)\\ \hspace{1em}\hspace{1em}+\left[\left({\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}\right)+{\mathsf{\lambda }}_2^2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}\left({\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}\right)\right]\left({\mathsf{\lambda }}_1{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}-{\mathsf{\lambda }}_1^{-3}{\mathsf{\lambda }}_2^{-2}{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)-\left({\mathsf{\epsilon }}^{-1}-{\mathsf{\epsilon }}_0^{-1}\right){\tilde{\mathrm{D}}}^2{\mathsf{\lambda }}_1^{-3}{\mathsf{\lambda }}_2^{-2}\\ {\mathsf{\sigma }}_2=\left({\mathsf{\mu }}_{\mathrm{R}1}+{\mathsf{\lambda }}_1^2{\mathsf{\mu }}_{\mathrm{R}2}\right)\left({\mathsf{\lambda }}_2-{\mathsf{\lambda }}_2^{-3}{\mathsf{\lambda }}_1^{-2}\right)\\ \hspace{1em}\hspace{1em}+\left[\left({\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}\right)+{\mathsf{\lambda }}_1^2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}\left({\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}\right)\right]\left({\mathsf{\lambda }}_2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}-{\mathsf{\lambda }}_2^{-3}{\mathsf{\lambda }}_1^{-2}{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)-\left({\mathsf{\epsilon }}^{-1}-{\mathsf{\epsilon }}_0^{-1}\right){\tilde{\mathrm{D}}}^2{\mathsf{\lambda }}_2^{-3}{\mathsf{\lambda }}_1^{-2}\\ \mathrm{Ẽ}=\left({\mathsf{\epsilon }}^{-1}-{\mathsf{\epsilon }}_0^{-1}\right){\tilde{\mathrm{D}}\mathsf{\lambda }}_1^{-2}{\mathsf{\lambda }}_2^{-2}\end{array}$$

(20b)

for Equation (18b).

Note that the equations for $\mathrm{Ẽ}$ (20a) (3rd equation) and (20b) (3rd equation) are the relation between the nominal electric field and the nominal displacement, to be used later.

Then, according to the methodology outlined in Appendix C and, after the pertinent algebra, one obtains for the kinetic restriction equations,

$$\begin{array}{c}\frac{{\mathrm{d}}^{\mathsf{\alpha }}{\mathsf{\xi }}_1}{{\mathrm{dt}}^{\mathsf{\alpha }}}=\frac{1}{{\mathsf{\tau }}_1^{\mathsf{\alpha }}}\left[\left({\mathsf{\lambda }}^2{\mathsf{\xi }}_1^{-3}{\mathsf{\xi }}_2^{-2}-{\mathsf{\lambda }}^{-4}{\mathsf{\xi }}_1^3{\mathsf{\xi }}_2^4\right)-\mathrm{k}\left({\mathsf{\lambda }}^{-2}{\mathsf{\xi }}_1{\mathsf{\xi }}_2^2-{\mathsf{\lambda }}^4{\mathsf{\xi }}_1^{-5}{\mathsf{\xi }}_2^{-4}\right)\right]\\ +\frac{1}{{\mathsf{\tau }}_2^{\mathsf{\alpha }}}\left[\left({\mathsf{\lambda }}^2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-3}-{\mathsf{\lambda }}^{-4}{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^3\right)-\mathrm{k}\left({\mathsf{\lambda }}^{-2}{\mathsf{\xi }}_1^2{\mathsf{\xi }}_2-{\mathsf{\lambda }}^4{\mathsf{\xi }}_1^{-4}{\mathsf{\xi }}_2^{-5}\right)\right]\end{array}$$

(21a)

$$\begin{array}{c}\frac{{\mathrm{d}}^{\mathsf{\beta }}{\mathsf{\xi }}_2}{{\mathrm{dt}}^{\mathsf{\beta }}}=\frac{1}{{\mathsf{\tau }}_1^{\mathsf{\beta }}}\left[\left({\mathsf{\lambda }}^2{\mathsf{\xi }}_1^{-3}{\mathsf{\xi }}_2^{-2}-{\mathsf{\lambda }}^{-4}{\mathsf{\xi }}_1^3{\mathsf{\xi }}_2^4\right)-\mathrm{k}\left({\mathsf{\lambda }}^{-2}{\mathsf{\xi }}_1{\mathsf{\xi }}_2^2-{\mathsf{\lambda }}^4{\mathsf{\xi }}_1^{-5}{\mathsf{\xi }}_2^{-4}\right)\right]\\ +\frac{1}{{\mathsf{\tau }}_2^{\mathsf{\beta }}}\left[\left({\mathsf{\lambda }}^2{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-3}-{\mathsf{\lambda }}^{-4}{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^3\right)-\mathrm{k}\left({\mathsf{\lambda }}^{-2}{\mathsf{\xi }}_1^2{\mathsf{\xi }}_2-{\mathsf{\lambda }}^4{\mathsf{\xi }}_1^{-4}{\mathsf{\xi }}_2^{-5}\right)\right]\end{array}$$

(21b)

where the initial conditions are ${\mathsf{\xi }}_1\left(0\right)={\mathsf{\xi }}_2\left(0\right)=1$, and $\mathrm{k}=\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}=\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}$.

It should be noted that there are many possibilities with respect to the choosing of the kinetic models to fulfil the thermodynamic requirements expressed by Equation (A14) in Appendix C. They are closely dependent on the properties of the electroelastic material under consideration.

5. Practical Application

For practical demonstration purposes, the acrylic elastomer VHB 4910 from 3M is the system to which the present approach is applied. Of course, other electroelastic materials should be considered.

5.1. Physicochemical Structure of VHB 4910

The chemical structure of VHB 4910 is formed by a two-block acrylic structure in which the length of each block and at least one of the substituted radicals in the lateral chain are commercially protected [43]. The average density of the material is 960 kg m⁻³, and the usual thickness of the tape is 1 mm.

In the following, it is considered that dielectric permittivity is not affected by stretching [40]. Despite the possible effect of the electric field on the morphological structure of the material, it is not taken into account in the present study.

5.2. Calorimetric Measurements

As reported in a previous paper [24], the glass transition temperature was calculated by differential scanning calorimetry (DSC) at 20 °Cmin⁻¹ using a TA Instruments Q-20 calorimeter. Measurements were carried out in a dry nitrogen atmosphere in a temperature range from −80 to 140 °C. The endothermic shift of the baseline was taken as the glass temperature (Tg = −37 °C). This temperature, which is in good agreement with previous results obtained by other authors [43,44,45], is about 60 degrees below that the temperature at which the experiments are carried out. Despite that, viscoelastic activity should be expected at this temperature. In fact, a requirement in order to have a good electroelastic response in the material is to show a glass transition temperature of about 50 °C to 60 °C lower than the temperature at which the experiments are carried out (typically, the room temperature). The material under study obviously accomplishes this requirement, but this is also valid for other rubber-like materials.

5.3. Dynamic Mechanical Measurements

As mentioned, the room temperature at which the future calculations are carried out is more than 60 degrees above the glass transition temperature. In these circumstances, the material is too soft to obtain reliable dynamic mechanical data.

Due to the above-mentioned problem, most authors try to use indirect methods to find truly reliable data. For example, Molberg et al. [46] use, as departing point, dynamic shear data from which the dynamic tensile modulus is obtained by assuming a constant Poisson´s ratio of 0.5, corresponding to an incompressible material. The data were obtained at 23 °C in a frequency range between 0.8 mHz and 80 Hz. The obtained results do not show an indication of a maximum in the curve representing loss modulus against the frequency. Moreover, the storage modulus is fitted to an exponential function of the frequency, a very unusual procedure in dynamic mechanical experiments. On the other hand, Suo et al. [4] use, in their calculations, data from the seminal papers by Wissler et al. [19] and Planté and Dubowski [20].

In [2], Lochmatter et al. fit a visco-hyperelastic film model by using a uniaxial tensile–creep–relaxation test sequence. They also assume a single relaxation time in the main equation, but they do not report the temperature of the experiment. As pointed out by Wang et al. [7], the model parameters according to the Lochmatter approach are obtained at very low stretch rates which led to the fact that they reflect basically the static (relaxed) viscoelastic modulus. For this reason, the relaxation strength, (Kp-Ks) in the terminology used by Lochmatter et al. [2] is very low. In fact, it is a well-known fact that the relaxation strength, which depends on the number of molecular moieties suffering relaxation processes, diminishes drastically with the temperature. Wang et al. [7] follow, in many respects, the main lines of Lochmatter´s paper using an ad hoc procedure based on uniaxial tensile strain at 25 °C with larger stretch rates, but they once more only consider a single relaxation time model. Moreover, the reported value of the relaxed modulus, Ep in Wang´s terminology, is larger than the lower value observed for the dynamic viscoelastic modulus (see Table 1 of Ref. [7]), which is unlikely, despite it can be due to the fitting procedure. When the data of Wang are used to fit an equation such as Equation (14), this fact is taken into account. On the other hand, when a single relaxation time is considered, the fitting procedure is relatively easy and the representative parameters of the viscoelastic model are obtained by using the nonlinear least-squares method [7]. However, when a more complex model equation such as Equation (14) is used for the fitting procedure in order to give an account of a distribution of relaxation times, the situation is more complicated. First of all, the starting data should be the dynamic shear modulus as appears in Equation (17).

Now, the two components of the dynamic shear modulus, that is, the storage and loss modulus ${\mathsf{\mu }}^{\prime }={\mathrm{Re} \mathsf{\mu }}^{*}, \mathsf{\mu } ″={\mathrm{Im} \mathsf{\mu }}^{*}$, respectively, need to be obtained. This problem has been classically addressed by using the Kramers–Kronig (KK) transforms [47,48], which are a consequence of the causality principle together with the linearity of the system, which is assumed in those that follow. A good account of the use of (KK) relations for the case of dielectric relaxation has been presented by Van Turnhout [49]. According to the theory, two alternatives are possible—the first one is to relate storage and loss modulus, and the second one, chose for convenience, is to relate the phase angle $\mathsf{\delta }\left(\mathsf{\omega }\right)={\tan }^{-1}\frac{ \mathsf{\mu }″}{{\mathsf{\mu }}^{\prime }}$ with the modulus |$\mathsf{\mu }$*| as follows [50]:

$$\mathsf{\delta }\left(\mathsf{\omega }\right)=\frac{2\mathsf{\omega }}{\mathsf{\pi }}{{\displaystyle \int }}_0^{\infty }\frac{\ln \left|{\mathsf{\mu }}^{*}\right|}{{\mathrm{x}}^2-{\mathsf{\omega }}^2}\mathrm{dx}$$

(22)

Partial integration of Equation (22) leads to

$$\mathsf{\delta }\left(\mathsf{\omega }\right)=\frac{1}{\mathsf{\pi }}{{\displaystyle \int }}_{-\infty }^{+\infty }\frac{\mathrm{dln}\left|{\mathsf{\mu }}^{*}\left(\mathrm{x}\right)\right|}{\mathrm{d}\ln \mathrm{x}}\ln \left|\frac{\mathrm{x}+\mathsf{\omega }}{\mathrm{x}-\mathsf{\omega }}\right|\mathrm{d}\ln \mathrm{x}$$

(23)

This equation has been proposed by Booij and Thoone [51]. These authors assume that the factor $\frac{\mathrm{dln}\left|{\mathsf{\mu }}^{*}\left(\mathrm{x}\right)\right|}{\mathrm{d}\ln \mathrm{x}}$ in (23) is nearly constant in a broad frequency band and after taking

$${{\displaystyle \int }}_{-\infty }^{+\infty }\ln \left|\frac{\mathrm{x}+\mathsf{\omega }}{\mathrm{x}-\mathsf{\omega }}\right|\mathrm{d}\ln \mathrm{x}=\frac{1}{2}{\mathsf{\pi }}^2$$

(24)

They obtain a simplified expression for $\mathsf{\delta }\left(\mathsf{\omega }\right)$ as

$$\mathsf{\delta }\left(\mathsf{\omega }\right)\cong \frac{\mathsf{\pi }}{2}\frac{\mathrm{dln}\left|{\mathsf{\mu }}^{*}\left(\mathrm{x}\right)\right|}{\mathrm{d}\ln \mathrm{x}}$$

(25)

However, this result is only a crude approximation for the phase angle.

Consequently, it is necessary to solve by numerical integration of Equation (23). For this purpose, the first step is to fit the empirical data for the viscoelastic shear modulus |${\mathsf{\mu }}^{*}$|, obtained from Equation (21a,b) and the data of [7] to a sigmoidal-type equation as follows:

$$\left|{\mathsf{\mu }}^{*}\right|=\left[{\mathsf{\mu }}_{\mathrm{R}}+\frac{{\mathsf{\mu }}_{\mathrm{U}}-{\mathsf{\mu }}_{\mathrm{R}}}{1+\exp \left(-1.6113\left(\mathrm{x}+3.7772\right)\right)}\right]·{10}^6 \mathrm{Pa}$$

(26)

where $\mathrm{x}=\ln \frac{\mathsf{\omega }}{2\mathsf{\pi }}=\ln \mathrm{f}$.

The relaxed and unrelaxed moduli in Equation (26) are approximately given by

$${\mathsf{\mu }}_{\mathrm{R}}=1.5·{10}^4 \mathrm{Pa};{ \mathsf{\mu }}_{\mathrm{U}}=6.25·{10}^4 \mathrm{Pa}$$

(27)

which are in sufficiently good agreement with the data of the Ref. [7].

It should be stressed that Equation (26) is only used to fit the data to a smooth function in order to obtain more data points for the forthcoming calculations and in any sense means a model for the viscoelastic data such as that has been previously chosen as Equation (17). In fact, Equation (26) does not yet predict the presence of two relaxation times corresponding to the high- and low-frequency side of the spectrum as appearing in the Huet model [28], represented by Equation (17).

After numerical integration of (23) using Equation (26), between the limits ${\mathrm{x}}_{\max }=\ln \mathrm{f}=-2{ \mathrm{and} \mathrm{x}}_{\min }=\ln \mathrm{f}=-7$, the phase angle is obtained in this convenient range of frequencies from which $\tan \mathsf{\delta }\left(\mathsf{\omega }\right)$ and, subsequently, the real and imaginary parts of the dynamic shear modulus (${\mathsf{\mu }}^{\prime }, \mathsf{\mu }″$) can be easily obtained according to

$${\mathsf{\mu }}^{\prime }=\left|{\mathsf{\mu }}^{*}\right|\cos \mathsf{\delta }\left(\mathsf{\omega }\right)$$

(28a)

$$\mathsf{\mu }″=\left|{\mathsf{\mu }}^{*}\right|\sin \mathsf{\delta }\left(\mathsf{\omega }\right)$$

(28b)

The obtained results are shown in

Table 1. Calculated values for viscoelastic parameters as a function of the frequency.

	Values from [7]			This Paper
ln f	tan δ	μ′ × 10−4 (Pa)	μ′′ × 10−4 (Pa)	tan δ	μ′ × 10-4 (Pa)	μ′′ × 10−4 (Pa)
−7				0.11	1.52	1.13
−6.5				0.18	1.54	0.28
−6	0.25	1.60	0.40	0.27	1.57	0.42
−5.5	0.38	1.70	0.65	0.39	1.66	0.65
−5	0.53	1.82	0.96	0.54	1.83	0.99
−4.5	0.68	2.156	1.46	0.68	2.17	1.48
−4	0.70	2.78	1.95	0.70	2.63	1.98
−3.5	0.54	3.75	2,03	0.55	3.65	2.12
−3	0.39	4.80	1.87	0.40	4.79	1.92
−2.5	0.28	5.51	1.54	0.28	5.51	1.54
−2	0.19	5.91	1.12	0.18	5.90	1.06
−1.5				0.10	6.10	0.61

, together with the data becoming from Ref. [7] after interpolation and using Equation (16) to obtain the shear modulus.

These data are used to fit the parameters appearing in Equation (17). The following parameters are obtained:

$$\mathsf{\alpha }=0.85, \mathsf{\beta }=0.81,{ \mathsf{\tau }}_1=6.77 \mathrm{s}. ,{ \mathsf{\tau }}_2=27.853\mathrm{s}.$$

(29)

For completeness, the analytical expressions for the storage and loss modulus of the dynamic modulus from Equation (17) are given in Appendix D.

5.4. Dielectric Measurements

As in the case of viscoelastic relaxation data, the frequency and temperature dependence of the dielectric permittivity should be taken into account. However, due to the fact that the viscoelastic relaxation times are five or more orders of magnitude larger than their dielectric counterparts, this effect is not taken into account in the present paper.

It is noted that the methodology of calculation should be the same as in the viscoelastic case. However, due to the fact that the frequency band of the dielectric data is larger than the viscoelastic one, the fitting procedure should be presumably easier.

6. Effect of the Electric Field without Mechanical Forces

In order to make application of the proposed methodology of calculation, let us consider now the case where the system is subjected to an external electric force field but not to mechanical forces. Then, Equation (19) is solved together with Equations (18a) or (18b) when ${\mathsf{\sigma }}_1={\mathsf{\sigma }}_2=0$ either for compliant or floating electrodes, which gives, respectively,

$$\left({\mathsf{\lambda }}^6-1\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^6{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}-{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)-\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^2-{\mathsf{\lambda }}^8\right)-\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^2{\mathsf{\xi }}_1^2{\mathsf{\xi }}_2^2-{\mathsf{\lambda }}^8{\mathsf{\xi }}_1^{-4}{\mathsf{\xi }}_2^{-4}\right)-\frac{\mathsf{\epsilon }}{{\mathsf{\mu }}_{\mathrm{R}1}}{\mathsf{\lambda }}^8{\mathrm{Ẽ}}^2=0$$

(30a)

for the case of compliant electrodes, and

$$\left({\mathsf{\lambda }}^6-1\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^6{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}-{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)-\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^2-{\mathsf{\lambda }}^8\right)-\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^2{\mathsf{\xi }}_1^2{\mathsf{\xi }}_2^2-{\mathsf{\lambda }}^8{\mathsf{\xi }}_1^{-4}{\mathsf{\xi }}_2^{-4}\right)+\left({\mathsf{\epsilon }}_0^{-1}-{\mathsf{\epsilon }}^{-1}\right)\frac{{\mathsf{\epsilon }}^2}{{\mathsf{\mu }}_{\mathrm{R}1}}{\mathsf{\lambda }}^8{\mathrm{Ẽ}}^2=0$$

(30b)

for the case of floating electrodes. In Equations (30a) and (30b), $\mathrm{Ẽ}$ is the nominal electric field.

In the case of the modified Mooney–Rivlin model, Equations (18a) and (b) for the relaxed and unrelaxed shear moduli satisfies

$$\mathsf{\mu }={\mathsf{\mu }}_1+{\mathsf{\mu }}_2$$

(31)

and using the relation

$$\mathrm{k}=\frac{{\mathsf{\mu }}_{\mathrm{U}2}}{{\mathsf{\mu }}_{\mathrm{U}1}}=\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}=0.12$$

(32)

one obtains ${\mathsf{\mu }}_{\mathrm{U}1}=5.5804,{ \mathsf{\mu }}_{\mathrm{R}1}=1.3393,{ \mathsf{\mu }}_{\mathrm{U}2}=0.6696,{ \mathsf{\mu }}_{\mathrm{R}2}=0.1607$, which will be used in the forthcoming calculations. As a consequence of the structure of Equation (30a,b), the use of the new figures for ${\mathsf{\mu }}_{\mathrm{U}1},{ \mathsf{\mu }}_{\mathrm{R}1},{ \mathsf{\mu }}_{\mathrm{U}2},{ \mathsf{\mu }}_{\mathrm{R}2}$ scarcely changes the obtained results with respect to the use of those of Equation (27).

The parameter $\mathrm{k}$ has been obtained from the data of Ref. [7] which, in time, has been taken from Figure 22 of Ref. [20] by using the classical procedure to fit the data to a Mooney–Rivlin equation (p. 103, [15]) for stretch ratios $\left(\mathsf{\lambda }\right)$ between 2 and 5.

7. Numerical Calculations

As the differential fractional algebraic equations (DFAE) given by (21a,b) do not have a closed solution, we approximate its solution by using a numerical method combining the fractional forward Euler method, for fractional differential equations applied to Equation (21a,b), and an iterative method to find the roots of the algebraic restriction given by Equation (30a,b). To explain the used methodology, an equation for a generic stretch $\mathsf{\xi }$ is considered. In this way, given a fractional initial value problem,

$$\frac{{\mathrm{d}}^{\mathsf{\alpha }}\mathsf{\xi }\left(\mathrm{t}\right)}{{\mathrm{dt}}^{\mathsf{\alpha }}}=\mathrm{f}\left(\mathrm{t},\mathsf{\xi }\left(\mathrm{t}\right)\right), \mathsf{\xi }\left(0\right)={\mathsf{\xi }}_0=1$$

(33)

By Theorem 3.1 of the Ref. [52], we can obtain the Volterra integral equation associated with the problem above, that is,

$$\mathsf{\xi }\left(\mathrm{t}\right)={ \mathsf{\xi }}_0+({\mathrm{J}}_0^{\mathsf{\alpha }}\mathrm{f})\left(\mathrm{t}, \mathsf{\xi }\left(\mathrm{t}\right)\right)$$

(34)

To obtain an approximate solution, first, we consider a discretisation of the fractional integral,${ \mathrm{J}}_0^{\mathsf{\alpha }}\mathrm{f}$, based on an equally spaced mesh for the time variable ${\left\{{\mathrm{t}}_{\mathrm{i}}=\mathrm{ih}\right\}}_{\mathrm{i}=0}^{\mathrm{N}},$ with $\mathrm{N}+1$ nodes defined in the interval $\left[0,\mathrm{T}\right], \mathrm{T}>0,$ and $\mathrm{h}:=\mathrm{T}/\mathrm{N}.$ The numerical approximation of the fractional integral of $\mathrm{f}\left(\mathrm{t}, \mathsf{\xi }\left(\mathrm{t}\right)\right)$ in ${\mathrm{t}}_{\mathrm{i}}$ is given by

$$\begin{array}{c}({\mathrm{J}}_0^{\mathsf{\alpha }}\mathrm{f})\left(\mathrm{t}, \mathsf{\xi }\left({\mathrm{t}}_{\mathrm{i}+1}\right)\right)=\frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)} {{\displaystyle \int }}_0^{{\mathrm{t}}_{\mathrm{i}+1}}{\left({\mathrm{t}}_{\mathrm{i}+1}-\mathrm{s}\right)}^{\mathsf{\alpha }-1} \mathrm{f}\left(\mathrm{s},\mathsf{\xi }\left(\mathrm{s}\right)\right)\mathrm{ds} =\frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)}{\displaystyle {\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{i}}}{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{j}}}^{{\mathrm{t}}_{\mathrm{j}+1}}{\left({\mathrm{t}}_{\mathrm{i}+1}-\mathrm{s}\right)}^{\mathsf{\alpha }-1} \mathrm{f}\left(\mathrm{s},\mathsf{\xi }\left(\mathrm{s}\right)\right)\mathrm{ds} \\ \approx \frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)}{\displaystyle {\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{i}}}\mathrm{f}\left({\mathrm{t}}_{\mathrm{j}},\mathsf{\xi }\left({\mathrm{t}}_{\mathrm{j}}\right)\right){{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{j}}}^{{\mathrm{t}}_{\mathrm{j}+1}}{\left({\mathrm{t}}_{\mathrm{i}+1}-\mathrm{s}\right)}^{\mathsf{\alpha }-1}\mathrm{ds}\\ =\frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }+1\right)} {\displaystyle {\displaystyle \sum }_{\mathrm{j}=0}^{\mathrm{i}}}\mathrm{f}\left({\mathrm{t}}_{\mathrm{j}},\mathsf{\xi }\left({\mathrm{t}}_{\mathrm{j}}\right)\right)\left[{\left(\mathrm{i}+1-\mathrm{j}\right)}^{\mathsf{\alpha }}-{\left(\mathrm{i}-\mathrm{j}\right)}^{\mathsf{\alpha }}\right]\end{array}$$

(35)

Observe that, in the first step, we have split the integral on the time interval $\left[0,{\mathrm{t}}_{\mathrm{i}+1}\right]$, into $\mathrm{i}$ integrals on the time intervals $\left[{\mathrm{t}}_{\mathrm{j}},{\mathrm{t}}_{\mathrm{j}+1}\right]$, in the second step, we have approximated the function $\mathrm{f}\left(\mathrm{s},\mathsf{\xi }\left(\mathrm{s}\right)\right)$ over the whole interval $\left[{\mathrm{t}}_{\mathrm{j}},{\mathrm{t}}_{\mathrm{j}+1}\right]$ just by $\mathrm{f}\left({\mathrm{t}}_{\mathrm{j}}, \mathsf{\xi }\left({\mathrm{t}}_{\mathrm{j}}\right)\right)$, and, finally in the last step, the remaining integral has been computed. With this approximation, we obtain the fractional forward Euler method, which is an extension of the classical Euler method as follows:

$$\mathsf{\xi }\left({\mathrm{t}}_{\mathrm{i}+1}\right)={ \mathsf{\xi }}_0+\frac{1}{\mathsf{\Gamma }\left(\mathsf{\alpha }+1\right)}{\displaystyle \sum }_{\mathrm{j}=0 }^{\mathrm{i}}\mathrm{f}\left({\mathrm{t}}_{\mathrm{j}},\mathsf{\xi }\left({\mathrm{t}}_{\mathrm{j}}\right)\right)\left[{\left(\mathrm{i}+1-\mathrm{j}\right)}^{\mathsf{\alpha }}-{\left(\mathrm{i}-\mathrm{j}\right)}^{\mathsf{\alpha }}\right]$$

(36)

According to our initial value problem, the fractional Euler method is given by

$${\mathsf{\xi }}_{\mathrm{i}+1}={ \mathsf{\xi }}_0+\frac{{\mathrm{h}}^{\mathsf{\alpha }}}{\mathsf{\Gamma }\left(\mathsf{\alpha }+1\right)}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{i}}\left[{\left(\mathrm{i}+1-\mathrm{j}\right)}^{\mathsf{\alpha } }-{\left(\mathrm{i}-\mathrm{j}\right)}^{\mathsf{\alpha }}\right] \left(\frac{1}{{\mathsf{\tau }}^{\mathsf{\alpha }}}\left({\mathsf{\lambda }}_{\mathrm{j}}^2{\mathsf{\xi }}_{\mathrm{j}}^{-3}-{\mathsf{\lambda }}_{\mathrm{j}}^{-4}{ \mathsf{\xi }}_{\mathrm{j}}^3\right)\right)$$

(37)

where $\mathrm{i}=0,\dots \mathrm{N}-1$, and ${\mathsf{\xi }}_{\mathrm{i}}$ and ${\mathsf{\lambda }}_{\mathrm{i}}$ are the approximations of $\mathsf{\xi }\left(\mathrm{t}\right)$ and $\mathsf{\lambda }\left(\mathrm{t}\right)$ in $\mathrm{t}={\mathrm{t}}_{\mathrm{i}}=\mathrm{ih}.$ Notice that in Equation (37), ${\mathsf{\xi }}_{\mathrm{i}+1}$ is obtained in terms of approximations constructed in previous nodes. The value of $\mathrm{h}$ is selected small enough to obtain a converged solution in each one of the models analysed.

These equations are valid for each one of the two $\mathsf{\xi }$ appearing in Equation (21a,b), and in this sense, the resulting equations are easily obtained from the former ones.

Once ${\mathsf{\xi }}_{\mathrm{i}}$ is obtained, the value of ${\mathsf{\lambda }}_{\mathrm{i}+1}$ is computed by finding the root of the equation for compliant electrodes (or their counterpart equation for floating electrodes),

$${\mathsf{\mu }}_{\mathrm{R}}\left({\mathsf{\lambda }}_{\mathrm{i}+1}^6-1\right)+\left({\mathsf{\mu }}_{\mathrm{U}}-{\mathsf{\mu }}_{\mathrm{R}}\right)\left({\mathsf{\lambda }}_{\mathrm{i}+1}^6{\mathsf{\xi }}_{\mathrm{i}+1}^{-2}-{\mathsf{\xi }}_{\mathrm{i}+1}^4\right)-{\mathsf{\epsilon }\mathsf{\lambda }}_{\mathrm{i}+1}^8{\mathrm{Ẽ}}^2=0$$

(38)

which is closer to ${\mathsf{\lambda }}_{\mathrm{i}}$. This root is computed numerically by using a combination of bisection, secant, and inverse quadratic interpolation methods [53].

This strategy is easily generalised when, as in our case, two parameters, ${\mathsf{\xi }}_1,{\mathsf{\xi }}_2$, appear in the restriction Equation (30a,b).

The following parameter is used in the calculations, $\mathsf{\epsilon }=39.843·{10}^{-12}{\mathrm{Fm}}^{-1}$, being ${\mathsf{\epsilon }}_0=8.854·{10}^{-12}{\mathrm{Fm}}^{-1}$ the permittivity of the evacuated space. The critical electric field is given by ${\mathrm{Ẽ}}_{\mathrm{c}}=0.6203\sqrt{\frac{{\mathsf{\mu }}_{{\mathrm{R}}_1}}{\mathsf{\epsilon }}}$ (see Appendix E). The corresponding results for the time evolution of $\mathsf{\lambda },{\mathsf{\xi }}_1,{\mathsf{\xi }}_2$ are shown in Figure 6 and Figure 7 for, respectively, compliant and floating electrodes and several electric fields.

The calculations have also been carried out for the case of a neo-Hookean model with a viscoelastic response in compliant mode governed by a single relation time ($\mathsf{\alpha }=1,\mathsf{\tau }=5.27 \mathrm{s}$), together with the following parameters: ${\mathsf{\mu }}_{\mathrm{R}}=1.5\times {10}^4\mathrm{Pa},{ \mathsf{\mu }}_{\mathrm{U}}=6.25\times {10}^4\mathrm{Pa}$ as in the Equation (27) and ${\mathrm{Ẽ}}_{\mathrm{c}}=0.6873\sqrt{\frac{{\mathsf{\mu }}_{\mathrm{R}1}}{\mathsf{\epsilon }}}$ as indicated in Appendix E and in agreement with Ref. [44]. The corresponding equations are

$$\left({\mathsf{\lambda }}^6-1\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}}-{\mathsf{\mu }}_{\mathrm{R}}}{{\mathsf{\mu }}_{\mathrm{R}}}\left({\mathsf{\lambda }}^6{\mathsf{\xi }}^{-2}-{\mathsf{\xi }}^4\right)-\frac{\mathsf{\epsilon }}{{\mathsf{\mu }}_{\mathrm{R}}}{\mathsf{\lambda }}^8{\mathrm{Ẽ}}^2=0$$

(39)

with the corresponding kinetic equation given by

$$\frac{\mathrm{d}\mathsf{\xi }}{\mathrm{dt}}=\frac{1}{\mathsf{\tau }}\left[{\mathsf{\lambda }}^2{\mathsf{\xi }}^{-3}-{\mathsf{\lambda }}^{-4}{\mathsf{\xi }}^3\right]$$

(40)

The results are shown in Figure 8.

For purpose of comparison, the results obtained by Suo et al. for a single relaxation time have been recovered by using our calculation program and taking into account the following parameters taken from Ref. [4] ${\mathsf{\mu }}_{\mathrm{U}}=8.78\times {10}^4 \mathrm{Pa},{ \mathsf{\mu }}_{\mathrm{R}}=2.28\times {10}^4 \mathrm{Pa}, \mathsf{\alpha }=1, \mathsf{\tau }=200 \mathrm{s}.$, together with a value of $\mathsf{\epsilon }=4.7$ for the relative permittivity, as reported by Wissler and Mazza in Ref. [16]. The agreement is excellent.

8. Stability and Bifurcation Analysis

In the linear theory of elasticity, solutions of traction boundary problems are unique, and the stress fields are also unique. This is a consequence of the linearity of the constitutive equations relating stress and strain. However, under finite deformations, the stress–strain relationships for rubbers are nonlinear, and unlike the linear theory, the uniqueness of the solutions is, a priori, not warranted. Concomitantly, the appearance of instabilities and possible bifurcation phenomena are expected. For this reason, stability analysis in order to check the possibility of the appearance of bifurcations is pertinent because nonlinear elasticity is exhibited by these electroelastic materials.

An analysis of the instability and possible bifurcation phenomena requires using Equation (18a) as starting point because in the case of the neo-Hookean model no instability and subsequent bifurcation, as expected, are predicted. The classical procedure outlined in [54,55] and based on the Hessian approach is followed here. Let us consider, as above, a slab whose thickness in the third direction is small, in comparison with the lateral dimensions subjected to two pairs of dead loads in these two directions. Moreover, a homogeneous electric field is applied in the third direction. Since the deformation and the field are basically homogeneous, the field equations are automatically fulfilled. For the case of equal dead loads on the four lateral surfaces, one has for each of Equation (18a)

$$\frac{\partial \mathrm{W}}{\partial {\mathsf{\lambda }}_1}=\frac{\partial \mathrm{W}}{\partial {\mathsf{\lambda }}_2}$$

(41)

After the pertinent algebraic calculations, the subsequent equation which can be factorised to give a symmetric solution defined by

$${\mathsf{\lambda }}_1={\mathsf{\lambda }}_2$$

(42)

that represents the bisector of the first quadrant in a ${\mathsf{\lambda }}_1 \mathrm{vs}.{ \mathsf{\lambda }}_2$ plot. Moreover, a nonsymmetric solution is found, which intersects with the symmetric solution after taking into account (19) for the values of $\mathsf{\lambda }$ solving the following equations for, respectively, compliant and floating electrodes

$$\left({\mathsf{\lambda }}^6+1\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^6{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}+{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)+\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left(3{\mathsf{\lambda }}^2-{\mathsf{\lambda }}^8\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left(3{\mathsf{\lambda }}^2{\mathsf{\xi }}_1^2{\mathsf{\xi }}_2^2-{\mathsf{\lambda }}^8{\mathsf{\xi }}_1^{-4}{\mathsf{\xi }}_2^{-4}\right)+{\mathsf{\epsilon }\mathsf{\mu }}_{{\mathrm{R}}_1}^{-1}{\mathsf{\lambda }}^8{\mathrm{Ẽ}}^2=0$$

(43a)

$$\left({\mathsf{\lambda }}^6+1\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}1}-{\mathsf{\mu }}_{\mathrm{R}1}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left({\mathsf{\lambda }}^6{\mathsf{\xi }}_1^{-2}{\mathsf{\xi }}_2^{-2}+{\mathsf{\xi }}_1^4{\mathsf{\xi }}_2^4\right)+\frac{{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left(3{\mathsf{\lambda }}^2-{\mathsf{\lambda }}^8\right)+\frac{{\mathsf{\mu }}_{\mathrm{U}2}-{\mathsf{\mu }}_{\mathrm{R}2}}{{\mathsf{\mu }}_{\mathrm{R}1}}\left(3{\mathsf{\lambda }}^2{\mathsf{\xi }}_1^2{\mathsf{\xi }}_2^2-{\mathsf{\lambda }}^8{\mathsf{\xi }}_1^{-4}{\mathsf{\xi }}_2^{-4}\right)-\left({\mathsf{\epsilon }}_0^{-1}-{\mathsf{\epsilon }}^{-1}\right){\mathsf{\epsilon }}^2{\mathsf{\mu }}_{{\mathrm{R}}_1}^{-1}{\mathsf{\lambda }}^8{\mathrm{Ẽ}}^2=0$$

(43b)

The time-dependent parameters $\mathsf{\lambda },{\mathsf{\xi }}_1,{\mathsf{\xi }}_2$ at which the bifurcation occurs are obtained by solving Equation (43a,b), together with Equation (21b). The parameters for the Equation (43a,b) used in the corresponding numerical calculations are the same as the used in the previous calculations. Figure 9 and Figure 10, respectively, show the effect of the electric field on the time evolution of the bifurcation stretch $\mathsf{\lambda }$ and ${\mathsf{\xi }}_1,{\mathsf{\xi }}_2$ for compliant and floating electrodes because we are dealing with a time-dependent bifurcation. Of course, other bifurcation modes and more complex bifurcation maps should be obtained by using, for example, incremental methods (see, for example, Ref. [56]), but for the present purposes, the approach used here is sufficient.

9. Discussion

The plot of the stretching $\mathsf{\lambda }$ and of the parameters ${\mathsf{\xi }}_1,{\mathsf{\xi }}_2$ of the dielectric elastomer under study for several applied electric fields in the compliant configuration is, respectively, shown in Figure 6 and Figure 7. In the case of compliant electrodes (Figure 6), starting from t = 0, the material relaxes with time, and subsequently, the stretch increases. For $\mathrm{Ẽ}<{\mathrm{Ẽ}}_{\mathrm{c}}$, the observed parameters tend to a new equilibrium state, but for $\mathrm{Ẽ}\ge {\mathrm{Ẽ}}_{\mathrm{c}}$, the elastomer becomes unstable at progressively higher electric fields. These results are in qualitative agreement with those observed in [4]. However, as can be seen in Figure 7, in the case of floating electrodes, the values of $\mathsf{\lambda },{ \mathsf{\xi }}_1,{\mathsf{\xi }}_2$ are lower than the unity decreasing with the time until an equilibrium value. This is due to the fact that the material sample is in compression under the effect of progressively higher electric fields.

In the case of a neo-Hookean material governed by a single relaxation time corresponding to Figure 8, the instability appears at lower times, as previously observed [4]. This fact suggests that the use of a more realistic model than the Mooney–Rivlin model is a more realistic strategy to achieve valuable results in the study of these materials.

The bifurcation analysis, which is shown in Figure 9 and Figure 10, indicates that the values of $\mathsf{\lambda },{ \mathsf{\xi }}_1,{\mathsf{\xi }}_2$ at the bifurcation increase with the time for the case of compliant, as well as for the floating electrodes. It should be noted that for ${\mathsf{\xi }}_1={\mathsf{\xi }}_2=0$ and in absence of electric field, a value of $\mathsf{\lambda }=2.897$ is obtained with $\frac{{\mathsf{\mu }}_{\mathrm{R}}}{{\mathsf{\mu }}_{\mathrm{U}}}=0.12$ from the following bifurcation equation:

$$\left({\mathsf{\lambda }}^6+1\right)+0.12\left(3{\mathsf{\lambda }}^2-{\mathsf{\lambda }}^8\right)=0$$

(44)

in agreement with the classical results. However, this value increases (diminishes) with the electric field for the case of compliant (floating) electrodes, in qualitative agreement with the results obtained in Ref. [54].

10. Conclusions

In this paper, as usual, an empirical model for the free energy density for a visco-hyperelastic material, VHB 4910, has been used as an example. In the present approach, no modelling is carried out on the possible electrical relaxation or conduction processes, because, as mentioned in the Introduction Section, when the viscoelastic relaxation processes take place, their electrical counterparts have occurred several decades before. Moreover, the treatment of the viscoelastic effects on the response of electroelastic materials, as outlined by Zhao and Suo [54], is modified in order to consider (a) the Mooney–Rivlin model instead of the neo-Hookean one, (b) a distribution of relaxation times instead of a single relaxation time, and (c) the asymmetric character [27] of the observed viscoelastic relaxations instead of a symmetric (CC) arc [26]. For this purpose, fractional derivatives are used in order to construct the respective constitutive equations, as explained in Section 3. The results obtained in the present paper suggest that the use of fractional derivatives for modelling the viscoelastic behaviour of electroelastic materials can give an account of the distribution of dielectric relaxation times. This seems to be a more realistic approach than those previously used to calculate the time evolution of the stretch $\mathsf{\lambda }$ and the viscoelastic parameters ${\mathsf{\xi }}_1,{\mathsf{\xi }}_2$ of these materials on the basis of a single relaxation time. Moreover, from the experimental point of view, the cases of compliant as well as floating electrodes are analysed. Both cases give very different results.

In order to compare our results with the former studies, the calculations are also carried out for the case of a neo-Hookean model with a viscoelastic model governed by a single relaxation time taking into account the parameters as in [4], and the agreement is very good.

On the other hand, the classical bifurcation analysis carried out for compliant as well as floating samples show results that are in agreement with those previously obtained for electroelastic, nonviscoelastic materials (Ref. [54]). It is to be stressed that the results obtained should potentially be useful in the design of electroelastic sensors and transducers as well as biomedical applications [57,58], prostheses, and robots, for example, in the case of cylindrical geometry for aortic prosthesis [59]. This is one of the main highlights of the present study, that is, to prevent the appearance of unexpected anomalous electroelastic responses.