Using a Fully Fractional Generalised Maxwell Model for Describing the Time Dependent Sinusoidal Creep of a Dielectric Elastomer Actuator

: Actuators made of dielectric elastomers are used in soft robotics for a variety of applications. However, due to their mechanical properties, they exhibit viscoelastic behaviour, especially in the initial phase of their performance, which can be observed in the ﬁrst cycles of dynamic excitation. A fully fractional generalised Maxwell model was derived and used for the ﬁrst time to capture the viscoelastic effect of dielectric elastomer actuators. The Laplace transform was used to derive the fully fractional generalised Maxwell model. The Laplace transform has proven to be very useful and practical in deriving fractional viscoelastic constitutive models. Using the global optimisation procedure called Pattern Search, the optimal parameters, as well as the number of branches of the fully fractional generalised Maxwell model, were derived from the experimental results. For the fully fractional generalised Maxwell model, the optimal number of branches was determined considering the derivation order of each element of the branch. The derived model can readily be implemented in the simulation of a dielectric elastomer actuator control, and can also easily be used for different viscoelastic materials.


Introduction
Dielectric elastomer actuators, also called DEAs, can have different structures.The first DEAs had parallel or cylindrical structures [1,2].From then on, many different structures were developed: spring roll actuators, helical actuators, and stacked actuators where the actuators are stacked on top of each other [3][4][5].All these actuators have the same basic structure and activation principle.Thus, regardless of whether the DEAs are parallel or cylindrical, their activation principle is the same.Their structure is like that of a capacitor.It has an upper and a lower electrode and an elastic dielectric.Some special rules apply to both the electrodes and the dielectric.The dielectric must be elastic.The electrodes must be conductive, and their mechanical structure must be the same or similar to that of the dielectric [6].If their mechanical properties are different, the DEA may lose its mobility and efficiency.Usually, conductive pastes of carbon are chosen for the electrodes.
The material properties of the selected elastomer VHB 4910 are viscoelastic, which means that the material has both viscous and elastic properties.When the DEA is subjected to a sinusoidal load, it exhibits sinusoidal creep behaviour.When a high voltage is applied, coupling occurs between the electrical and mechanical forces.Suo et al. were the first researchers to describe the electromechanical coupling in detail in a DEA [7].In the work of Gu et al. [8], the viscoelastic behaviour was captured using continuum mechanics and nonequilibrium thermodynamics, where the Helmholtz free energy along with the Gent model was used to derive the constitutive equations.The work shows a good match between the experimental and simulation results.However, extensive knowledge in the above areas is required, and the model cannot be integrated easily into the control simulations.The Prandtl-Ishlinskii model and the modified Prandtl-Ishlinskii model are used in the work of Zuo et al. [9,10].The work is based on the fourth order polynomial to describe the asymmetric behaviour, and fixed weights with thresholds of the play operators are used to describe the rate dependence, while the second order derivative of the input voltage is introduced into the fourth order polynomial to describe the peak-to-peak shifts that depend on the frequencies.This is a complex process that results in a good match between the experimental data and the calculated responses.The Prandtl-Ishlinskii model can be used as a feedforward compensator.Wissler et al. [11] used the Prony series theory to model the time dependent viscoelastic behaviour of the elastomer used in soft actuators.The drawback of this method is the inability to use it in the control simulations.A standard linear solid rheological model which satisfies the thermodynamic consistency was also used to capture rate-dependant behaviour of the soft materials in the work of Zhao et al. [12].This model could capture the rate dependent mechanical behaviour of soft materials.However, the introduction of the internal variables in the modelling process significantly increases the complexity of the proposed constitutive model and the number of governing equations.
Fractional calculus has proved to be a powerful mathematical tool in the approach used to modelling the time dependent mechanical behaviour of viscoelastic materials.Its advantages lie in the remarkable reduction in model parameters where fractional orders of derivatives are used.This can be especially seen in the field of fractional viscoelasticity, including soils [13], polymers [14,15] and construction materials [16].Xu et al. [17] used fractional constitutive models of a fractional Kelvin-Voigt model, a fractional Maxwell model and a fractional Poynting-Thomson model, which is a springpot connected in series with two parallel springpots.In the work of Barretta et al. [18], hereditariness and nonlocality of bending problems were presented with the help of fractional operators.Time-dependent hereditary behaviour, which is typical of viscoelastic materials, has been modelled with the help of a springpot, a fractional Kelvin-Voigt and a fractional Maxwell model.
To capture the sinusoidal creep behaviour, the fully fractional generalised Maxwell model is derived in Section 2. The generalised fractional Maxwell model has already been used in the work of Luo et al. [19] to determine the storage and loss modulus, but not with the Laplace transform, nor to derive the governing equation of the DEA.Our model is the first to describe the material behaviour of the DEA using the Laplace transform.The fully fractional generalised Maxwell model was used because it does not require complex knowledge of continuum mechanics and thermodynamics.It can readily be derived using the Laplace transform to obtain the transfer function of the system.Using the lsim function in the Matlab software, the excitation voltage with three different frequencies can be implemented in the transfer function of the system, and the response of the DEA can be calculated easily.An optimisation procedure using the Pattern Search global optimisation solver was used to derive the material parameters based on the experimental results.Section 3 shows the results of the optimisation procedures.Section 4 summarises the conclusions.

