Solving viscoelastic problems with a Laplace transform approach supplanted by ARX models, suggesting a way to upgrade Finite Element or spectral codes

Finite Element codes used for solving the mechanical equilibrium equations in transient problems associated to (time-dependent) viscoelastic media generally relies on time-discretized versions of the selected constitutive law. Recent concerns about the use of non-integer differential equations to describe viscoelasticity or well-founded ideas based upon the use of a behavior’s law directly derived from Dynamic Mechanical Analysis experiments in frequency domain, could make the Laplace domain approach particularly attractive if embedded in a time discretized scheme. Based upon the inversion of Laplace transforms, this paper shows that this aim is not only possible but also gives rise to a simple algorithm having good performances in terms of computation times and precision. Such an approach, which fully relies on the Laplace-defined Behavior or Transfer Function (LBTF) can be promoted if it uses AutoRegressive with eXogeneous input parametric models perfectly substitutable to the real LTBF. They avoid the hitherto prohibitive pitfall of having to store all past data in the computer’s memory while maintaining an equal computation precision.


Introduction
To be firmly characterized, the ViscoElastic (VE) behavior of solid materials requires a mathematical model associated to this rheology, as well as a metrological approach to correctly identify the involved physical parameters.For a Representative Elementary Volume (REV) of viscoelastic media, the model relates in a biunivocal correspondence both the stress and the strain applying onto this REV as a function of time.In the case of linear viscoelasticity, these variables are linked through a simple convolution product in time domain which kernel describes the VE behavior (relaxation or retardation depending on whether the strain or the stress are considered for the excitation).Through the convolution theorem, using Laplace or Laplace-Carson transforms of the behavior's model -one should call it the łtransfer functionž of the material-is a standard approach to compute the response (or output) of the material volume to a given time-dependent solicitation (or input) (Tschoegl 1989).This is particularly true since the development of numerical algorithms that perform inverse Laplace transform precisely (Davies and Martin 1979;de Hoog et al. 1982;Dingfelder and Weideman 2014;den Iseger 2006;Stehfest 1970;Talbot 1979).The transfer function formally corresponds to the output obtained for a Dirac delta distribution (pulse) considered as input.Inherent to experiments in solid mechanics, the case of the impulse response in quasi-static conditions1 cannot be obtained as in other scientific fields (the flash method for instance to measure the thermal diffusivity of materials (Jannot and Degiovanni 2018).Therefore, creep, relaxation (Dooling et al. 1997;Emri and Tschoegl 1993) or even indentation tests (Uluutku et al. 2022) are generally preferred to access a retardance or relaxance (Tschoegl 1961), corresponding respectively to an admittance or to an impedance in the electrical analogy.One interest of the Laplace formulation of VE behaviors is for the harmonic excitation which corresponds to the Dynamic Mechanical Analysis (DMA) experimental technique.In that case, the harmonic modelling corresponds to a reduction of the one-sided Laplace transformed problem where the Laplace variable ( notation used in this article) is restricted to the imaginary axis  = , where  is the harmonic pulsation of the input excitation.But although many DMA devices are produced by a few worldwide manufacturers, they can still be hardly used for metrological precise characterizations.Depending on the conducted loading test (bending, torsion, traction, shearing tests), it is still a research concern to obtain a perfect match of DMA-measured mechanical parameters and those reached with quasi-static experiments for other type of excitations (ideally with multisequence loading).The generally small specimen sizes, very weak considered forces and displacements in linear VE regime, add to this difficulty.These experiments are however widely used for qualitative purposes, when discrimination between different materials is looked for, on damping properties typically.DMA is considered to be a very precise tool to investigate phase transitions (Menard and Menard 2020) for example but, in this case, experiments are rather conducted at a given frequency for a sweep in temperature (measurement of a glass transition typically).
Anyway, it would be of great benefit, especially for viscoelastic heterogeneous materials (Agbossou et al. 1993;Shaterzadeh et al. 1998), to make the technique enters in a fruitful interaction with simulation tools (André et al. 2021;Fuentes et al. 2017;Gallican and Brenner 2019;Trumel et al. 2019).The perspective that could be draw, is the following.Assuming that DMA could characterize more efficiently the behavior's law in terms of transfer function, and that this transfer function could be implemented directly in a FEM mechanical code through the behavior's law module, then simulations at the scale of a structure could be directly performed.The problem to circumvent is that numerical codes will necessarily proceeds with a time-step algorithm (an integration scheme in time) which appears incompatible with the use of a Laplace transform defined behavior's law.The present paper shows how this can be achieved with a different philosophy than that implemented in previous works on heterogeneous materials (Lévesque et al. 2007;Rekik and Brenner 2011) or on homogeneous materials by using the classical collocation method (Schapery 1962) to identify the relaxation kernel from its expression in Laplace domain.The approach leads to input-output convolution type series, common for Linear Time Invariant (LTI) systems, which suggests the introduction of AutoRegressive with eXogeneous input (ARX) models, well known in the field of Control Theory and Automation.
In Section 2, we recall the equations and vocabulary associated to the VE problem and define the central behavior's law considered for this study, which can be reduced to either the simple analogical Standard Linear Solid (SLS) behavior nor asymptotically, to a non-integer or fractional (NIF) order operator.In Section 3, the incremental approach is built by making use of a Laplace-defined behavior's law and of its response to a unitary ramp excitation.This solution is shown precise but involves all the history of the mechanical variables which makes it inefficient to use with meshed structures.ARX models' structure is shown here to be very powerful to avoid this drawback.Section 4 illustrates how it works and the achieved performances with the three types of VE behavior considered in this study.

