ON AN INTEGRATED DYNAMIC CHARACTERIZATION OF VISCOELASTIC MATERIALS BY FRACTIONAL DERIVATIVE AND GHM MODELS

Abstract The passive vibration control of mechanical systems under unwanted vibrations can be accomplished in a very effective way by using devices incorporating viscoelastic materials. The design of such devices requires a broad knowledge of the dynamic properties of the employed viscoelastic material, usually supplied by adequate mathematical models. Among the available mathematical models, the fractional derivative (FD) model and the Golla-Hughes-McTavish (GHM) model, along with either the Williams-Landel-Ferry (WLF) equation or the Arrhenius equation, are now very prominent. The current work investigates the use of these models in a wide and integrated dynamic characterization of a typical and thermorheologically simple viscoelastic material. It focuses on experimental data collected from 0.1 to 100 Hz and –40 °C to 50 °C, which are simultaneously manipulated to raise both the frequency and the temperature dependencies of the material. In fitting the models, a hybrid approach combining techniques of genetic algorithms and nonlinear optimization is adopted. The ensuing results are evaluated by means of objective function values, comparative experimental-predicted data plots, and the Akaike’s Information Criterion (AIC). It is shown that the four-parameter fractional derivative model presents excellent curve fitting results. As for the GHM model, its modified version is the most adequate, although a higher number of terms is required for a satisfactory goodness-of-fit. None the less the fractional derivative model stands out.


INTRODUCTION
Effective actions of vibration control with viscoelastic devices generally require a previous and wide knowledge of the dynamic behavior of the employed viscoelastic materials, particularly their so-called dynamic properties: the dynamic elasticity modulus and the corresponding loss factor.Among other factors, these properties are dependent on frequency and temperature, and this dependence can be quite pronounced (Nashif et al., 1985).
A usual and established way of describing the dynamic behavior of viscoelastic materials, combining the previously mentioned properties, is the representation by complex moduli (Snowdon, 1968).In a complex modulus description,

DYNAMIC CHARACTERIZATION
All viscoelastic solids have elastic and dissipative properties.Such properties may be well represented, in the linear region, by complex functions.These complex functions are called 'complex moduli of elasticity' (Snowdon, 1968;Pritz, 1998).Thus, there are the complex bulk, shear and Young moduli, and also the complex Poisson's ratio.
The complex Young modulus  � [Pa] represents the dynamic relation between the uniaxial longitudinal stress σ and the associated uniaxial longitudinal deformation ε.Among other factors, it is clearly dependent on frequency and temperature (Nashif et al., 1985).Therefore, it can be expressed as where  is frequency [rad/s], T is temperature [K],   is the real Young modulus [Pa], associated with the stored energy, and   is the imaginary Young modulus [Pa], associated with the dissipated energy, and i is the imaginary unit.
Defining the loss factor   as the ratio of   to   (it is, therefore, dimensionless), the modulus can  � be rewritten as (Ω, ) =   (Ω, )[1 +   (Ω, )] (2) There is a similar expression for the complex shear modulus G [Pa], which is DERIVATIVE AND GHM MODELS Wagner Barbosa de Medeiros Júnior et al.
Latin American Journal of Solids and Structures, 2019, 16(2), e164 3/19 (Ω, ) =   (Ω, )[1 +   (Ω, )] (3) Experimental data concerning the so-called dynamic properties R E and E  is obtained from tests in a range of frequencies and temperatures of interest, generating curves like those in Figure 1.The experimental data can be assembled by the frequency-temperature superposition principle (Nashif et al., 1985).This principle states that the various curves of dynamic properties obtained over a frequency range in distinct temperatures (Fig. 1) can be superimposed, in a temperature of reference, by means of shifts in frequency.Thus, as seen in Figure 2, single master curves are formed in a standardized plot called reduced frequency nomogram (Jones, 2001).In mathematical terms, the superposition principle establishes that and where Ω  is the so-called reduced frequency [rad/s],  0 is the temperature of reference [K],  is the material density [kg/m 3 ], and  0 is the material density at the reference temperature [kg/m 3 ].
The above expressions show that, except for the factor ( 0  0 /) (usually negligible), the values of the dynamic properties at a frequency Ω and a temperature T are equal to the values of the same properties at a reduced frequency Ω  and the reference temperature  0 .
The reduced frequency Ω  is given by where   (), which is dimensionless, is known as the 'shift factor'.The shift factor, as the name implies, accounts for the shifts in frequency required to take the curves of dynamic properties from their specific measurement temperatures to the reference temperature.
The shift factor can be modeled in several ways.The most employed models are: the WLF (Williams-Landel-Ferry) equation, in which where  1 and  2 are specific material parameters to be raised from experimental data; and the Arrhenius equation, in which where   is the so-called 'activation temperature' [K], also a material parameter.It is considered that, from a practical point of view, any of these equations can be employed (Mead, 2000).
For the sake of usability, mathematical models can be fitted to the master curves presented in a reduced frequency nomogram.Two of those models, the fractional derivative model and the GHM model, stand out.

