Franck Renaud
email: franck.renaud@supmeca.fr

Jean-Luc Dion
email: jeanluc.dion@supmeca.fr

Gaël Chevallier
email: gael.chevallier@supmeca.fr

Imad Tawfiq

Rémi Lemaire

A new identification method of viscoelastic behavior : Application to the Generalized Maxwell Model

Keywords: Viscoelasticity, Identification, Generalized Maxwell, Fractional Derivatives, Transfer Function, Pole-Zero, Damping

This paper focuses on the Generalized Maxwell Model (GMM) identification. The formulation of the transfer function of the GMM is defined, as well as its asymptotes. To compare identification methods of the parameters of the GMM, a test transfer function and two quality indicators are defined. Then, three graphical methods are described, the enclosing curve method, the CRONE method and an original one. But the results of graphical methods are not good enough. Thus, two optimization recursive processes are described to improve the results of graphical methods. The first one is based on an unconstrained nonlinear optimization algorithm and the second one is original and allows constraining identified parameters. This new process uses the asymptotes of the modulus and the phase of the transfer function of the GMM. The result of the graphical method optimized with the new process is very accurate and fast.

Introduction

Many mechanical systems are damped with viscoelastic materials. This helps to avoid instabilities and to limit the levels of vibration. Although the viscoelastic behavior of materials is of great importance in order to obtain accurate results, the assumption of purely elastic materials is very commonplace for frequency analysis with Finite Element (FE) models. In order to carry out realistic Complex Eigenvalue Analysis (CEA) in dynamics, one needs to model viscoelasticity. Linear viscoelasticity has been described by many authors, for citing just a few of them: Ferry [START_REF] Ferry | Viscoelastic Properties Of Polymers[END_REF], Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF], Caputo and Mainardi [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF], Lakes [START_REF] Lakes | Viscoelastic Solids[END_REF], Chevalier and Vinh [START_REF] Chevalier | Mechanics of Viscoelastic Materials and Wave Dispersion[END_REF], or Balmès and Leclère [START_REF] Balmès | Viscoelastic vibration toolbox[END_REF]. Linear viscoelasticity assumes the existence of a relation between stress, σ i j , of a material and its strain, ε kl , history. Let us call h i jkl (t), the relaxation function which is also an element of the complex stiffness tensor. The Fourier transform is denoted with a hat.

               σ i j (t) = t -∞ h i jkl (t -τ) ε kl (τ) dτ σi j (ω) = ĥi jkl (ω) εkl (ω) (1)
Linear viscoelasticity is defined by the Relations [START_REF] Ferry | Viscoelastic Properties Of Polymers[END_REF]. In the present paper, identification will be conducted on ĥi jkl (ω) by considering a rheological model between σi j and εkl . Linear viscoelasticity of rheological models is written in terms of force and displacement: F(ω) = H(ω)X(ω), with H(ω), the complex stiffness of a viscoelastic material, equivalent to ĥi jkl (ω).

Nomenclature σ i j

The stress and σi j its Fourier transform. ε kl

The strain and εkl its Fourier transform. h i jkl

The relaxation function and ĥi jkl its Fourier transform. Tests on viscoelastic materials led by many authors exhibit strong stiffness frequency dependence on both modulus and phase (see Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF], Soula et al. [START_REF] Soula | Transient responses of polymers and elastomers deduced from harmonic responses[END_REF][START_REF] Soula | Measurements of isothermal complex moduli of viscoelastic materials over a large range of frequencies[END_REF]). For example, Figure 1 shows the stiffness of a multi layer material with layers of glue, rubber and steel. The non-Zero value of the stiffness phase of viscoelastic materials is well known, because it generates damping. But the stiffness modulus is also frequency dependent, hence, simple models, with constant modulus, such as constant complex modulus or modal damping are not sufficient enough for carrying out CEA. In order to show how six of these viscoelastic models behave, they are listed with their complex impedance formulas in Table 1. All these models are linear ones. The Figure 2 shows the modulus and the phase versus frequency of each model. The parameters used to produce this figure are given in Table 3 in Appendix A.

The constant complex modulus model is non-causal so it is only suitable in the frequency domain, Gaul et al. [START_REF] Gaul | Damping description involving fractional operators[END_REF]. The constant complex modulus is not a relevant model since its modulus is constant. The models of Maxwell and Voigt are efficient only on a small frequency range, but they are unrealistic, respectively, at low and high frequency, where their modulus is infinitely small and high. Moreover, a ninety-degree-of-phase material would tend towards a viscous fluid. 2

Linear Models

Symbol

Frequency representation Z(ω)

Constant complex modulus

K * K(1 + jη) Maxwell ( jωKC)/(K + jωC) Voigt K + jωC Generalized Maxwell K 0 + N i=1 jωK i C i K i + jωC i Generalized Fractional Calculus Maxwell γ 1 γ N K 0 + N i=1 ( jω) γ i K i G i K i + ( jω) γ i G i
Table 1: Rheologic representation of some common linear viscoelastic models

We would like to point out that Generalized Mawell Model is classically composed of Maxwell cells in parallel (see Koeller [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF]). With such definition this model is not able to display reversible creep (see Caputo and Mainardi [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF]). As this paper deals only with viscoelastic solids, Generalized Maxwell Model (GMM) and Generalized Fractional Maxwell Model would refer to a spring in parallel with respectively Maxwell cells and Fractional Calculus Maxwell cells (see Koeller [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF]). Thus the GMM defined here is the same as that used by Chevalier and Vinh [START_REF] Chevalier | Mechanics of Viscoelastic Materials and Wave Dispersion[END_REF] and the same as the Maxwell representation given by Caputo and Mainardi [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF], without the first dashpot. The fractional derivatives allow describing behaviours between purely elastic and purely viscous. Oldham and Spanier [START_REF] Oldham | The Fractional Calculus[END_REF] and Podlubny [START_REF] Podlubny | Fractional Differential Equations[END_REF] have described the mathematical properties of fractional derivatives. Podlubny [START_REF] Podlubny | Fractional Differential Equations[END_REF] has given some examples for viscoelasticity. Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF], Koeller [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF], Bagley and Torvik [START_REF] Bagley | On the fractional calculus model of viscoelastic behavior[END_REF], Gaul et al. [START_REF] Gaul | Damping description involving fractional operators[END_REF], Chevalier and Vinh [START_REF] Chevalier | Mechanics of Viscoelastic Materials and Wave Dispersion[END_REF] have studied the application of fractional derivatives in viscoelasticity. Lion [START_REF] Lion | Thermomechanically consistent formulations of the standard linear solid using fractional derivatives[END_REF][START_REF] Lion | On the thermodynamics of fractional damping elements[END_REF] has investigated rheological models incorporating "fractional damping elements", called spring-pots, from the point of view of thermodynamics. Heymans [START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF] has extended the fractional calculus to non-linear viscoelasticity. Koeller [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF] has defined the springpot as a rheological element whose stress is proportional to the fractional derivative of the strain.