Generic equations
The mathematical problem characterizing the viscoelastic behavior is of convolution nature through the well-known Boltzmann superposition principle (Christensen 1982).The stress for example is given by the convolution product where  () figures the relaxance kernel (Tschoegl 1989) and  () the test function i.e. the excitation (strain loading path) imposed on the system.Alternatively, in the second equality,  () figures

André & Noûs
Solving viscoelastic problems with ARX models the relaxation function or relaxation modulus with  () = ∫  0  () d and  () = d d +  (0) () as a result of  being discontinuous at  = 0 but with a definite limit.The lower integral bound initiates at  = 0 as for any causal problem.This expression suggests that the full-time history dependence of variable  (since the experiment starts) contributes to the actual value of the stress.The VE problem considered here in 1-D is referred to as a SISO model (Single Input-Single Output) with the vocabulary of control theory and referring to the role it plays in the historical development of the parametric models used below (Gevers 2006).To this mathematical formulation corresponds an alternative form in terms of the Ordinary Differential Equation (2), through a linearity property.With  ( ) denoting the -th derivative of function  , one can have the general equation which produces a given model structure for the relaxance kernel  (), readily obtained through Laplace transforming of Equation ( 2).Laplace-Carson transform is in general invoked in the mechanical community but is just a commodity practice to access directly the relaxation modulusÐ or alternatively creep function, considering a permutation of the independent variables  and  in Equation ( 1)Ðwhen a Heaviside function is considered for  (): the measured stress is then the image of the time-dependent relaxation modulus.However, it is not necessary for the reasoning and even blurs the full generality achievable with the transfer function concept.In the Laplace domain, Equation (2) transforms into where function H () describes the Laplace Behavior or Transfer Function (LBTF).Then, from the convolution theorem, and the selection of the excitation function, the solution is produced according to the alternative form of Equation (1).H () can also be defined directly through a specific mathematical function and one can think here to the empirical Davidson-Cole relaxation model or to the more recently promoted models with non-integer exponents leading to fractional differential operators substituted to integer ones in Equation (2) (Bagley and Torvik 1986;Schiessel et al. 1995).This aspect is mentioned here to draw attention of the reader to the large application that can be made from the results presented here.When ε () = 1, H () is named the transfer function of the system.Using a metaphorical expression from biology, the transfer function characterizes the łfull DNAž of the system (material) in terms of frequency information.It corresponds to the response of the system to the impulse (Dirac Delta) function H () = σ () and is obviously never realized experimentally in mechanics due to an impossibility related to inertial effectsÐhowever, it is a quasi-systematic approach for thermal systems excited with laser pulses for example, as exemplified in (Corbin and Turriff 2012).When such a transfer function is known from the engineer, then it allows computing the response of the system to any input excitation.If  () =  (), the multiplicative identity of a convolution product is the realization of the transfer function and then, in any other known specific strain loading path  ℓ (), we have fully determined by the impulse response   ().