FRACTIONAL DERIVATIVE (FD) MODEL
The one-dimensional fractional derivative constitutive equation at a given temperature is (Bagley, 1979;Pritz, 2003) where   �   �,   (dimensionless),  0 [Pa],   [Pa], and   (dimensionless) are material parameters, with m from 1 to M and n from 1 to N. As to    and    , they are fractional order derivatives, the Riemann-Liouville (RL) definition of which is where  is the fractional derivative order, such that 0 <  < 1, Γ is the gamma function, and  is an integration variable (Bagley and Torvik, 1983;Nashif et al., 1985).
As pointed out by Bagley andTorvik (1983, 1986) and Pritz (2003), the fractional calculus approach has been shown to be an efficient and robust tool for describing the dynamic behavior of real viscoelastic materials, particularly elastomers used in noise and vibration control designs.That has also been observed by Lopes et al. (2004), Espíndola et al. (2005) andMedeiros Jr (2010), in instances particularly familiar to the current authors.
Similarly to integer order derivatives, it can be shown that the Fourier transform (denoted by F) of the above defined fractional order derivative of a certain function f(t) is given by Latin American Journal of Solids and Structures, 2019, 16(2), e164 5/19 Considering Eq. ( 9) with only four material parameters, it follows that in which  1 =  1 , based on thermodynamic considerations (Bagley and Torvik, 1986).The preceding equation can be represented by an arrangement of mechanical elements as shown in Cinielo et al. (2017).Fourier transforming both sides of Eq. ( 12) gives An alternative expression for (Ω) is where  1 =  and  1 = , for simplicity, and  1 =  ∞  for convenience.This expression is known as the fourparameter fractional derivative model for the complex Young modulus.
The inclusion of temperature dependence in the above equation, via the shift factor T  , gives 4 GHM MODEL

Standard GHM (SGHM) model
The standard (original) GHM model can be introduced as follows: consider a simple structure under longitudinal stress and assume the existence of a material relaxation function E(t) [Pa] given by () =   ()+H(t) (16) where () is a unit step function,   = lim  →∞ () is the elastic equilibrium modulus of the material [Pa], and H(t) is the time dependent relaxation function [Pa] (Moschen, 2006;Christensen, 1982).
It is known that the resultant stress for a sudden deformation applied at t = 0 can be written as either or, given Eq. ( 16), Besides linearity and causality, it is also assumed above that the material is relaxed.
As the GHM model is developed in the Laplace domain, it is worth noting that, in the current context, the Laplace transform (denoted by L) of a function f(t) can be given by (Papoulis, 1962) Thus, Laplace transforming Eqs. ( 17) and ( 18) provides, respectively, and From Eqs (20) and ( 21), it follows that For the sake of simplicity, the terms () and (), which are the complex Young modulus and the complex dissipation function in the Laplace domain, can be denoted as ( ) E s and ( ) h s , in such a way that The complex dissipation function ℎ() has several forms and, therefore, it must be chosen by the analyst (Golla and Hughes, 1985;Barbosa, 2000).Usually, a function obtained from the Biot series (Biot, 1955) is employed in the following form This series is truncated to the second term, with  2 >  1 >  0 , to produce what is called a 'dissipation pole pair' (Golla and Hughes, 1985), commonly written as where Applying the restriction  � =  � ̂, established by Golla and Hughes (1985) when the modeling of vibrating systems is considered, the expression for the complex Young modulus in the Laplace domain is obtained, which is Assuming that the Fourier transform can be treated as a special case of the Laplace transform, with  =  (see Papoulis, 1962), the complex Young modulus in the frequency domain, for a given temperature, can be written, from Eq. ( 26), as Equation ( 27) is the standard GHM counterpart of Eq. ( 14), which describes the complex Young modulus by the use of fractional calculus.
Including the temperature dependence in Eq. ( 27), as done in the previous section, gives

