# Review of "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, Laurence Brassart

► To cite this version:

Stéphane André, Camille Noûs, Laurence Brassart. Review of "Solving viscoelastic problems with a Laplace transform approach supplanted by ARX models, suggesting a way to upgrade Finite Element or spectral codes". 2023. hal-04174179

HAL Id: hal-04174179

https://hal.science/hal-04174179

Submitted on 14 Dec 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

Open Review

OAI hal-04174179

Reviewed Article

DOI 10.46298/jtcam.10304

Licence

CC BY 4.0

©The Authors

# Review of “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" />, Camille NOÛS<fnref n="2" />, and Laurence BRASSART<fnref n="3" />,<fnref n="4" />

## Review of version 1

Permalink: hal-03845394v1

### Reviewer 1 (Anonymous)

**Reviewer**

This paper deals with the use of linear viscoelastic laws expressed in the Laplace domain in structural analysis codes using time stepping schemes. Two ways are proposed: The first one relies on a numerical algorithm which gives an efficient estimate for the convolution product used to compute the inverse Laplace transform which is named STEP-LTBF. The second relies on a new way (to the reviewer knowledge) to express the visco-elastic constitutive relation and give the model named ARX. Both methods are used to simulate the time response of various viscoelastic laws including fractional models. To be published this paper should be modified to address the following points.

**Reviewer**

In view of the conclusions, it appears that the numerical scheme STEPLTBF is in fact only used to identify the parameters of the ARX model. If this is the case, its interest is less than presented in this paper because there is no problem of computation time for this step of calibration.

**Authors**

This is a true observation and it is exact that regarding the ARX model identification, the computation time of the Step-LBTF approach is not an argument. ARX model identifications could have been performed from any precise solution obtained by one or another mathematical method.

These are the few elements which led us to favour the method.

1. The Laplace perspective: We develop this approach specifically because the original starting point was to proceed from VE behavior's law defined directly in Laplace (or frequency) domain. This was stated in the introduction.

2. According to this perspective, we had to account for cases where the kernel relaxation has no formal expression in time domain. This was actually a strong motivation for us because we show in other works (cited paper Andre et al. 2003) that the fractional admittance model referred to as Oustaloup Model in the article, is very efficient and parsimonious to represent VE behavior with "fractal dynamics" or as homogenized model for heterogeneous materials with constituents being VE (work not published yet but Figures 5 and 6 next, will illustrate this point).

3. We show it offers a nice derivation for a 'precise' incremental scheme. Incremental schemes are obvious for (inherent to) LTI systems from a mathematical stand point. Actually, it is possible to

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

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

<footnote n="3">³ Department of Engineering Science, University of Oxford, Oxford, United Kingdom</footnote>

<footnote n="4">⁴ Editor</footnote>
describe the loading path in terms of successive crenels and to apply the same type of reasoning
for an incremental scheme, but with the introduction of a mid-point rule to stay in the same level
of precision. Here the choice was to invoke Duhamel's theorem to work with test functions of
ramp-type from which an incremental algorithm is naturally derived.

4. Being precise and with low computation times here in the homogenized case, the step-LBTF approach by itself may be of interest when performing parameter estimation through optimization procedures applied onto real data.

5. Last but not least, the key idea to develop an ARX strategy is to obtain a so-called system input-output relationship, where the contribution to the current output (stress state increment in our step-LBTF approach) appears as a linear combinations of past inputs (strain rate increment) weighted by a the urss-response sequence. Maybe naively, we therefore think that introducing ARX models in viscoelasticity should be done (quasi-epistemologically) through this entry point. The remark of the referee and present answer led us to make this clearer in the new version of the text to be more informative regarding this approach. With more important consequences on the new proposed version, we consider that this remark may originate from an apparent excessive insistence around the step-LBTF approach in the initial version while the central point remains the introduction of ARX models. Therefore the mention to it has been suppressed from the title.

**Reviewer** The proposed ARX model seems very promising to identify complex viscoelastic behavior laws with few parameters and requiring only to store few variables. But, to really discuss the advantages of this model in a vision of use by an engineer, it would be necessary to validate it for complex 3D loads for which the loading directions are not constant.

