Identifying Mechanical Properties of Viscoelastic Materials in Time Domain Using the Fractional Zener Model

The present paper aims at presenting a methodology for characterizing viscoelastic materials in time domain, taking into account the fractional Zener constitutive model and the influence of temperature through Williams, Landel, and Ferry’s model. To that effect, a set of points obtained experimentally through uniaxial tensile tests with different constant strain rates is considered. The approach is based on the minimization of the quadratic relative distance between the experimental stress-strain curves and the corresponding ones given by the theoretical model. In order to avoid the local minima in the process of optimization, a hybrid technique based on genetic algorithms and non-linear programming techniques is used. The methodology is applied in the characterization of two different commercial viscoelastic materials. The results indicate that the proposed methodology is effective in identifying thermorheologically simple viscoelastic materials.


INTRODUCTION
In recent decades, we have witnessed a huge increase in the use of polymer materials in applications within structural engineering. The good acceptance of mechanical components made from this type of material is, to a certain extent, due to the fact that they can be easily molded into basically all shapes and also to its excellent general performance in corrosive environments, and to its being a mechanical energy dissipater.
Polymers are highly viscoelastic materials (VEMs), and their mechanical properties vary greatly depending on the load application speed. This is the reason why it is important to have access to models that make it possible to predict with accuracy the mechanical behavior of this type of viscoelastic material during its life time or to predict its behavior in a wide range of frequencies.
The mathematical modeling of the mechanical behavior of VEMs must therefore be performed not only in frequency domain but also in time domain. Generally speaking, regarding the frequency domain approach, the idea is to impose a harmonic loading with variable frequency and known range, and use some dynamic model of the system in order to obtain the storage modulus and the loss factor in a predetermined frequency interval (Moreira et al., 2010;Martinez-Agirre and Elejabarrieta, 2010). Regarding time domain, the idea is to impose a stress history in order to obtain a strain history or, reversely, impose a strain history and obtain a stress history and, subsequently, apply a model of the system (Liu and Subhash, 2006;Plaseied and Fatemi, 2008;Chae et al., 2010).
One way of mathematically representing the models of linear viscoelastic behavior is through integer order differential equations. The main advantage of using such tool lies in its ease of storing and handling the terms, almost always trivial ones, in their derivatives and integrals. A disadvantage is that, in order to represent the behavior of the material properly, a high number of terms is usually needed, generating a higher computational cost (Brinson and Brinson, 2008).
Another way of representing such materials is by means of models based on fractional calculus (FC). These models prove to be attractive for describing the viscoelastic behavior of many materials when added to engineering structures aiming at increasing the amount of dissipated energy (Bagley and Torvik, 1983). Models based on FC present advantages over the classical linear viscoelasticity (based on integer order derivatives) since a smaller number of terms is required for a good representation of the material's behavior -four or five parameters, for example -which, a priori, render calculation and numerical simulation easier (Pritz, 1996).
An improvement on the classical models of linear viscoelasticity may occur when replacing Newton's viscous dampers by Scott-Blair's fractional dampers (Mainardi, 2010). Therewith, differential equations of the mechanical type written in terms of integer order derivatives are replaced by equations involving non-integer order differentiations (not necessarily of a fractional order). Among several possible models (Mainardi, 2010), the fractional Zener model (FZ) proves to be very effective in predicting material behavior in a large frequency band (Pritz, 2003).
During the process of characterizing VEMs, the influence of temperature in the behavior of those materials becomes clear. In this case, a material can be defined as thermorheologically simple when all its relaxation times are affected by temperature in the same way, thus allowing the application of the Time-Temperature Superposition Principle (TTSP).
When applying the TTSP at a given temperature, the viscoelastic behavior may be associated to a behavior at another temperature through a simple shift in the time scale. Such shift allows predicting the material properties over a long period of time, from data collected in tests performed in comparatively short periods of time.
The correspondence between time and temperature and its effects in viscoelastic behaviors may be mathematically represented by some models. According to Tschoegl et al. (2002), the models that stand out are: the Williams-Landel-Ferry model (WLF) (Willians et al., 1955), the Bueche model, the Bestul-Chang model, the Goldstein model, and the Adam-Gibbs model. Except for Bueche's model, the equation forms of all the other time-temperature superposition models are similar to that of the WLF model or reducible to it, differing only in what concerns the interpretation of the equation parameters (Tschoegl et al., 2002).
Specifically in VEMs identification processes in time domain, some works are worth mentioning. Welch et al. (1999) analyzes the quasi-static viscoelastic response of polymeric materials using constitutive models based on FC in time domain. The author uses fractional models of five and seven parameters for the stress relaxation functions, which are capable of modeling materials in which relaxation data manifest at glass level. Park (2001) characterizes a viscoelastic fluid by comparing experimental data to the generalized Maxwell model and to the fractional derivative model. The analysis is performed both in time and frequency domains. Jiménez et al. (2002) report a stress relaxation study in samples of PMMA and PTFE polymers (Methylmethacrylate and Polytetrafluorethylene) indicating that there is not only a relaxation time, as predicted in Maxwell's classic model, but two relaxation time distributions. This is approached using the fractional Maxwell model. Moschen (2006) performs a numerical simulation of the approximate method of the inverse Laplace transform and the deconvolution method, which are used in the numerical calculation of the stress-relaxation function. The author uses the fractional standard solid model (or fractional Zener), with four parameters, applied to VEMs. Pacheco (2015) presents a methodology for characterizing viscoelastic materials through an inverse problem of identification. To that end, the authors uses experimental data in time domain, extracted from stress-strain curves at different temperatures and strain rates. The methodology developed allows characterizing -in time domain -models of VEMs, based on Prony series, considering the influence of temperature and pressure over the mechanical behavior of those materials.
Thus, the main aim of the present work is to propose and implement a methodology for identifying mechanical properties of VEMs in time domain, using the fractional Zener model and taking into account the influence of temperature via WLF model. The methodology is based on the search for the material parameters that minimize the quadratic relative distance between the stress values experimentally measured in tests with constant strain rates and the values obtained from the constitutive model.
The present article includes a short theoretical review of concepts that prove to be important for the development of the proposed methodology (integrals and fractional derivatives, Mittag-Leffler functions, fractional Zener constitutive model, and temperature influence via WLF model). Subsequently, the formulation of the identification problem and the computational structure used are presented. The validation of the proposed methodology is performed through pseudo-tests of a real VEM (EAR C-2003), which was identified in frequency domain and, through a time-frequency interconversion process, the material parameters were obtained in time domain. Next, the results obtained in the characterization of the material CYCOLAC BDT5510 by the proposed model are presented. After that, a graphic comparison of the modulus relaxation function of that material with the results obtained through a model based on Prony series is presented. The article closes with the final conclusions.