Principles of the DEA
The structure of the DEA is similar to that of a capacitor.It has an upper and lower conductive electrode and an elastic dielectric.The structure of the DEA can be seen in Figure 1.Usually, elastic materials are chosen for dielectrics, such as VHB tapes.VHB 4910 tape was used in this study [20].A high voltage DC is required to activate the DEAs.When the voltage is applied, the electrical force also known as the Maxwell force is generated between the upper and lower conductive electrodes, causing the DEA to contract vertically and expand longitudinally, as elastomer is known to be incompressible.
Fractal Fract.2022, 6, 720 3 between the upper and lower conductive electrodes, causing the DEA to contract v cally and expand longitudinally, as elastomer is known to be incompressible.The Maxwell force can be calculated as The relative permittivity was chosen as 4.7 according to [8].A conductive paste chosen from BareConductive ® [21].

Derivation of the Fully Fractional Generalised Maxwell Model
To derive the fully fractional generalised Maxwell model, one needs to use fracti derivatives.Fractional derivatives are derivatives that are not restricted to positive gers, but can be any real number.There are three definitions of fractional derivat namely, Riemann-Liouville, Caputo, and Gruenwald-Letnikov [22].The Gruenwaldnikov definition is used, namely time limits integer order of derivation/integration 0, derivation 0, integration because it can handle the fractional derivation and integration easily.It is also suitabl the numerical calculation of fractional derivatives.
The generalised Maxwell model is used to describe constitutive models with vis lastic effects.It consists of a spring connected in parallel to the branches of the Max elements.The Maxwell element consists of a spring connected in series with a das [15,16].To convert this into a fully fractional generalised Maxwell model, all element The Maxwell force can be calculated as ε 0 − absolute permittivity ε r = 4.7 − relative permittivity E − electric field V − voltage l 1 , l 2 , l 3 − dimensions of the DEA The relative permittivity was chosen as 4.7 according to [8].A conductive paste was chosen from BareConductive ® [21].