**Authors** The term “suggested” has been retained in the title instead of “promising”, as suggested by referee 1 because indeed, we do not present results which are the output of a numerical code having the ARX models embedded into it.

Actually the implementation in a FFT solver is done (under AMITEX developed by L. Gélébart at the CEA, French Atomic Energy and Alternative Energy Commission) and a new text is in view to discuss numerical aspects and explore various case of composite materials under more complex loads. We can give here very preliminary results which obviously cannot claim for the moment to be published. A particulate composite is considered with viscoelastic matrix and elastic spherical beads (volume fraction 30%, log distribution in size around 0.26 for a VER of side length equal to 1). The mechanical solicitation is in shear (not a complex loading then).

<figure><img src="image_2.png" /><figcaption>**Figure 1** Left: Microstructure of a Particulate Composite REV for a 3D application. Right: Stress output (Left axis) to the ramp+dwell macroscopic strain applied on the REV (right axis). Curves for 32³, 128³ and 256³ voxels. Behavior's law is DLR model with 50 modes. Figures kindly provided by Julien Boisse (LEMTA laboratory) who implemented the ARX approach in the AMITEX software.</figcaption></figure>

We compare

* the semi-analytical model (Label sa) derived from the DLR approach (50 modes) according to the scheme of Sorvari and Hämäläinen (2010). This incremental scheme is analytically exact.

* the step-LBTF approach making use of the urss-response (Label urss)

* the ARX substitution (Label ARX130)

By changing the number of voxels used to describe the microstructure, we observe the obvious
influence on the precision but more interestingly the CPU times evolution. For the 128³ discretization, the ARX130 is the fastest one 1420 seconds (24 processor Intel Xeon Gold 5118 CPU at 2.30 GHz) even if the sa algorithm works with analytical expressions (×10 ARX CPU time) and with the step-approach being then highly penalized (×50 ARX CPU time). The initial project of the manuscript was to couple the ARX models presentation and performances with an example of implementation in a 3D code but it appears during the writing process and implementation of the computations that it will bring the arguments and discussion much too far. Firstly, although rather straightforward, there are many aspects to discuss regarding the implementation in a 3-D code. Secondly, it was already a sum for a single paper to present:

* on one hand, the computations made to obtain and validate an approach based on Laplace transforms for the expression of constitutive VE laws (the step-LBTF approach) and,

* on the other hand, its substitution for ARX models which philosophy had to be introduced first and applicability to a wide range of behavior's law demonstrated.

It was then mainly a matter of choice at one stage.

Besides, we can describe how the application of the present approach to the 3D cases must be implemented in a numerical code and bring arguments to show that no technical problems are to be expected.

* In a 3D perspective, two ARX models must be associated to the VE behavior description in terms of $\kappa$ and $\mu$ moduli which is made independently in the exact same way as shown in the present paper for the pure 1D behavior.

* As mentioned in the text (and more precisely in the new version), one drawback of using ARX models is that they are identified for a defined and constant time step. Simulations in a numerical code must be carried out with the same time step.

* To mimic the behavior of the incremental process that will happen in a numerical code and to prove the consistency of ARX models, we have performed the following steps:

1. For a given input strain path $\epsilon^{\text{true}}(t)$ (Case 0: “smooth” here), the stress output is computed ‘exactly’ by the step-LBTF approach (we denote it $\sigma^{\text{exact}}$)

2. For a non-optimal selected ARX($na, nb, nk$) model (identified on a non-optimal calibration experiment for example or with parameters $na, nb$ being too low), a corrected strain path $\epsilon^{\text{corr}}(t)$ is determined incrementally by requiring for any time step $t_i$ a minimization of the gap $\epsilon = |\sigma^{\text{exact}}(t_i, \epsilon^{\text{true}}) - \sigma^{\text{ARX}}(t_i, \epsilon^{\text{corr}})|$ This minimization process figures out the iterative one of a numerical code to obtain a convergence of a criterium based on the verification of the equilibrium equation $\overline{\text{div}\sigma} = 0$