André & Noûs
Solving viscoelastic problems with ARX models 2.2 Description of the considered VE behavior's laws and loading sequences For the computations presented later, we will consider as generic viscoelastic transfer function, the one furnished by the Dynamic of Linear Relaxations (DLR) viscoelastic model (Cunat 2001).This model originates from a TIP framework (Thermodynamics of Irreversible Processes) combined to a modal approach of the dissipative processes.To say it shortly, it ends in a finite set of ODEs identical to Equation (2) but where coefficients   and   are in auto-recursive geometric progression, or equivalently, when the Laplace counterpart of Equation (3) exhibits a recursive Pole and Zero Distribution.This model has been shown efficient to describe and measure associated material constants of HDPE, a thermoplastic, semi-crystalline polymer (Blaise et al. 2016).Furthermore, it has the proper mathematical structure used in the collocation method and therefore, Section 4.2 provides some sort of comparison between the latter and the use of ARX models.This model summarizes to • DLR constitutive equation • DLR constitutive equation in Laplace Domain (LBTF) • DLR impulse response in time (Transfer function) with ś   , the unrelaxed (Young or glassy) modulus, ś (  ,   ), the spectrum of relaxation times and weights defined as ensuring  =1   = 1 with the recursion  +1 /  =  and  +1 /  = ; also  max is the material parameter corresponding to the maximum relaxation time,  and  are the number of decades and of relaxation modes of the spectrum, respectively, ś   (), the relaxed state defined simply here as   () =  r  () where  r is the relaxed modulus (rubber modulus for polymers).Starting from this generic VE behavior, three other transfer functions are considered.They all participate to a validation of the incremental approach described next and to the understanding of its underlying subtilities.1.The first one is the Standard Linear Solid (SLS) or 3-parameter Voigt model and corresponds to the reduction of the DLR model to one single relaxation mode.The transfer function and other formulations of the SLS behavior's law are directly obtained from Equation ( 6) by discarding the summation operator.2. The second one is obtained when considering asymptotically the case of an infinite number of modes  → ∞ which we have shown earlier to be related to a specific fractional model (André et al. 2003), i.e. a model based on non-integer derivative operators.It is obtained directly in Laplace domain without any analytical expression of its original and will be described later in Section 4.3.1.This case is considered here in view of the abundant literature available on these NIF (Non-Integer or Fractional) models since the 90's as they allow for the description of a plethoric set of VE behaviors.3. The third one is also a fractional viscoelastic model.Because the previous one stems from the DLR behavior and remains rather confidential, another fractional law will also be quickly investigated in the same Section 4.3.Based on Rabotnov's suggestion in 1948 to use a fraction-exponential operator, it was widely considered since then (Koeller 1984;Sevostianov et al. 2015).