Generalized Maxwell Model (GMM) and Generalized Fractional Calculus Maxwell Model (GFCMM) (see Table 1) seem to be very relevant. Firstly, their modulus is constant at low frequency, secondly, it increases on a frequency range and finally it is also constant at high frequency (see Figure 2). Moreover their phase is firstly nil at low frequency, secondly it is non-Zero on a frequency range and finally it is also nil at high frequency (see Figure 2). GFCMM use non-integer derivative orders. The use of a GFCMM in state space model requires that these derivative orders are fractions of integer. Both GMM and GFCMM allow one to describe behaviors between purely elastic and purely viscous materials (see Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF], Soula et al. [START_REF] Soula | Transient responses of polymers and elastomers deduced from harmonic responses[END_REF]). Figure 2 plots a transfer function of a GFCMM with one Pole-Zero couple with a derivation order of γ = 1/2. The GMM plotted in blue in Figure 3 has 3 couples of spring and dashpot Figure 2: Complex stiffness of the viscoelastic models described in Table 1. The parameters of these models are detailed in Table 3 in Appendix A.

in series. It allows describing the same behavior as the Pole-Zero Fractional Derivatives model. By increasing the number of spring and dashpot in series, the blue curve will tend towards the black one. Although the two models are relevant for viscoelasticity, the present paper focuses on GMM.

Given a viscoelastic material, one needs to identify the parameters of the GMM from experiments. Common tests can be sorted in 3 kinds of methods. The first ones are based on the Oberst's beam method. It consists in comparing the frequency response function of a structure with and without the viscoelastic material to be tested (see Oberst and Frankenfeld [17], Zhang and Richards [START_REF] Zhang | Parameter identification of analytical and experimental rubber isolators represented by maxwell models[END_REF], Castello et al. [START_REF] Castello | Constitutive parameter estimation of a viscoelastic model with internal variables[END_REF], de Lima et al. [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF]). Although this works well for low damping materials, this only allows one to measure values taken at the frequency of the first Eigenmodes of the beam. The second kind of tests consists in studying the creep and relaxation curves of materials (see Chen [START_REF] Chen | Determining a Prony Series for a viscoelastic material from time varying strain data[END_REF]). This is very efficient to get values at low frequency, when the material takes time to respond to the excitation. But to get high frequency values, one needs to assess a perfect unit step function when excitating the material, which is technically hard. The last kind of test is the Dynamic Mechanical Analysis. This consists in forcing oscillations of a material at a frequency and measuring its strain and its stress far from the resonance frequencies of the test bench. By sweeping in frequency, one is able to build the transfer function between stress and strain on a wide frequency range. This last kind of test is the most suitable for acquiring knowledge on the behavior of a material versus frequency.

Given an experimental transfer function characterizing a complex stiffness H(ω), the identification of the parameters of the GMM can be carried out with some graphical methods. Two of them were described by Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF] and Oustaloup [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF]. These methods are led on a Pole-Zero Formulations which are equivalent to the GMM. For example, the equivalence between the classical Standard Linear Solid 1 and a Pole-Zero Formulation with one Pole-Zero couple is demonstrated by Equations [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF].

                     Z(ω) = K 0 + jωK 1 C 1 K 1 + jωC 1 = K 0 1 + jω(C 1 /K 1 + C 1 /K 0 ) (1 + jωC 1 /K 1 ) = K 0 1 + jω/ω z,1 1 + jω/ω p,1 ω z,1 = C 1 K 1 + C 1 K 0 -1 = K 0 K 1 C 1 (K 0 + K 1 )
and

ω p,1 = K 1 C 1 (2)
1 Standard Linear Solid is also called Zener model. It is a simple case of Generalized Maxwell Model where N = 1.

The graphical identification methods use the characteristics of the asymptotes of Pole-Zero Formulations. The present paper presents an algorithm based on the principles of these graphical methods. This algorithm allows one to identify parameters from both the modulus and the phase curves with more efficiency than the classical graphical methods thanks to the optimization algorithm based on asymptotes.

Theoretical aspects

The Generalized Maxwell Model

Experimental transfer functions, characterizing complex stiffnesses, are well described by their modulus and phase. They are usually plotted versus the logarithm of the frequency. Moreover, the modulus is usually represented in logarithmic scale. In such scale, models with a product of terms, like Pole-Zero Formulations (PZF given by Equation ( 3)), become a sum of terms. This is why PZF are very suitable for the identification of parameters on transfer functions. The rheological formulation of GMM already presented in Table 1, is reminded in Equation (4).

Z(ω) = K 0 N i=1 1 + jω/ω z,i 1 + jω/ω p,i (3)
Z(ω) = K 0 + N i=1 jωK i C i K i + jωC i (4)
Equation ( 4) is very difficult to use in logarithmic scale. Fortunately, Relations (5) allow one to compute GMM parameters from the parameters of PZF (see Dion [START_REF] Dion | Modélisation et identification du comportement dynamique de liaisons hydro-élastiques[END_REF], Dion and Vialard [START_REF] Dion | Identification or rubber shock absorber mounts[END_REF]). The proof of these relations is given in Appendix B.

                     K i = K 0 N h=1 ω p,h ω z,h ω p,i -ω z,h ω p,i + ω p,h (δ ih -1) C i = K i ω p,i (5)