3. Report the Recalculated Strain $\epsilon^{\text{corr}}$ and appreciate the consistency of the ARX modelling. In Figure 2, this process was made with an ARX(1, 3, 0) computed through equation (5.6) (now Equation (10)) i.e. the polynomial form. Two level of precision were considered $\epsilon = 10^{-3}$ and $\epsilon = 10^{-5}$ which do not impact the result. It is shown that for this “basic” ARX model, the correction lies around 3% of the signal variation.

In Figure 3, the same process was reconducted with a more precise ARX model ARX(3, 5, 1) and it was computed with the state-space representation (Equation (11)) to show that both approaches can be used and will behave similarly. Because this ARX model is more refined than the previous one, the reconstructed strain path is nearly the original one. The gap plotted in the figure (right axis) is now below 1.3% the signal variation.

According to us, this proves the consistency of the ARX approach. In the FFT solver we use, the iterations performed within a same time-step to calculate the new voxel strain through a set of convergence criteria lead to the good result. The numerical study we are currently conducting mainly consider the relationship between accuracy of the result and the CPU time associated to the simulations, which depends on the number of iterations within a convergence step and therefore the rationale in selecting whether a raw nor a refined ARX model.

**Reviewer** It would also be necessary to confront it with classical methods of collocation which make it possible to identify this type of behavior by Prony series, already implemented in structural analysis codes such as Abaqus, and, perhaps, discuss it according to the number of terms of the series.

**Authors** This is a very interesting point to discuss and obviously our sentence in the introduction with
<figure><img src="image_3.png" /><figcaption>**Figure 2** ARX 130 – Left axis: Recalculated strain path to “adjust” the ARX output response to the exact one (Stress). Right axis: gap between the recalculated and true strains.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3** ARX 351 – Left axis: Recalculated strain path to “adjust” the ARX output response to the exact one (Stress). Right axis: gap between the recalculated and true strains.</figcaption></figure>

associated references could have been developed more. We wrote: "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 & Brenner, 2011) or on homogeneous materials by using the classical collocation method (Schapery, 1962) when trying to identify the relaxation kernel from its expression in Laplace domain."

This reference to the difference in approach between the collocation method and arx was probably too sober. Following the referee's remark, we add some elements of discussion in the new version of the manuscript in order to mention the connections and at same time, distinctions between the collocation method and the arx approach which, from our point of view, doesn't really call for a confrontation. And we hope that the comments below will help refine our perspective on this issue.

* It must be said first (and this comment was added in the text) that we do not consider the ARX approach as a model for viscoelasticity. Although mathematically in link with it through the system approach in case of Linear Time Invariant property, the arx models are not considered as, and have not to be substituted to, a pertinent physical approach of behavior's law which is the only one able to involve clear physical concepts related to the dynamics of the microstructure.

* ARX models could be considered as "proxy" models in a numerical perspective, which is what is done here.

* Confrontation in terms of philosophy

1. What is in common? The collocation method and present ARX approach both offer a way to start from a Laplace or Laplace-Carson formal expression of the behavior's law (more precisely of a respondance) and then use algorithms in time to compute the response to a given excitation.

2. What is different? The philosophy of the collocation method is providing a surrogate behavior's law of fixed mathematical structure which parameters are the output of an inverse problem (Parameter estimation problem). The drawbacks and correlative appropriate precautions have been widely discussed (Levesque, Rekik, Tschoegl, references were initially given in the paper)
As far as we know they are quasi exclusively considered to fit relaxation moduli curves. The parameters identified by a collocation scheme may lead to physical interpretations or analysis. On the contrary, parameters of ARX models do not authorize any physical comments. They capture the LTI behavior on an adaptable horizon of time - which depends on the required precision - acting like models compressing the number of degree of freedoms.

* **Practical confrontation**