André & Noûs
Solving viscoelastic problems with ARX models Irrespective of the łfractionalž behavior considered, it will next be referred to as a NIF model.Their natural definition in Laplace domain will legitimate even more the tool raised by this study in view of an introduction into numerical codes.
For the simulations presented later, we will also consider a set of different loading paths  ℓ () presented in Table 1 along with their known Laplace domain expression.This set will allow a critical appraisal for the methods described later.Case 0 corresponds to a smooth function in time.Case 1 corresponds to the crenel excitation and Cases 2 to 4 correspond to more complex sequences made of ramps.The associated graphs of these loading paths (inputs on the system) will be shown later along with the calculated responses (outputs).
stab : time at which a dwell is imposed to the strain  = ( +  + ) stab Table 1 Considered test cases for the loading strain path (H () : Heaviside function) 3 Towards an incremental time-step approach of the problem 3.1 The step-approach reached from the unit ramp response FE codes in solid mechanics and more recent FFT-based spectral solvers are based on a timediscretization approach for elasto-visco(-plastic) materials submitted to transient excitation.
The equilibrium equations are solved at specific locations of the spatial discretized domain.
For each increment in time step, a call to the behavior's law computation is required which is always specified in time discretized version, with different possible schemes (Sorvari and Hämäläinen 2010).The objective of the present study is to allow the behavior law of the material be defined directly in the Laplace domain.We describe in this section how to compute precisely and incrementally the temporal response of a viscoelastic material with its Laplace-defined Behavior or Transfer Function (LBTF).By incremental we precisely mean running a loop over a time-discretized vector and calculating the output stress according to  ( + Δ) =  () + Δ  .For the different test-cases considered, this incremental computation referred to as step-LBTF, needs to be validated.To this aim, it will be compared in the whole time interval [0;   ] with: 1. the computation of the response directly in time-domain through Equation ( 5), with the convolution product between the transfer function (impulse response   ())Ðwhen known mathematicallyÐand the input excitation (strain loading path   ()).The Convolution built-in Matlab function is used to establish this output response and the result will be referred to as   *   ) () in the curve labels of the coming plots.2. a full inverse Laplace Transform approach, meaning that the computation will be based on the Laplace product of Equation (3b) and inverted over the whole discretized vector of the time domain under consideration [0;   ].We will use the Laplace inversion based on De Hoog's algorithm (de Hoog et al. 1982).These computations will be labelled in the figures as full-LBTF.
Note that if nothing proscribes the use of De Hoog's algorithm embedded in a time loop to calculate incrementally the output stress, this leads to largely increased CPU times (our tests give a multiplication factor of 50 for a time vector of dimension 7000) and with no significant enhancement of precision.This certainly precludes any application of the Laplace inverse in the incremental schemes in times of FEM or Spectral solver codes.The step-LBTF approach originates from Duhamel's theorem, see (Özisik 1993) or (Sneddon 1951, Theorem 33, page 164) which allows various different form of the generic convolution response according to Equation (1).Going one step further, it is possible to express this convolution in terms of the response to ramp functions.More precisely, we compute the elementary time response to a Unitary Ramp (ramp of unit slope over [0; Δ]) which then remains Stucked at a Stationary value equal to 1 × Δ.This input is described as: We named it hereafter URSS-Response and noted it  0 with subscript 0 referring to the time origin when this excitation takes place.It can be computed once for all in a preliminary stage using for example the Laplace inverse and stored as a vector of the discretized values  0 (  ) over all time steps.De Hoog's algorithm will also be used for this purpose.The loading path is described through the discretized slope variations Δ  =   −   −1 at time   , see Figure 1, with   = ( ℓ ( +1 ) −  ℓ (  ))/Δ.Hence, it is approximated as a piecewise-linear function:  ℓ (  ) =  Δ  (  − ( − 1)Δ).t From the Laplace transform application, it can be shown that By the time translation invariance property, the incremental scheme  (  + Δ) =  (  ) + Δ   follows with Δ   being obviously the cumulative dot product of the backward URSS-response originating recursively from time   with the vector of slope increments up to time   .A first ramp of slope Δ 1 =  1 acts on the VE response over all subsequent time steps.A second ramp driven by the slope change pulse Δ 2 acts on the VE response delayed from one time-step, and so