Latin American Journal of Solids and Structures 14 (2017) 131-152
Considering that, in the present study, the structural system is initially found at rest, there is no need to treat information occurring for a time t < 0 . So, Riemann-Liouville's definitions (Mainardi, 2010) prove to be the most adequate for the present work. For this purpose, is consid- , a is a positive real number and the Euler gamma function G (Li and Zeng, 2015). A Riemann-Liouville fractional integral to the left of order a is defined by (1) and the Riemann-Liouville fractional integral to the right of order a is defined by (2) Hence, the left and right Riemann-Liouville fractional derivatives with order a and n being an integer positive value such that n n a -£ £ 1 , are defined as respectively.

Mittag-Leffler Functions
The interest in functions of the Mittag-Leffler (ML) type arises from the close connection between those functions with differential and integral equations of fractional order, and integral equations of the Abel type . In this case, as the nature of the exponential function is the solution of differential equations of integer order and constant coefficients, the ML function has an analogous role in the solution of differential equations of non-integer order.
The ML function of order n , denoted by ( ) E z n and with n > 0 , is defined by a series representation, convergent over the whole complex plane, as (Haubold et al., 2009) According to Mainardi (2010), for the convergence of the power series in Eq. 5, the parameter n may be complex, provided An extension to the ML function is obtained by replacing the additive constant 1 of the Euler gamma function argument (Eq. 5) by an arbitrary complex parameter m . Thus, the representation of the ML function with two parameters, n and m , takes the form In the present work, the ML function is of the ( )  , n Î Ì 0 1 order, and argument z is always real and negative. In this case, this argument is conveniently redefined as a t z where  a t Î is a positive fixed parameter and characteristic of the analyzed system. Thus, from Eqs. 5 and 6, one obtains Some important limiting properties of the ML function for a parameter may be highlighted from Eq. 8. Analyzing asymptotic behaviors for a very short time and for a very long time, one obtains , the ML function presents a very fast decrease to ( ) + a t 0 t  , and when ( ) a t t  +¥ , a very slow decrease occurs (Mainardi, 2010).
There are several integral relations associated to the ML functions that can be established by applying the formulae of the Gamma and Beta functions and other techniques. Here, the following integral associated to the ML function (Shukla and Prajapati, 2007) stands out: Considering the special case in which b = 1 e w = -1 , Eq. 11 leads to is the ML function with two parameters (Eq. 6).
Keeping in mind its application to a VEM constitutive model, Eq. 12 is rearranged considering the upper limit variable of integration t as , where k t is any given time. Thus, perform- Equating this result with the R.H.S. of Eq. 12, one obtains an important result for the present work:

Fractional Zener Constitutive Model for VEMs
Constitutive models are used as idealizations of the mechanical behavior of materials. Modifications of classical models of viscoelasticity may be obtained, for example, by proposing a relationship between stress and strain in the material based on non-integer order derivatives (often named 'fractional derivatives' in the literature). Integer order derivatives depend only on the local behavior of the function, whereas fractional derivatives depend on the whole history of that function. This becomes evident in Eqs. 1-4. Among several possibilities of constitutive models, one can highlight here the fractional Zener model (FZ) ( Fig. 1) which, in frequency domain, proved to be efficient in predicting VEMs' behavior (Pritz, 2003). Considering a VEM element based on the mechanical model in Fig. 1, the differential equation that relates normal (or shear) stress, ( ) for one-dimensional problems is given by (Mainardi and Spada, 2011) where 1 R and R 2 are the stiffness modulus of the elastic elements, h is the pertinent viscosity coefficient and n is the differentiation order of Scott-Blair's damper (Mainardi, 2010).
One of the classical tests for characterizing a VEM in time domain is the relaxation test in which the material is subjected to a constant strain e 0 and the stress gradually decreases with time.
For a viscoelastic solid material, stress decreases until it reaches a point of equilibrium, whereas for viscoelastic fluids, stress tends to zero. In the present test, there is a rearrangement of the structure of the material under constant strain, dissipating part of the inner energy as heat, while part of it can be recovered when stress is removed.
In the relaxation test, a constant deformation e 0 is applied, and stress ( ) The relaxation modulus function of the FZ model can be obtained by solving Eq. 16, considering ( ) 0 t e e = constant, and is given by (Mainardi and Spada, 2011) ( ) Using the asymptotic expansion presented in Eq. 10, two important material parameters of that VEM stand out. They are: the instantaneous relaxation modulus, E 0 , given by VEMs can also be defined as materials with memory. In other words, the state of stress at a given point of the material, for a certain moment in time, does not depend only on the state of Latin American Journal of Solids and Structures 14 (2017) 131-152 strain at that same moment, but also on the whole history of strains to which that material point has been subjected. The reverse, on the other hand, is also true. This relationship can be described using Boltzmann's Superposition Principle (Christensen, 1982). In this case, for a specific test of uniaxial stress or pure shear, the material behavior can be represented taking into account that, at each instant, the strain is increased by e D in intervals of time x D . Thus, the expression obtained in terms of the relaxation modulus function ( ) x is the strain rate at each instant x and 0 t x £ < .
In a specific case where the material is loaded at a constant strain rate, that is, the stress represented by Eq. 20 can be rewritten as Replacing Eq. 17 in Eq. 22, and using the integral shown in Eq. 14, one has the stress in a viscoelastic solid described according to FZ model, at any k t instant, strained at a constant rate  e placed as ( )

Temperature Influence on Mechanical Behavior
It is consensual in the literature (Brinson and Brinson, 2008;Lakes, 2004) that the strong influence of temperature onto the viscoelastic behavior of polymers can be described by the Time-Temperature Superposition Principle (TTSP), which proposes a change in the time scale of viscoelastic response. In this case, the material can be defined as thermorheologically simple. This means that each record of the relaxation process (each numerical value of relaxation stress) at a reference temperature T 0 , can be shifted according to (Lakes and Capodagli, 2008) ( ) ( ) where ( ) T T a is a time-temperature shift factor that makes it possible to obtain the experimental data in very large time scales (very long test periods) or very short ones. A mathematical model for the time-temperature shifting factor -widely used in other workswas proposed by Williams et al. (1955), here referred to as the WLF model, and is given by where C 1 and C 2 are constant specific characteristics of each material, and T 0 represents the reference temperature. Although such equation is recommended for materials operating above the glass transition region (Tschoegl et al., 2002), in this work it is used without such restriction in the present work.