This confrontation is done actually. Although it was not presented as such, such confrontation can be analyzed precisely through the results and discussion made in section 4.2 when comparing the ARX models to the case of the DLR model. The DLR model can be considered as nothing else than a collocation method but where the time constants of the Prony series are selected in an exhaustive manner (a “great” number, spaced equally in logarithmic scale) and the weight associated to each term directly defined in a recursive manner from these time constants (according to a thermodynamic principle). Full details on that matter can be found in the cited paper Andre et al (2003, Rheologica Acta). In other words, and from our point of view in having applied this DLR model to real or synthetic data in an inversion process (see for example Blaise et al. MTDM 2016 and Andre et al. MTDM, 2020), this DLR model is an optimized way for applying the collocation method.

To go a little bit further in the discussion:

* This justifies the introduction in this article of the fractional 'Oustaloup' model because the DLR model (a regular Ladder-type collocation scheme) leads to this latter in the asymptotic limit of an infinite number of terms in the Prony series and preserving the recursivity. For the Oustaloup model, there is no analytical expression of the relaxation kernel in time domain and therefore the recommended methodology (precise computation with LBTF approach and ARX model identification) seems to us of potential interest.

* The collocation method initially “suffers” from the risk to introduce a bias in representing the VE behavior because in general, no sensitivity analysis and control of the identifiability of the parameters is made when applying this method with a series of 3, 4, 5, or whatever terms (we do not check the Abaqus way of combining with this problem). Since the proposition of Schapery, great efforts have been made especially in the 90’s, 2000’s to provide good inversion algorithms for the collocation method (Rekik and Brenner 2011 for example) and find appropriate representation of the behavior’s law in time domain, avoiding the pitfall’s associated to this ill-posed problem. This being said (and well known), the collocation scheme already found its place in any incremental scheme in time.

* To conclude, in the theory of LTI systems, the close connections between these apparently disconnected/different approaches (maybe because introduced successively and outside precisely a global viewpoint?) are ‘soluble’ in the frequency approach and the representation of a transfer function with appropriate numbers and scheme of distributed zeros and poles of its rational form. It is clear by the way that the collocation method was recognized as a path to solve the same question as we addressed here. We present here an alternative path.

If the DLR model, as a possible collocation scheme, can be proved efficient to model some VE behavior, then this model and its incremental version have to be used for sure in numerical computations. But this is not always possible. We consider below the case of a particulate composite material *with a VE matrix* represented by the SLS with one relaxation time. The volume fraction of elastic spherical particles is 50% but we consider a log distribution in size of the particle diameters to introduce a huge spectrum of length scales. The behavior of the heterogeneous material is computed by a FFT solver (CRAFT Virtual DMA), on a 512³ voxels REV, and the frequency mode of this simulation allows the direct identification of the LBTF (Laplace behavior or transfer function).

The two figures below show the conservative and loss shear moduli versus frequency.

* in Figure 4, the DLR model applied as effective homogeneous RVE model is not proved very efficient in this case (50 modes were used and additional modes i.e. “number of terms in the series” do not change things).

* in Figure 5, the fractional “Oustaloup” model presented in this article is able to perfectly represent the REV effective behavior and as a consequence, because its corresponding time description is not available, would be efficiently replaced by an ARX model for the modelling of a larger
structure and the prediction of a temporal behavior.

<figure><img src="image_5.png" /><figcaption>**Figure 4** DLR model applied to represent the effective behavior of a Particulate VE composite (3-parameter model, 2 hyper-parameters)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5** $Y_n$ fractional model applied to represent the effective behavior of a Particulate VE composite (5-parameter model)</figcaption></figure>

The alternative path of using a collocation method to handle Laplace transformed definition of VE behavior's law is well evoked in the introduction. In this regard, we have been careful to speak of a different philosophy. But we found this remark very valuable and to avoid some kind of "dead angle" in our text, we added a few sentences in the new version to this end, precisely when the DLR model is presented.

**Reviewer** Typos and minor corrections: It should be said in the text after equation (2) that $f(i)$ is the $i$th derivative of $f$ with respect to time.

**Authors** OK, done

**Reviewer** There is a problem in the last term of equation (3), I think that the numerator and the denominator are inversed.

**Authors** Absolutely, thanks. modified

**Reviewer** Page 11 there is a problem in the number of the section, probably 3.2 instead of 3.

**Authors** Yes, modified