Additional GHM Terms
In order to achieve satisfactory curve fitting results for the GHM model (regarding experimental data), the strategy of expanding the model by increasing the number of dissipation pole pairs in the dissipation function ( ) h s is commonly employed.Each new pair adds its specific contribution to the behavior of the dynamic property curves.
Thus, the complex Young modulus in Laplace domain is expressed by where k is the number of dissipation pole pairs or, simply, GHM terms.
In the frequency domain, it gives However, as new GHM terms are added, new problems also appear, since the terms are not orthogonal functions, and it is very difficult to obtain initial estimates for the additional parameters.To circumvent these problems, Gibson and McTavish (1995) introduced a parameter reduction, relating parameters  ̂ and  ̂ by where 0 < u < 1.
Also, the values of ˆk  are restricted to the geometric series where r is the geometric growth ratio of  ̂.
Given the modifications introduced by Eqs ( 31) and ( 32), instead of (3k+1) parameters, the expanded GHM model contains (k+4) parameters.It should be stressed that the above discussion refers to a condition of constant temperature.The temperature dependence can be introduced, as explained in a previous section.

Alternative GHM (AGHM) Model
An alternative form of the GHM model was presented by Friswell et al. (1997), in which the restriction  � =  � ̂, established by Golla and Hughes (1985) and pointed out before, was no longer necessary.Deriving again the expression for the complex Young modulus in the frequency domain, without the above restriction, it follows that (see Eqs. ( 25) and ( 26)) Rearranging terms in Eq. ( 33) and introducing new parameters produces where This form can also be expanded to contain additional terms, in a similar fashion to that previously presented for Eq. ( 30).As for the temperature dependence, it can be also included in Eq. ( 34) to provide

Modified GHM (MGHM) Model
A modified GHM model was originally introduced by Martin (2011) and further reported by Martin and Inman (2013), resorting again to the restriction  � =  � ̂ and multiplying parameter  ̂ by an auxiliary parameter  � , either in the numerator or the denominator of Eq. ( 27), to create a novel parameter  � =  �  ̂.Choosing to introduce this parameter into the numerator, the result is Latin American Journal of Solids and Structures, 2019, 16(2), e164 8/19 As before, the MGHM model can also be expanded to contain additional terms, as shown in Eq. ( 30).Regarding the temperature dependence, it can be included in Eq. ( 36) to provide As it is of concern to the current paper, the complex Young modulus equation for the MGHM model is presented below including the temperature dependence and additional terms.It is

INFLUENCE OF MODEL PARAMETERS
The parameters of Eq. ( 14), for the fractional derivative model; of Eq. ( 27), for the standard GHM model; of Eq. ( 34), for the alternative GHM model; and of Eq. ( 36), for the modified GHM model can be related to typical features of the master curves.These features, shown in a pictorial fashion in Fig. (3), are the asymptotic behavior of the real modulus at lower and upper frequencies, the slope of the real modulus in the transition region, and the location of a fixed maximum loss factor.Such relationships highlight the influence of those parameters on the curves and can be accomplished by analyzing, in each case, the corresponding equation regarding the feature of concern (Bagley and Torvik, 1986;Lopes, 1998;Medeiros Jr. et al., 2011).Table 1 compiles the relevant features and the related model parameters, as far as the complex Young modulus is concerned.In that table, t  is the transition frequency (Jones, 2001).Parameter analyses can provide some directions regarding the search regions for the parameter values.The FD and the SGHM models have well defined search regions for their parameters.As for the AGHM model, it shows a broader search region than the SGHM model, while the MGHM model has a loosely defined search region, although it regains the control over the slope of the real modulus in the transition region, regarding the AGHM model.Furthermore, it is possible to verify that some parameter variations in the GHM models can cause abnormal behaviors on the dynamic curves (Medeiros Jr., 2010;Medeiros Jr. et al., 2011).