Derivation of the Fully Fractional Generalised Maxwell Model
To derive the fully fractional generalised Maxwell model, one needs to use fractional derivatives.Fractional derivatives are derivatives that are not restricted to positive integers, but can be any real number.There are three definitions of fractional derivatives, namely, Riemann-Liouville, Caputo, and Gruenwald-Letnikov [22].The Gruenwald-Letnikov definition is used, namely a, t − time limits m − integer order of derivation/integration p− = p > 0, derivation p < 0, integration Because it can handle the fractional derivation and integration easily.It is also suitable for the numerical calculation of fractional derivatives.
The generalised Maxwell model is used to describe constitutive models with viscoelastic effects.It consists of a spring connected in parallel to the branches of the Maxwell elements.The Maxwell element consists of a spring connected in series with a dashpot [15,16].
To convert this into a fully fractional generalised Maxwell model, all elements are replaced by the so-called springpot element [23,24].The springpot element has two parameters, one of which represents the material properties and the other the derivative order.The derivative order is bounded between 0 and 1, since there is no known physical interpretation above 1 [24].Figure 2 shows both models.
Fractal Fract.2022, 6, 720 4 of 1 parameters, one of which represents the material properties and the other the derivative order.The derivative order is bounded between 0 and 1, since there is no known physica interpretation above 1 [24].Figure 2 shows both models.The Laplace transform with Laplace operator is used to derive the governing equa tion of motion for the DEA using the fully fractional generalised Maxwell model.The La place transform turns the derivative into a multiplication and the integration into a divi sion.The fully fractional generalised Maxwell model describes the material behaviour The electric force is calculated using Equation (1).Since the weight is used to hold the DEA in the stretched position, its inertia must be included in the governing equation o motion.All branches of the fully fractional generalised Maxwell model are subject to the same displacement, as shown in Figure 3.    parameters, one of which represents the material properties and the other the derivative order.The derivative order is bounded between 0 and 1, since there is no known physical interpretation above 1 [24].Figure 2 shows both models.The Laplace transform with Laplace operator is used to derive the governing equation of motion for the DEA using the fully fractional generalised Maxwell model.The Laplace transform turns the derivative into a multiplication and the integration into a division.The fully fractional generalised Maxwell model describes the material behaviour.The electric force is calculated using Equation (1).Since the weight is used to hold the DEA in the stretched position, its inertia must be included in the governing equation of motion.All branches of the fully fractional generalised Maxwell model are subject to the same displacement, as shown in Figure 3.Each branch has two displacements because it has two springpots.The electric force and the force generated by the weight are distributed jointly between the branches of the fully fractional generalised Maxwell model.The governing equation of motion is calculated as Displacement of each branch with the fractional Maxwell element is calculated as The force in the first branch is calculated as The force in each branch containing two springpots is the same in each springpot, which is calculated as Expressing x i,1 and x i,2 from Equation ( 6) with the fractional integration, one gets Inserting the results from Equation (7) into Equation ( 4), one gets Force F i should be expressed from Equation (8).Laplace transformation is used, since force F i is part of the fractional integration with different orders of integration.Performing Laplace transformation on Equation (8) and expressing F i , one gets Performing the Laplace transform on Equation (3) and inserting Equation ( 9), one gets Rearranging Equation (10) to get the transfer function, one gets Equation ( 11) represents the transfer function of the governing equations of motion for the fully fractional generalised Maxwell model with n branches, where the influence of the electric force and the force of the weight of the mass 316 g are combined in the equation.The add-on FOMCON is required to implement Equation (11) in Matlab [25].Equation (11) is written in the relation force-displacement.To convert Equation ( 11) into the relation stress-strain, one needs to use and

Experiments
The elastomer VHB 4910 was used to set up the experiment.The elastomer was cut to the initial dimensions of 49 mm × 50 mm × 1 mm, with only a 10 mm wide area used to expand the elastomer to the initial dimension, while the rest was used for clamping.The initial dimensions for the sinusoidal force excitation were 49 mm × 60 mm × 0.16 mm.The active area to which the conductive paste was applied was 40 mm × 60 mm.
Three different frequencies with an amplitude of 6 kV DC were used for the sinusoidal voltage excitation.The three frequencies were F 1 = 1 13 Hz, F 2 = 1 7 Hz, and F 3 = 1 5 Hz.These frequencies were chosen so that they were not multiples of each other.In this way, the frequencies are not associated with a common factor.The experiment lasted 155 s to capture the sinusoidal creep behaviour of the DEA.The displacement of the DEA was measured using a Wenglor laser sensor [26].The experimental setup is shown in Figure 4.