André & Noûs
Solving viscoelastic problems with ARX models on.This makes the contribution to the current output Δ   to be the URSS-response sequence weighted by the current input Δ +1 and past input Δ  , Δ  −1 , . . .
which is some kind of standard primer to derive an ARX model.As first element of comparison between those three approaches, Figure 2 shows the results in the case where the relaxation spectrum is limited to one mode (SLS rheological element).The considered excitation is made of three ramps (Case 3) and is shown in the same figure (right axis).The agreement is very good but one can notice in the insert of Figure 2 that the step-LBTF approach is even more precise around discontinuities when compared to the full-LBTF.As second element of comparison, we report in Table 2 the indicative CPU times for the methods and relative discrepancies quantified with the following standard quality metrics: RMS-error and fit level (in %) defined as where  and  ref will be defined along the text but refer always to model simulations outputs taken alternatively as the łmeasurementž signal  and the model response signal  ref .
⟨⟩ refers to the mean of signal , when averaged over the whole number of points (time steps) considered.
The Step-LBTF simulation is taken in  here a more refined VE behavior (DLR model) where relaxation runs over a spectrum made of  = 50 modes.The simulation holds for the difficult case of the crenel excitation (Case 1), with Δ = 0.001 and  = 7000 time steps.This time step corresponds to a choice which ensures a negligible error to be made for the material response computed by the three different approaches and for the full set of excitation types.It is also convenient to produce representative CPU times to compare the methods.But it must be mentioned that it has no influence on the results, the 'full' or 'step' Laplace formulations being both based on the Laplace inverse of an exact function, Equation (3b) or Equation ( 9).De Hoog's algorithm behaves equally for any selected sampled time (with the above-mentioned restriction, that it performs badly near strong discontinuities in the signal to restore when the full-time vector is passed to it).This effect can be seen in the insert of Figure 2.
The following comments can be made: • Firstly, the convolution approach based on a Matlab algorithm appears to be the fastest when considering the same time step (row 1) as for the other methods.But this approach never supports the comparison in terms of precision.Generally, time steps more than 100 times smaller (row 2*) are needed to obtain a solution within the same order of magnitude in precision as for the present Step-LBTF approach (row 4).Therefore, in the next sections, results obtained with this convolution routine will be given for Δ = 10 −5 s to serve as a reference for the different approaches.
• Secondly, approaches based on Laplace inversion can be compared when looking at row 3 (full-LBTF) and row 4 (Step-LBTF incremental version).The precision appears of same order, except near discontinuities as mentioned above.For the crenel excitation (Case 1), strong perturbations over a few time steps around the crenel discontinuity are produced.As outliers due to a numerical artifact, they have been excluded from the RMS and fitting error computations so as to reflect more exactly the expected performances reached in the case of smooth variations (like for Case 0).Regarding the CPU times, they also appear of same order but it must be recalled that the incremental scheme offers a computation by more than 50 times faster when compared to the repeated application of DeHoog's algorithm at each time step.Hence the superiority of the present approach, both in terms of precision and CPU times, keeping in mind that these simulations are for the homogeneous case.At this stage, an incremental formulation has been derived from the system input-output relationship given in Equation (11), assuming zero initial conditions.It produces VE behavior computations with good accuracy and small CPU times.In itself, this step-LBTF method has certain advantages.It can be used when the inverse Laplace transform of a relaxation kernel is not known analytically.It also allows to identify the parameters of the behavior's law by inverse method applied to experimental data, for which the łexcitationž signal in deformation would be described as a piecewise linear function.However, used as is, there is little chance that this approach will find a useful application by being implemented in a numerical code for applications to heterogeneous materials, because it formally requires to store in memory all the values of the variables of previous times; the six components of both the stress and strain tensors stored for each spatial node point and for each time step may represent a too large amount of data.
Therefore, based on the literature available in control theory for LTI systems, we found the ARX models an efficient manner to push forward this idea of using an incremental scheme initially based on the present step-LBTF approach but allowing for a drastic reduction of the information to be stored in the computer's memory.

The step approach achieved with ARX parametric models
ARX AutoRegressive models with eXogeneous inputs (Ljung 1998) are parametric models which are shown in this paper to be very suited for being substituted to Laplace transforms and at the same time, providing an incremental approach of time domain computations.When the convolution problem requires the all input history, autoregressive models provide the model output at a given time as a linear combination of the only  previous time steps of the output history and  time steps of the input history.Retardation or delay can be considered using  time steps.The structure of the model is given as where   = Δ for  ⩾ 1,  ARX () = 0 for  ≤ 0 and   () denotes the ideal (strain) excitation imposed on the (material REV) system.In command systems vocabulary, the stress figures the output and   the input or command signal.Coefficients   ,  are the observer Markov

André & Noûs
Solving viscoelastic problems with ARX models parameters.They have to be identified for a selected triplet (, , )Ðthus the common notation ARX(, , ) used hereafterÐduring a so-called calibration experiment.Providing this latter to a specially designed optimization procedure2 allows to obtain the optimal set of parameter vectors   ,   (dimensions and values) in the sense of a minimization of the residuals, while avoiding any over-adjustment.ARX model identifications are performed in 'simulation' mode, that is considering a cost function to minimize which is based on the residuals between the 'measured' and 'model' signals.We must insist again on the fact that these ARX approaches are not strictly speaking physical models describing the viscoelasticity phenomenon in connection with microstructural mechanisms.It is obvious from the mathematical structure of Equation (13), when compared to the VE behavioral models that are investigated in this study, that the ARX structure and parameters encode in a compact form the structure of the VE relaxation kernel, its physical/mechanical parameters and the time step used to produce the calibration experiment.This latter will be made here from a first simulation obtained in a direct way for a given strain excitation  ℓ (selected from the test cases referenced in Table 1).Use was made of the step-LBTF approach to this aim.
The second phase is a validation experiment which aims at quantifying the robustness of the identified ARX model i.e. its intrinsic character.This means ruling on the fact that parameters   and   are now able, through the ARX(, , ) model structure, to reproduce the real behavior in all situations.We will use the available test cases responses for various input excitations as data sets for the validation experiment.
It should be noted here that this ARX approach is deployed for a computing strategy but can be used by experimentalists.In that case a calibration experiment really obtained from a metrological set-up is necessary and the effect of the measurement noise has to be considered.
Practically, the identification of the observer Markov parameters (polynomial structure of Equation ( 13)) will allow the computation of the system Markov parameters which correspond to a classical state-space representation of the system in control theory.This means that, rather than Equation (13), we implemented the first-order matrix difference equation with  () the state of the system,  () the output stress observable,  () the input strain.Matrix  is  × ,  is  × 1 because only one input signal is considered,  is 1 ×  because a single output signal is considered and  is 1×1.These matrices will be fully determined from identified observer Markov parameters and directly given with Matlab System Identification Toolbox.The formulation of Equation ( 14) is preferred in our codes because Equation (13) leads to a reconstruction of the validation experiment which depends from a deadbeat observer gain (not given by the Matlab System Identification Toolbox) and which is inherent to the compression mechanism that limits this input-output model to a small number of parameters (Phan and Longman 1996;Wu 2022).Additionally, we observe that this formulation produces smaller computation times than with Equation (13).Although both are of the same order, ARX formulations reduce them by a factor 10 compared with the Step-LBTF approach, see last row of Table 2.