EXPERIMENTAL AND CURVE FITTING METHODS
In order to compare the adequacy of the above models, the dynamic properties were raised with the aid of a specific equipment called NETZSCH DMA 242 C. A typical viscoelastic material -commercially known as ISODAMP C-1002, manufactured by EAR Aearo Technologies LLC (3M Group) -was employed.Previous experiences in the characterization of viscoelastic materials, drawn from Jones (1992Jones ( , 2001)), Lopes (1998), Lopes et al. (2004), and Espíndola et al. (2005), led to the selection of 10 frequencies, namely: 0.1Hz, 0.2Hz, 0.5Hz, 1Hz, 2Hz, 5Hz, 10 Hz, 20 Hz, 50Hz, and 100Hz.
The temperatures were also chosen based on previous experiences, on the manufacturer's nomogram, and on preliminary tests of dynamic characterization performed with the complete experimental set-up.These temperatures, in a total of 7, were the following: -40 °C, -20 °C, -10 °C, 0 °C, 10 °C, 20 °C, and 50 °C.The smallest interval around 0 °C was due to the fact that the transition region of the material lies within that range.At each temperature, a prior 30 minute stabilization time was respected.
In order to avoid nonlinear effects, all measurements were made at low strain amplitudes.Fig. 4 shows a sample of viscoelastic material in the corresponding clamping device of the equipment, just above its thermal chamber.More details regarding the experiment can be found in Balbino et al. (2013).
Having obtained the experimental dynamic properties, the relative error to fit the models was defined as the difference between the experimental complex Young modulus ( e E ) and the predicted (computed) complex Young modulus ( E ), divided by the absolute value of the experimental complex Young modulus ( e E ), such that where Ω  was a given frequency, and   was a given temperature.The objective function f(x) for the curve fitting (minimization) process was the overall sum of the products of the relative errors by their corresponding complex conjugates, such that where p and q were the number of frequencies and the number of temperatures, respectively.
Latin American Journal of Solids and Structures, 2019, 16(2), e164 10/19 As defined above, the objective function f(x) employs all the experimental complex Young modulus data (in both frequency and temperature) simultaneously in order to determine the model parameters of the viscoelastic material, including those parameters associated with the shift factor   () (Lopes et al., 2004).Thus, the necessary elements for a complete dynamic characterization of the material are obtained in a single procedure, from which the reduced frequency nomogram can be directly built.
The minimization was performed by the use of a hybrid strategy, which combined, in sequence, a genetic algorithm and a nonlinear optimization method, as implemented in the ga and fmincon proprietary MATLAB® routines, respectively.This choice proved to be particularly valuable when the GHM models were handled.
The design vectors for the FD, SGHM, AGHM and MGHM models were as the Arrhenius equation (Eq.( 8)) was chosen to model the shift factor   (), in line with an observation made by Jones (1992) regarding the investigated viscoelastic material.
The MGHM model was also employed with one and two sets of additional parameters, denoted MGHM2 and MGHM3 models, respectively.The corresponding design vectors were Thus, the minimization problem was to find the design vector x for each model in order to minimize the objective function f(x).Loose lower and upper bounds were also applied, based on the remarks made in section 5.
The curve fitting results were assessed in terms of the final values of the objective function and of those values divided by the number of degrees of freedom (the number of data points minus the number of fitted parameters) (Motulsky and Ransnas, 1987).Plots of experimental-predicted real Young modulus data and of the corresponding relative residuals were also generated, since the real Young modulus is a key dynamic property when calculating the stiffness of vibration isolators and dynamic neutralizers (Espíndola et al., 2010).The relative residuals were computed as the difference between the experimental and the predicted real modulus over the experimental modulus.
The Akaike's Information Criterion (AIC) was also employed in the analysis.Although not widely known, the AIC proved to be a very valuable and consistent tool.Its theoretical basis can be found in Burnham and Anderson (2002).The main concepts and equations -with the corresponding interpretation -are summarized below, following the lines of Motulsky and Christopoulos (2003).
The AIC is defined by where N is the number of data points, SS is the sum of squares (given, in that case, by the value of the objective function), and K is the number of fitted parameters plus one As pointed out by Motulsky and Christopoulos (2003), the AIC takes into account the sum of squares and the number of fitted parameters in evaluating the goodness-of-fit.It is also highlighted that the value of the AIC lies in comparing models, so that what really matters is the difference between AIC values of models of concern.The individual AIC values are computed and the model with the smallest AIC value is the most likely to be correct.
Given its better adequacy, it is recommended that the second-order (corrected) AIC values should always be used.It is calculated by and its use is the same as that explained above.When two particular models A and B are compared in terms of their AICc values, the probability Pr of choosing the correct model is computed by where Δ is the difference between AICc values (smaller minus greater).It can then be stated that there is a (Pr x 100%) probability that model A is correct, and a [(1 -Pr) x 100% ] probability that model B is correct.
Taking the above definitions into account, the probability that one model is correct can be divided by the probability the other model is correct in order to give a quantity called evidence ratio (ER).It is given by 50) and allows one to state, within the limits of the corresponding experimental design, that one model is (ER) times more likely to be correct than the other model.