Optimisation
To obtain the material parameters of the fully fractional generalised Maxwell model, the experimental results were optimised with the model.For the optimisation procedure,

Optimisation
To obtain the material parameters of the fully fractional generalised Maxwell model, the experimental results were optimised with the model.For the optimisation procedure, the global optimisation solver Pattern Search was chosen, which is also known as a direct search method in Matlab software.The flow of the optimisation procedure is shown in Figure 5. First, the user must specify the number of branches for the fully fractional generalised Maxwell model and the initial parameters for the model.

Optimisation
To obtain the material parameters of the fully fractional generalised Maxwell model, the experimental results were optimised with the model.For the optimisation procedure, the global optimisation solver Pattern Search was chosen, which is also known as a direct search method in Matlab software.The flow of the optimisation procedure is shown in Figure 5. First, the user must specify the number of branches for the fully fractional generalised Maxwell model and the initial parameters for the model.The electric and mechanical forces are calculated as the sum of the Maxwell force given by Equation ( 1) and the weight of the dead mass.In Equation (1), V is replaced with the calculation of the appropriate voltage regarding the frequency in use.Using the Matlab function lsim, the response of the transfer function of the fully fractional generalised Maxwell model to the electric and mechanical forces can be calculated easily.The responses for all frequencies were calculated and compared with the experimental results using the least squares method.The efficiency of the fully fractional generalised Maxwell model with a different number of branches was calculated using the method of the coefficient of determination, also known as the  method, calculated as  The electric and mechanical forces are calculated as the sum of the Maxwell force given by Equation (1) and the weight of the dead mass.In Equation (1), V is replaced with the calculation of the appropriate voltage regarding the frequency in use.Using the Matlab function lsim, the response of the transfer function the fully fractional generalised Maxwell model to the electric and mechanical forces can be calculated easily.The responses for all frequencies were calculated and compared with the experimental results using the least squares method.The efficiency of the fully fractional generalised Maxwell model with a different number of branches was calculated using the method of the coefficient of determination, also known as the R 2 method, calculated as

Start of optimisation
Total sum of squares.
Residual sum of squares.
After the optimisation procedure was completed, the best results were recorded and the final value of R 2 was calculated for each frequency.At the end, the average value of R 2 was calculated as R 2 mean for the whole frequency range.The number of branches was chosen to be between 1 and 5. Finally, each frequency was optimised individually for the model, to compare the results when only one frequency was optimised to the results where the whole frequency range was optimised.It was examined which the essential parameters of the model were for each frequency.Table 1 shows all symbols used in the research.

Results
After the optimisation procedure was set up using Matlab software and the Pattern Search global optimisation algorithm, the R 2 mean values shown in Table 2 were obtained.As can be seen, the number of branches increased the accuracy of the model, but not significantly.It can also be seen that adding more than three branches did not affect the efficiency of the model.The initial parameters used in the optimisation, as well as the optimal parameters obtained from the optimisation, are shown in Table 3.    From Table 3, in some branches the order of derivation does not change from the initial values.The derivation orders of 1 represent dashpot elements, and if both elements of the branch have the order of 1, they can be combined into one dashpot element.The same applies for spring elements if the order of derivation equals 0.
Figure 6 shows the experimental and calculated results from the optimised parameters of the fully fractional generalised Maxwell model for one to five branches.Adding more than two branches did not improve the efficiency of the model.The best matching between the data is seen for the middle frequency of 1/7 Hz, where the matching was 88%.For the frequency of 1/5 Hz, the worst matching was achieved of only 12.9%.
The optimisation of the model for an individual frequency was performed after the optimisation for the whole frequency range was carried out.From Table 2, it is seen that two additional branches of the fully fractional Maxwell elements are the most optimised.In Table 4, the initial and optimised parameters of only two additional branches are shown for each individual frequency for the fully fractional generalised Maxwell model.Each frequency demands its own material parameters, as well as different topology of the fully fractional generalised Maxwell model.Figure 7 shows matching between the experimental and calculated results if each frequency was optimised individually for the model.
Figure 6 shows the experimental and calculated results from the optimised param ters of the fully fractional generalised Maxwell model for one to five branches.Add more than two branches did not improve the efficiency of the model.The best match between the data is seen for the middle frequency of 1/7 Hz, where the matching was 88 For the frequency of 1/5 Hz, the worst matching was achieved of only 12.9%.

