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

Stéphane André, Camille Noûs

► To cite this version:

Stéphane André, Camille Noûs. Solving viscoelastic problems with a Laplace transform approach supplanted by ARX models, suggesting a way to upgrade Finite Element or spectral codes. Journal of Theoretical, Computational and Applied Mechanics, 2023, pp.1-18. 10.46298/jtcam.10304. hal-03845394v2

HAL Id: hal-03845394

https://hal.univ-lorraine.fr/hal-03845394v2

Submitted on 10 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.10304

HAL hal-03845394v2

**History**

Received Nov 15, 2022

Accepted Jul 20, 2023

Published Oct 6, 2023

**Associate Editor**

Laurence BRASSART

**Reviewers**

Anonymous

Anonymous

**Open Review**

OAI hal-04174179

**Supplementary Material**

Software
SWHID [SoftwareHeritage]

**Licence**

CC BY 4.0
©The Authors

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

Stéphane ANDRÉ<fnref n="1" />,<fnref n="2" /> and Camille NOUS<fnref n="2" />

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.

**Keywords:** Laplace transform, ARX models, iterative algorithm, viscoelasticity, fractional relaxation kernels

## 1 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 conditions<fnref n="1" /> 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.

<footnote n="1">¹ Université de Lorraine, CNRS, LEMTA, F-54000 Nancy, France</footnote>

<footnote n="2">² Cogitamus Laboratory, F-75005 Paris, France</footnote>

<footnote n="1">¹ If inertial effects are taken into account, VE materials can be characterized through complex moduli measured with the vibrational response of beam-like structures where the excitation can be produced in forms of an impulsion (hammer blow).</footnote>
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 (p notation used in this article) is restricted to the imaginary axis $p = j\omega$, where $\omega$ 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 (Ag-bossou 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.

# 2 Mathematical viscoelastic problem

## 2.1 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

$$ \sigma(t) = \int_{0}^{t} H(t - \tau)\varepsilon(\tau) \, d\tau = \int_{0}^{t} G(t - \tau)\dot{\varepsilon}(\tau) \, d\tau \quad (1) $$

where $H(t)$ figures the relaxance kernel (Tschoegl 1989) and $\varepsilon(t)$ the test function i.e. the excitation (strain loading path) imposed on the system. Alternatively, in the second equality, $G(t)$ figures
the relaxation function or relaxation modulus with $G(t) = \int_{0}^{t} H(\tau) d\tau$ and $H(t) = \frac{dG}{dt} + G(0)\delta(t)$ as a result of $G$ being discontinuous at $t = 0$ but with a definite limit. The lower integral bound initiates at $t = 0$ as for any causal problem. This expression suggests that the full-time history dependence of variable $\epsilon$ (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 $f^{(i)}$ denoting the $i$-th derivative of function $f$, one can have the general equation

$$ \sum_{i=0}^{n} \alpha_{i} \sigma^{(i)}(t) = \sum_{j=0}^{m} \beta_{j} \varepsilon^{(j)}(t) \qquad (2) $$

which produces a given model structure for the relaxance kernel $H(t)$, 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 $\sigma$ and $\epsilon$ in Equation (1)—when a Heaviside function is considered for $\epsilon(t)$: 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

$$ \mathcal{L}\{H(t)\} = \bar{H}(p) = \frac{\bar{\sigma}(p)}{\bar{\varepsilon}(p)} = \frac{\sum_{j=0}^{m} \beta_j p^j}{\sum_{i=0}^{n} \alpha_i p^i} \qquad (3a) $$

where function $\bar{H}(p)$ 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

$$ \sigma(t) = \mathcal{L}^{-1}\{\bar{H}(p)\bar{\varepsilon}(p)\}, \qquad (3b) $$

the alternative form of Equation (1). $\bar{H}(p)$ 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 $\bar{\varepsilon}(p) = 1$, $\bar{H}(p)$ 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 $\bar{H}(p) = \bar{\sigma}^{\delta}(p)$ 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 $\epsilon(\tau) = \delta(t)$, the multiplicative identity of a convolution product

$$ \sigma^{\delta}(t) = \int_{0}^{t} H(t - \tau)\delta(\tau) d\tau = H(t) \qquad (4) $$

is the realization of the transfer function and then, in any other known specific strain loading path $\epsilon^{\ell p}(t)$, we have

$$ \sigma(t) = (\sigma^{\delta} * \epsilon^{\ell p})(t) = \int_{0}^{t} \sigma^{\delta}(t - \tau) \epsilon^{\ell p}(\tau) d\tau \qquad (5) $$

fully determined by the impulse response $\sigma^{\delta}(t)$.
## 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 $\alpha_i$ and $\beta_j$ 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

$$ \dot{\sigma}(t) = \sum_{j=1}^{N} \dot{\sigma}_j = \sum_{j=1}^{N} \left( p_j E^u \dot{\epsilon} - \frac{\sigma_j - p_j \sigma^r}{\tau_j} \right) \qquad (6a) $$

* DLR constitutive equation in Laplace Domain (LBTF)

$$ \tilde{H}(p) = \frac{\bar{\sigma}(p)}{\bar{\epsilon}(p)} = \sum_{j=1}^{N} p_j \left( \frac{\tau_j E^u p + E^r}{1 + \tau_j p} \right) \qquad (6b) $$

* DLR impulse response in time (Transfer function)

$$ \sigma^{\delta}(t) = \mathcal{L}^{-1}\{\tilde{H}(p)\}(t) = \sum_{j=1}^{N} p_j E^u \delta(t) + \sum_{j=1}^{N} p_j (E^r - E^u) \frac{e^{-t/\tau_j}}{\tau_j} \qquad (6c) $$

with

- $E^u$, the unrelaxed (Young or glassy) modulus,

- $(\tau_j, p_j)$, the spectrum of relaxation times and weights defined as

$$ \tau_j = \tau_{\max} 10^{-(\frac{N-j}{N-1})d} \quad \text{and} \quad p_j = \frac{\sqrt{\tau_j}}{\sum_{j=1}^{N} \sqrt{\tau_j}} \qquad (7) $$

ensuring $\sum_{j=1}^{N} p_j = 1$ with the recursion $p_{j+1}/p_j = \alpha$ and $\tau_{j+1}/\tau_j = \beta$; also $\tau_{\max}$ is the material parameter corresponding to the maximum relaxation time, $d$ and $N$ are the number of decades and of relaxation modes of the spectrum, respectively,

- $\sigma^r(t)$, the relaxed state defined simply here as $\sigma^r(t) = E^r \epsilon(t)$ where $E^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 $N \to \infty$ 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).
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 $\varepsilon^{lp}(t)$ 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).