Given a transfer function of GMM, Z(ω) = |Z(ω)| exp( jϕ(ω)), the modulus and the phase of the associated PZF are defined in Equation [START_REF] Balmès | Viscoelastic vibration toolbox[END_REF].

                             |Z(ω)| = K 0 N i=1 |Z(ω)| i = K 0 N i=1 1 + ω/ω z,i 2 1 + ω/ω p,i 2 ϕ(ω) = N i=1 ϕ i (ω) = N i=1 tan -1 ω ω z,i -tan -1 ω ω p,i (6)
It is worth noticing that PZF can be generalized to fractional calculus. Then, Equation (3) becomes Equation [START_REF] Soula | Transient responses of polymers and elastomers deduced from harmonic responses[END_REF]. The modulus and phase of this Fractional Calculus Pole-Zero Formulation (FCPZF) is given by Equation [START_REF] Soula | Measurements of isothermal complex moduli of viscoelastic materials over a large range of frequencies[END_REF]. Thus, PZF is a particular case of the FCPZF, where the derivation order is γ i = 1. One could expect FCPZF to be equivalent to GFCMM of Table 1, but the link between these two formulations has not been established yet.

Z(ω) = K 0 N i=1 1 + jω/ω z,i γ i 1 + jω/ω p,i γ i (7)                                              |Z(ω)| = K 0 N i=1 1 + 2 ω/ω z,i γ cos (γπ/2) + ω/ω z,i 2γ 1 + 2 ω/ω p,i γ cos (γπ/2) + ω/ω p,i 2γ ϕ(ω) = N i=1 tan -1 ω/ω z,i γ sin (γπ/2) 1 + ω/ω z,i γ cos (γπ/2) -tan -1          ω/ω p,i γ sin (γπ/2) 1 + ω/ω p,i γ cos (γπ/2)          (8)
2.2. Asymptotes of the PZF Relations ( 6) are non-linear compared with frequency and make the identification difficult to carry out. Fortunately, it is possible to approach the behavior of PZF using the asymptotes of the modulus and the phase (see Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF], Oustaloup [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF] and Dion [START_REF] Dion | Modélisation et identification du comportement dynamique de liaisons hydro-élastiques[END_REF][START_REF] Dion | Identification or rubber shock absorber mounts[END_REF]). The authors want to point out that thanks to Relation (6

): log(|Z(ω)|) = log(K 0 ) + N i=1 log(|Z(ω)| i ) and ϕ(ω) = N i=1 ϕ i (ω).
In statics, when the angular frequency is nil, log (|Z(ω = 0)|) = log (K 0 ) and ϕ (ω = 0) = 0. In dynamics, the behavior of PZF is the superposition of the behaviors of Pole-Zero couples. Each of these couples adds a non-linear contribution compared with frequency to the modulus and to the phase (see Equations ( 9)).

                     log |Z(ω)| i = 1 2 log        1 + ω ω z,i 2        - 1 2 log        1 + ω ω p,i 2        ϕ i (ω) = tan -1 ω ω z,i -tan -1 ω ω p,i (9)
At this point, we would like to define χ = log(ω), the logarithm of the angular frequency. Assuming an angular frequency far from the Zero and the Pole and taking the limit, Relations (9) tend towards the asymptotes given by Equations [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF]. These asymptotes are shown in Figure 3. The curves generated by a Pole-Zero couple, as well as their asymptotes, are symmetrical around the medium angular frequency,

χ c,i = log( √ ω z,i ω p,i ) = (χ z,i + χ p,i )/2. log(|Z(χ)| i ) -→ log(|Z asy (χ)| i ) =          0 if χ < χ z,i < χ p,i log(ω/ω z,i ) = χ -χ z,i if χ z,i ≤ χ ≤ χ p,i log(ω p,i /ω z,i ) = χ p,i -χ z,i if χ z,i < χ p,i < χ ϕ i (χ) -→ ϕ asy i (χ) =          0 if χ < χ z,i < χ p,i π/2 if χ z,i ≤ χ ≤ χ p,i 0 if χ z,i < χ p,i < χ (10)
It is worth noticing that the asymptote of the modulus is the primitive of the asymptote of the phase divided by π/2 (see Equations [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF] and Figure 3). Considering the analogy between the asymptote of the phase and the derivative of the asymptote of the modulus with respect to χ, we would like to define the stiffening, S , of a modulus and the area, A, under a phase curve. Given two angular frequencies, ω a and ω b , S and A are defined thanks to Relations [START_REF] Oldham | The Fractional Calculus[END_REF]:

∀ χ a = log(ω a ) ∈ R and χ b = log(ω b ) ∈ R,                    S (|Z|, χ a , χ b ) = χ b χ a ∂ ∂χ log(|Z(χ)|) dχ = log(|Z(χ b )|) -log(Z|(χ a )|) The stiffening A(ϕ, χ a , χ b ) = χ b χ a ϕ(χ) dχ
The area under the phase curve (11) 6), ( 9) and ( 10)). The asymptotes are plotted in dashed line.

The analogy between the phase and the derivative of the modulus allow one to get a useful relation between the stiffening and the area of the asymptotes of PZF (see Equation ( 12)). Moreover, given two angular frequencies logarithms χ a < χ z,i and χ b > χ p,i , Equations [START_REF] Podlubny | Fractional Differential Equations[END_REF] becomes S (|Z asy (χ

)| i , χ a , χ b ) = 2/π×A(ϕ asy i (χ), χ a , χ b ) = (χ p,i -χ z,i ). ∀ χ a ∈ R and χ b ∈ R, S (|Z asy (χ)| i , χ a , χ b ) = 2 π A(ϕ asy i (χ), χ a , χ b ) (12)
The asymptotes present the advantage to be easy to manipulate and to be simultaneously representative of modulus and phase. Asymptotes allow one to perform simple graphical identification.

Existing methods

The Enclosing Curve Method