Discussion
The fully fractional generalised Maxwell model was used for the first time to describe the material behaviour of a DEA.Table 2 shows that the number of branches contributes to the effectiveness of the model.However, the contribution of each additional branch was small.Increasing the number of branches up to two increased the effectiveness.Surprisingly, adding more branches did not improve the effectiveness of the model.It can be seen from Table 3 that when adding more than two branches, the optimised parameters were chosen such that the model could be reduced to the fully fractional generalised Maxwell model with only two branches.This was possible because the two α parameters within a branch were 1, which represents dashpots.Dashpots in series can be reduced to one dashpot.The reduced and rearranged model can be seen in Figure 8.
order of derivative equal to 0.2 which means it has 80% characteristics of a spring and only 20% characteristics of a damper.The other springpot has the order of derivative equa to 0.52 which means that nearly half of its characteristics are those of a spring, and hal those of a damper.
The fully fractional generalised Maxwell model best describes the frequency of 1/7 Hz, where the match between the experimental and modelled response was 0.88, which was a good match.However, at the highest frequency of 1/5 Hz, the match between the data was the lowest and was only 0.12.The average match between the data over all three frequency ranges was 0.533, which is less than the methods used in the work of Gu et al [8] and Zuo et al. [9] where agreement between data was more than 0.9.If only data from frequency of 1/7 Hz were compared to the data from the work of Gu and Zuo, then ou method proves very efficient since it is much easier to implement it and use it in the con trol.Concluding points can be itemized as follow:

•
The number of fully fractional Maxwell elements slightly affected the effectiveness of the model.

•
Adding more than two branches did not increase the effectiveness of the model.

•
The fully fractional Maxwell model was reduced to the model seen in Figure 8.

•
The middle frequency of 1/7 Hz had the best agreement of 0.88 between data.

•
Optimising each frequency individually drastically improved the overall agreement between data to 0.745.

•
Optimising each frequency individually has a drawback since each frequency requires its own material parameters.

•
Topology optimisation cannot be included into the Pattern Search algorithm.
The reduced and rearranged model has two springpots.The first springpot has the order of derivative equal to 0.2 which means it has 80% characteristics of a spring and only 20% characteristics of a damper.The other springpot has the order of derivative equal to 0.52 which means that nearly half of its characteristics are those of a spring, and half those of a damper.
The fully fractional generalised Maxwell model best describes the frequency of 1/7 Hz, where the match between the experimental and modelled response was 0.88, which was a good match.However, at the highest frequency of 1/5 Hz, the match between the data was the lowest and was only 0.12.The average match between the data over all three frequency ranges was 0.533, which is less than the methods used in the work of Gu et al. [8] and Zuo et al. [9] where agreement between data was more than 0.9.If only data from frequency of 1/7 Hz were compared to the data from the work of Gu and Zuo, then our method proves very efficient since it is much easier to implement it and use it in the control.
The initial and optimised parameters for the optimisation of the individual frequency with the fully fractional generalised Maxwell model are listed in Table 3.The match between the experimental and calculated responses was improved.The average match increased to 0.745, which is a good match and can be easily compared with the work of Gu and Zuo.However, the material parameters of the fully fractional generalised Maxwell model were different for each frequency, which was a drawback, since one would like to have universal material parameters for the entire frequency range.
Finding the optimal number of branches for the model cannot be included in the optimisation, because changing the number of branches changed the number of lower and upper bounds on the model.This is something that cannot be included in the optimisation solver for the Pattern Search.This can only be done by observing the material parameters of the individual branches within the model manually.