<table><caption>**Table 1** Considered test cases for the loading strain path ($H(t)$: Heaviside function)</caption><thead><tr><th>ε<sup>lp</sup>(t)</th><th>Time domain</th><th>Laplace domain</th></tr></thead><tbody><tr><td>Case 0: "smooth"</td><td><sup>2</sup>&frasl;<sub>&radic;3</sub>e<sup>-t/2</sup> sin(&radic;3t/2)</td><td><sup>1</sup>&frasl;<sub>s<sup>2</sup> + s + 1</sub></td></tr><tr><td>Case 1: Crenel</td><td><table><tr><td>ε<sub>0</sub></td><td>0 ≤ t ≤ t<sub>c</sub></td></tr><tr><td>0</td><td>t > t<sub>c</sub></td></tr></table><br>t<sub>c</sub> = 2s and ε<sub>0</sub> = 1 in the simulations</td><td><sup>ε<sub>0</sub></sup>&frasl;<sub>s</sub>(1 - e<sup>-st<sub>c</sub></sup>)</td></tr><tr><td>Case 2: 2 successive ramps</td><td>at + b(t - t<sub>r1</sub>)H(t - t<sub>r1</sub>)<br>t<sub>r1</sub>: time of first slope change in ramp<br>b: slope change; 2nd ramp of slope a + b</td><td><sup>1</sup>&frasl;<sub>s<sup>2</sup></sub>(a + be<sup>-st<sub>r1</sub></sup>)</td></tr><tr><td>Case 3: 3 successive ramps</td><td>at + b(t - t<sub>r1</sub>)H(t - t<sub>r1</sub>) + c(t - t<sub>r2</sub>)H(t - t<sub>r2</sub>)<br>t<sub>r2</sub>: time of second slope change in ramp<br>c: new slope change</td><td><sup>1</sup>&frasl;<sub>s<sup>2</sup></sub>(a + be<sup>-st<sub>r1</sub></sup> + ce<sup>-st<sub>r2</sub></sup>)</td></tr><tr><td>Case 4: 3 successive ramps<br>and stabilization</td><td>at + b(t - t<sub>r1</sub>)H(t - t<sub>r1</sub>)<br>+c(t - t<sub>r2</sub>)H(t - t<sub>r2</sub>)<br>-d(t - t<sub>stab</sub>)H(t - t<sub>stab</sub>)<br>t<sub>stab</sub>: time at which a dwell is imposed to the strain<br>d = (a + b + c)t<sub>stab</sub></td><td><sup>1</sup>&frasl;<sub>s<sup>2</sup></sub>(a + be<sup>-st<sub>r1</sub></sup> + ce<sup>-st<sub>r2</sub></sup> - de<sup>-st<sub>stab</sub></sup>)</td></tr></tbody></table>

# 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 time-discretization 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 $\sigma(t + \Delta t) = \sigma(t) + \Delta\sigma^t$. 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; t_f]$ with:

1. the computation of the response directly in time-domain through Equation (5), with the convolution product between the transfer function (impulse response $\sigma^\delta(t)$)—when known mathematically—and the input excitation (strain loading path $\varepsilon^{lp}(t)$). The Convolution built-in Matlab function is used to establish this output response and the result will be referred to as $\sigma^\delta * \varepsilon^{lp}(t)$ 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; t_f]. 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; Δt]) which then remains Stucked at a Stationary value equal to 1 × Δt. This input is described as:

$$ \varepsilon^{\ell p}(t) = \begin{cases} 0 & t < 0 \\ t & 0 \le t \le \Delta t \\ \Delta t & t > \Delta t. \end{cases} \qquad (8) $$