Vinh [START_REF] Vinh | Sur le passage du régime harmonique au régime transitoire viscoelastique[END_REF] has described a graphical method based on asymptotes (see Figures 4 and7). This consists in plotting 2 enclosing curves above and below the experimental modulus curve, here : |H(ω)| + = |H(ω)| × 10 0.06 and |H(ω)| -= |H(ω)| /10 0.06 . Then the 2 enclosing curves are linked alternatively by 2 types of curves: constant ones ( f = a) and linear one with a unitary slope ( f = log(ω) + b). The Zeros are given by the angular frequency where the linear segments cross the bottom curve and the Poles are given by the angular frequency where the linear segments cross the top curve. The closer to modulus are the enclosing curves, the better the model fits the modulus curve and the more it requires Pole-Zero couples. The parameters obtained with the enclosing curves method are listed in Table 4 in Appendix C. Each linear segment is in fact the asymptotic behavior of a Pole-Zero couple. This old method shows how powerful it is to use asymptotic curves in the logarithmic scale instead of the non-linear modulus and phase relations in the linear scale. Oustaloup [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF] has proposed the same kind of method for the phase curve. His method consists in equalizing the area under the experimental curve, ψ(χ), and the area between the model curve, ϕ(χ). [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF]. The areas under the asymptotes curve and under the test curve are equal. The width α i give the distance between poles and zeros.

The CRONE Method

The CRONE Method is described by Oustaloup [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF] in his book. First, the number of Pole-Zero couples is defined to have at least one Pole-Zero couple per decade in the frequency domain, here we have 8 couples for 5 decades. The frequency domain is divided in frequency subdomains, χ ∈ χ a,i ... χ a,i+1 , one for each Pole-Zero couple, (see Figure 5). The subdomains have the same length, i . In each subdomain, the area under the test curve is approached by the area of a rectangle whose height is the first phase value of the subdomain, ψ i . It has to be equal to the area under the asymptotic curve of the model phase. Zeros are arbitrary placed on the left of subdomains, thus, χ z,i+1 = i + χ z,i . Poles are placed in such a way that χ p,i = χ z,i + α i , with α i = ψ i i × 2/π. This method generates edge effects with a decreasing phase near the edge. To avoid it, the first Zero and the last Pole are corrected by an arbitrary value around 0.08 : χ corrected z,1 = χ z,1 -0.08 and χ corrected p,end = χ p,end + 0.08. This helps to increase the phase value of the first and the last Pole-Zero couple. Then in Figure 5 the asymptotic area of the first subdomain begins on the left of it. As it is shown in Figure 5, this method is not very accurate for non smooth phase curves. The parameters of the method used to obtain Figures 5 and8 are listed in Table 5 in Appendix C.

A test transfer function

In order to benchmark the methods presented in this paper, it is useful to define a common transfer function with modulus and phase, H(χ) = |H(χ)| exp( jψ(χ)). For this purpose we have chosen and arbitrary Pole-Zero function with big variations in order to test the robustness of the methods. The logarithm of the Zeros and Poles are presented in Table 2, and K 0 = 1000.

It is worth noticing that a Pole-Zero couple with the Zero bigger than the Pole is not causal. Such a couple generates a decreasing modulus with a negative phase, what would mean that the material respond before the excitation. One can expect an identification method of viscoelastic properties to always give causal Pole-Zero couples, i.e. with a Zero lower than the Pole. In order to test the robustness of the methods with the test transfer function, its 4th Pole-Zero couple is not causal (see Table 2). Nevertheless, the transfer function, which is the sum of all couples behavior, is still globally causal with an increasing modulus and a positive phase (Figure 6). It is worth noticing that for all methods, experimental curves are usually interpolated exponentially in the frequency domain for the modulus and linearly for the phase. This helps to treat values between experimentally measured points. The present test function is defined with 1000 points equally spaced in logarithmic scale between 10 2 and 10 6 rad/s. The subscript k will be used to refer 2).

to one of these points.

Two fitting indicators

To compare the methods presented in this paper, two curves fitting indicators are defined, one for the modulus and one for the phase. The purpose is to quantify how good a model Z(χ) = |Z(χ)| exp( jϕ(χ)), is able to fit an experimental transfer function H(χ) = |H(χ)| exp( jψ(χ)). To be equivalent, these indicators must respect the analogy between the stiffening of the modulus asymptote and the area under the asymptote of the phase see Relation [START_REF] Podlubny | Fractional Differential Equations[END_REF].

δ |Z| = χ max χ min ∂ ∂χ log (|H(χ)|) - ∂ ∂χ log (|Z(χ)|) 2 dχ χ max χ min ∂ ∂χ log (|H(χ)|) 2 dχ ( 13
)
δ ϕ = χ max χ min (ψ(χ) -ϕ(χ)) 2 dχ χ max χ min ψ(χ) 2 dχ ( 14
)
Taking advantage of the definition of the stiffening as an integral in Equation ( 11), the indicators are given in Equations ( 13) and [START_REF] Lion | Thermomechanically consistent formulations of the standard linear solid using fractional derivatives[END_REF]. They are normalized to be dimensionless. A model curve that fits really well an experimental curve will exhibit a low indicator. These indicators only allow evaluating the identification of Poles and Zeros, but do not evaluate the identification of the static modulus, K 0 . Thus, even if a GMM has low indicators, an offset can occur between the modulus of the GMM and the modulus of the experimental curve. Indeed, the offset between the two modulus curves is driven by the static modulus: K 0 . This static modulus is chosen after the identification of poles and zeros. One can either set it before the identification and deals the offset with the first zero (see Section 4.2.2).

As the indicators of these two methods show (Figures 7 and8), the enclosing curves method gives better results than the CRONE method. But these results can be improved, furthermore, none of both methods allow identifying parameters from both modulus and phase curves.

The Asymptotical Variations Method

Figure 9 represents the identification process. This begins by an initialization of the parameters of the PZF. Three initialization methods are presented in this paper and they are detailed in Section 3 for the two firsts and in Section 4.1 for the last one. Given the initial parameters of the PZF, optimization processes allow one to improve the fitting of experimental curves. The indicators are used to quantify the quality of the identified parameters.

Initialization of the PZF parameters with a new graphical method

Given an experimental transfer function with modulus and phase, one can expect to use both of them for the identification of GMM parameters. The previous methods have already shown that these two kinds of curves content the same information. The next methods presented in this paper are described for the modulus and for the phase, then to manage the two kinds of value, the Poles and the Zeros will be the averaged with a weighting coefficient r :