**Reviewer** I can't understand the relationship between equations (5.7) and (5.6).
**Authors** These equations are now equations (10) and (11) in the corrected version of the manuscript. The relationship is not so easy to establish and is described in a clear manner (from our point of view) in the cited paper of Phan and Longman. We also add the reference of Wu (Equivalence of SS-based MPC and ARX-based MPC). Both the polynomial (equation 10) and State-Space (equation 11) formulations can be seen as equivalent. They are part of a very classical framework in control Theory. We used both of them in our own computations (as shown previously in Figure 1).

They do not behave in exactly the same manner of course. Both versions of these ARX formulations were included in the FFT solver code AMITEX but this is part of the ongoing work to examine the effect of one or the other formulation in terms of number of iterations, precision...

We slightly modify the text to cover the points discussed above. We leave the interested reader in the mathematical aspects of going from one representation to the other, with the selected references.

**Reviewer** In the appendix, it will be better to give the algorithm with mathematical equations and not with the python syntax.

**Authors** Actually, we suggest to suppress this appendix. Firstly because it is not really an algorithm and it directly follows the computation process raised by the Step-LBTF approach that we detailed a little bit more by adding two equations. But the main reason is that in between the submission and now, we develop an executable (and executive) paper in the form of a digital notebook (a LiveScript in Matlab denomination). This means that we will attach an open-source digital document to this article, providing both the useful explanations outlined in the article and lines of an executable code. This will allow the reader of the article to see how the computations were made.

**Reviewer 2 (Anonymous)**

**Reviewer** General appreciation.

Finite Element codes as well as FFT solvers are based on time-discretization approach for viscoelastic materials under transient excitation. As a result, defining viscoelastic laws in terms of their Laplace Transforms (in the frequency domain as provided by DMA) is up to now not possible. The paper proposes a promising and efficient strategy in order to face this problem. It is successfully tested for a material point and different viscoelastic laws of increasing complexity. In a first part, a way to compute incrementally the temporal response of a viscoelastic material with its Laplace-defined Behavioral Transfer Function (Step-LBTF) is advanced and validated by comparison to other more "global" approaches (Full-LBTF, convolution approach). However, the whole history of the variables is still required which is considered as too heavy at each point for future structural calculations. This leads the Authors to propose the second original contribution of the paper consisting in substituting ARX parametric models to Laplace Transforms in the incremental time-step approach. The good performances of the strategy in terms of CPU times and accuracy are illustrated for the Standard Linear Solid (SLS) behaviour, the Dynamic of Linear Relaxations (DLR) model and fractional models (Non-integer Outstaloup model and Rabotnov model). Special attention is paid to the calibration and validation procedures based on four different loading histories. At last, the paper is very well written.

As a result, I recommend the publication of the paper in the International Journal of Theoretical, Computational and Applied Mechanics (JTCAM) provided the following corrections and remarks/questions be addressed.

**Authors** Thanks first for this refereeing work. The main modifications have been made in section 3.

**Reviewer** All the simulations are done for a unique material point. The strategy is not yet implemented in a FE code for structural calculations. I propose to change the title by adding the word "promising": "a promising way to upgrade Finite Element or spectral codes." It will be more accurate and will avoid any confusion.

**Authors** Yes we agree, maybe a "?" was missing in this proposition. Actually, the implementation in a FFT solver is done and studies are on the way to make clear the gain offered by ARX models. The initial project of the paper was indeed to couple ARX models presentation and performances

<figure><img src="image_7.png" /></figure>
with an example of implementation in a 3D code but it appears during the writing process and
implementation of the computations that it will bring the arguments and discussion much too
far. Firstly, although rather straightforward, there are many aspects to discuss regarding the
implementation in a 3D code. Secondly, it was already a sum for a single paper to present:

* on one hand, the computations made to obtain and validate an approach based on Laplace transforms for the expression of constitutive VE laws (the step-LBTF approach) and,

* on the other hand, its substitution for ARX models which philosophy had to be introduced first and applicability to a wide range of behavior's law demonstrated.
It was then mainly a matter of choice at one stage.