Simulation results on three VE behavior
In this section, we validate the step approach reached by substituting ARX models to the Laplace inversion of behavior's law formalized in Laplace domain.Three VE behaviors will be considered to assess the methodology which correspond to the three previously introduced behaviors: SLS, DLR and NIF models.

SLS model
The response to the SLS model with parameters   = 1000 MPa,  r = 100 MPa,  SLS = 1 s is computed using Full-LBTF, convolution and Step-LBTF approaches.Figure 2 illustrates how well these approaches compare in test case 3 (three successive ramps).We use then the Step-LBTF approach as the reference output  ref to identify an ARX model.The basic ARX(1,3,0) is shown in Table 3 to produce very high level of precision.Of course, for a given calibration experiment (rows of Table 3), the highest performances are obtained when the case considered for the validation experiment (columns) is the same (diagonal of Table 3).It is also clear that case 4 used as calibration experiment (Calib4SLS) offers the best matching conditions for all other cases considered for the validation.Hence the identified ARX model used for the incremental strategy will be retained from this test case.For the validation test case Nr3 (Valid3SLS), Figure 2 shows that the ARX model is perfectly able to compute the expected response incrementally using only three historical values of the input and Table 2 reports an excellent fitting performance of 99.987 %.Indeed, the RMS error appears more than ten times smaller than with the convolution product approach (row 2) for a CPU times reduces by a factor larger than 3000.The same applies for all other validation cases.

DLR model
To go beyond the simple SLS case, we consider now the more complex DLR viscoelastic behavior.It relies on Equation (6a) with a relaxation spectrum described over  = 6 decades, with a discrete set of  = 50 relaxation modes and parameters defined like for the SLS simulation:   = 1000 MPa,  r = 100 MPa,  max = 1 s.The Step-LBTF approach is chosen again to identify the two ARX(1,3,0) and ARX(2,4,0) models considered here.Table 4 gives the fitting error obtained after identification of ARX(2,4,0) model.It is still clear that Test Case 4 offers a signal with sufficiently broad information in frequency domain to cover correctly other types of excitation.Table 4 also illustrates that for this VE behavior, fitting performances, if still excellent, decrease when compared to the simple SLS behavior.That is also the reason why the simple ARX(1,3,0) structure now fails in describing accurately the response in various other situations.Figure 3 compares all approaches in test case 1: the input excitation is a crenel of duration   = 3.2 s and amplitude  0 = 1 (Right axis).Figure 4  The sudden input step induces an instantaneous response.A relaxation stage is then observed and fully completed here because   >  max before it restarts oppositely at   when the applied excitation ceases.It can be observed, especially at discontinuities, that for this more complex Table 5 Comparison of LBTF and ARX approaches for different excitations with Δ = 10 −3 s except for the convolution Δ = 10 −5 s.DLR behavior's law.behavior's law involving a spectrum of relaxation times, the ARX(2,4,0) has superior abilities than ARX(1,3,0).This is clearer from the values of errors reported in Table 4 and the graph of residuals in Figure 4.Because the ARX mathematical structure does not change and involves very few time steps, the CPU times are of same order, see Table 5, much less than other approaches, especially the Step-LBTF which requires the full history of the input command.For the smooth excitation (case 0), Figure 5 reports the stress output for all approaches and Table 5 reports the achieved performances.In that case, both ARX identified models give similar performances.The RMS errors of 2.6 and 1.5 MPa for respectively the ARX(1,3,0) and ARX(2,4,0), when compared to the maximum signal amplitude of 250 MPa correspond to about 1 % and 0.6 %.They were 16.5 and 2.5 % for the crenel excitation.ARX(2,4,0) can then be favored in case of expected discontinuities in the strain input path, for nearly equal and small CPU times in any case.The residuals between the reference Step-LBTF signal and all other approaches (not shown here) behave similarly as those shown in Figure 4 for the crenel excitation, except where discontinuities take place.Through these results obtained for the DLR viscoelasticity model, we compare in a way the ability of this ARX approach to substitute for a collocation method.Of course, if the true behavior of a material is perfectly described using collocation, then the incremental scheme for numerical computation follows.Regardless of the CPU time issue, to be considered only from the perspective of heterogeneous viscoelastic microstructures, it will produce very accurate simulations.For some materials, other strategies to model viscoelasticity can be favoured such as the use of fractional derivatives (Friedrich et al. 1999) and lead naturally to consider such non integer or fractional (NIF) models in the next section.