ω z,i = (1 -r) × ω modulus z,i + r × ω phase z,i
and ω p,i = (1r) × ω modulus p,i + r × ω phase p,i . Here, the weight r will always be taken to 0.5. Moreover, in order to prevent edge effects, the test curves are extrapolated with a constant modulus and a null phase on about half a decade before and after the frequency domain.

Unlike the previous method already described by Oustaloup, the area under the phase is evaluated by the Trapezes Method (see Figure 10). This helps to give better results when the transfer function is only known in a few points. Then Pole-Zero couples are placed according to the need of stiffening and phase. This is realized by dividing the frequency domain in N frequency subdomains, one for each Pole Zero couple. Each of these subdomains, χ ∈ χ a,i ..χ b,i , has the same experimental stiffening

α i = χ b,i χ a,i d dχ [log(|H(χ)|)]dχ = log(|H(χ b,i )|) -log(|H(χ a,i )|
), (see the modulus curve in Figure 10) or the same area of phase, α i π/2 = χ b,i χ a,i ψ(χ)dχ, (see the phase integral curve in Figure 10). Thus, if high stiffening or high phase occurs in the vicinity of a frequency, then couples will concentrate around this frequency. In a subdomain, the medium frequency of a couple χ c,i = (χ z,i + χ p,i )/2, is placed with respect to this principle. The medium frequency divides the subdomain in two parts of equal stiffening or area of phase. In Figure 10, the stiffening or phase integral apportionment is represented by horizontal lines and the frequency domain apportionment by vertical lines, continuous lines separating subdomains and dashed lines denoting medium frequencies.

As in the previous method, the ratio α i is used to place the Zeros and the Poles, but here it is defined for both modulus and phase. The factor α i is the value calculated in a subdomain and it affects only the Pole Zero couple of this subdomain. Because of the frequency domain apportionment, here α i is the same for all subdomains. Unfortunately, there are some cases where curves with big variations can produce very smooth fitted model. To avoid these cases, α i is weighted by a coefficient ν i which depends on the length of the subdomains. A short subdomain means that the phase or the slope of the stiffness is big, then a big ν i coefficient is required, Equation ( 15):

                         total = χ max -χ min i = χ b,i -χ a,i ν i = total -i total N N -1 =⇒                χ z,i = χ c,i - ν i α i 2 χ p,i = χ c,i + ν i α i 2 (15)
Notice that N i=1 ν i = N, and that when all the subdomains have the same length, the weighting coefficients are ν i = 1. The results of this method are shown in Figure 10. The blue curves are produced by a fitting on the test modulus curve, the Pole Zeros coefficients of such a model are shown in Table 6 in Appendix C. Table 7 lists the Pole Zeros coefficients of the model produced by a fitting on the test phase curve (red curves). The static modulus, K 0 , in Equation ( 3) is taken as the first value of the test modulus curve. Notice that all Zeros are lower than their associated Poles. Finally, the average Pole Zeros coefficients are shown in Table 8, they lead to the model plotted in Figure 11.

This method allows improving the placement of the couples. It is better than the CRONE method but it is worse than the enclosing curve method which intrinsically places well couples. Nevertheless, this method uses both the modulus and the phase. The values of the convergence indicators are lower for a fitting on the experimental modulus : 6.3% and 1.7% than for a fitting on the experimental phase: 12.31% and 3.8%. The averaged parameters score intermediate indicators: 8.87% and 2.57%. These indicators value are not so good, this means that the model is not close enough to the test curves especially for the test curve of phase (see Figure 11). The results of this method can be optimized and they will be used as an initialization for the optimization processes proposed in the following sections.

Optimization of identified parameters

Given a GMM with initial values for Poles and Zeros, some ways to optimize the fitting of experimental curves are described in this section. For this purpose, iterative methods are required. One can use the classical Newton-Raphson 8 and Figure 11) are obtained by averaging the ones from the modulus curve (see Table 6) and from the phase curve (see Table 7). method to find the parameters leading to the least mean square value. This method does not work very well because of the existence of several local minima.

To overcome this difficulty, two optimizations processes are presented here. First, a Matlab function is used, it is based on the Nelder-Mead Simplex Method. Second, a new method based on the modulus and phase asymptotes is described. These asymptotes are used to translate the differences between the model and the experimental curves into step of variations.

An unconstrained nonlinear optimization algorithm

"fminsearch" is a Matlab function based on the Nelder-Mead Simplex Method which convergence properties have been discussed by Lagarias and al. [25]. This method finds the minimum of a scalar-valued nonlinear function W of n real variables, the Zeros and the Poles. For the present identification problem, this function will be defined as the weighted sum of the 2 indicators previously defined in Equations ( 13) and ( 14) (see Equation ( 16)). This allows fitting both modulus and phase in the same time:

W = (1 -r) × δ |Z| + r × δ ϕ (16)
Although in general cases, taking the weighting coefficient r 0.5 gives better results, in the present paper, r has been taken to 0.5 for all methods in order to compare them. This was already the case in Section 4.1. The simplex method will be used to minimize W. The calculation time used by this method seems to be strongly dependent of the initial parameters. These initial parameters where found by the method described in Section 4.1 and are listed in Table 8 in Appendix C.

According to the indicators : 0.11% and 0.03%, this method gives very good results (Figure 12). The parameters of this method are reported in Table 9 in Appendix D. Since it is an unconstrained method, non causal Pole-Zero couples can occur like the fourth couple of Table 9 for which its Zero is greater than its Pole. Despite its good results, this method remains unconstrained and can lead to non causal solutions. Figure 12: The GMM obtained after the optimization of parameters given by the method described in Section 4.1 (see Table 8 in Appendix C). The optimization is proceed thanks to a Matlab function based on the Nelder-Mead Simplex Method. The optimized parameters are detailed in Table 9 in Appendix D. The indicators of such parameters are δ |Z| = 0.11% and δ ϕ = 0.03%. Notice that the fourth Pole is lower than the fourth Zero, thus it is non-causal.

Optimization with the Asymptotical Variations Method

The method described here allow constraining Pole Zero couples to be causal, i.e. with the Zero smaller than the Pole. The principle of this method is to find each local difference between the experimental curves and the model curves in terms of stiffening or phase integral. These local differences are attributed to all different couples in function of their local weight. Local lack and excess of stiffening and integral of phase can be seen as a gradient integrated on a variable step. The asymptotes provide the simple relations for the translation of this variable step into parameter variations.