Conclusions
The match between experimental and calculated results was lower for the whole frequency range than in the work of others.Observing only the middle frequency, the method would be easily compared with the work of others.However, the fully fractional generalised Maxwell model can be derived and implemented easily, and the responses of the model can be determined quickly.Implementation of the model in simulation control is straightforward.The proposed method can easily be used on other materials with viscoelastic behaviour.In future work, topology optimisation could be included into the optimisation procedure.

Figure 1 .
Figure 1.Structure of the parallel DEA: (a) initial state and (b) activated state.

Figure 1 .
Figure 1.Structure of the parallel DEA: (a) initial state and (b) activated state.

Figure 2 .
Figure 2. Two types of models used: (a) Fully fractional model and (b) Generalised model.

Figure 3 .
Figure 3. Displacement of the elements.

Figure 2 .
Figure 2. Two types of models used: (a) Fully fractional model and (b) Generalised model.

c i, 1 ,
c i,2 − material properties of springpot α i,1 , α i,2 − order of derivation k i − spring constants m − mass of weight n − number of branches with the Maxwell element Fel − electrical force applied Fg − weight of the mass The Laplace transform with Laplace operator is used to derive the governing equation of motion for the DEA using the fully fractional generalised Maxwell model.The Laplace transform turns the derivative into a multiplication and the integration into a division.The fully fractional generalised Maxwell model describes the material behaviour.The electric force is calculated using Equation (1).Since the weight is used to hold the DEA in the stretched position, its inertia must be included in the governing equation of motion.All branches of the fully fractional generalised Maxwell model are subject to the same displacement, as shown in Figure3.

Figure 2 .
Figure 2. Two types of models used: (a) Fully fractional model and (b) Generalised model.

Figure 3 .
Figure 3. Displacement of the elements.Figure 3. Displacement of the elements.

Figure 3 .
Figure 3. Displacement of the elements.Figure 3. Displacement of the elements.
al Fract.2022, 6, 720 7 of 14 capture the sinusoidal creep behaviour of the DEA.The displacement of the DEA was measured using a Wenglor laser sensor [26].The experimental setup is shown in Figure 4.

Figure 4 .
Figure 4. Experimental set up of the DEA.

Figure 4 .
Figure 4. Experimental set up of the DEA.

Figure 4 .
Figure 4. Experimental set up of the DEA.

Figure 6 .
Figure 6.Experimental vs. calculated responses of the DEA for different numbers of branches (n the fully fractional generalised Maxwell model.

Figure 6 .
Figure 6.Experimental vs. calculated responses of the DEA for different numbers of branches (n) in the fully fractional generalised Maxwell model.

Figure 7 .
Figure 7. Responses of individually optimised frequencies with the generalised fractional Maxw model.

Figure 7 .
Figure 7. Responses of individually optimised frequencies with the generalised fractional Maxwell model.

Figure 8 .
Figure 8. Reduced and rearranged optimised fully fractional generalised Maxwell model with two additional branches.

Figure 8 .
Figure 8. Reduced and rearranged optimised fully fractional generalised Maxwell model with two additional branches.

Table 1 .
Nomenclature table with SI units.

Table 2 .
R 2 mean value of the fully fractional generalised Maxwell model for one to five branches.

Table 3 .
Initial and optimised parameters of the fully fractional generalised Maxwell model for one to five branches.

Table 4 .
Initial and optimised parameters of the fully fractional generalised Maxwell model for each individual frequency.

Table 4 .
Initial and optimised parameters of the fully fractional generalised Maxwell model for e individual frequency.