NIF models
Finally, the so-called fractional models are considered i.e. models involving a frequency LBTF with non-integer exponents in terms of the Laplace (or frequency) variable.These later imply a non-integer differentiation operator in temporal domain.We will consider two of them, having distinct conceptual basis.

Non-integer Oustaloup model
The first model was proposed initially by Sabatier et al. (2015) to represent recursive schemes of the impulse response of a LTI system whose mathematical structure is made of Prony series with recursive factors.We proved that such a model can be related to the DLR approach in the asymptotic limit of infinite modes of relaxation (André et al. 2003).The formal expression of the behavior law in Laplace domain is known directly as with  high = 1/ min ,  low = 1/ max defining the bounds of a spectrum of relaxation times covered by the Laplace variable range.The values   = 1000 MPa and  r = 100 MPa are the same as those taken for the previous SLS and DLR behaviors.The relaxation component of the transfer function Ȳ depends on the last three parameters appearing in Equation (15).We consider the case of  = 0.62,  min = 0.01 s, and  max = 16 s.

Solving viscoelastic problems with ARX models
This transfer function is rather elaborated, see (André et al. 2003) for details.In the limit  high → ∞, it reduces to the well-known Davidson-Cole approach in frequency domain.In temporal domain, this transfer function implies non-integer derivatives.Therefore, the simulation presented in this section with such a transfer function and its łreductionž to ARX models will illustrate how non-integer models of viscoelasticity, formulated in Laplace domain, can be used easily in incremental time domain computations.Tables 6 and 7 give the  fit errors obtained for the different possible combinations of the set of calibration and validation simulations for both ARX(1,3,0) and ARX(2,4,0) respectively.The ARX(2,4,0) model appears of course superior to the ARX(1,3,0).It is seen when comparing the results of line 2 versus column 2 in both tables that it   Step-LBTF.The residuals RMS error given in the 6th and last sub-line corresponds to test case 1 and 4 respectively and can be compared by looking at the signal response variation shown in Figures 6 and 7. CPU times for ARX models are greatly diminished, by an order of magnitude of 10 in average between the Full-LBTF and ARX(2,4,0) model.Looking at lines 3 and 4 (same test case 4 but calibration of ARX models from case 1 or 4) show again that the calibration performed with the reference crenel excitation gives results as good as for the proper test case.