We named it hereafter URSS-Response and noted it $\sigma_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

$$ \sigma_0^-(t) = \mathcal{L}^{-1} \left\{ \bar{H}(p) \frac{1}{p^2} (1 - e^{-p\Delta t}) \right\} \qquad (9) $$

and stored as a vector of the discretized values $\sigma_0^-(t_k)$ 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 $\Delta S_k = S_k - S_{k-1}$ at time $t_k$, see Figure 1, with $S_k = (\varepsilon^{\ell p}(t_{k+1}) - \varepsilon^{\ell p}(t_k))/\Delta t$. Hence, it is approximated as a piecewise-linear function: $\varepsilon^{\ell p}(t_k) = \sum_k \Delta S_k R(t_k - (k-1)\Delta t)$.

<figure><img src="image_2.png" /><figcaption>**Figure 1** Discretization of the input $\varepsilon^{\ell p}$ in terms of slope variations.</figcaption></figure>

From the Laplace transform application, it can be shown that

$$ \sigma(t) = \mathcal{L}^{-1} \left\{ \bar{H}(p) \frac{1}{p^2} (1 - e^{-p\Delta t}) \Delta S_k e^{-p\Delta t} \right\} = \sigma_0^- * \ddot{\varepsilon}^{\ell p}(t) \qquad (10) $$

By the time translation invariance property, the incremental scheme $\sigma(t_k + \Delta t) = \sigma(t_k) + \Delta\sigma^{t_k}$ follows with $\Delta\sigma^{t_k}$ being obviously the cumulative dot product of the backward URSS-response originating recursively from time $t_k$ with the vector of slope increments up to time $t_k$. A first ramp of slope $\Delta S_1 = S_1$ acts on the VE response over all subsequent time steps. A second ramp driven by the slope change pulse $\Delta S_2$ acts on the VE response delayed from one time-step, and so
on. This makes the contribution to the current output $\Delta\sigma^{t_k}$ to be the URSS-response sequence weighted by the current input $\Delta S_{k+1}$ and past input $\Delta S_k, \Delta S_{k-1}, \dots$

$$ \Delta\sigma^k = \sigma_0^{-,k+1}\Delta S_1 + \sigma_0^{-,k}\Delta S_2 + \dots + \sigma_0^{-,1}\Delta S_{k+1} \quad (11) $$

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

$$ e_{\text{RMS}} = \sqrt{\frac{1}{n} L_2^2(y - y_{\text{ref}})} \quad \text{and} \quad e_{\text{fit}} = 100 \left( 1 - \frac{L_2^2(y - y_{\text{ref}})}{L_2^2(y - \langle y \rangle)} \right) \quad (12) $$

where $y$ and $y_{\text{ref}}$ will be defined along the text but refer always to model simulations outputs taken alternatively as the “measurement” signal $y$ and the model response signal $y_{\text{ref}}$. $\langle y \rangle$ refers to the mean of signal $y$, when averaged over the whole number of points (time steps) considered. The Step-LBTF simulation is taken in Table 2 as the output reference $y_{\text{ref}}$. We have considered

<table><caption>**Table 2** CPU times, RMS error and Fitting error for the different numerical approaches (DLR model, Case 1 excitation)</caption><thead><tr><th>Method</th><th>Δt (s)</th><th>CPU time (s)</th><th>e<sub>RMS</sub> (MPa)</th><th>e<sub>fit</sub> (%)</th></tr></thead><tbody><tr><td>Convolution</td><td>0.001</td><td>0.0019</td><td>40.88</td><td>86.4</td></tr><tr><td>Convolution (*)</td><td>10<sup>-5</sup></td><td>20.5</td><td>3.4</td><td>99.90</td></tr><tr><td>Full-LBTF (De Hoog)</td><td>0.001</td><td>0.085</td><td>0.45</td><td>99.998</td></tr><tr><td>Step-LBTF (incremental)</td><td>0.001</td><td>0.065</td><td>-</td><td>-</td></tr><tr><td>Step-ARX(1,3,0)</td><td>0.001</td><td>0.0065</td><td>0.28</td><td>99.987</td></tr></tbody></table>

here a more refined VE behavior (DLR model) where relaxation runs over a spectrum made of $N = 50$ modes. The simulation holds for the difficult case of the crenel excitation (Case 1), with $\Delta t = 0.001$s and $N = 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 $\Delta t = 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
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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.

## 3.2 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 na previous time steps of the output history and nb time steps of the input history. Retardation or delay can be considered using nk time steps. The structure of the model is given as

$$ \sigma_{\text{ARX}}(t) = - \sum_{i=1}^{na} a_i \sigma_{\text{ARX}}(t - \tau_i) + \sum_{j=1}^{nb} b_j \varepsilon_{id}(t - \tau_{j+nk}) \quad (13) $$