"Promising" is an adjective that makes us really uncomfortable (for reasons related to the ideology governing the politics of science nowadays). Maybe "suggesting a way" to upgrade FE codes is more fairly measured. After all, it is primarily a retrospective evaluation that determines whether a scientific result has proven to be relevant or not.

**Reviewer** In the abstract, please detail the abbreviation ARX for non-specialist Readers.

**Authors** It has been done

**Reviewer** Page 10: Except for the convolution approach, the time step seems to be set to 0.001 s. This time step leads to accurate results with respect to the two quality metrics employed. However, a deeper discussion regarding the time step influence on results and this choice could be added.

**Authors** The Laplace approaches (Full-LBTF and Step-LBTF responses) rely on exact analytical formulas and the return to time-domain relies on an algorithm (De Hoog) which is not sensitive to the time sampling. It proceeds unequivocally for a given sampled time. Figures ?? and ?? below illustrate the results obtained with both approaches where the time step is set to $\Delta t = 0.01$ s instead of 0.001 s for the case of excitations of type 0 (smooth excitation) and 1 (crenel). The convolution approach (Matlab algorithm) with $\Delta t = 10^{-5}$ s serves as a reference. It can be seen that both approaches give the same result.

<figure><img src="image_8.png" /><figcaption>**Figure 6** Step-LBTF and Full-LBTF for $\Delta t = 0.01$ s (Smooth excitation)</figcaption></figure>

Figure ?? corresponds to the same case but with $\Delta t = 0.2$ s (35 points over the time range [0; 7]). It can be seen that the input being exactly (even if less) discretized, there is no induced error.

As a conclusion, using Laplace transforms and of course the same precise inverse algorithm, there is no influence of the time step over the precision of the approaches. The difference lies in the fact that, on one hand, the full-LBTF with same time step introduces errors near discontinuities (figure 3) and, on the other hand, the application of DeHoog's algorithm for each time step multiplies then the CPU times (by a factor of 50 for a time step of 0.001 s), which clearly undermines its use in the incremental schemes of numerical codes.

The time step of 0.001 s was only chosen in order to produce representative CPU times over a large number of time steps (7000) and to compare test cases having this same time step, the crenel excitation requiring it near discontinuities.

The paragraphs of Section 3.1 discussing these aspects have been enriched with more details that probably were lacking in the original version. We also recomputed all the values in Table 1 for the case of the DLR model (instead of the SLS one) which is more representative of a true
<figure><img src="image_9.png" /><figcaption>**Figure 7** Zoom of the previous figure.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8** Step-LBTF and Full-LBTF for $\Delta t = 0.01$ s (Crenel excitation)</figcaption></figure>

VE behavior. This, we believe, makes the specificity of the different approaches and the gain provided by the StepLBTF, more evident.

**Reviewer** In the rest of the paper, and for all comparisons, the time step seems to remain equal to 0.001 s for every strategy (except convolution approach)/models. I fully understand this choice. However, what is the sensibility of the ARX models to the time step for a given viscoelastic law? How does such a sensitivity evolve when moving to one model to another? I think that a comment could help a non-specialist Reader and would value the work and the potential of the advanced strategy.

**Authors** Now regarding the ARX approach, there is no discussion to have because it is a model adapted to a given time step to produce the best predictive polynomial behavior of a given signal (no matter how the signal is obtained). Of course, as said in the text, this implies a constraint: the numerical simulations based on ARX have to use the same time step.

Going a little bit further:

* An ARX model, whatever its degree of refinement (the tuple na,nb,nk), is obtained by identification from calibration data with a given known sampling time step (whatever it is). The one chosen for this article (0.001, as explained above) also makes it convenient to provide representative execution times (more precisely a ratio between the different methods). As this value ensures a precise computation of the exact response with the Step-LBTF approach, it naturally secures good metrics for the parameters identification avoiding possible biased effects if the degree of the ARX model (the tuple na,nb,nk) is changed and then make the comparison between one and the other ARX model pertinent.