Fractional Rabotnov model
The second considered model relies on the use of a relaxation kernel noted ∋ *  (, ) which was proposed by Rabotnov and Scott-Blair (Sevostianov et al. 2015).Like the previous one, it is based on a very tractable analytical expression in Laplace domain which leads to the LBTF It does not correspond to the same mathematical structure as given by Equation ( 15).In the above expression,  denotes the non-integer exponent,  the inverse of a (relaxation) time, and the Rabotnov kernel is defined in temporal domain by   dynamics like relaxation phenomena in disordered systems under various physical approaches of 'slow-dynamics' type.For example in (Weron and Kotulski 1996) the Fractal Time Random Walk model of dynamical phenomena in complex condensed matter systems is shown to produce the Mittag-Leffler function, as the inverse Fourier transform of the Cole-Cole function.When random walks are used to approach anomalous diffusion or relaxation in complex systems, the introduction of statistically distributed waiting times before each walk step has been also shown to produce dynamics described by ML functions (Metzler and Klafter 2000).They appear as the formal solution of Fokker-Planck equations used in a subdiffusion (i.e.non Brownian) phenomenon (Metzler and Klafter 2002).
Concluding the study with this example presents the interest of being able to produce the temporal expression of the impulse response and then a convolution product just as it was possible for the SLS and DLR models but not with the NIF 'Oustaloup' model.We simply have in that case for the 'Rabotnov' impulse response: As a recall, this leads for any loading path  ℓ () to compute the convolution product  () = (  *  ℓ ) ().Within a full identical strategy, we just show in Figure 9 the performance of ARX(1,3,0) model in this new case of wide interest in the field of relaxation phenomena.ARX(1,3,0) parameters were identified from the simulated outputs to  ℓ 4 (Case 4 of the multisequence loading path).Simulations were performed for the same values of   = 1000 MPa,  r = 100 MPa as before, with  = 1 s −1 and the four  values 0.2, 0.4, 0.6, 0.8 considered in Figure 8.The identified ARX(1,3,0) models gave respectively a fit error of 98.05 %, 98.3 %, 99 % and 99.58 % and were then used to compute the output to  ℓ 3 .Figure 9 reports the four sets of curves obtained to compare the Step-LBTF and the ARX responses.Even with this low parameterized ARX structure, all approaches (Full-LBTF, Convolution, Step-LBTF, ARX) compare again very well in all cases (error fit above 99.99 %).

Conclusion and perspectives
We presented a strategy for solving the quasi-static mechanical problem under any kind of excitation for a linear viscoelastic material which constitutive law is prescribed in Laplace domain but usefully substituted by parametric models.For the simple problem of a homogeneous body, ARX models could provide an interesting alternative to the use of the exact models in the case where very short computation times are required.For example, an UMAT subroutine in Abaqus could work as a generic behavior's law with the ARX parameters being introduced by the user according to the ARX structure he decided to use.ARX models have been shown here to perform remarkably well for all types of responses in a mechanical VE problem, in the linearity domain w.r.t. the input signals.But the more promising interest lies probably in the use of ARX parametric models as proxy models for behavior's laws defined in terms of their Laplace transform (in the frequency domain as produced by DMA for example) that cannot be readily used in incremental schemes of numerical codes.Therefore, it could provide straightforward alternatives in the case of heterogeneous media with viscoelastic phases to produce the homogenized response of the REV or to simulate a complex structure.

Figure 1
Figure 1 Discretization of the input  ℓ in terms of slope variations.

Figure 2
Figure 2 ARX and LBTF approaches compared in the SLS behavior case, with excitation made of three ramps (Case 3).Calibration of ARX(1,3,0) model based on Case 4.

Figure 4
Figure 4 Residuals corresponding to Figure 3. Step-LBTF łsignalž is used as reference for the calibration of ARX models.

Figure 5
Figure 5 ARX and LBTF approaches compared in the DLR behavior case, with excitation made of the smooth excitation  ℓ =  () (Case 0).Calibration of ARX models on Case 4.

Figure 6 (Figure 7 Figure 8
Figure 6 ARX and LBTF approaches compared in the fractional behavior case (NIF-Y n ), with the crenel excitation  ℓ 1
Table 2 as the output reference  ref .We have considered

Table 2
CPU times, RMS error and Fitting error for the different numerical approaches (DLR model, Case 1 excitation)

Table 3
Fitting error for the various test cases (0 to 4) and for the different sets of Calibration/Validation numerical experiments for the SLS behavior's law plots the residuals obtained for the different methods (always considering the Step-LBTF approach for  ref ).
Table 4  fit errors (%) obtained when using the different test cases for calibration and validation experiments to identify the ARX(2,4,0) model substitutable to the DLR behavior's law.

Table 7
Identification of ARX(2,4,0) model:  fit errors (in %) obtained when using the different test cases for calibration and validation experiments.is preferable to consider the test case 1 (the crenel) as calibration test to qualify the ARX model.It covers better the other validation cases although test case 4 is superior if the input signal is free of any discontinuity.Table 8 gathers CPU times,  fit and RMS errors for different test cases with the specified ARX calibration model.The  fit and RMS errors are established with reference taken from the

Table 8
Comparison of LBTF (NIF-Y n ) and ARX approaches for different excitations (time step Δ = 10 −3 s) and ARX calibration models.