RESULTS
The results of the current investigation are presented in this section with the corresponding comments.
Figure 5 displays the wicket plot for the experimental data.As shown by Jones (1992Jones ( , 2001)), adequate experimental data should form an "inverted" U curve in a wicket plot, which is satisfactorily seen in that figure.Tables 3 to 8 contain the values of the fitted parameters, regarding the investigated models.Figures 12 to 17 show the associated plots of relative residuals (each relative residual is the difference between the experimental and the predicted values over the experimental value), whereas the histograms for these relative residuals, regarding their relative frequencies, are displayed in Fig. 18.It can be observed that there is a general trend for the SGHM, the AGHM, and the MGHM models to underestimate the real Young modulus.As additional terms are included in the MGHM model, the relative residuals become more balanced, such as the case of the FD model.Table 9 contains the AICc values for the models of concern.It is observed that the FD model has the lowest value among all the models and that, when compared to the SGHM and the AGHM models, the MGHM model has the lowest value.As additional terms are considered in the MGHM model, the AICc value becomes progressively lower.These observations are fully in line with the results in Table 2 and Figs  As for the probabilities, they are displayed in Table 10, in which each model on the left is compared to a distinct model on the top.From these results, which follow from those in Table 9, it can be stated that the FD model is the most likely to be correct, with 100% probability in all the pertinent comparisons.The MGHM model, on its turn, is more likely to be correct than the SGHM and the AGHM models, also with 100% probability.
As additional terms are included in the MGHM model, a progressive improvement is clearly observed.In fact, it can be stated that, in the context of a wide frequency and temperature dynamic characterization, not only is the inclusion of additional terms desirable but highly recommended.

CONCLUSIONS
A comparative study involving fractional derivative and GHM models was performed, regarding a broad dynamic characterization of a typical viscoelastic material, in terms of both frequency and temperature.Distinctive use of the Akaike's Information Criterion was made.
For a wide frequency and temperature representation of the viscoelastic material, it is shown that the fractional derivative model is the most likely to be correct.It stands out even when additional terms are considered in the MGHM model, which is the most likely to be correct among its counterparts.Additional terms are highly recommended when the MGHM is to be employed in the current context.

Figure 1 :
Figure 1: Experimental plots of dynamic properties.

Figure 3 :
Figure 3: Typical features of dynamic property curves.
Parameter values -MGHM2 model (one set of additional parameters).
Parameter values -MGHM3 model (two sets of additional parameters).

Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :
Figures 6 to 11 portrait the plots of experimental-predicted real Young modulus data.As pointed out previously, the real modulus is a key dynamic property for designing viscoelastic devices for vibration control.It can be visually observed that the FD model is clearly the best in fitting the experimental data, when compared to the SGHM, the AGHM, and the MGHM models.When additional terms are included in the MGHM model, the corresponding plots show progressive improvement (as possibly expected), but the FD model still seems to be the best.

Figure 18 :
Figure 18: Histograms of relative residuals for real modulus data.

Table 1 :
Typical features and model parameters.

Table 2 :
Final values of the objective function.