To compute the stiffening and the phase integral, the experimental curves are interpolated in logarithmic scale thanks to linear functions. After the interpolation, curves reach a rate of 100 points per decade. The subscript k will be used to denote these points. The interval between two successive points is a subdomain where integration is done with the Trapezes Method. In Figure 13, Figure 15 and Figure 17, the interpolated modulus and phase are represented with a wide line.

While identifying the modulus, the relevant quantity to look at is the difference between the modulus slope of the experimental curve,

χ k+1 χ k d dχ [log(|Z(χ)|)]dχ = log(|Z(χ k+1 )|) -log(|Z(χ k )|).
The local stiffening difference, δS k , is defined by Equation [START_REF] Oberst | Damping of the bending vibrations of thin laminated metal beams connected through adherent layer[END_REF]. In the upper part of Figure 13, the model lack and excess of local stiffening are represented in painted area.

δS k = log |Z(χ k+1 )||H(χ k )| |Z(χ k )||H(χ k+1 )| (17)
Relation [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF], describing the behavior of asymptotes of the model is used to convert the local differences of stiffening into parameter variations steps. These variations steps are shared between all couples in agreement with their local weight in each subdomain. The local weight of a couple, I |Z| k,i is defined as the ratio of its local stiffening contribution to the model's local stiffening (Equation ( 18)). The local stiffening contribution of the i th couple is log(|Z(χ k+1 )| i )log(|Z(χ k )| i ), this is the difference between 2 successive points of the modulus logarithm generated by 9)).

I |Z| k,i = log(|Z(χ k+1 )| i ) -log(|Z(χ k )| i ) log(|Z(χ k+1 )|) -log(|Z(χ k )|) (18)
Figure 13 shows, in its upper part, the contribution of each particular couple to the derivative of the modulus, 14. Notice that for each couple, the peak of influence occurs in the vicinity of its medium frequency, χ c,i . These medium frequencies are represented in the lower part of Figure 14 with the local stiffening difference, δS k (Equation ( 17)). In each subdomain, this difference already shared between couples according to their influence must be shared again between the Pole and the Zero of a couple. The sum of stiffening difference occuring in subdomains on the left of the medium frequency is attributed to the Zero and the sum of the ones of the right is attributed to the Pole (see Equation [START_REF] Castello | Constitutive parameter estimation of a viscoelastic model with internal variables[END_REF]). Thus at each iteration step denoted by the superscript u, the parameters are changed according to the lack and excess of stiffening.

                   χ modulus z,i = χ u z,i - χ k ≤χ c,i I |Z|,u k,i δS u k χ modulus p,i = χ u p,i + χ c,i <χ k I |Z|,u k,i δS u k (19)
The way to identify the phase curve is analogous. The relevant quantity to look at is the difference between the phase of the experimental curve, ψ(χ), in continuous line in Figure 15 and the phase of the model curve, ϕ(χ), in dashed line in Figure 15. Integrating this phase between 2 points, ω ∈ [ω k ..ω k+1 ], provides the local area of phase :

χ k+1 χ k ϕ(χ)dχ = (ϕ(χ k+1 ) + ϕ(χ k ))/2.
The local area of phase difference, δA k , is defined by Equation [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF]. The model lack and excess of area of phase are represented in painted area.

δA k = ψ(χ k+1 ) + ψ(χ k ) -ϕ(χ k+1 ) -ϕ(χ k ) 2 (20)
Relation [START_REF] Koeller | Applications of fractional calculus to the theory of viscoelasticity[END_REF], describing the behavior of asymptotes of the model will be used to convert the local differences into parameter variations step. These variations steps are shared between all couples in agreement with their local k,i is defined as the ratio of its local area of phase contribution to the model's local area of phase (Equation ( 21)). The local area of phase contribution of the i th couple is (ϕ i (χ k+1 ) + ϕ i (χ k ))/2, this is the integral of phase between two successive points of the phase generated by the i th couple (Equation ( 9)).

I ϕ k,i = ϕ i (χ k+1 ) + ϕ i (χ k ) ϕ(χ k+1 ) + ϕ(χ k ) (21)
Figure 15 shows, in its upper part, the phase contribution of each particular couple ϕ i (χ) and, in the lower part, the model phase, ϕ(χ), as the sum of all couple contributions. The influence curves of Pole Zero couple are shown in the upper part of Figure 16. Notice that for each couple, the peak of influence occurs in the vicinity of its medium frequency, χ c,i . These medium frequencies are represented in the lower part of Figure 16 with the local area of phase difference, δA k (see Equation [START_REF] De Lima | An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring[END_REF]). In each subdomain, this difference already shared between couples according to their weight must be shared again between the Pole and the Zero of a couple. The sum of area of phase difference occuring in subdomains on the left of the medium frequency is attributed to the Zero and the sum of the ones of the right is attributed to the Pole (see Equation [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF]). Thus at each iteration step denoted by the superscript u, the parameters are changed according to the lack and excess of area of phase.

                       χ phase z,i = χ u z,i - 2 π χ k ≤χ c,i I ϕ,u k,i δA u k χ phase p,i = χ u p,i + 2 π χ c,i <χ k I ϕ,u k,i δA u k ( 22
)
The sharing of differences between couples in function of their weight has the advantage to only affects couples close to a local deviation, thus it does not disturb couples without any weight. The sharing of differences between Pole and Zero, in function of the position of a subdomain to the medium frequency, allows two different kinds of parameters behavior during iteration. First, the medium frequency of a Pole Zero couple can move in frequency without any change of the relative position of the Pole to the Zero. It seems that the whole couple is moving in frequency. The second behavior has a standstill medium frequency while the coefficient Pole and the Zero are changing their value. This behavior induces the growth or the decrease of the whole couple. Figure 17: The GMM obtained after the optimization of parameters given by the method described in Section 4.1 (see Table 8 in Appendix C). The optimization is proceed thanks to the method, described in Section 4.2.2, using asymptotes to translate model deviations from the test curves into parameters steps. The optimized parameters are detailed in Table 10 in Appendix D. The indicators of such parameters are δ |Z| = 0.77% and δ ϕ = 0.15%.