Formulation of the Identification Problem
The main aim of the present work is to formulate and computationally implement an identification process of the mechanical properties of a VEM in time domain considering the FZ constitutive model and the influence of temperature by WLF model. To that end, a family of points experimentally obtained through a set of traction uniaxial tests (or pure shear) at different temperatures, with different constant strain rates and zero initial strain is analyzed.
For that purpose, one establishes the total number of temperatures as NTemp, the total number of strain rates as NSr, and the number of points sampled from each curve as NPts.
Starting from Eq. 26 and using the experimental data, it is possible to formulate an adequate standard optimization problem for obtaining the properties of VEM. In this case, a relative distance function is defined, ( ) At this point, the aim is to write the problem of identifying the material parameters as an optimization problem. For this purpose, the objective function is defined through a measure of the Latin American Journal of Solids and Structures 14 (2017) 131-152 quadratic relative distance between the curves of the FZ theoretical model and those obtained experimentally. According to Rektorys (1980), the distance between those two functions can be defined through the norm of the difference between them. The present work chose to use the L 2 -norm of that difference, producing Furthermore, when several curves are simultaneously adjusted, the total quadratic relative distance, D 2 , is divided by the total number of curves, NC, producing As discussed above, the present aim is to adjust a theoretical model of material behavior (in this case, the FZ model) to a set of experimental curves. A scalar value providing a measure of the relative distance between those curves (the experimental and theoretical ones) is given by Eq. 36, and this is defined as the objective function of the optimization problem. Thus, it is possible to explicit the standard optimization problem as follows:

Implementation of the Mittag-Leffler Function
In Eq. 26, it is possible to observe that the argument of the ML function is given by In other words, having defined the material parameters, a t and n , and the temperature and therefore T a , the parameter c can be obtained for any time t . However, the evaluation of this function is significantly costly, for it is produced after obtaining a summation (Eq. 9), which is strongly dependent on the desired accuracy. In order to reduce the time spent on this operation, the present work chose to perform an approximation of the ML function values in the interval in which all experimental curves are found. In the optimization process, for a given configuration, argument c of the ML function is dependent just on time t and on factor T a (which depends on the temperature of each curve). In this sense,  Gorenflo et al ( , 2003. Such points and the respective values of the ML function are tabulated and stored.
Once the minimization process starts, for the ic-th curve, the k-th point must be evaluated via Eq. 26. Therefore, keeping in mind the argument of the ML function for the k-th point, one can resort to the previously obtained table and, having located the interval to which the argument belongs, one obtains the corresponding value of the ML function for the argument evaluated through a linear interpolation. By performing such procedure, one observes a significant decrease in the computational time spent calculating the Mittag-Leffler function every time this recurrence becomes necessary.

Computational Implementation
Initial results of the identification process based only on non-linear programming techniques (NLP) indicated the existence of several local minima points. Thus, obtaining the optimal properties strongly depends on the initially arbitrated point in the optimization process. In order to avoid such problem, the identification methodology proposed in the present work is anchored on a hybrid optimization technique. Such computational process, showed in Fig. 2, was implemented in a set of routines in the MatLab ® software. First of all, the experimental data are read, followed by the definition of the parameters for executing the routine through genetic algorithms (GA), and the respective optimization process. Subsequently, the best point obtained from the optimization process by GA is used as entry datum for the NLP routine and its execution. As the optimization technique by GA is not deterministic, and having in mind the aim of increasing the probability of finding the global optimum point, the hybrid optimization process (GA and NLP) is performed 10 times, which is an arbitrated number. The optima material parameters are taken considering the set of values that produces the smallest quadratic relative distance D 2 . In this optimization process, functions ga.m and fmincon.m (belonging to the MatLab ® toolbox) were used.
The parameters used in the optimization process by GA are: population size is 100; number of generations is 1000; stop tolerance is 1.0E-8 (the algorithm stops when the average relative shift along generations is smaller or equal to the set value); and mutation rate is 0.02. The optimization process by NLP is based on the quadratic sequential programming method where the Hessian matrix of the objective function is updated through the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method. The optimization parameters by NLP are: maximal number of evaluations of the function is 1000, maximal number of iterations allowed is 600, and stop tolerance is 1.

Validation of the Proposed Approach
For validating the proposed methodology, a reference material named EAR C-2003 was used. According to the manufacturer (EAR-Aearo Technologies), such material belongs to the ISODAMP vinyl and ISOLOSS urethane family. No further information is given about the composition of the material, except for the nomogram and its mechanical properties in frequency domain.  From the nomogram provided by the manufacturer (Fig. 3) several curves are constructed, with constant temperatures for each curve. Having those curves, an adjustment is made by the ordinary least squares method in order to obtain the material parameters in frequency domain. In this case, considering the constitutive model of fractional Zener, the expression for the complex Young modulus presented by Pritz (2003) is

Latin
where R W is the reduced frequency (given by Thus, when the material parameters of VEM in frequency domain are known, it is possible to perform an interconversion to time domain and obtain the properties of the model in that space. Therefore, the Fourier Transform is applied to each term of the differential equation of the FZ model (Eq. 16) and, considering the influence of temperature, the constitutive equation of the described model in frequency domain, R W , is given by    -] 8.086e -10 6.030e -10 4.526e -7 5.664e -7 --- Observing Tab. 1, one notices that the identification using all combinations of temperatures and rates presents better convergence of the parameters when compared to the reference material. As expected, the test performed at a fast strain rate ( C e =  0.100 (mm/mm)/s) and cold temperature  Fig. 4. This function considers the history of all temperatures and strain rates, and represents the master curve for a reference temperature (in this case, 41.9 o C) obtained through the optimization process. It is possible to observe that the best result obtained presents very low level of error when compared to the reference material curve, for the curves almost overlap.
At this point, it is worth making an amendment related to the graphic representation of the relaxation modulus presented in the present work.

RESULTS AND DISCUSSIONS
This section presents the results obtained in the process of identifying the CYCOLAC BDT5510 material. This material is a polymer of the Acrylonitrile Butadiene Styrene type (ABS) manufactured by SABIC®. The experiments were performed according to the ASTM D638 norm, considering four different temperatures (T = - and, at each temperature, the material was subjected to traction uniaxial with a fixed strain rate equal to 0.0833 (mm/mm)/s. The difference between this section and the previous one is that, here, the experimental data do not result from 'pseudo-tests'. Fig. 5 presents the experimental data (+) and the corresponding stress-strain deformation curves considering the FZ constitutive model (Eq. 27), obtained from the identification process, for a strain rate of 0.0833 (mm/mm)/s. Tab. 2 presents the results obtained in the identification process, considering the analysis of each temperature individually under a single strain rate, as well as the simultaneous analysis of all temperatures involved. In all analyses, the reference temperature is 23 o C.    Fig. 7 shows the graphical representation of the relaxation modulus functions of the results presented in Tab. 2.
The results obtained previously are confronted with results obtained by a similar identification process, but using another constitutive model. In that case, the methodology used is presented by Pacheco et al. (2015) and refers to VEMs characterization, in time domain, but using the Prony series, via Weichert's constitutive model, which are based on constitutive equations of integer order. In the present model, the relaxation modulus is given by In Fig. 8 are presented the relaxation modulus functions obtained in an identification process for the CYCOLAC BDT5510 material by Weichert model (with 4, 8, and 16 terms) and the fractional Zener model. Analyzing the curves only within the time window available for the test (solid lines), it is possible to note a relative convergence between two material models. That may be due to a little number of tests with a unique strain rate and few temperatures which carries little information about the material behavior and a very short window of time (less than one second).

CONCLUSIONS
The present work proposes a methodology for characterizing VEM in time domain, on the basis of the Fractional Zener constitutive model under the influence of temperature, here considered via WLF model. The methodology developed allows characterizing viscoelastic materials with a simple thermorheologically behavior. To this end, experimental data extracted from stress-strain curves at different strain rates and temperatures were used as a starting point. Having these curves in mind, a reverse identification problem was proposed aiming at reducing the quadratic distance between the experimental stress-strain curves and the corresponding ones given by the theoretical model. The formulation was implemented in the Matlab ® software and a hybrid optimization technique (GA and NLP) was used.
The evaluation of the model using a real VEM (EAR C-2003), but with experimental pseudotests, under the results obtained via simultaneous rates and temperatures, presents a good level of convergence. This becomes evident in the graphical representation of the relaxation modulus functions, when a superposition of curves occurs from the results obtained via the proposed methodology and from the results obtained by the interconversion between frequency and time domains. Furthermore, it was observed that the test performed at a low strain rate and high temperature provides a good convergence for the equilibrium relaxation modulus, E ∞ . Moreover, tests with high strain rate and low temperatures result in good estimates of instantaneous relaxation modulus, E 0 .
Although this approach enables the identification of the relaxation module with few pure tensile tests, it is clear that these tests must be conducted in a suitable and wide range of temperatures and strain rates which enable the acquisition of information for better curve fit throughout your domain. The results indicate that the methodology is adequate for characterizing VEMs in time domain, provided they are thermorheologically simple materials.