Viscoelastic Effects on the Response of Electroelastic Materials

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.


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.

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 Polymers 2021, 13, 2198 5 of 28 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.
Polymers 2021, 13, x FOR PEER REVIEW 5 of 28 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.

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, ER, EU, and ηE).
The differential equation governing this model is given by where σ and ϵ are, respectively, the stress and the strain, τ = η /(E − E ) is the characteristic relaxation time, η (η > 0) is the viscosity coefficient associated with the

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 where σ and are, respectively, the stress and the strain, is the viscosity coefficient associated with the dashpot, and E R , E U are, respectively, the relaxed and unrelaxed tensile modulus. In Polymers 2021, 13, 2198 6 of 28 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: 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 = 0 H(t), where H(t) is the unit Heaviside function (step function), is given by This seems to be the type of response depicted in Figure 2b of [4]. Equation (8) may be considered as corresponding to a 'fractional dashpot' where α = 0, 1 refer, respectively, to the ideal plasticity and the Newtonian liquid. In our case, being 0 < α < 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) where τ = , being η the pseudo-viscosity, and f(t) denotes the fractional derivative of order 0 < α ≤ 1 of a function f(t).
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 f(t) 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 where Γ(α) is the Euler's Gamma function. Then, the Caputo derivative is defined by where n = −[−α], (being [−α], the greatest integer less than or equal to (−α ), and f ( ) (t) denotes de n-th ordinary classical derivative of f(t). Intuitively, in order to obtain the Caputo fractional derivative of f(t) of order α, first, it is necessary to differentiate n times the function f(t) and then to integrate the exceeding part. This latter term, consisting of calculating an integral (thus involving the evaluation of the function f(t) 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 f(t  Then, the corresponding relaxation modulus in the time domain is given by and the dynamic modulus in the frequency domain is obtained from the Laplace transform of both sides of (4) and is defined by where s = jω, E(s) = σ(s)/ 0 , and E * = σ(s)/ (s) 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 where α ∈ (0, 1). 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.
Polymers 2021, 13, 2198 7 of 28 Since s represents in the time domain the first-order derivative with respect to the time, s α can be interpreted in the present context as a fractional derivative of order α. The introduction of the constant α 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: where σ 2 and ε 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 η E is a pseudo-viscosity. Equation (8) may be considered as corresponding to a 'fractional dashpot' where α = 0, 1 refer, respectively, to the ideal plasticity and the Newtonian liquid. In our case, being 0 < α < 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) 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 f(t) 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 where Γ(α) is the Euler's Gamma function. Then, the Caputo derivative is defined by where n = −[−α], (being [−α], the greatest integer less than or equal to (−α), and f (n) (t) denotes de n-th ordinary classical derivative of f(t). Intuitively, in order to obtain the Caputo fractional derivative of f(t) of order α, first, it is necessary to differentiate n times the function f(t) and then to integrate the exceeding part. This latter term, consisting of calculating an integral (thus involving the evaluation of the function f(t) 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 f(t). 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 f(t). Furthermore, it is easy to check, by means of the application of integration by parts the formula to (J α 0 f)(t), that the Caputo derivative becomes the ordinary derivative when α is a positive integer, thus retaining the physical meaning of the classical derivative. In our case, we consider that 0 < α ≤ 1; thus, the Caputo derivative of this specific order is given by 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 timetemperature superposition [38].
It is noteworthy that the dielectric counterpart of the mechanical model (Appendix A) consists of a capacitor with capacitance C 1 in parallel with a series coupling of a capacitor with capacitance C 2 and a CPE, as represented in Figure 3.
It is noteworthy that the dielectric counterpart of the mechanical model (Appendix A) consists of a capacitor with capacitance C in parallel with a series coupling of a capacitor with capacitance C and a CPE, as represented in Figure 3.

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 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 where 0 < β < α < 1, τ = τ, τ = τγ / , and γ > 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

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.
It is noteworthy that the dielectric counterpart of the mechanical model (Appendix A) consists of a capacitor with capacitance C in parallel with a series coupling of a capacitor with capacitance C and a CPE, as represented in Figure 3.

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 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 where 0 < β < α < 1, τ = τ, τ = τγ / , and γ > 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  In a similar way to how is conducted in Equation (6), one can write 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 where 0 < β < α < 1, τ 1 = τ, τ 2 = τγ −1/β , and γ > 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 where ν * is the viscoelastic Poisson ratio. For incompressible materials, the Poisson ratio is commonly taken as ν * = 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: In these conditions, in those that follow, a similar expression to that of Equation (14) is taken for the shear modulus, that is, The use of the shear modulus µ is justified because from the Hooke's law and the relation between the stretching λ and the deformation for relatively small deformations, the classical result µ = Nk B T is obtained, where N is the number of chains of the network, k 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 µ Ui and µ 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 ω −s , (0 < 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), α, β and τ 1,2 are, respectively, the fractional-order of the ξ 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 In the present approach, it is assumed that the permittivity ε 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, µ Ri and µ Ui with 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, ξ ij , i, 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 D being the nominal dielectric displacement. The relationship between the nominal dielectric displacement and the nominal electric field is given by D = ελ 2 1 λ 2 2Ẽ , where ε is the dielectric permittivity, andẼ 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.
Polymers 2021, 13, x FOR PEER REVIEW 10 of 28 and fourth terms in this Figure correspond to the bottom branch of the before mentioned scheme. Moreover, ξ , i, 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 D being the nominal dielectric displacement. The relationship between the nominal dielectric displacement and the nominal electric field is given by D = ελ λ E , where ε is the dielectric permittivity, and E 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 (18) 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 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 where ε 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 Then, the following equations for the stresses and the electric field together with two kinetic equations are found (see Appendix C).
for Equation (17), and for Equation (18b). Note that the equations forẼ (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, where the initial conditions are ξ 1 (0) = ξ 2 (0) = 1, and . 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.

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.

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 −3 , 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.

Calorimetric Measurements
As reported in a previous paper [24], the glass transition temperature was calculated by differential scanning calorimetry (DSC) at 20 • Cmin −1 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.

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 tensilecreep-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 µ = Re µ * , µ = Im µ * , 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 δ(ω) = tan −1 µ µ with the modulus |µ*| as follows [50]: Partial integration of Equation (22) leads to This equation has been proposed by Booij and Thoone [51]. These authors assume that the factor (23) is nearly constant in a broad frequency band and after taking They obtain a simplified expression for δ(ω) as 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 |µ * |, obtained from Equation (21a,b) and the data of [7] to a sigmoidal-type equation as follows: where x = ln ω 2π = ln f. The relaxed and unrelaxed moduli in Equation (26) are approximately given by µ R = 1.5·10 4 Pa; µ U = 6.25·10 4 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 x max = ln f = −2 and x min = ln f = −7, the phase angle is obtained in this convenient range of frequencies from which tan δ(ω) and, subsequently, the real and imaginary parts of the dynamic shear modulus ( µ , µ ) can be easily obtained according to The obtained results are shown in Table 1, together with the data becoming from Ref. [7] after interpolation and using Equation (16) to obtain the shear modulus. Table 1. Calculated values for viscoelastic parameters as a function of the frequency.

Values from [7]
This Paper These data are used to fit the parameters appearing in Equation (17). The following parameters are obtained: For completeness, the analytical expressions for the storage and loss modulus of the dynamic modulus from Equation (17) are given in Appendix D.

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.

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 σ 1 = σ 2 = 0 either for compliant or floating electrodes, which gives, respectively, for the case of compliant electrodes, and for the case of floating electrodes. In Equations (30a) and (30b),Ẽ 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 µ = µ 1 + µ 2 (31) and using the relation one obtains µ U1 = 5.5804, µ R1 = 1.3393, µ U2 = 0.6696, µ R2 = 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 µ U1 , µ R1 , µ U2 , µ R2 scarcely changes the obtained results with respect to the use of those of Equation (27).
The parameter 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 (λ) between 2 and 5.

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 ξ is considered. In this way, given a fractional initial value problem, By Theorem 3.1 of the Ref. [52], we can obtain the Volterra integral equation associated with the problem above, that is, To obtain an approximate solution, first, we consider a discretisation of the fractional integral, J α 0 f, based on an equally spaced mesh for the time variable {t i = ih} N i=0 , with N + 1 nodes defined in the interval [0, T], T > 0, and h := T/N. The numerical approximation of the fractional integral of f(t, ξ(t)) in t i is given by Observe that, in the first step, we have split the integral on the time interval [0, t i+1 ], into i integrals on the time intervals t j , t j+1 , in the second step, we have approximated the function f(s, ξ(s)) over the whole interval t j , t j+1 just by f t j , ξ t j , 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: According to our initial value problem, the fractional Euler method is given by where i = 0, . . . N − 1, and ξ i and λ i are the approximations of ξ(t) and λ(t) in t = t i = ih. Notice that in Equation (37), ξ i+1 is obtained in terms of approximations constructed in previous nodes. The value of 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 ξ appearing in Equation (21a,b), and in this sense, the resulting equations are easily obtained from the former ones.
Once ξ i is obtained, the value of λ i+1 is computed by finding the root of the equation for compliant electrodes (or their counterpart equation for floating electrodes), which is closer to λ 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, ξ 1 , ξ 2 , appear in the restriction Equation (30a,b).
The following parameter is used in the calculations, ε = 39.843·10 −12 Fm −1 , being ε 0 = 8.854·10 −12 Fm −1 the permittivity of the evacuated space. The critical electric field is given byẼ c = 0.6203 . The corresponding results for the time evolution of λ, ξ 1 , ξ 2 are shown in Figures 6 and 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 (α = 1, τ = 5.27 s), together with the following parameters: µ R = 1.5 × 10 4 Pa, µ U = 6.25 × 10 4 Pa as in the Equation (27) andẼ c = 0.6873 µ R1 ε as indicated in Appendix E and in agreement with Ref. [44]. The corresponding equations are with the corresponding kinetic equation given by 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] µ U = 8.78 × 10 4 Pa, µ R = 2.28 × 10 4 Pa, α = 1, τ = 200 s., together with a value of ε = 4.7 for the relative permittivity, as reported by Wissler and Mazza in Ref. [16]. The agreement is excellent.      For purpose of comparison, the results obtained by Suo et al. for a time have been recovered by using our calculation program and taking following parameters taken from Ref. [4] μ = 8.78 × 10 Pa, μ = 2.28 × 200 s., together with a value of ε = 4.7 for the relative permittivity, as rep and Mazza in Ref. [16]. The agreement is excellent.

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) After the pertinent algebraic calculations, the subsequent equation which can be factorised to give a symmetric solution defined by that represents the bisector of the first quadrant in a λ 1 vs. λ 2 plot. Moreover, a nonsymmetric solution is found, which intersects with the symmetric solution after taking into account (19) for the values of λ solving the following equations for, respectively, compliant and floating electrodes The time-dependent parameters λ, ξ 1 , ξ 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. Figures 9 and 10, respectively, show the effect of the electric field on the time evolution of the bifurcation stretch λ and ξ 1 , ξ 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. Polymers 2021, 13, x FOR PEER REVIEW 21 of 28

Discussion
The plot of the stretching λ and of the parameters ξ 1 , ξ 2 of the dielectric elastomer under study for several applied electric fields in the compliant configuration is, respectively, shown in Figures 6 and 7. In the case of compliant electrodes (Figure 6), starting from t = 0, the material relaxes with time, and subsequently, the stretch increases. ForẼ <Ẽ c , the observed parameters tend to a new equilibrium state, but forẼ ≥Ẽ 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 λ, ξ 1 , ξ 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 Figures 9 and 10, indicates that the values of λ, ξ 1 , ξ 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 ξ 1 = ξ 2 = 0 and in absence of electric field, a value of λ = 2.897 is obtained with µ R µ U = 0.12 from the following bifurcation equation: 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].

Conclusions
In this paper, as usual, an empirical model for the free energy density for a viscohyperelastic 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 λ and the viscoelastic parameters ξ 1 , ξ 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.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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