The parameters found according to modulus and phase are averaged as in the previous method: ω z,i = (1r) × ω modulus z,i + r × ω phase z,i , with r = 0.5 and the static modulus, K 0 , defined in Equation ( 3) is taken as the first value of the test modulus curve. The method described here for the modulus allows fitting its derivative. Thus, although the global form of the modulus is respected, an offset always occurs between the model and the test function. To remove this offset, the first Zero is modified according to the mean distance between the 2 curves :

χ u+1 Z 0 = χ u Z 0 -1 N V k=1 log(|H(χ k )|) -log(|Z u (χ k )|).
Indeed, this means distance can be seen as a stiffening. Moreover, the causality is verified at each iteration. If the Zero of one couple is bigger than its Pole then their values are forced, for example, χ p,i -χ z,i = 10 -6 centred on the medium frequency. Notice that such a couple is superfluous for the model since it does not produce any phase or stiffening. Finally, the problem of edge effects have already been treated by the extrapolation of tests curves. If for some reasons the medium frequency of a couple would step out the frequency domain, a virtual medium frequency would be defined inside the frequency domain in order to allow the computation of Equations ( 19) and [START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF].

The iterative process stops when the maximum parameter step, δ u ZP max = max(|χ u+1 z,i -χ u z,i |, |χ u+1 p,i -χ u p,i |) is smaller than a defined value, ε = 10 -3 for example, or after a certain amount of iteration, 125 for example. For robustness, the maximum parameter step is imposed to be smaller than the previous one : δ u ZP max ≤ δ u-1 ZP max . If it is not the case all parameter steps are normalized by multiplying per δ u-1 ZP max /δ u ZP max . During the first iteration, this step limitation is initialized to a certain value, 0.25 for example.

Conclusion and perspectives

Generalized Maxwell Model (GMM) can be expressed with a Pole Zero Formulation (PZF). In logarithmic scale this formulation is a superposition of Pole Zero couple behavior. These behaviors are non-linear against frequency and it is easier to manipulate the asymptotes to approach the value of stiffening and integral of phase. The new identification method presented in this paper is split in 2 parts, initialization of parameters and optimization.

The initialization places the medium frequencies of Pole-Zero couples in the frequency domain such that they are equally distant according to stiffening and integral of phase. Thus, the couples are placed where they are most needed. The stiffening and area of phase are translated into parameter steps: distance between Pole and Zero of each couple, thanks to the relations of asymptotes. This gives the initialized parameters of the GMM. Then they are optimized by considering each local deviation between the model and the experimental curve. Each local deviation is translated into parameter step thanks to asymptotes. The steps are shared between couples in function of their local weight. Then they are shared between Poles and Zeros in function of the position of the local deviation compared to the medium frequency of the couple. Iterating until the parameters steps become too small provides the optimized parameters. This new method is very efficient for the identification of a Pole Zero Formulation according to both modulus and phase curves. It is neither restricted to viscoelasticity nor to the GMM. A very analogous version of this method could be applied to identify the parameters of a Fractional Calculus Pole Zero Formulation (FCPZF). Indeed the only difference would be the asymptotic value. Given a derivation order, γ i and considering Equation ( 8), the stiffening and the integral of phase of one Pole-Zero couple's asymptotes become for a FCPZF, respectively, S i = γ i × (χ p,i -χ z,i ) and A i = γ i × π/2 × (χ p,i -χ z,i ).

The advantage of this method is the possibility to constrain the parameters. In the present case of viscoelasticity, for the seek of causality, Zeros were constrained to be smaller than their Pole. Although this method is very quick and efficient, some improvements would still be possible. It is possible that the Zero of a couple would become close but different to its Pole. In such a situation the contribution of this couple is very low inducing very low step variation because of the weighting parameters. That means that it will be difficult to modify this kind of couple. This problem can be solved by adding a constant to the weighting coefficient, I |Z| k,i and I ϕ k,i . This constant would be equal for all couples, for example: ζ = 1/N. Other improvements can be done when, for some reasons, a Zero tends to be equal to a Pole. In this case the behavior of the Zero is exactly cancelled by the behavior of the Pole. In that case, the number of Pole Zero couples to identify can be reduced. One other possibility is to suppress this couple and to replace it by another one located in the frequency region with the highest deviations.

B. The parameters of GMM in function of the parameters of PZF

The proof of Relations ( 5) is given below. It allows linking a Pole-Zero Formulation to a Generalized Maxwell Model. Starting from a Pole-Zero Formulation:

Z( jω) = K 0 N i=1 1 + jω/ω z,i 1 + jω/ω p,i = K 0 jω        N h=1 ω p,h ω z,h               1 jω N i=1 jω + ω z,i jω + ω p,i       
Z( jω) can be expressed as a sum of fractions:

Z( jω) = K 0 jω        N h=1 ω p,h ω z,h               A 0 jω + N i=1 A i jω + ω p,i       
Considering the following relation:

A 0 jω + N i=1 A i jω + ω p,i = 1 jω N i=1 jω + ω z,i jω + ω p,i
The coefficient A 0 is obtained by multiplying the previous relation by jω and taking jω = 0. The coefficients A i are obtained by multiplying the previous relation by ( jω + ω p,i ), and taking jω = -ω p,i :

A 0 = N h=1 ω z,h ω p,h
and

A i = N h=1 ω z,h -ω p,i ω p,h (1 -δ ih ) -ω p,i ∀i ∈ [1..N] (23)
By replacing Equations ( 23) into the Pole-Zero Formulation, one obtains :

Z( jω) = K 0 jω        N h=1 ω p,h ω z,h               A 0 jω + N i=1 A i jω + ω p,i        = K 0 + K 0        N h=1 ω p,h ω z,h        N i=1 jω A i jω + ω p,i
Then by comparing the above relation with the Generalized Maxwell Model, one is able to identify the coefficients K stk , K i and C i , see Equation [START_REF] Dion | Identification or rubber shock absorber mounts[END_REF]:

Z(ω) = K stk + N i=1 jωK i C i K i + jωC i
Thus, Equations (5) are retrieved:

                                 K stk = K 0 K i = K 0 N h=1 ω p,h ω z,h ω z,h -ω p,i ω p,h (1 -δ ih ) -ω p,i C i = K i ω p,i (24)