where $\tau_i = i\Delta t$ for $i \ge 1$, $\sigma_{\text{ARX}}(t) = 0$ for $t \le 0$ and $\varepsilon_{id}(t)$ denotes the ideal (strain) excitation imposed on the (material REV) system. In command systems vocabulary, the stress figures the output and $\varepsilon_{id}$ the input or command signal. Coefficients $a_i, b_j$ are the observer Markov
parameters. They have to be identified for a selected triplet (na, nb, nk)—thus the common notation ARX(na, nb, nk) used hereafter—during a so-called calibration experiment. Providing this latter to a specially designed optimization procedure<fnref n="2" /> allows to obtain the optimal set of parameter vectors $a_i, b_j$ (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 $\epsilon^{\ell p}$ (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 $a_i$ and $b_j$ are now able, through the ARX(na, nb, nk) 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

$$ 
\begin{aligned}
x(k + 1) &= Ax(k) + B\varepsilon(k) \\
\sigma(k) &= Cx(k) + D\varepsilon(k)
\end{aligned} 
\qquad (14) $$

with $x(k)$ the state of the system, $\sigma(k)$ the output stress observable, $\varepsilon(k)$ the input strain. Matrix $A$ is $n \times n$, $B$ is $n \times 1$ because only one input signal is considered, $C$ is $1 \times n$ because a single output signal is considered and $D$ is $1 \times 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.

# 4 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.

## 4.1 SLS model

The response to the SLS model with parameters $E^u = 1000$ MPa, $E^r = 100$ MPa, $\tau_{SLS} = 1$ s is computed using Full-LBTF, convolution and Step-LBTF approaches. Figure 2 illustrates how well

<footnote n="2">² The functions made available by Ljung from the Identification Toolbox of Matlab will be used for this purpose (Ljung 1998).</footnote>
these approaches compare in test case 3 (three successive ramps). We use then the Step-LBTF approach as the reference output $y_{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.

<table><caption>**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</caption><thead><tr><th></th><th></th><th>ValidoSLS</th><th>Valid1SLS</th><th>validation<br>Valid2SLS</th><th>Valid3SLS</th><th>Valid4SLS</th></tr></thead><tbody><tr><td rowspan="5">calibration</td><td>CaliboSLS</td><td>99.68</td><td>96.89</td><td>96.5</td><td>98.81</td><td>95.14</td></tr><tr><td>Calib1SLS</td><td>99.68</td><td>99.97</td><td>99.85</td><td>99.93</td><td>99.84</td></tr><tr><td>Calib2SLS</td><td>98.95</td><td>98.26</td><td>99.99</td><td>99.89</td><td>99.84</td></tr><tr><td>Calib3SLS</td><td>98.9</td><td>95.91</td><td>99.72</td><td>99.98</td><td>99.63</td></tr><tr><td>Calib3SLS</td><td>99</td><td>96.49</td><td>99.83</td><td>99.98</td><td>99.98</td></tr></tbody></table>

## 4.2 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 $d = 6$ decades, with a discrete set of $N = 50$ relaxation modes and parameters defined like for the SLS simulation: $E^u = 1000$ MPa, $E^r = 100$ MPa, $\tau_{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 $t_c = 3.2$ s and amplitude $\epsilon_0 = 1$ (Right axis). Figure 4 plots the residuals obtained for the different methods (always considering the Step-LBTF approach for $y_{ref}$).

<table><caption>**Table 4** $\epsilon_{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.</caption><thead><tr><th></th><th></th><th colspan="4">validation</th></tr><tr><th></th><th></th><th>ValidoDLR</th><th>Valid1DLR</th><th>Valid3DLR</th><th>Valid4DLR</th></tr></thead><tbody><tr><td rowspan="4">calibration</td><td>CaliboDLR</td><td>99.36</td><td>89.56</td><td>96.89</td><td>91.09</td></tr><tr><td>Calib1DLR</td><td>97.7</td><td>96.8</td><td>99.15</td><td>96.13</td></tr><tr><td>Calib3DLR</td><td>97.99</td><td>87.43</td><td>99.93</td><td>99.54</td></tr><tr><td>Calib4DLR</td><td>98.14</td><td>86.39</td><td>99.85</td><td>99.82</td></tr></tbody></table>

The sudden input step induces an instantaneous response. A relaxation stage is then observed and fully completed here because $t_c > \tau_{max}$ before it restarts oppositely at $t_c$ when the applied excitation ceases. It can be observed, especially at discontinuities, that for this more complex
<table><caption>**Table 5** Comparison of LBTF and ARX approaches for different excitations with $\Delta t = 10^{-3}$ s except for the convolution $\Delta t = 10^{-5}$ s. DLR behavior's law.</caption><thead><tr><th>Method</th><th></th><th>Full-LBTF</th><th>Convolution</th><th>Step-LBTF</th><th>ARS(1,3,0)</th><th>ARX(2,4,0)</th></tr></thead><tbody><tr><td rowspan="3">ε<sup>lp</sup><br>Case 0</td><td>CPU</td><td>0.08</td><td>21.8</td><td>0.19</td><td>0.006</td><td>0.011</td></tr><tr><td>e<sub>fit</sub> (%)</td><td>99.96</td><td>98.96</td><td>—</td><td>99.88</td><td>99.96</td></tr><tr><td>e<sub>RMS</sub>(MPa)</td><td>1.53</td><td>1.41</td><td>—</td><td>2.6</td><td>1.52</td></tr><tr><td rowspan="3">ε<sup>lp</sup><br>Case 1</td><td>CPU</td><td>0.074</td><td>1.062</td><td>0.161</td><td>0.005</td><td>0.008</td></tr><tr><td>e<sub>fit</sub> (%)</td><td>98.68</td><td>99.94</td><td>—</td><td>112</td><td>95.24</td></tr><tr><td>e<sub>RMS</sub>(MPa)</td><td>22.98</td><td>4.6</td><td>—</td><td>291</td><td>43.7</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3** Case 1: Crenel strain excitation and DLR viscoelastic behavior response (Valid1DLR). Comparison of LBTF approaches and ARX(1,3,0) and ARX(2,4,0) parametric models identified on test case 4 (Calib4DLR). Curve resulting from the convolution $\sigma^{\delta} \star \epsilon^{lp}$ obtained for a time step $\Delta t/100$.</figcaption></figure>

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.

<figure><img src="image_5.png" /><figcaption>**Figure 4** Residuals corresponding to Figure 3. Step-LBTF “signal” is used as reference for the calibration of ARX models.</figcaption></figure>

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
<figure><img src="image_6.png" /><figcaption>**Figure 5** ARX and LBTF approaches compared in the DLR behavior case, with excitation made of the smooth excitation $\epsilon^{\ell p} = f(t)$ (Case 0). Calibration of ARX models on Case 4.</figcaption></figure>

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.

## 4.3 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.

### 4.3.1 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

$$ \bar{Y}_n(p) = \frac{\bar{\sigma}(p)}{\bar{\varepsilon}(p)} = E^u + (E^r - E^u) \frac{\left(1 + \frac{p}{\omega_{\text{high}}}\right)^{n-1}}{\left(1 + \frac{p}{\omega_{\text{low}}}\right)^n} \quad (15) $$

with $\omega_{\text{high}} = 1/\tau_{\text{min}}$, $\omega_{\text{low}} = 1/\tau_{\text{max}}$ defining the bounds of a spectrum of relaxation times covered by the Laplace variable range. The values $E^u = 1000$ MPa and $E^r = 100$ MPa are the same as those taken for the previous SLS and DLR behaviors. The relaxation component of the transfer function $\bar{Y}_n$ depends on the last three parameters appearing in Equation (15). We consider the case of $n = 0.62$, $\tau_{\text{min}} = 0.01$ s, and $\tau_{\text{max}} = 16$ s.
This transfer function is rather elaborated, see (André et al. 2003) for details. In the limit $\omega_{high} \to \infty$, 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 $e_{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

<table><caption>**Table 6** Identification of ARX(1,3,0) model: $e_{fit}$ errors (%) in test cases for calibration and validation experiments.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="4">validation</th></tr><tr><th>Valid0Yn</th><th>Valid1Yn</th><th>Valid3Yn</th><th>Valid4Yn</th></tr></thead><tbody><tr><td rowspan="4">calibration</td><td>CaliboYn</td><td>99.11</td><td>36.77</td><td>99.31</td><td>97.08</td></tr><tr><td>Calib1Yn</td><td>96.36</td><td>95.13</td><td>99.67</td><td>95.01</td></tr><tr><td>Calib3Yn</td><td>95.44</td><td>61.2</td><td>99.87</td><td>99.31</td></tr><tr><td>Calib4Yn</td><td>97.67</td><td>37.12</td><td>99.7</td><td>99.31</td></tr></tbody></table>

<table><caption>**Table 7** Identification of ARX(2,4,0) model: $e_{fit}$ errors (in %) obtained when using the different test cases for calibration and validation experiments.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="4">validation</th></tr><tr><th>ValidoYn</th><th>Valid1Yn</th><th>Valid3Yn</th><th>Valid4Yn</th></tr></thead><tbody><tr><td rowspan="4">calibration</td><td>Calib0Yn</td><td>99.78</td><td>91.37</td><td>99.39</td><td>97.52</td></tr><tr><td>Calib1Yn</td><td>98.53</td><td>98.5</td><td>99.38</td><td>99.16</td></tr><tr><td>Calib3Yn</td><td>99.02</td><td>90.93</td><td>99.99</td><td>99.0</td></tr><tr><td>Calib4Yn</td><td>99.28</td><td>88.54</td><td>99.96</td><td>99.93</td></tr></tbody></table>

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, $e_{fit}$ and RMS errors for different test cases with the specified ARX calibration model. The $e_{fit}$ and RMS errors are established with reference taken from the 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.

### 4.3.2 Fractional Rabotnov model

The second considered model relies on the use of a relaxation kernel noted $\exists_{\alpha}^*(\beta, t)$ 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

$$ \bar{Y}_{\alpha}(p) = \frac{\bar{\sigma}(p)}{\bar{\epsilon}(p)} = E^u + (E^r - E^u) \mathcal{L}[\exists_{\alpha}^*(\beta, t)] = E^u + \frac{E^r - E^u}{p^{1-\alpha} + \beta}. \quad (16) $$

It does not correspond to the same mathematical structure as given by Equation (15). In the above expression, $\alpha$ denotes the non-integer exponent, $\beta$ the inverse of a (relaxation) time, and the Rabotnov kernel is defined in temporal domain by

$$ \exists_{\alpha}^*(\beta, t) = E_{\alpha+1, \alpha+1}(\beta t^{\alpha+1}) = t^{\alpha} \sum_{n=0}^{\infty} \frac{(\beta t^{1+\alpha})^n}{\Gamma[(n+1)(1+\alpha)]} \quad (17) $$
<table><caption>**Table 8** Comparison of LBTF (NIF-Y<sub>n</sub>) and ARX approaches for different excitations (time step Δt = 10<sup>-3</sup> s) and ARX calibration models.</caption><thead><tr><th>Method</th><th></th><th>Full-LBTF</th><th>Step-LBTF</th><th>ARS(1,3,0)</th><th>ARX(2,4,0)</th></tr></thead><tbody><tr><td rowspan="3">ε<sup>ℓp</sup><br>Case 0<br>Calib0</td><td>CPU</td><td>0.2034</td><td>0.1325</td><td>0.0054</td><td>0.01</td></tr><tr><td>e<sub>fit</sub> (%)</td><td>99.99</td><td>—</td><td>99.64</td><td>99.9995</td></tr><tr><td>e<sub>RMS</sub> (MPa)</td><td>1.7</td><td>—</td><td>3.37</td><td>0.4</td></tr><tr><td rowspan="3">ε<sup>ℓp</sup><br>Case 1<br>Calib1</td><td>CPU</td><td>0.085</td><td>0.126</td><td>0.005</td><td>0.0076</td></tr><tr><td>e<sub>fit</sub> (%)</td><td>99.6</td><td>—</td><td>99.65</td><td>99.87</td></tr><tr><td>e<sub>RMS</sub> (MPa)</td><td>26.4</td><td>—</td><td>23.9</td><td>14.7</td></tr><tr><td rowspan="3">ε<sup>ℓp</sup><br>Case 4<br>Calib1</td><td>CPU</td><td>0.08</td><td>0.25</td><td>0.006</td><td>0.009</td></tr><tr><td>e<sub>fit</sub> (%)</td><td>99.99</td><td>—</td><td>99.97</td><td>99.99</td></tr><tr><td>e<sub>RMS</sub> (MPa)</td><td>0.54</td><td>—</td><td>5.8</td><td>3.3</td></tr><tr><td rowspan="3">ε<sup>ℓp</sup><br>Case 4<br>Calib4</td><td>CPU</td><td>0.079</td><td>0.23</td><td>0.0065</td><td>0.009</td></tr><tr><td>e<sub>fit</sub> (%)</td><td>99.99</td><td>—</td><td>99.99</td><td>100</td></tr><tr><td>e<sub>RMS</sub> (MPa)</td><td>0.54</td><td>—</td><td>2.36</td><td>0.25</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6** ARX and LBTF approaches compared in the fractional behavior case (NIF-Y<sub>n</sub>), with the crenel excitation ε<sub>1</sub><sup>ℓp</sup> (Case 1). Calibration of ARX models on Case 1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7** ARX and LBTF approaches compared for the NIF behavior case (NIF-Y<sub>n</sub>), with excitation made of the three ramps+dwell excitation ε<sub>4</sub><sup>ℓp</sup> (Case 4). Calibration of ARX models on Case 4.</figcaption></figure>
This kernel is expressed through the 2-parameters Mittag-Leffler function $E_{\alpha+1,\alpha+1}$ and $\Gamma$ denotes the Gamma function. Figure 8 provides an example of the Mittag-Leffler function calculated from Equation (17) and by performing De Hoog's inverse Laplace algorithm on the second term of the rhs of Equation (16). Mittag-Leffler functions are evident functions arising from irreversible

<figure><img src="image_9.png" /><figcaption>**Figure 8** Rabotnov Kernel – Comparison of expression $\exists_{\alpha}^{*}(\beta, t)$ defined as a Mittag-Leffler function in Equation (17) and the Laplace inverse of $\frac{1}{\rho^{1-\alpha+\beta}}$: $\alpha = 0.8$ to $0.2$ and $\beta = 1$ i.e. characteristic relaxation time of $1$ s.</figcaption></figure>

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:

$$ \sigma^{\delta}(t) = \mathcal{L}^{-1}\bar{H}(p)(t) = E^{u}\delta(t) + (E^{r} - E^{u})t^{-\alpha}E_{1-\alpha,1-\alpha}(-\frac{1}{\tau}t^{1-\alpha}) \quad (18) $$

As a recall, this leads for any loading path $\epsilon^{\ell p}(t)$ to compute the convolution product $\sigma(t) = (\sigma^{\delta} * \epsilon^{\ell p})(t)$. 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 $\epsilon_{4}^{\ell p}$ (Case 4 of the multisequence loading path). Simulations were performed for the same values of $E^{u} = 1000$ MPa, $E^{r} = 100$ MPa as before, with $\beta = 1 s^{-1}$ and the four $\alpha$ 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 $\epsilon_{3}^{\ell p}$. 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%).

# 5 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
<figure><img src="image_10.png" /><figcaption>**Figure 9** Fractional Rabotnov VE kernel with excitation $\epsilon_3^{fp}$ made of three ramps (Case 3). ARX(1,3,0) models calibrated on Case 4.</figcaption></figure>

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.

## References

Agbossou, A., A. Bergeret, K. Benzarti, and N. Alberola (1993). Modelling of the viscoelastic behaviour of amorphous thermoplastic/glass beads composites based on the evaluation of the complex Poisson's ratio of the polymer matrix. *Journal of Materials Science* 28(7):1963–1972. [DOI], [HAL].

André, S., J. Boisse, and C. Noûs (2021). An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite. *Journal of Theoretical, Computational and Applied Mechanics*. [DOI], [OA].

André, S., Y. Meshaka, and C. Cunat (2003). Rheological constitutive equation of solids: a link between models based on irreversible thermodynamics and on fractional order derivative equations. *Rheologica Acta* 42(6):500–515. [DOI], [HAL].

Bagley, R. L. and P. J. Torvik (1986). On the fractional calculus model of viscoelastic behavior. *Journal of Rheology* 30(1):133–155. [DOI].

Blaise, A., S. André, P. Delobelle, Y. Meshaka, and C. Cunat (2016). Advantages of a 3-parameter reduced constitutive model for the measurement of polymers elastic modulus using tensile tests. *Mechanics of Time-Dependent Materials* 20(4):553–577. [DOI], [HAL].

Christensen, R. (1982). *Theory of Viscoelasticity. An Introduction*. 2nd ed. Elsevier. [DOI].

Corbin, S. F. and D. M. Turriff (2012). Thermal diffusivity by the laser flash technique. *Characterization of Materials*. Wiley, pp 1–10. [DOI].

Cunat, C. (2001). The DNLR approach and relaxation phenomena. Part I – Historical account and DNLR formalism. *Mechanics of Time-Dependent Materials* 5(1):39–65. [DOI].

Davies, B. and B. Martin (1979). Numerical inversion of the Laplace transform: a survey and comparison of methods. *Journal of Computational Physics* 33(1):1–32. [DOI].

De Hoog, F. R., J. H. Knight, and A. N. Stokes (1982). An improved method for numerical inversion
of Laplace transforms. *SIAM Journal on Scientific and Statistical Computing* 3(3):357-366. [DOI].

den Iseger, P. (2006). Numerical Transform Inversion using Gaussian quadrature. *Probability in the Engineering and Informational Sciences* 20(1):1-44. [DOI].

Dingfelder, B. and J. A. C. Weideman (2014). An improved Talbot method for numerical Laplace transform inversion. *Numerical Algorithms* 68(1):167-183. [DOI], [ARXIV].

Dooling, P., P. Buckley, and S. Hinduja (1997). An intermediate model method for obtaining a discrete relaxation spectrum from creep data. *Rheologica Acta* 36(4):472-482. [DOI].

Emri, I. and N. W. Tschoegl (1993). Generating line spectra from experimental responses. Part I: Relaxation modulus and creep compliance. *Rheologica Acta* 32(3):311-322. [DOI].

Friedrich, C., H. Schiessel, and A. Blumen (1999). Constitutive behavior modeling and fractional derivatives. *Rheology Series*. Vol. 8. Advances in the Flow and Rheology of Non-Newtonian Fluids. Elsevier, pp 429-466. [DOI].

Fuentes, F., L. Demkowicz, and A. Wilder (2017). Using a DPG method to validate DMA experimental calibration of viscoelastic materials. *Computer Methods in Applied Mechanics and Engineering* 325:748-765. [DOI], [OA].

Gallican, V. and R. Brenner (2019). Homogenization estimates for the effective response of fractional viscoelastic particulate composites. *Continuum Mechanics and Thermodynamics* 31(3):823-840. [DOI], [HAL].

Gevers, M. (2006). System Identification without Lennart Ljung: what would have been different? *Forever Ljung in System Identification*. Ed. by T. Glad and G. Hendeby, pp 61-85. [HAL].

Jannot, Y. and A. Degiovanni (2018). *Thermal Properties Measurement of Materials*. Wiley. [DOI].

Koeller, R. C. (1984). Applications of fractional calculus to the theory of viscoelasticity. *Journal of Applied Mechanics* 51(2):299-307. [DOI].

Ljung, L. (1998). System Identification. *Signal Analysis and Prediction*. Ed. by A. Procházka, J. Uhlíř, P. W. J. Rayner, and N. G. Kingsbury. Applied and Numerical Harmonic Analysis. Birkhäuser Boston. Chap. 11, pp 163-173. [DOI].

Lévesque, M., M. Gilchrist, N. Bouleau, K. Derrien, and D. Baptiste (2007). Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media. *Computational Mechanics* 40(4):771-789. [DOI], [HAL].

Menard, K. and N. Menard (2020). *Dynamic Mechanical Analysis*. 3rd ed. CRC Press. [DOI].

Metzler, R. and J. Klafter (2000). Kramers' escape problem with anomalous kinetics: non-exponential decay of the survival probability. *Chemical Physics Letters* 321(3-4):238-242. [DOI].

Metzler, R. and J. Klafter (2002). From stretched exponential to inverse power-law: fractional dynamics, Cole-Cole relaxation processes, and beyond. *Journal of Non-Crystalline Solids* 305(1-3):81-87. [DOI].

Phan, M. Q. and R. W. Longman (1996). Relationship between state-space and input-output models via observer Markov parameters. *Dynamics and control of structures in space III* 19:185-200. [HAL].

Rekik, A. and R. Brenner (2011). Optimization of the collocation inversion method for the linear viscoelastic homogenization. *Mechanics Research Communications* 38(4):305-308. [DOI], [HAL].

Sabatier, J., P. Lanusse, P. Melchior, and A. Oustaloup (2015). Fractional order differentiation and robust control design. *Intelligent systems, control and automation: science and engineering* 77:13-18. [DOI].

Schapery, R. A. (1962). Approximation methods of transform inversion for viscoelastic stress analysis. *Fourth US National Congress on Applied Mechanics* (University of California, Berkeley, USA, June 18-21, 1962), pp 1075-1084. [HAL].

Schiessel, H., R. Metzler, A. Blumen, and T. F. Nonnenmacher (1995). Generalized viscoelastic models: their fractional equations with solutions. *Journal of Physics A: Mathematical and General* 28(23):6567-6584. [DOI].

Sevostianov, I., V. Levin, and E. Radi (2015). Effective properties of linear viscoelastic microcracked materials: Application of Maxwell homogenization scheme. *Mechanics of Materials* 84:28-43. [DOI], [OA].

Shaterzadeh, M., C. Gauthier, J. F. Gerard, C. Mai, and J. Perez (1998). Dynamic mechanical prop-
erties of spherical inclusions in polymer composite: A self-consistent approach considering
morphology. *Polymer Composites* 19(6):655–666. [DOI].

Sneddon, I. N. (1951). Chapter 5: The conduction of heat in solids. *Fourier Transforms*. McGraw-Hill. ISBN: 9780486685225.

Sorvari, J. and J. Hämäläinen (2010). Time integration in linear viscoelasticity—a comparative study. *Mechanics of Time-Dependent Materials* 14(3):307–328. [DOI].

Stehfest, H. (1970). Algorithm 368: Numerical inversion of Laplace transforms. *Communications of the ACM* 13(10):624. [DOI].

Talbot, A. (1979). Accurate numerical inversion of Laplace Transforms. *IMA Journal of Applied Mathematics* 23(1):97–120. [DOI].

Trumel, H., F. Rabette, F. Willot, R. Brenner, E. Ongari, M. Biessy, and D. Picart (2019). Understanding the thermomechanical behavior of a TATB-based explosive via microstructure-level simulations. Part I: Microcracking and viscoelasticity. *Europyro 44th International Pyrotechnics Seminar* (Tours, France, June 3–7, 2019). [HAL].

Tschoegl, N. (1961). The characterization of linear viscoelastic behaviour by respondance functions. Stress circuit theory. *Kolloid-Zeitschrift* 174(2):113–123. [DOI].

Tschoegl, N. (1989). *The Phenomenological Theory of Linear Viscoelastic Behavior*. Springer. [DOI].

Uluutku, B., M. R. McCraw, and S. D. Solares (2022). Direct measurement of storage and loss behavior in AFM force-distance experiments using the modified Fourier transformation. *Journal of Applied Physics* 131(16):165101. [DOI], [OA].

Weron, K. and M. Kotulski (1996). On the Cole-Cole relaxation function and related Mittag-Leffler distribution. *Physica A: Statistical Mechanics and its Applications* 232(1-2):180–188. [DOI].

Wu, L. (2022). Equivalence of SS-based MPC and ARX-based MPC. [ARXIV].

Özisik, N. (1993). The use of Duhamel’s theorem. *Heat Conduction*. 2nd ed. Wiley, pp 195–213.
ISBN: 9780471532568.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

**Authors' contributions** The author confirms sole responsibility for carrying the study and writing the manuscript.

**Supplementary Material** Explanations and Matlab scripts have been made available in the form of an executable companion paper as a live script (or notebook) at the permalink: swh:1:dir:3d4ad405846eb4dcf7ca35f9e72e85c6e817657c; origin=https://github.com/stolasandre/JTCAM-11178-AndreNous-ARX-models-for-Linear-ViscoElasticity;visit=swh:1:snp:9b0c5be3f62cd6fa3156cbc5ea38238f863816e3;anchor=swh:1:rev:dfcca448c484c4b52bb53a006c65f5doc6e8420c.

**Acknowledgements** The authors wish to express their posthumous gratitude to José RAGOT, professor at ENSG (Ecole Nationale Supérieure de Géologie) in Nancy, a specialist in control theory, who passed away in 2022, for the exchanges we were able to have together on some issues related to this article.

**Funding** This work was supported by the Energy, Materials, Processing, Products (EMPP) research department of Lorraine University and the ICEEL Carnot Institute.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.