* It has to be recalled (and we insist more on this in the new version) that through the mathematical ARX model structure, the identified parameters encode both the time step used to produce the calibration data (even if they are of experimental origin) along with the physical parameters used in the “true” model (i.e. the moduli, the time constants, the non-integer exponents). It really behaves like a black-box model. To be clear, the parameters of the same ARX model structure (ARX130 for example) will be different if the calibration data are produced with a different time
step (or the experimental data recorded with a different sampling rate). But the behavior will be equivalently described (for a given input signal).

* If calibration data convey a rich frequency spectrum (like when the input signal corresponds to a continuous but not differentiable function) then a smaller time step is required. The time step of 0.001 is although a compromise to be able to capture the effect of strong discontinuities in the input signal like for case 1 (the crenel excitation).

**Reviewer** Page 17: Case 1 appears to be the best one for the calibration of the ARX model for the first fractional model. Why? Indeed, in other cases, it was the fourth case.

**Authors** This comment is in connection with the last point we discussed. It is important to apply the ARX model identification on calibration data that offers a rich informational content in terms of frequency (frequency spectrum of the input signal). Therefore, it is not surprising to observe that case 1 (crenel) or case 3 (successive ramps and dwell) are the best situations to use for calibrating the ARX model. To qualify which of these two cases is the most favorable could be possible thanks to fitting criteria but actually the small differences we observe do not make real sense. It is actually the user who has to select the ARX model according to the objective behind.

It is important to have in mind here that the smooth excitation was considered here for validating the ARX models mainly because it will correspond to the real behavior of a single voxel among many others involved in the description of a heterogeneous material. One cannot expect obviously a material voxel to behave in discontinuous manner in strain or stress.

**Reviewer** Minor corrections:

* The numbering of equations is not classical. Equations (1) to (4) then (5.XX) in the rest of the paper. To be homogenized.

* The notations for the ARX models must be homogenized throughout the paper (text, figures and tables): sometimes one has ARX(1,3,0), sometimes ARX130. The remark is also valid for ARX(2,4,0) and ARX240.

* Figures 2, 3, 6 and 7: please precise in the legends that strain excitation corresponds to the red line in the Figure, by for ex. "(Case X, red line)", or point the line by an arrow as in Figure 5.

* After Eq. (5.3): "All of them will participate to a validation of the incremental approach described next and to the understanding its underlying subtleties". The word "to" was missing, I guess.

* Table 1, last line, second column: a sign "-" to be suppressed at the end of the equation.

* Page 11: check the number of the subsection: 3.2 and not 3, I guess.

* Table3: "different sets" in the legend.

* After Equation (5.8): prefer $E^u = 1000$ MPa, $E^r = 100$ MPa than the formulation $E^u = 1000$, $E^r = 100$ in MPa. $\tau_{min} = 1 \times 10^{-2}$ s: please check the typo.

* Equation (5.10): define $E_{\alpha+1,\alpha+1}$ and $\Gamma$.

* 3 lines later: what is rhs?

**Authors** All minor corrections have been considered to produce the new version. rhs is the common abbreviation used in math for right-hand side of an equation.

# Review of version 2

Permalink: hal-03845394v2

## Reviewer 1 (Anonymous)

**Reviewer** The authors have very clearly answered the reviewers questions and made the necessary corrections and clarifications. Then, the paper can be published at JTCAM.

## Reviewer 2 (Anonymous)

**Reviewer** The Authors have precisely answered my questions and performed relevant corrections in the revised version. I recommend the publication of the paper in JTCAM.
Editor's assessment (Laurence BRASSARD)

The paper proposes a novel numerical strategy for time-integrating viscoelastic constitutive equations specified in terms of transfer functions in the Laplace domain. First, an incremental approach based on Laplace transform (step-LBTF approach) is presented to integrate the constitutive law over one time-step. Second, the step-LBTF approach is substituted by an ARX (AutoRegressive with exogeneous input) model, which has few parameters and much lower storage requirement. The accuracy and efficiency of the methods are tested for a material point for several constitutive models. The ARX method has potential to improve the precision and efficiency of finite element or spectral codes, which will be addressed in a forthcoming work. Both reviewers gave positive comments on the first version of the manuscript and mainly requested some clarifications. All comments were well addressed in the revised version and both Reviewers recommended the revised manuscript for publication.

**Open Access** This review 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.

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