ω

  Angular frequency. χ = log(ω) Logarithm of the angular frequency. H(ω) = |H(ω)| × exp( jψ(ω)) Measured impedance of the material. A test function in this paper. Z(ω) = |Z(ω)| × exp( jϕ(ω)) Impedance of the identified model. |Z u (χ k )| Modulus of the model, during the u th iteration taken at the angular frequency 10 χ k of the k th point of the modulus curve. With k ∈ [1..V], and u ∈ [1..U]. ϕ u (χ k ) Phase of the model, during the u th iteration taken at the angular frequency 10

Figure 1 :

 1 Figure 1: Complex shear stiffness of a multilayer viscoelastic material at ambiant temperature in the [200-3500Hz] frequency range. The top graph represents the normalized modulus H( f )/H( f = 200Hz) (dimensionless) and the bottom graph represents the phase ψ (in degrees).

Figure 3 :

 3 Figure 3: The modulus and the phase of a PZF with only one one Pole-Zero couple (see Equations (6), (9) and (10)). The asymptotes are plotted in dashed line.

Figure 4 :

 4 Figure 4: Enclosing Curves Method allows identifying Poles and Zeros thanks to the enclosing curves of the modulus. The asymptotes are plotted between the two enclosing curves.

Figure 5 :

 5 Figure 5: CRONE's graphical method[START_REF] Oustaloup | La commande CRONE : commande robuste d'ordre non entier[END_REF]. The areas under the asymptotes curve and under the test curve are equal. The width α i give the distance between poles and zeros.

Figure 6 :

 6 Figure 6: The modulus and the phase of the test transfer function used to compare identification methods. Particular attention should be payed on the fourth Pole Zero couple which has been arbitrary taken to be non-causal (see Table2).

Figure 7 :

 7 Figure 7: Enclosing Curves Method. Pole Zero couples are detailed in Table 4 in Appendix C. The indicators of such parameters are δ |Z| = 4.5% and δ ϕ = 1.25%

Figure 8 :

 8 Figure 8: CRONE's graphical method [22]. Pole Zero couples are detailed in Table 5 in Appendix C. The indicators of such parameters are δ |Z| = 23.48% and δ ϕ = 8.24%.

Figure 9 :

 9 Figure 9: The iterative identification process with three different methods for initialization and two methods for optimization.

  Figure10: A new graphical method using asymptotes to identify parameters on both modulus and phase curve. The parameters are detailed in Appendix C. The final parameters (see Table8and Figure11) are obtained by averaging the ones from the modulus curve (see Table6) and from the phase curve (see Table7).

Figure 11 :

 11 Figure 11: The final results of the graphical method described in Section 4.1. The parameters are in detailed in Table 8 in Appendix C. They have been averaged from the ones plotted on Figure 10. The indicators of such parameters are δ |Z| = 8.87% and δ ϕ = 2.57%

  d dχ [log(|H(χ)|)], in dashed line in Figure 13 and the modulus slope of the model, d dχ [log(|Z(χ)|)], in dashed-dotted line in Figure 13. Integrating this modulus slope in between 2 points, ω ∈ [ω k ..ω k+1 ], provides the local stiffening :

Figure 13 :

 13 Figure 13: Modulus of the test function and its derivative compared to the model's ones with the couple contributions during the 1st iteration

Figure 14 :

 14 Figure 14: Sharing between the Poles and the Zeros of the modulus derivative difference between the model and the test function

  d dχ [log(|Z(χ k )| i )], and, in the lower part, the derivative of the modulus, d dχ [log(|Z(χ k )|)], as the sum of all couple contributions. The influence curves of each Pole Zero couple are shown in the upper part of Figure

Figure 15 :

 15 Figure 15: Phase of the test function compared to the model's phase with the couple contributions during the 1st iteration

Figure 16 :

 16 Figure 16: Sharing between the Poles and the Zeros of the phase difference between the model and the test function

  χ k of the k th point of the phase curve. With k ∈ [1..V], and u ∈ [1..U]. The Zero of the i th Pole-Zero couple during the u th iteration and χ u z,i = log ω u z,i its logarithm. With i ∈ [1..N].

	η	Damping rate.
	K 0	Static stiffness taken at ω = 0, i.e. t = +∞.
	K i	Stiffness of the i th spring.
	C i	Damping of the i th dashpot.
	G i	Coefficient of the i th spring-pot.
	γ i	Fractional derivative order of the i th Pole-Zero couple.
	ω u z,i
	ω u p,i	A Pole and χ u p,i = log ω u p,i its logarithm.
	ω u c,i
	δ i j	Kronecker's symbol.
	A	Area under the phase curve.
	S	Stiffening of the impedance modulus.
	δ |Z|	Modulus convergence indicator.
	δ ϕ	Phase convergence indicator.
	r	Weighting coefficient.
	Acronyms
	GMM	Generalized Maxwell Model
	PZF	Pole-Zero Formulation
	GFCMM Generalized Fractional Calculus Maxwell Model
	CEA	Complex Eigenvalue Analysis

= ω u z,i ω u p,i the medium frequency of a Pole-Zero couple and χ u c,i = (χ u z,i + χ u p,i )/2 its logarithm.

Table 2 :

 2 Poles and Zeros list of the test function and K stk1 = 1000

		1	2	3	4	5	6
	χ z	2.52 3.10 3.20 4.10 4.40 4.80
	χ p	2.95 3.50 3.75 3.80 4.60 5.25
	χ p -χ z	0.43 0.40 0.55 -0.30 0.20 0.45

This method gives good results (see Figure17), and the indicators prove it : 0.77% and 0.15%, moreover the causality is respected for all Pole Zero couples (see Table10in Appendix D). Finally, this method is a little bit faster than the simplex method one.

Acknowledgments

The authors wish to thanks Thierry Pasquet and Remi Lemaire from the CSB division, NVH department of Robert Bosch France for their technical and financial help.

A. Rheological models

Model

Generalized Maxwell Model Generalized Fractional Calculus Maxwell with K 0 = 7.746 and γ 1 = 2/3 K 1 = 100 and C 1 = 0.02