Antoine Riaud

Michael Baudoin

Jean-Louis Thomas

Olivier Bou Matar

come

SAW synthesis with IDTs array and the inverse filter: toward a versatile SAW toolbox for microfluidics and biological applications

SAW synthesis with IDTs array and the inverse filter: toward a versatile SAW toolbox for microfluidics and biological applications.

Antoine Riaud, Michael Baudoin, Jean-Louis Thomas, Olivier Bou Matar .

Abstract-Surface acoustic waves (SAWs) are versatile tools to manipulate fluids at small scales for microfluidics and biological applications. A non-exhaustive list of operations that can be performed with SAW includes sessile droplet displacement, atomization, division and merging but also the actuation of fluids embedded in microchannels or the manipulation of suspended particles. However, each of these operations requires a specific design of the wave generation system, the so-called interdigitated transducers (IDTs). Depending on the application, it might indeed be necessary to generate focused or plane, propagating or standing, aligned or shifted waves. Furthermore, the possibilities offered by more complex wave-fields such as acoustical vortices for particle tweezing and liquid twisting cannot be explored with classical IDTs. In this paper, we show that the inverse filter technique coupled with an interdigitated transducers array (IDTA) enables to synthesize all classical wave-fields used in microfluidics and biological applications with a single multifunction platform. It also enables to generate swirling SAWs, whose potential for the on-chip synthesis of tailored acoustical vortices has been demonstrated lately. The possibilities offered by this platform is illustrated by performing successively many operations on sessile droplets with the same system. Index Terms-Acoustofluidics, SAW, Microfluidics, Inverse Filter

I. INTRODUCTION

S INCE the seminal work of Shiokawa in 1989 on the atomization of droplets [START_REF] Shiokawa | Liquid streaming and droplet formation caused by leaky rayleigh waves[END_REF], the potential of SAWs for fluid actuation at microscale has been widely explored in the literature: Progressive waves synthesized by straight IDTs enable sessile droplets displacement [START_REF] Wixforth | Acoustic manipulation of small droplets[END_REF], [START_REF] Renaudin | Saw nanopump for handling droplets in view of biological applications[END_REF], [START_REF] Brunet | Droplet displacement and oscillations induced by ultrasonic surface acoustic waves: a quantitative study[END_REF], [START_REF] Baudoin | Low energy droplet actuation via modulated surface acoustic waves[END_REF] or fluid pumping in a micro-channel [START_REF] Girardo | Polydimethylsiloxane-linbo 3 surface acoustic waves micropumps for fluid control into microchannels[END_REF], [START_REF] Cecchini | Acoustic-counterflow microfluidics by surface acoustic waves[END_REF]. Acoustic field generated by pairs of IDTs can be combined to synthesize stationary waves and manipulate collectively particles [START_REF] Shilton | Particle concentration and mixing in microdrops driven by focused surface acoustic waves[END_REF], [START_REF] Shi | Focusing microparticles in a microfluidic channel with standing surface acoustic waves (ssaw)[END_REF], [START_REF] Shi | Acoustic tweezers: patterning cells and microparticles using standing surfac acoustic waves (ssaw)[END_REF], [START_REF] Raghavan | Particle concentration via acoustically-microcentrifugation[END_REF], [START_REF] Tran | Fast acoustic tweezers fo the two-dimensional manipulation of individual particles in microfluidic channels[END_REF], [START_REF] Ding | On-chip manipulation of single microparticles, cells and organisms using surface acoustic waves[END_REF], [START_REF] Guo | Rf-activated standing surface acoustic wave for on-chip particle manipulation[END_REF], [START_REF] Collins | The particle valve: on demand particle trapping, filtering, and release from microfabricated polydimethylsiloxane membrane using surface acoustic waves[END_REF], or cells [START_REF] Shi | Acoustic tweezers: pattening cells and microparticles using standing surface acoustic waves (ssaw)[END_REF], [START_REF] Franke | Surface acoustic wave actuated cell sorting (sawacs)[END_REF], [START_REF] Hartmann | A novel tool for dynamics cell adhesion studies ) the de-adhesion number investigator dani[END_REF], [START_REF] Sivanantha | Characterization of adhesive properties of red blood cells using surface acoustic wave induced flows for rapid diagnosis[END_REF], [START_REF] Ding | Cell seperation using tilted-angle standing surface acoustic waves[END_REF], [START_REF] Li | Acoustic separation of circulating tumor cells[END_REF], [START_REF] Guo | Controlling cell-cell interactions using surface acoustic waves[END_REF] in droplets or embedded in a microfluidic chamber. Focused waves synthesized by concentric IDTs are suitable for the fusion of droplets [START_REF] Babetta | Droplet translocation by focused surface acoustic waves[END_REF] or for high power applications such as droplet atomization [START_REF] Tan | Interfacial jetting phenomena induced by focused surface vibrations[END_REF]. Finally, more complex operations such as mixing at low Reynolds number through chaotic advection or droplet division may require some fancy combination of plane waves with either rotating [START_REF] Frommelt | Microfluidic mixing via acoustically driven chaotic advection[END_REF] or shifted [START_REF] Collignon | Planar microfluidic drop splitting and merging[END_REF] wavefronts. So far, it has thus been demonstrated that SAWs enable to perform many basic operations required in microfluidics and biological applications [START_REF] Friend | Microscale acoustofluidics: Microfluidics driven via acoustics and ultrasonics[END_REF], [START_REF] Yeo | Surface acoustic wave microfluidics[END_REF].

Nevertheless each of these operations requires a specific optimized design, which is not compatible with a multi-function actuation platform required for the development of many labson-chips. Moreover, classical IDTs do not allow the synthesis of swirling SAWs [START_REF] Riaud | Anisotropic swirling surface acoustic waves synthesis from inverse filtering for on-chip generation of acoustical vortices[END_REF] envisioned for on-chip synthesis of tailored acoustical vortices [START_REF] Riaud | Taming the degeneration of bessel beams at anisotropic-isotropic interface: toward three-dimensional control of confined vortical waves[END_REF] and consequently threedimensional single particle manipulation [START_REF] Baresch | Observation of a singlebeam gradient force acoustical trap for elastic particles: acoustical tweezers[END_REF] and vorticity control [START_REF] Riaud | Cyclones and attractive streaming generated by acoustical vortices[END_REF]. To overcome these shortcomings, we developed an array of 32 optimized IDTs driven by a programmable electronics that enables independent control of each transducer. To synthesize the targeted wave-field, the inverse filter technique [START_REF] Tanter | Time reversal and the inverse filter[END_REF] originally developed for bulk waves has been adapted for SAWs [START_REF] Riaud | Anisotropic swirling surface acoustic waves synthesis from inverse filtering for on-chip generation of acoustical vortices[END_REF]. In this paper, we show that it is possible to synthesize plane waves in different directions, focused waves with focal point located at arbitrary position and swirling SAWs with the same device. The potential of this system is illustrated on sessile droplets by showing successively droplet controlled displacement, division, fusion and nebulization with the same platform.

II. METHODS

Rayleigh surface acoustic waves (R-SAWs) are synthesized at the surface of a X-cut, 1.05 mm thick, Niobate Lithium (LiNb03) piezoelectric substrate by an array of 32 unidirectional interdigitated transducers (SPUDT IDTs) (see Fig. 1). The X-cut was chosen for its good electro-mechanical coefficients along Z (K 2 = 5, 9%) and Y (K 2 = 3, 1%) directions and its weaker anisotropy compared to the 128 • Y-cut generally used for unidirectional IDTs. The IDTs array (IDTA) was fabricated using a lift-off process: (i) the substrate is coated with a AZ15-10 photoresist sacrificial layer patterned with conventional photolithography technique, (ii) a titanium (Ti) layer of 20 nm and a gold (Au) layer of 200 nm are deposited by evaporation (titanium is used for its good adherence on LiNbO3 and gold for its good electrical conductivity), (iii) the sacrificial layer is washed out by a developer.

The IDTA has been designed for a working frequency of 11.9 MHz, with a wavelength adapted to each wave propagation direction of the anisotropic LiNbO3 crystal. Compared to previous work by the authors [START_REF] Riaud | Anisotropic swirling surface acoustic waves synthesis from inverse filtering for on-chip generation of acoustical vortices[END_REF], the design has been further optimized by placing the IDTs on a slowness curve and slightly curving the IDTs to promote diffraction. These two modifications improve the illumination by each transducer of the central zone (radius of 5 mm) of the substrate that we call the "acoustical scene" and where microfluidic operations are performed. An optimal spatial coverage is indeed essential for the synthesis of a wide variety of acoustic wave-fields.

Each of the IDTs is excited independently with a dedicated programmable electronics that enables the synthesis of wave packets at carrying frequencies up to 12 MHz. Impedance matching for each transducers was achieved with external electronic components (inductances). Finally, the inverse filter method [START_REF] Tanter | Time reversal and the inverse filter[END_REF] was used to determine the optimal input signal for each IDT to synthesize a targeted wave-field. Indeed, the inverse filter is a very general method to synthesize a specific wave-field in a linear medium given a set of independent programmable sources. This process can be basically decomposed into three distinct stages. First, the signal emitted by each transducer (impulse response) is measured in the acoustics scene. In practice this response is recorded in a number of control points whose distance cannot exceed λ/2 according to Nyquist-Shannon sampling principle. This allows to define, in the Fourier space, a transfer matrix H ij (ω) (called the propagation operator) between the Fourier transform of the entrance signal emitted by transducer j, E j (ω), and the response signal at control point i, S i (ω): S i = H ij E j (with Einstein notations). In the present experiments, the surface vibrations of the substrate at control point i (typically of the order of a few nanometers) are measured by a homemade polarized Michelson interferometer whose principle is given in reference [START_REF] Riaud | Anisotropic swirling surface acoustic waves synthesis from inverse filtering for on-chip generation of acoustical vortices[END_REF]. Then, the targeted output signal S is defined and the transfer matrix H is inverted to compute the optimal input signal E = H -1 S. Finally, the optimal signal is synthesized by each transducer. If e j (t), s i (t) and h ij (t) are the inverse Fourier transform of E j (ω), S i (ω) and H ij (ω), it is worth noting that the output time signal s(t) is the convolution product of h ij (t) and e j (t):

s i (t) = h ij (t) * e j (t).
While the inverse filter method is simple in principle, some complexity arises when implementing it. Indeed, the propagation operator is generally ill-conditioned since small errors in the measurements produce very large errors in the reconstructed results. Then, the number of control points is not necessarily the same as the number of sources (transducers) and thus the propagation operator is not necessarily a square matrix. So, the pseudo-inverse of the propagation operator is obtained through singular value decomposition. Finally, the inverse filter technique had been initially developed to generate acoustical wave-fields in 3D media. In this case, the target field is a surface and has a smaller dimensionality (2D) than the propagative medium (3D), whereas for SAWs the target field has the same dimensionality as the propagative medium (both 2D). So the control points are not independent and the wave-field must fulfill the dispersion relation. This requires some refinements in the method (see Ref. [START_REF] Riaud | Anisotropic swirling surface acoustic waves synthesis from inverse filtering for on-chip generation of acoustical vortices[END_REF]). To enable droplet manipulation, the central zone was treated with a hydrophobic self assembled monolayer. Otherwise, the droplet would spread on the gold layer at the center of the Niobate Lithium substrate since it is a perfectly wetting medium. Nevertheless, the adherence of alkane-thiol molecules on the central gold layer was not optimal leading to contact angle below 90 • (θ ≈ 70 • ), large hysteresis and thus pining of the contact line. Despite these bad properties of the hydrophobic layer, we nevertheless succeeded to control the drop displacement with the IDTA as we will see in the next section. The droplet dynamics was recorded with a Flea 3 camera on a MZ1 macroscope at 75 frames per second. A sketch of the experimental setup is shown on Fig. 2.

III. RESULTS

A. Complex wave-fields synthesis

To demonstrate the potential of the SAW toolbox for fluid samples actuation, we first investigated the possibilities offered by the IDTA and the inverse filter technique to synthesize the main surface acoustic wave-fields used in the literature: plane, focalized and swirling SAWs. While our setup is optimized for the synthesis of focalized waves and swirling SAWs (since the IDTs are disposed radially around the slowness curve), we have shown that it is possible to synthesize waves with plane wavefronts in the desired direction with peak to peak amplitudes larger than 5.5 nm, that is to say well above the values classically used for droplets displacement [START_REF] Brunet | Droplet displacement and oscillations induced by ultrasonic surface acoustic waves: a quantitative study[END_REF], [START_REF] Baudoin | Low energy droplet actuation via modulated surface acoustic waves[END_REF]. We provide on Fig. 3 two examples of plane progressive waves synthesized in different directions. These results show that despite the curvature of the IDTs, the wavefront aberration is weak in the acoustical scene and compensated by the inverse filter.

Then, we synthesized focalized waves with different positions of the focal point in the acoustical scene and apodization around a preferential direction φ o (see Fig. 4). It is important to note that contrarily to previous attempts to synthesize focalized waves for microfluidic applications with concentric IDTs [START_REF] Friend | Microscale acoustofluidics: Microfluidics driven via acoustics and ultrasonics[END_REF], the anisotropy of the piezoelectric medium is considered here in the definition of the targeted wave-field:

F(r, θ) = 1 2π π -π h(φ -φ 0 , σ)e ikr(φ)r cos(φ-θ) dφ, with h(φ -φ 0 , σ) = exp - φ -φ 0 2 4π 2 σ 2 ,
(r, θ) the polar coordinates, k r (φ) = ω/c R (φ), the radial wave vector, c R (φ), the phase speed of R-SAWs in the φ direction (defining the anisotropy of the medium, σ the aperture (here σ 2 = 0.2) and the symbol denotes the shorter angular distance (which can be computed as follows φ-φ 0 = arccos(cos(φφ 0 ))). The function h then represents an apodization over an aperture σ.

Consideration of the anisotropy is indeed essential to ensure real focalization of the acoustic wave. Excellent results are obtained (see Fig. 4) with maximum normal amplitude of the R-SAW of up to 100 nm. Naturally, small discrepancies are observed compared to the targeted wave-field (the phase does not exactly follows the slowness curve) owing (i) to the finite numbers of transducers used for the field synthesis and (ii) to the finite number of measurement points obtained with the interferometer and used for the reconstruction of the experimental acoustic field. Finally, anisotropic swirling SAWs of topological order 0 and 2 have been synthesized with the new transducers array presented in this paper (see Fig. 5). The targeted wave-field is defined according to the formula introduced in ref. [START_REF] Riaud | Anisotropic swirling surface acoustic waves synthesis from inverse filtering for on-chip generation of acoustical vortices[END_REF]:

W l (r, θ) = 1 2πi l +π -π e ilφ+ikr(φ)r cos(φ-θ) dφ, (1)
with l the topological order of the swirling SAW. It is interesting to note that swirling SAWs of order 0 are nothing but focalized waves with no angular apodization. These waves are not appropriate for particles tweezing or vorticity control since they do not have phase singularity (and thus dark spot) at their center to trap particles and they do not carry angular momentum contrarily to higher order swirling SAWs. Nevertheless, they can be of much practical use for applications where high intensity is required at the focal point.

B. Droplet manipulation with the platform

These specific wave-fields have been used to perform operations on water droplets of initial volume 2 µl. Since it is not possible to synthesize continuous wave-fields with the programmable electronics, some burst of duration 25 µs, carrying frequency 11.9 MHz and variable repetition rates of a few kHz have been used for droplet actuation (for each operation, the exact sequence used for actuation is described on the corresponding figure). Since droplet displacement is the result of cumulative nonlinear effects (radiation pressure and acoustic streaming), and the characteristic hydrodynamic times associated with droplets motion are slow compared to the repetition rates, the forcing is essentially seen as a continuous forcing. Nevertheless, the repetition rates of several kHz are compatible with droplets high order inertio-capillary vibrations [START_REF] Rayleigh | On the capillary phenomena of jets[END_REF], [START_REF] Qi | Interfacial destabilization and aomization driven by surface acoustic waves[END_REF], [START_REF] Baudoin | Low energy droplet actuation via modulated surface acoustic waves[END_REF], [START_REF] Blamey | Microscale capillary wave turbulence excited by high frequency vibration[END_REF], which may explain why spurious atomization of the drop has been observed. Indeed, both the frequency of the excitation and the acoustical powers used in the present experiments are compatible with droplet nebulization [START_REF] Qi | Interfacial destabilization and aomization driven by surface acoustic waves[END_REF] (larger power than usual were used in the present experiments to overcome the retention of the contact line). Nevertheless, this shortcoming could be simply overcome with appropriate hydrophobic treatment of the surface with low hysteresis [START_REF] Crudden | Ultra stable self-assembled monolayers of n-heterocyclic carbenes on gold[END_REF]. Droplet displacement (Fig. 6, Movie S1), fusion (Fig. 7, movie S3) and atomization (Fig. 8, second part of movie S3) were obtained respectively with focalized waves, swirling SAWs of second and zero order. The characteristic of the wave-field are summarized on the corresponding figure. For the droplet division (Fig. 9, movie S2), we used a different method than the one proposed in ref. [START_REF] Collignon | Planar microfluidic drop splitting and merging[END_REF]: we alternatively synthesized some burst of focalized waves with two different focal points as presented in the previous subsection. With this specific wave-field, we were able to separate droplets even at these high contact line hysteresis.

IV. DISCUSSION

In this paper, we have shown that the combination of interdigitated transducers array (IDTA) and the inverse filter technique enables the synthesis of the most common surface acoustic wave-fields used for actuation of fluids at the microscale: plane waves in different directions, anisotropic focalized SAWs and anistropic swirling SAWs. As a proof of concept, we show that this SAW toolbox enables to perform many basic operations required in droplet based digital microfluidics: droplet displacement, division, fusion and atomization. Since it is virtually possible to synthesize any acoustic wave-field compatible with the anisotropy of the substrate, it might be easily demonstrated that it is possible to perform other operations previously demonstrated in the literature such as micro-mixing, collective particle manipulation with standing waves, or jetting. In the same way, it might be possible to demonstrate the versatility of this platform for fluid samples manipulations in microchannels.

The most thrilling perspectives of this SAW toolbox nevertheless lie in the development of new operations that cannot be achieved with other techniques, such as selective 3D manipulation of one or a few particles or cells (displacement, rotation), precise vorticity control independently of the boundary conditions with swirling SAWs, or collective manipulation of several sessile droplets. This last operation is indeed not possible with a limited number of transducers since (i) there is a shadow zone behind each droplet that prevents the displacement of a second droplet situated behind it due to the strong absorption of Rayleigh waves by liquid samples, and (ii) it would be necessary to localize the droplets position in real time. This problems can be overcome with IDTA. Indeed, the combination of the acoustic waves emitted by each transducers would enable better spatial coverage of the substrate and thus the elimination of shadows zone. Moreover the IDTs can be used not only as actuators but also as sensors of the acoustical echo as demonstrated by Alzuaga et al. [START_REF] Alzuaga | A large scale x-y positioning and localisation system of liquid droplet using saw on linbo3[END_REF], [START_REF] Bennes | Detection and high-precision positioning of liquid droplets using saw systems[END_REF] and Renaudin et al. [START_REF] Renaudin | Saw nanopump for handling droplets in view of biological applications[END_REF] to determine the location of the drop without extra visualization setup.

This work opens prospects toward the design of potentially the most versatile toolbox for active control of small fluid samples, for microfluidics and biological applications.

  Manuscript received December 25, 2015; revised XXX; accepted XXX. This work was supported by Agence Nationale de la Recherche ANR-12-BS09-0021-01 and ANR-12-BS09-0021-02 and Region Nord Pas de Calais. The authors are with Institut d'Electronique, de Microélectronique et de Nanotechnologie (IEMN), LIA LICS, Université de Lille and EC Lille, UMR CNRS 8520, 59652 Villeneuve d'Ascq, France and Sorbonne Universités, UPMC Univ Paris 06, CNRS UMR 7588, Institut des NanoSciences de Paris (INSP), F-75005 Paris, France. Corresponding author : Michael Baudoin (email : michael.baudoin@univ-lille1.fr)

Fig. 1 .

 1 Fig.1. Optimized array of 32 interdigitated unidirectional transducers used for the synthesis of various surface acoustic wave-fields (focused, plane, swirling) and the manipulation (displacement, division, fusion, atomization) of droplets. The central black zone is a gold layer used as a mirror for the measurement of the transducer response with the Michelson interferometer.

Fig. 2 .

 2 Fig. 2. Sketch of the experimental setup used for the synthesis of various wave-fields and the actuation of droplets.

Fig. 3 .

 3 Fig. 3. Phase of plane progressive SAWs synthesized with the IDTA in two different directions (a and b). The maximum peak to peak amplitude in the two directions are respectively 6 and 10 nm.

Fig. 4 .

 4 Fig. 4. Phase and amplitude of focalized propagative anisotropic SAWs converging at two different focal points (a and b). The maximum peak to peak amplitude for the two different focal points are respectively 60 nm and 100 nm. Colorbar: peak to peak amplitude ×10 -2 nm.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Phase and amplitude of anisotropic swirling SAWs of topological order 0 and 2 (a and b). The maximum peak to peak amplitude of the W 0 and W 2 swirling SAWs are respectively 108 nm and 80 nm. Colorbar: peak to peak amplitude ×10 -2 nm.

Fig. 7 .

 7 Fig. 7. Stack of images showing the fusion of two droplets of 2 µl. The fusion is obtained with a centered swirling SAWs of topological order 2. The time elapsed between two successive images is 27 ms. See also movie S3.

Fig. 8 .

 8 Fig. 8. Stack of images showing the atomization (nebulization) of a 2 µl droplet. It is obtained with a high intensity centered swirling SAW of order 0 (annular wave). The time elapsed between two successive images is 13 ms. See also second part of movie S3.

Fig. 9 .

 9 Fig. 9. Stack of images showing the asymmetric division of a 2 µl droplet into two daughter droplets of different volumes. It is obtained with waves focalized at two focal points situated on both sides of the drops and with angular apodization. The time elapsed between two successive images is 27 ms. See also movie S2.

ACKNOWLEDGMENT

The authors would like to thank Pr. James R. Friend and Pr. David. W. Greve. for their invitation to IEEE International Ultrasonic Symposium 2015 in Taipei, that enabled our group to present these results.


Rafael Penas

Etienne Balmes
email: balmes@sdtools.com

Arnaud Gaudin

S Laflamme

A unified non-linear system model view of hyperelasticity, viscoelasticity and hysteresis exhibited by rubber

Keywords:

Rubber modeling is an old subject and so many models exist that it is difficult to have a clear vision of what exists and is more appropriate. Rather than attempting a standard review, this paper proposes classification using the traditional system modeling strategy, where raw measurements are either processed to obtain non-parametric models, or used to identify parametric models, whose accuracy can be controlled by order selection or by numerical implementation considerations. A full test campaign, including multi-step relaxation, low speed triangular and sine tests, on a large deformation compression sample is used to illustrate the need to model and combine the base behaviors known as hyperelasticity, viscoelasticity, and rate independent hysteresis. The equivalence between linear viscoelasticity and linear time invariant systems is used to clarify the link between order selection and accuracy of a generalized Maxwell model. Rate independent hysteresis is analyzed using a convolution product like the one used for viscoelastic transients by introducing a relaxation modulus. Measurements of the hysteretic relaxation modulus are used to propose strategies to measure the asymptotic hyperelastic modulus and discriminate between different hysteretic model forms. A parallel between Iwan and Maxwell models is detailed, and non-parametric models are used to show that the two overlap in the low frequency small deformation regime. Regularized rate independent hysteresis and non-linear viscoelasticity are finally shown to lead to a similar view allowing a transition between the rate independent and linear relaxation models. The instantaneous ratio analytic force and displacement signals, or instant complex modulus, is introduced as novel non-parametric estimation of sine measurements and shown to be a powerful tool to analyze and validate the fact that a force rate relaxation with non-linear relaxation frequencies is most appropriate to represent the non-linear coupling of all three effects.

Introduction

System models provide a relation between inputs and outputs. Identification is the usual name describing an inverse problem where a class of models is chosen, and the associated values are selected to minimize a distance between test results and the model predictions. It is classical [START_REF] Ljung | System Identification: Theory for the user[END_REF] to distinguish non-parametric identification where the model is described by a curve whose values can be derived from test and parametric identification, where the curve is described in functional form using a finite number of parameters. Parametric models are typically decomposed in a number of evolution equations describing how system model states change for given inputs, and observation equations relating the outputs with states. The number of states involved in the evolution equation corresponds to a model order.

Parametric forms have the advantage of being high resolution, meaning that they can be used to build non-parametric models with any degree of accuracy (frequency or time resolution). The inverse transformation from non-parametric to parametric models is called identification and is a numerical process that requires a choice of the parametric form and an optimization process to estimate the parameters. Its results notably depend on a proper choice of representative non-parametric data, motivating the need for notable work in selecting appropriate tests and data extraction. For the considered application of rubber traction tests, detailed in Section 2, instant modulus will be shown to be relevant for both low speed triangular and sine testing. The use of instantaneous ratio of analytic signals to characterize the non-linear behavior during a period, thus considering the system to be linear time dependent system rather than non-linear is a notable contribution to the analysis of the considered test.

It is also useful to distinguish two forms of parametric models: selected order and order independent. The first ones are simpler to implement as they are represented by an assembly of basic components and with order selection driven by accuracy requirements. The order independent models are normally based on physical assumptions on the behavior. They thus demand less data and less parameters for identification. Furthermore, having a predefined type of behavior, they tend to allow data extrapolation. Fractional derivatives [START_REF] Koeller | Application of fractional calculus to the theory of viscoelasticity[END_REF][START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF][START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF][START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF] and triboelastic [START_REF] Coveney | A triboelastic model for the cyclic mechanical behavior of filled vulcanizates[END_REF][START_REF] Segalman | Using series-series iwan-type models for understanding joint dynamics[END_REF] models are examples of order independent viscoelastic and hysteretic models that will be detailed. These models are however generally not suited for transient integration so that a transformation to a selected order form is necessary. Such transformations will be analyzed in terms of discretization and relation between order selection and accuracy will be illustrated. Despite the advantages, if such model is not consistent with data, bias is induced and thus problematic.

Material models characterize the relation between stress outputs and the history of strain inputs. For scalar or 0D bushing models, one can also consider the system model to give a force/state map [START_REF] Al-Hadid | Application of the force-state mapping approach to the identification of non-linear systems[END_REF][START_REF] Witters | Black-box model identification for a continuously variable, electro-hydraulic semi-active damper[END_REF]. At the material level, both scalar uniaxial and tensorial multivariate [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF] strain/stress inputs/outputs have been considered, even though the focus will here be placed on the scalar case.

While materials for flexible mounts or bushings are used for the present illustration, the same models are considered in other fields: biologic tissues [START_REF] Franchini | Viscoelastic characterization of human descending thoracic aortas under cyclic load[END_REF], composite materials [START_REF] Cappelli | Multi-scale identification of the elastic properties variability for composite materials through a hybrid optimisation strategy[END_REF], and even meta-materials [START_REF] Dykstra | Viscoelastic snapping metamaterials[END_REF]. These materials are known to exhibit different types of behavior, and in particular

• hyperelasticity, which assumes stress to depend non-linearly on the current value of strain;

• viscoelasticity, which describes a dynamic dependence of stress on strain history; • rate-independent hysteresis, also called friction in a number of settings, characterizes the dependence of stress on the strain path.

Each effect has been extensively studied independently, combined using both series and parallel combinations of constitutive laws, while including non-linearity representations in a number of ways. The contribution of this paper is thus to present a unified perspective on how to select non-parametric test data for identification, parameter forms for identification and model order selection strategies for a good accuracy.

For the traction test that will be detailed in Section 2, Fig. 1 illustrates how relaxation, constant speed and sine tests, can be combined in a single non-parametric map giving hyperelastic, hysteretic and viscoelastic complex moduli (ratio between stress and strain variations) as a function of static operating strain. In this map, the solid red curve illustrates the hyperelastic behavior. The dashed green ones correspond to low speed tests and illustrate the rate-independent hysteresis (with largest difference from hyperelastic curve visible right after turning points), detailed in Section 5.2. Finally, the blue to yellow maps give the average complex modulus dependence on frequency (blue to yellow), pre-strain (black vertical lines) and amplitude (distance to black vertical lines). This non-parametric representation readily illustrates that hyperelasticity is a base branch with an increased stiffness, or parallel branch, for hysteretic behavior and further increase for viscoelasticity. After separate presentations of each effect, Section 6 will detail how this data can be used to validate combined non-linear material models. Dynamic dependence on strain history is first assumed to be linear in Section 4. This is called linear viscoelasticity in the materials community and Linear Time Invariant (LTI) systems in the control community. The direct output of tests are non parametric representations that characterize LTI systems either through frequency dependent transfer functions, called complex modulus in the viscoelastic material literature [START_REF] Nashif | Vibration Damping[END_REF][START_REF] Lakes | Viscoelastic Materials[END_REF] 𝜎(𝜔) = 𝐸(𝜔)𝜀(𝜔) or 𝐹 (𝜔) = 𝐾(𝜔)𝑥(𝜔) [START_REF] Ljung | System Identification: Theory for the user[END_REF] or their inverse Fourier transform, known as a time domain impulse response, which corresponds to the time derivative of the relaxation function [START_REF] Salençon | Viscoélasticité, Presse des Ponts et Chaussés[END_REF], as will be further detailed in Section 4. These transfers verify conditions on stability and causality enforcing constraints on the amplitude and phase of the transfer [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF].

Parametric LTI models are also often used in system dynamics with multiple forms of interest: rational fractions, pole-zero representations, state space models for which bidirectional transforms exist. In viscoelasticity, classical parametric forms are the generalized Maxwell, which is a subset of rational fractions (see [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF] for example), and fractional derivative models (see [START_REF] Koeller | Application of fractional calculus to the theory of viscoelasticity[END_REF]). While both are often considered as separate models, Section 4.3 will sustain the argument that finite horizon approximations of the fractional derivative model (such as the Grünwald approximation [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF][START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF]) have an explicit expression as rational fractions. The fractional derivative model can thus be seen as an order independent parametrization of a subclass of Maxwell models with an identification process including order selection controlling the distance between the two forms.

The literature on rate-independent behavior is very rich from the mathematical point of view [START_REF] Brokate | Hysteresis and Phase Transitions[END_REF]. Base rules ensuring physically realistic behavior are Madelung rules on hysteresis loop closure and Masing's law [START_REF] Masing | Eigenspannungen und verfestigung beim messing[END_REF][START_REF] Festjens | Identification and Modeling of Jointed Structures for Dynamic Analysis[END_REF] on initial and subsequent loops. Following a path similar to viscoelasticity, Section 5 introduces a non-parametric formulation of hysteresis as a convolution problem associated with an hysteretic relaxation modulus. This formulation induces automatic verification of Madelung rules. Selected order parametric models are parallel combination of classical Jenkins cells combining a spring and a friction element. Combined in an Iwan series [START_REF] Segalman | Using series-series iwan-type models for understanding joint dynamics[END_REF] this is a hysteretic equivalent of the Maxwell model. Finally, as for fractional derivatives, some authors [START_REF] Coveney | A triboelastic model for the cyclic mechanical behavior of filled vulcanizates[END_REF][START_REF] Segalman | Using series-series iwan-type models for understanding joint dynamics[END_REF] have proposed order independent parametric forms that will be discussed in Section 5.3.

Fig. 1 illustrates that all three behavior hyperelasticity, hysteresis and viscoelasticity are coupled and associated strategies are described in Section 6. Assuming purely hysteretic models is not realistic and a number of authors have proposed regularization techniques [START_REF] Festjens | Identification and Modeling of Jointed Structures for Dynamic Analysis[END_REF][START_REF] Ismail | The hysteresis bouc-wen model, a survey[END_REF]. Assuming linear viscoelasticity is not realistic either, and Heymans [START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF] sustains that a structural temperature depending on physical temperature, history of anelastic deformation, and other environmental effects should be used. This idea is shown to correspond to the use of non-linear shifts of viscoelastic relaxation times, which was proposed in [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF] under the name reduced time. The unified result is thus a tribo-viscoelastic model with an underlying base hyperelastic behavior and parallel stress, rather than strain, relaxation contributions. For low frequency cells, relaxation frequencies notably increase for significant velocities, which can be seen as an increased structural temperature, leads to an effect that is difficult to distinguish from rate-independent hysteresis, and effectively models the so-called Payne effect.

Description of experimental test case

The traction/compression sample show in Fig. 2, developed and tested by Vibracoustic NE, is used throughout this paper. The test article shape ensures a constant deformation on the central section of the specimen. The relation between displacement of the machine and strain in the central section is obtained using a polynomial spline generated by image correlation. The first component of the nominal (1st Piola-Kirchhoff) stress is obtained by the ratio between the measured force and the reference surface

𝑃 11 𝑡𝑒𝑠𝑡 =  𝑡𝑒𝑠𝑡 ∕𝐴.
A number of tests were performed using the profiles illustrated in Fig. 3. The preexisting process used multi-step for hyperelastic characterization, triangular excitation to characterize rate-independent hysteresis and sine testing for viscoelastic complex modulus testing. The present paper is the result of rethinking how to best exploit the data of these tests and propositions for revising the procedure will be given in the conclusion. The tests were performed on a MTS Landmark 200 Hz Elastomer Test System machine. The acquisition bandwidths are 1 Hz for the multi-step relaxation tests, 3 Hz for triangles at 10 mm/min and 50 times the sine frequency for sine testing.

Hyperelasticity or non-linear static gain

Rubber bushings often undergo very large deformations. Linear elasticity (or constant static gain in LTI systems terminology) may be not suited. The origins of this non-linearity are both material and geometric, but here the distinction will not be made, as the specimen is supposed to capture only material effects. The non-parametric hyperelastic estimation corresponds to the isolation of an elastic backbone curve.

Hyperelasticity can be characterized by multi-step relaxation tests, illustrated in Fig. 4. At the end of each step, one waits a long time (60 min, in this case) to obtain a reasonable estimate of the long term/low frequency elastic behavior. One can thus extract an asymptotic behavior shown as red square markers in the figure which corresponds to the end point of relaxation experiments shown as black dotted lines. One can however note that the points on the upward and downward part of the experiment do not perfectly coincide and when points are repeated there is some spread. There is thus, either a long term or a path effect, which can be represented by either hysteresis or extremely slow relaxation. Note that the relaxation transient can also be analyzed using a viscoelastic assumption leading to consistent results.

The figure also illustrates the possibility to extract the constant part of sine tests (harmonic 0) shown as blue circles slowly converging towards the hyperelastic behavior, and finally the much more detailed hyperelastic behavior shown as a solid red curve that was extracted from triangular or quasi-static test as will be detailed next.

To obtain detailed non-parametric hyperelastic model, the novel proposition of this work is to replace the traditional stress/ deformation curve shown in Fig. 5 left, by the instantaneous slope or tangent modulus illustrated in Fig. 5 right. In this representation it appears that some distance away from the turning point, all curves converge towards a single hyperelastic backbone shown as solid lines while close to the turning point as materialized by dashed lines the modulus is higher. Section 5.1 will introduce the idea of hysteretic relaxation.

The oscillations present in the solid hyperelastic modulus are here due to signal processing difficulties. The test used was not designed to extract the hyperelastic modulus. Thus, a rather low 0.3 Hz bandwidth was used and an order four cascaded second order section filter with a padding fraction of 80% was found to yield reasonable results. But in future tests, acquisition should be redesigned with the objective of accurately measuring the instantaneous tangent modulus.

To recover the single stress/strain hyperelastic curve shown in Fig. 4, the hyperelastic moduli found in different parts of the test are averaged and a piece-wise cubic interpolation is built. For harmonic sines tests, the linear assumption is that the time responses contain harmonics 0 (static operating condition) and 1 of excitation frequency 𝜔, leading to a signal of the form

𝑢(𝑡) = 𝑢 0 + ℜ ( 𝑢 1 (𝜔) exp(𝑖𝜔𝑡) ) = 𝑢 0 + ℜ(𝑢 1 (𝜔)) cos(𝜔𝑡) -ℑ(𝑢 1 (𝜔)) sin(𝜔𝑡) (2)
If the system is linear time invariant, one can classically obtain the harmonic 1 coefficients 𝑢 1 (𝜔) through the mean over one period

ℜ(𝑢 1 (𝜔)) = ∫ 2𝜋∕𝜔 0
𝑢(𝑡) cos(𝜔𝑡), . . . Here however, hyperelastic contribution is non-linear, and the classic approach would be to consider harmonic balance, where the signal would be composed of harmonics

𝑢(𝑡) = ∑ 𝑘 ℜ ( 𝑢 𝑘 (𝜔) exp(𝑖𝑘𝜔𝑡) ) .
Here, rather than focusing on the first harmonic of stress to obtain complex modulus estimate, a new approach will seek to define an instantaneous complex modulus as follows. First, the signals are assumed to be analytic, as done in the Hilbert transform [START_REF] Feldman | Hilbert Transform Applications in Mechanical Vibration[END_REF]. One thus defines a zero mean signal

ū(𝑡) = 𝑢(𝑡) - 𝜔 2𝜋 ∫ 2𝜋∕𝜔 0 𝑢(𝑡)𝑑𝑡 (3)
and builds a sliding window of 𝑁 points computing the instantaneous complex amplitude of a signal by solving the linear least squares problem

min 𝑢 1 (𝑡 1 ) ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ⎡ ⎢ ⎢ ⎣ cos(𝑡 1 ) -sin(𝑡 1 ) ⋮ ⋮ cos(𝑡 𝑁 ) -sin(𝑡 𝑁 ) ⎤ ⎥ ⎥ ⎦ { ℜ(𝑢 1 (𝜔)) ℑ(𝑢 1 (𝜔)) } - ⎧ ⎪ ⎨ ⎪ ⎩ 𝑢(𝑡 1 ) ⋮ 𝑢(𝑡 𝑁 ) ⎫ ⎪ ⎬ ⎪ ⎭ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ (4)
The instantaneous complex modulus is then obtained as the ratio 𝐸(𝜔, 𝑡 1 ) = 𝜎(𝜔, 𝑡 1 )∕𝜀(𝜔, 𝑡 1 ) of instantaneous stress and strain amplitudes. This can be seen as analyzing the problem as linear time varying rather than non-linear. Fig. 6, using direct sine measurements stabilized for 20 periods, clearly illustrates the cycles shown left cannot be considered as resulting from an LTI system where the stress/strain trajectories would be ellipses, and the instantaneous stiffness would be constant. The instantaneous complex stiffness illustrates a few useful signal features: the overall slope is decreasing with the hyperelastic modulus (shown by the black dotted line), while on the right the ratio with the hyperelastic modulus is more constant, as [START_REF] Franchini | Viscoelastic characterization of human descending thoracic aortas under cyclic load[END_REF] points out for the case of dynamic stiffness. The upward ⊳ and downward ⊲ parts of the cycle are notably different which will be related to non-linear effects in Section 6.1. The modulus increase with frequency, not shown in the figure, will be related to viscoelasticity in Section 4.

Finally, the instantaneous modulus is clearly lower at higher cycle amplitudes (the solid blue curves are above the dashed red ones). This is classically known as Payne effect [START_REF] Mullins | Softening of rubber by deformation[END_REF][START_REF] Fletcher | Nonlinearity in the dynamic properties of vulcanized rubber compounds[END_REF] and is actually much more visible as an overall decrease of the instant modulus for all values of pre-strain. This effect is quite notable for carbon black filled rubbers and has been extensively studied [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF][START_REF] Bhuiyan | A rheology model of high damping rubber bearings for seismic analysis: Identification of nonlinear viscosity[END_REF][START_REF] Cantournet | Mullins effect and cyclic stress softening of filled elastomers by internal sliding and friction thermodynamics model[END_REF]. Section 6.2 will use the analysis of how models reproduce these features to discriminate their validity.

The window length has a significant influence on the detailed shape. In the presence of noise, short windows have strong fluctuation. Here, the impact is mostly visible for low amplitude cycles. In the future, using Kalman estimation of the analytic signals [START_REF] Dion | Tracking and removing modulated sinusoidal components: A solution based on the kurtosis and the extended Kalman filter[END_REF] would seem the best approach. The idea of filtering measured signals using harmonics of the excitation was tested but tends to induce notable distortion. For test/analysis correlation that will be detailed in Section 6.2, using the same fraction of the period was found to be important.

Linear viscoelasticity

As stated in the introduction, linear viscoelasticity is the transcription to material behavior of the system dynamics assumption of a linear time invariant system where stress and strain are linked. LTI systems have many equivalent representations. In particular, one distinguishes

• non-parametric forms: the transfer function which is an analytic function of frequency, its inverse transform which is a time domain impulse response and this can be exploited using convolution products as detailed in Section 4.1 • rational selected order forms: continuous time differential equations, their representation in the form of rational fractions characterized by their poles and zeros, the partial fraction expansion as a sum of first order rational fractions with associated gains, or finally the discrete time recursion equations used for numerical integration are detailed in Section 4.2 • fractional derivatives often considered in the viscoelasticity literature are analyzed as order independent forms in Section 4.3.

The theory of viscoelasticity is applied to macroscopic 0D models, where the variables of interest are displacement 𝑥 and force 𝐹 , but as the test was for material characterization, scaling discussed in Section 2 is used to show results as strain 𝜀 and stress 𝜎.

Non-parametric transfer: complex and relaxation modulus

Three equivalent representations are considered: a transfer function in the frequency domain [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF], an impulse response associated with convolution (6) or a differential equation in the time domain [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF]. In material analysis, the transfer is called complex modulus and one has

𝜎(𝜔) = 𝐸(𝜔)𝜀(𝜔) = 𝐸 ′ (𝜔)(1 + 𝑖𝜂(𝜔))𝜀(𝜔) = |𝐸(𝜔)| exp(𝑖𝛿)𝜀(𝜔).
(

) 5
where the complex modulus may also be separated in its real part 𝐸 ′ and its complex part 𝐸 ′′ , respectively named storage and loss modulus, and the loss factor 𝜂 = 𝐸 ′′ ∕𝐸 ′ = tan(𝛿) is used here preferably to the transfer phase 𝛿. Estimation of linear transfers is a classical process [START_REF] Heylen | Modal Analysis Theory and Testing[END_REF], which can be performed using broadband signals (impulse response) or stepped sine testing which is the classical approach in Dynamic Mechanical Analysis (DMA) tests. Estimation of transfers at a set of frequencies is called non-parametric in the sense that it only gives the estimate of the complex modulus at a selected set of frequencies.

The inverse Fourier transform of the complex modulus is the relaxation modulus

𝐸(𝑡) = 𝑅 ′ (𝑡) = ℱ -1 {𝐸(𝜔)} = 1 2𝜋𝑖 ∫ +∞ -∞ 𝐸(𝜔) exp(𝑖𝜔𝑡)d𝜔 (6)
which can be measured directly through relaxation testing. The relations between stress and strain can be obtained through convolution of either strain with relaxation modulus or strain rate with relaxation function

𝜎(𝑡) = ∫ +∞ 0 𝐸(𝑟)𝜀(𝑡 -𝑟)d𝑟 = ∫ +∞ 0 𝑅(𝑟) ε(𝑡 -𝑟)d𝑟 (7)
To ease handling of initial conditions, classical texts on viscoelasticity [START_REF] Salençon | Viscoélasticité, Presse des Ponts et Chaussés[END_REF] use the Laplace-Carson transform and the relaxation function 𝑅(𝑡), but the Laplace transform can be used for 𝐸(𝑡).

The direct use of such curves is not very easy in transient simulation since it requires alternating between time and frequency domain. This is a very strong reason to prefer parametric models. Any system identification textbook [START_REF] Ljung | System Identification: Theory for the user[END_REF] will describe parametric models in terms of differential equations in time or rational fraction representations of transfers. But since the viscoelasticity term comes from material modeling, the process is called rheologic modeling as will be detailed next.

Selected order base representation: rational complex modulus

Rational fraction representations of transfers are, in the frequency domain, classically characterized by poles 𝑝 𝑖 and zeros 𝑧 𝑖 or separated in a sum of first order rational fractions with constant numerators called residues and first order denominators associated with poles.

𝐾(𝑠) = 𝐹 (𝑠) 𝑥(𝑠) = 𝐾 0 𝑁 ∏ 𝑖=1 1 -𝑠 ∕𝑧 𝑖 1 -𝑠 ∕𝑝 𝑖 = 𝐾 ∞ ( 1 + 𝑁 ∑ 𝑖=1 -𝑔 𝑖 𝑠 ∕𝜔 𝑖 + 1 ) = 𝐾 ∞ ( 𝑔 0 + 𝑁 ∑ 𝑖=1 𝑔 𝑖 𝑠 𝑠 + 𝜔 𝑖 ) . ( 8
)
For stability, all poles must have negative real parts. As most of the literature, real poles and positive gains will be considered here, but second order rational fractions are discussed in [START_REF] Mctavish | Modeling of linear viscoelastic space structures[END_REF] allowing the use of complex conjugate poles.

Real poles are associated with a characteristic time 𝜏 𝑖 = 2𝜋∕𝜔 𝑖 . The low frequency asymptotic gain is noted 𝐾 0 , the high frequency asymptote is noted 𝐾 ∞ and the residue associated with each pole is expressed as a fraction of the high frequency gain 𝐾 𝑖 = 𝐾 ∞ 𝑔 𝑖 . Note that this is the SDT [START_REF]Structural Dynamics Toolbox (for Use with MATLAB[END_REF] convention, where load fractions are such that ∑ 𝑁 𝑐𝑒𝑙𝑙 𝑖=0 𝑔 𝑖 = 1, corresponding to the use the high frequency modulus as reference. But viscoelastic materials in Abaqus [START_REF]ABAQUS/Standard, User's Manual[END_REF] have a different convention where ĝ𝑖 = 𝑔 𝑖 ∕𝑔 0 , using the low frequency modulus as reference.

These two equations imply that there are two different ways to see the problem: the middle form considers a high frequency total force 𝐾 ∞ 𝑥(𝑠) and removes relaxation of the displacement, and the right form considers a base elastic contribution 𝐾 0 𝑥(𝑠) and adds relaxations forces proportional to velocity.

In mechanics, this standard LTI model is known as the generalized Maxwell model illustrated in Fig. 7, also known as Prony series in the time domain, where the total load is the sum of a series of forces

𝐹 = 𝑁 ∑ 𝑖=0 𝐹 𝑖 , (9)
each associated with a first order relaxation equation which contributes with a gain increase and dissipation maximum centered at the pole frequency. Each branch is described by a time domain differential equation, equivalent to the LTI frequency domain formulation [START_REF] Al-Hadid | Application of the force-state mapping approach to the identification of non-linear systems[END_REF], which corresponds to the equilibrium of the internal cell point 𝐹 𝑖 = 𝑐 𝑖 ẋ𝑖 = 𝐾 𝑖 (𝑥 -𝑥 𝑖 ). The differential equation can be stated using displacements

𝑐 𝑖 𝐾 𝑖 ẋ𝑖 = ẋ𝑖 𝜔 𝑖 = (𝑥 -𝑥 𝑖 ) (10)
but it can also be written in terms of relaxation loads to be removed from the high frequency asymptotic force, as in the middle form of the equation ( 8),

Ḟ 𝑖 𝜔 𝑖 + 𝐹 𝑖 = -𝑔 𝑖 𝐹 ∞ (𝑥), (11)
or written as loads added to low frequency asymptotic force, as in the last form of (8), leading to a load rate relaxation

Ḟ 𝑖 + 𝜔 𝑖 𝐹 𝑖 = 𝑔 𝑖 Ḟ ∞ (𝑥). (12)
While all these forms are equivalent in the linear case, Section 6 will show notable differences when coupled with non-linearity, with the force rate relaxation, given by equation ( 12), being clearly more appropriate.

Having selected rational models, the next step is to choose the order, which refers to the equivalent notions of number of Maxwell cells, number of poles in the model and number of internal states. The relation between 𝐾 𝑖 , 𝑐 𝑖 , poles and zeros is for example given in the appendix of [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF].

The term selected order is used here to emphasize the fact that identification of a Maxwell model should be viewed as a discretization problem and thus related to an accuracy objective. The modulus being an increasing function of frequency, the poles 𝑝 𝑖 = -2𝜋∕𝜏 𝑖 are real and negative, and the rational fraction gains 𝐾 ∞ 𝑔 𝑖 are positive. Renaud & al. [17] give a clear overview both graphical and numerical methods used for complex modulus approximation methods. Of particular interest, is the discussion that poles can be placed arbitrarily, and accuracy controlled by the number of poles per decade.

To illustrate the discretization argument, it is chosen to discretize the storage modulus interval and place poles accordingly. In Fig. 8, the order 4 model, illustrates the fact that low order models are less accurate on storage modulus and loss factors. Even though the 30% error on loss may appear poor for the order 4 model, the need for more accuracy is debatable considering all other unknowns. In the end, our design rule is to use more than one pole per decade.

For numerical implementations, the differential equations must be discretized, in other words go from continuous Laplace forms to the discrete z-transform, with all the associated trade-off on accuracy and stability. Using integration over a time step, [START_REF] Simo | Computational Inelasticity[END_REF] obtains a recursion equation for the load relaxation which present a z-domain pole exp(-𝜔 𝑖 𝛥𝑡) which is always stable, but may have numerical conditioning problems if 𝜔 𝑖 𝛥𝑡 is too close to zero.

𝐹 𝑖 𝑛 = exp ( -𝜔 𝑖 𝛥𝑡 ) 𝐹 𝑖 𝑛-1 -𝜔 𝑖 𝑔 𝑖 𝛥𝑡 exp ( - 𝜔 𝑖 𝛥𝑡 2 ) 𝐹 ∞ 𝑛 + 𝐹 ∞ 𝑛-1 2 (13)
If one prefers to integrate the force rate, the recursion is given by

𝐹 𝑖 𝑛 = exp ( -𝜔 𝑖 𝛥𝑡 ) 𝐹 𝑖 𝑛-1 + 𝑔 𝑖 exp ( - 𝜔 𝑖 𝛥𝑡 2 ) (𝐹 ∞ 𝑛 -𝐹 ∞ 𝑛-1 ). ( 14
)

Order independent parametric models : fractional derivatives

Fractional derivatives [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF][START_REF] Lakes | Viscoelastic Materials[END_REF] are known to be a good alternative to Maxwell or rational models for viscoelastic materials. The main advantage is the use of only 4 parameters, or slightly more [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF]. The base time domain expression is

𝐹 + 1 𝜔 𝛼 𝑐 d 𝛼 𝐹 d𝑡 𝛼 = 𝐾 0 𝑥 + 1 𝜔 𝛼 𝑐 𝐾 ∞ d 𝛼 𝑥 d𝑡 𝛼 . ( 15
)
with the low frequency asymptotic stiffness 𝐾 0 , the high frequency stiffness 𝐾 ∞ , the peak damping frequency 𝜔 𝑐 and the fractional derivative value 𝛼. The Laplace transform of this equation is

( 1 + ( 𝑠 𝜔 𝑐 ) 𝛼 ) 𝐹 (𝑠) = ( 𝐾 0 + ( 𝑠 𝜔 𝑐 ) 𝛼 𝐾 ∞ ) 𝑥(𝑠), ( 16
)
leading to simple expression for the complex stiffness,

𝐾(𝑠) = 𝐾 ∞ + 𝐾 0 -𝐾 ∞ 1 + ( 𝑠 ∕𝜔 𝑐 ) 𝛼 , ( 17
)
where maximum loss is now also controlled by 𝛼 and not only dependent on the modulus ratio 𝐾 ∞ ∕𝐾 0 . The advantage of using few parameters is balanced by the need for convolutions. Indeed, the Riemann-Liouville definition of fractional derivatives is expressed as

d 𝛼 𝑓 (𝑡) d𝑡 𝛼 = 1 𝛤 (1 -𝛼) d d𝑡 ∫ 𝑡 0 𝑓 (𝑟) 𝑡 -𝑟 d𝑟 (18)
Since convolutions are not acceptable in practice, approximations are used. To illustrate that this corresponds to selecting a fixed order rational model, one will consider Grünwald's approximation [3]

d 𝛼 𝑓 (𝑡) d𝑡 𝛼 ≈ ( 𝑡 𝑁 ) -𝛼 𝑁 ∑ 𝑗=0 𝐴 𝑗+1 𝑓 (𝑡 -𝑗𝛥𝑡) , with 𝐴 𝑗 = 𝛤 (𝑗 -𝛼 -1) 𝛤 (-𝛼)𝛤 (𝑗) , ( 19
)
where 𝛤 is the gamma (continuum factorial) function. Considering fixed values for 𝛥𝑡, one obtains a summation over values with a fixed delay. Equation ( 15) can thus be expressed as a z-domain rational expression

( 1 + 𝑐 𝑁 ∑ 𝑗=0 𝐴 𝑗+1 z -𝑗 ) 𝐹 (z) = ( 𝐾 0 + 𝐾 ∞ 𝑐 𝑁 ∑ 𝑗=0 𝐴 𝑗+1 z -𝑗 ) 𝑥(z), (20)
with 𝑐 = (𝜔 𝑐 𝛥𝑡) -𝛼 . This expression can be rewritten as the displacement to force transfer

𝐹 𝑥 (z) = 𝐾 0 + 𝑁 ∑ 𝑗=0 ( 𝐾 ∞ 𝑐𝐴 𝑗+1 ) z -𝑗 1 + 𝑁 ∑ 𝑗=0 ( 𝑐𝐴 𝑗+1 ) z -𝑗 , ( 21
)
and classical discrete to continuous approximations, such as Tustin's transform z = 1+𝑠𝛥𝑡∕2 1-𝑠𝛥𝑡∕2 can be used to recover the rational form [START_REF] Al-Hadid | Application of the force-state mapping approach to the identification of non-linear systems[END_REF] with its alternate descriptions using poles and zeros or gains and characteristic times. It is worth noting that the number of time steps used in Grünwald's approximation corresponds to the number of poles or model order or number of Maxwell branches. Fig. 9 compares the response of the rational/Maxwell model corresponding to the Grünwald approximation and the analytic response. The approximation is valid in a band between 1∕𝑁𝛥𝑡 and 1∕𝛥𝑡, which gives some control of the accuracy. The conclusion to be drawn is that Grünwald's formula is not ideal and that identification strategies [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF] used for test would probably give better result. But in the end, numerical implementations use rational fractions. The real and only advantage of fractional derivatives is linked to separating the parameter identification and order selection phases.

Rate-independent hysteresis

Hysteresis is a system model assumption, where stress is still assumed dependent of strain history. Viscoelasticity, which assumed linear system dependence on strain history fits this definition. One will thus restrict the use of the term hysteresis to effects that only depend on the history of position but not velocity, hence, the full term rate-independent hysteresis. For the applications on the present work, hysteresis will always be associated to the force/displacement, or stress/strain. As for viscoelasticity in Section 4.1, one will discuss a non-parametric form. The usual order dependent Iwan series model is then discussed in Section 5.2, before the discussion of a few order independent forms in Section 5.3.

A non-parametric model: hysteretic relaxation

In classical analyses of hysteresis [START_REF] Brokate | Hysteresis and Phase Transitions[END_REF], Madelung's rules are among those listed as necessary to ensure physically realistic behavior. These state that every loop closes itself where it started, no matter how many internal loops were made. In more precise terms, and using the illustration in Fig. 10. • any curve 𝛤 1 emanating from a turning point 𝐴 of the input-output graph is uniquely determined by the coordinates of 𝐴;

• if any point 𝐵 on the curve 𝛤 1 becomes a new turning point, then the curve 𝛤 2 originating at 𝐵 will eventually lead back to 𝐴; • if the curve 𝛤 2 is continued beyond the point 𝐴, then it coincides with the continuation of the curve 𝛤 which led to 𝐴 before the 𝛤 1 -𝛤 2 cycle was traversed.

The quasi-static test result shown, drawn from experiments detailed in Section 2, demonstrate a clear verification of the Madelung rules despite the presence of some noise and not perfect cycle repetition. The other feature of interest is the hyperelastic change of global slope around zero strain, as seen in Fig. 5. One will thus seek to introduce a non-parametric model that allows the automatic verification of these rules while being compatible with non-linear hyperelasticity.

To follow a process similar to the one used in Section 4.1 for viscoelasticity, the first proposition made here is to state that a system that verifies Madelung's rules must have a force slope, or tangent stiffness, that only depends on distance to the turning point

d𝐹 d𝑥 (𝑥) = 𝐾 𝑓 ( | | 𝑥 -𝑥 𝑡𝑢𝑟𝑛 | | ) (22)
The function 𝐾 𝑓 (𝑥) will be called hysteretic relaxation stiffness or modulus since the force is obtained through a convolution product

𝐹 (𝑥) = 𝐹 (𝑥 𝑡𝑢𝑟𝑛 ) + ∫ |𝑥-𝑥𝑡𝑢𝑟𝑛| 0 𝐾 𝑓 (𝑟)𝑑𝑟 (23)
which is similar to the relaxation modulus convolution (7) of viscoelasticity. Note however that this is not a dynamic stiffness since it depends on distance to turning point and not velocity. If force evolution is expressed by [START_REF] Feldman | Hilbert Transform Applications in Mechanical Vibration[END_REF], one can readily see that Madelung's rules are verified since the loop closure can be obtained from

0 = 𝐹 (𝐴) -𝐹 (𝐵) + 𝐹 (𝐵) -𝐹 (𝐴) = ∫ 𝑥 𝐵 -𝑥 𝐴 0 𝐾 𝑓 (|𝑥|)𝑑𝑥 + ∫ 𝑥 𝐴 -𝑥 𝐵 0 𝐾 𝑓 (|𝑥|)𝑑𝑥 (24)
Fig. 11 illustrates the different segments described during a force/displacement cycle (using the continuous STS model described in Section 5.3). When ignoring coupling with hyperelasticity, the figure illustrates that the downward and upward parts of the cycle have the same slope as consistent with [START_REF] Mullins | Softening of rubber by deformation[END_REF] but not with test shown in Fig. 5. Proper strategies for coupling with hyperelasticity will thus be discussed in the next section.

For the first branch starting from an initial null state, the same range of values is found, but with a horizontal scaling factor. Masing's law [START_REF] Masing | Eigenspannungen und verfestigung beim messing[END_REF][START_REF] Festjens | Identification and Modeling of Jointed Structures for Dynamic Analysis[END_REF], provides a classic procedure to deal with this branch. The statements are

• the first descending curve presents the same shape as the first ascending curve with an aspect ratio of two. Which considering the proposed formulation translates into the fact that when starting zero state, the convolution (23) must use 2𝑑𝑟 rather than 𝑑𝑟 for integration. • subsequent curves present the same shapes as the last one, which is readily found here.

Despite the fact that the turning point definition is simple for scalar models, or even for three dimensional models in small strains, such description is more complex for finite strains, where elasto-plastic formulations are not definitively established [START_REF] Hashiguchi | Multiplicative hyperelastic-based plasticity for finite elastoplastic deformation/sliding: A comprehensive review[END_REF], while the extension to three dimensional finite strains for viscoelasticity is a problem with established solution [START_REF] Simo | Computational Inelasticity[END_REF].

Selected order base representation: Iwan series model

As for viscoelasticity and rational transfer functions presented in Section 4.2, a discretization of the relaxation modulus is obtained using a series of first order friction elements. Using rheologic terminology, the equivalent of a Maxwell cell is a Jenkins cell where the damper is replaced by a friction element as shown in Fig. 12. Internal friction states 𝑥 𝑖 are introduced and the convolution ( 22) is replaced by a set of non-linear evolution equations on displacement rate, distinguishing sticking and sliding states as

ẋ𝑖 = 0, if ‖ ‖ 𝑥 -𝑥 𝑖‖ ‖ < 𝐹 𝑖 𝑓 ∕𝐾 𝑖 sticking state ẋ𝑖 = ẋ, if ‖ ‖ 𝑥 -𝑥 𝑖‖ ‖ = 𝐹 𝑖 𝑓 ∕𝐾 𝑖 sliding state (25)
and more realistic friction models use elasticity 𝐾 0 (potentially hyperelasticity), and a series Jenkins cells leading to a global force written as

𝐹 (𝑥) = 𝐾 0 𝑥 + 𝑁 𝑐𝑒𝑙𝑙 ∑ 𝑖=1 𝐾 𝑖 (𝑥 -𝑥 𝑖 ) = 𝑁 𝑐𝑒𝑙𝑙 ∑ 𝑖=0 𝐹 𝑖 (𝑥, 𝑥 𝑖 ) (26)
Rather than considering sliding points, the problem may be seen as load saturation, as for the Dahl model, which is more readily analyzed in a force rate formulation of the evolution equation

Ḟ 𝑖 = 𝑔 𝑖 Ḟ ∞ , if 𝐹 𝑖 sign( ẋ) < 𝐹 𝑖 𝑓 = 𝑔 𝑖 𝐾 ∞ 𝑥 𝑖 𝑓 sticking state Ḟ 𝑖 = 0, if 𝐹 𝑖 sign( ẋ) = 𝐹 𝑖 𝑓 sliding/saturated state (27)
Fig. 12 illustrates two different discretizations of the order independent STS model that will be described in the next section. Starting from a turning point, the tangent stiffness begins at a high asymptote ∑ 𝑁 𝑐𝑒𝑙𝑙 𝑖=0 𝐾 𝑖 and tends to the lower limit 𝐾 0 . Assuming the Jenkins cells to be ordered by increasing saturation forces, the hysteretic relaxation stiffness for 𝑘 sliding cells is

𝐾 𝑘 𝑏 = ∑ 𝑖≤𝑘 𝐾 𝑖 ( 28
)
with stiffness changes occurring at break points given by

𝑥 𝑖 𝑓 -𝑥 𝑡𝑢𝑟𝑛 = 𝐹 𝑖 𝑓 𝐾 𝑖 ( 29
)
In Fig. 12 the hysteretic relaxation stiffness is represented using a log scale to emphasize the parallel with the complex modulus representation of Fig. 7. Choices in the discretization strategy will be linked to order and positioning of points on the hysteretic relaxation curve. Here the 3 cell model does not fit the full data very well because the curve is set to be above the fitting model. Obviously increasing the order gives a better control on the accuracy so that a non-parametric hysteretic relaxation curve may be fitted to the desired accuracy, just as complex moduli accuracy could be controlled using more states.

When accounting for hyperelasticity, which is needed as clearly visible in Figs. 5 and 10 right, non-linear coupling needs to be considered. Fig. 13 left clearly shows that the hysteretic stiffness is much higher for a turning point at -50% strain than at 50%. Considering 𝐾 𝑖 constant is thus not appropriate. The ratio of hysteretic by hyperelastic stiffness shown right seems much more constant despite experimental limitations stiffness estimation. Considering that the hysteretic stiffness depends on the local hyperelastic stiffness thus seems appropriate. In terms of non-parametric models, the non-parametric formulation [START_REF] Feldman | Hilbert Transform Applications in Mechanical Vibration[END_REF] can be revised to use a gain depending on distance to turning point When coupled with hyperelasticity, a notable difference between the distance dependent relaxation [START_REF] Bhuiyan | A rheology model of high damping rubber bearings for seismic analysis: Identification of nonlinear viscosity[END_REF] or load saturation (27) models is that the breakpoints positions do not depend on static strain in the first case and are more widely spaced for lower hyperelastic stiffness in the second. The latter was chosen as a more realistic representation, as Fig. 13 shows that the relaxation distance is higher when going down and thus coming from a lower modulus. Fig. 14 illustrates such model, where a voluntarily small order with 3 Maxwell and 3 Iwan cells in parallel, without the non-linear transition that will be detailed later, in a force rate formulation, is considered in a triangular test.

Ḟ ℎ = 𝑔 ℎ𝑦𝑠𝑡 ( | | 𝑥 -𝑥 𝑡𝑢𝑟𝑛 | | ) Ḟ ∞ (𝑥) (30)
At turning point, velocity changes very rapidly and the hysteretic modulus is discontinuous. It starts from a high value that relaxes rapidly because of the 3 viscoelastic cells. This relaxation is clearly something that was hard to see in the considered experiment because of too coarse time sampling and the procedure should be corrected.

After viscoelastic relaxation, hysteretic relaxation is seen as three sudden drops in the hysteretic stiffness as sliding/force saturation starts. The spacing of drops on the strain axis is visibly different for the segments from -5% to 5% and from 5% to 15%. This indicates that a stress saturation model was used.

In conclusion, the proposed non-parametric hysteretic modulus characterization indicates that force rate integration with a force saturation is the most appropriate. As will be discussed in Section 6.1, this is consistent with work on regularization of the Iwan model as well as approaches considering non-linear viscosity.

Order independent parametric forms

While order dependent identification may represent tests to the desired accuracy, it may be desirable for design phases to use an order independent parametric form thus separating the issues of identification and order selection. As the Iwan model has the limitations of a selected order model, introducing order independent forms with few parameters is of clear interest for identification and modeling accuracy trade-offs. Coveney & al. [6] extended the standard triboelastic solid, shown in Fig. 15, to an infinite series with all K𝑖 and 𝐹 𝑖 equal, which will be detailed here. It is worth noting, that Segalman [START_REF] Segalman | Using series-series iwan-type models for understanding joint dynamics[END_REF] developed a four parameter model to represent bolted connections, that expresses the force using a convolution and is thus close to hysteretic modulus formulation detailed here.

The exact relation between the series STS model with K𝑗 parameters and the Iwan series with 𝐾 𝑗 is readily found using the stiffness decrease as each slider starts sliding

𝐾 𝑖 𝑏 = 𝐾 0 + ( 𝑖 ∑ 𝑗=1 ( K𝑗 ) -1 ) -1 = 𝑖 ∑ 𝑗=0 𝐾 𝑗 (31)
and the distance to the turning points, where this stiffness change occurs

𝑥 𝑖 𝑓 = 𝑥 𝑖-1 𝑓 + F 𝑖 𝑓 𝐾 𝑖 𝑏 = 𝐹 𝑖 𝑓 𝐾 𝑖 (32)
To obtain an order independent model, the model is extended to an infinite number of cells, using F 𝑖 𝑓 = 𝐹 𝐶 𝑓 → 0 and K𝑖 = 𝐾 𝐶 → ∞, for 𝑖 > 2, while keeping the product 𝐹 𝐶 𝑓 𝐾 𝐶 = 𝐶 constant. The hysteretic force is then given by

𝐹 -𝐹 𝑡𝑢𝑟𝑛 = ⎛ ⎜ ⎜ ⎜ ⎝ √ √ √ √ √ 𝑁𝐶 𝑡 |𝑥 -𝑥 𝑡𝑢𝑟𝑛 | + ( 𝑁𝐶 2 K0 ) 2 - 𝑁𝐶 2 K1 ⎞ ⎟ ⎟ ⎟ ⎠ sign(𝑥 -𝑥 𝑡𝑢𝑟𝑛 ), (33)
where 𝑁 is a variable that equals 2 when the level of force has already been reached before and 1 otherwise. This rule is quite difficult to implement, as there may be an infinite number of turning points to be stored, leading to a complicate management of internal states. Even though, this is a model with a very small number of parameters and might be useful for identification.

Ignoring the first loading conditions, one has for the hysteretic relaxation curve

𝐾(𝑥 -𝑥 𝑡𝑢𝑟𝑛 ) = 𝐾 0 + 𝐶 2 ( 𝐶(𝑥 -𝑥 𝑡𝑢𝑟𝑛 ) + 𝐶 2 4 K1 2 ) -1∕2 (34)
For this order independent hysteretic relaxation, an equivalent Iwan series model can be obtained selecting breakpoints and the asymptotic behavior, using a maximum amplitude, placing forces using [START_REF] Bourgeteau | Modélisation Numérique Des Articulations En Caoutchouc de la Liaison Au Sol Automobile En Simulation Multi-Corps Transitoire[END_REF] and using [START_REF]ABAQUS/Standard, User's Manual[END_REF] and [START_REF] Simo | Computational Inelasticity[END_REF] to obtain the Iwan parameters.

Accuracy is clearly directly linked to the choice of number and position of breakpoints. This was used for Fig. 12.

Combining non-linear effects

Intermediates between viscoelasticity and rate-independent hysteresis

Assuming rate-independence for hysteresis or exactly linear behavior for viscoelasticity are obviously idealizations and one should not expect reality to match either behavior. For this reason, several authors were interested in developing intermediate models. For instance, [START_REF] Bourgeteau | Modélisation Numérique Des Articulations En Caoutchouc de la Liaison Au Sol Automobile En Simulation Multi-Corps Transitoire[END_REF] developed the STVS model, where the Jenkins cells were also tied to a non-linear damper, [START_REF] Sjöberg | On Dynamic Properties of Rubber Isolators[END_REF] developed a model that combines non-linear elasticity, non-linear hysteresis and non-linear viscoelasticity, [START_REF] Bhuiyan | A rheology model of high damping rubber bearings for seismic analysis: Identification of nonlinear viscosity[END_REF] developed a model where these three branches were also present, [START_REF] Coveney | Rate-dependent modeling of a highly filled vulcanizate[END_REF] introduced viscoelastic dependence on the STS model, among other works.

On the hysteretic side, friction behavior is irregular, and a number of authors have proposed regularizations strategies. For example, the Dahl [START_REF] Aguirre | Parametric identification of the dahl model for large scale MR dampers[END_REF] model uses an evolution equation given by

Ḟ 𝑖 = ( 1 - 𝐹 𝑖 𝐹 𝑖 𝑓 sign( ẋ𝑖 ) ) 𝛼 𝐾 0 ẋ𝑖 , ( 35
)
with a shape parameter 𝛼 controlling the transition distance between the low force stiffness 𝐾 0 and the saturated friction force 𝐹 𝑖 𝑓 . The Bouc-Wen model [START_REF] Ismail | The hysteresis bouc-wen model, a survey[END_REF] is another differential model on force.

On the viscoelastic side, Heymans [START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF] convincingly argues that the fractional derivative model is representative of viscoelastic behavior of polymers and that the underlying characteristic times depend on a structural temperature which evolves with temperature but also strain history.

Temperature is the most classic environmental factor studied for viscoelastic behavior. Many materials verify the frequency/ temperature equivalence principle [START_REF] Nashif | Vibration Damping[END_REF][START_REF] Lakes | Viscoelastic Materials[END_REF] which states that the dependence on frequency and temperature

𝐸(𝜔, 𝑇 ) = 𝐸(𝜔 𝑟𝑒𝑑 ) = 𝐸(𝛼(𝑇 )𝜔) (36)
can be described using a single reduced frequency variable 𝜔 𝑟𝑒𝑑 which is the product of frequency and a temperature dependent scale factor 𝛼(𝑇 ). This assumption is consistent with the physical notion that frequency and temperature are equivalent in terms of particle excitation. Classical parametric expressions of the shift factor are the Williams-Landel-Ferry, and Arrhenius models [START_REF] Lakes | Viscoelastic Materials[END_REF]. The reduced frequency 𝜔 𝑟𝑒𝑑 = 𝛼(𝑇 )𝜔 can directly be used in rational models, since the temperature dependence of poles in the rational expression [START_REF] Al-Hadid | Application of the force-state mapping approach to the identification of non-linear systems[END_REF] is simply given by 𝜔 𝑖 (𝑇 ) = 𝛼(𝑇 )𝜔 𝑖 (𝑇 0 ) [START_REF] Aguirre | Parametric identification of the dahl model for large scale MR dampers[END_REF] where the reference temperature 𝑇 0 is selected such that 𝛼(𝑇 0 ) = 1.

Using a similar idea of viscous poles having a non-linear dependence on the material state, [START_REF] Höfer | Modelling of frequency-and amplitude-dependent material properties of filler-reinforced rubber[END_REF] proposed reduced relaxation times. This formulation can be rewritten as the non-linear stress relaxation equation

Ḟ 𝑖 + 𝜔 𝑖 (1 + 𝛽 ‖ ẋ‖)𝐹 𝑖 = 𝑔 𝑖 𝑔 0 Ḟ 0 (𝑥) (38)
with a velocity dependent shift factor 𝛼(‖ ẋ‖) = (1 + 𝛽 ‖ ẋ‖) related to the structural temperature discussed in [START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF]. This non-linear shift provides a saturation for a force that is linear viscoelastic at low velocities.

For 𝛽 ‖ ẋ‖ ≫ 1 and assuming a linear base branch 𝐹 0 = 𝐾 0 𝑥, the asymptotic value is a friction force

𝐹 𝑖 ≈ 𝑔 𝑖 𝐾 ∞ 𝛽𝜔 𝑖 ẋ ‖ ẋ‖ ( 39
)
leading to the expression using hysteretic constants 𝛽 = 1 𝜔 𝑖 𝑥 𝑖 𝑓 [START_REF] Vincent | Du Modèle Matériau à la Mécanique Des Systèmes : Étude Dynamique d'une Liaison Souple En Silicone Chargée de Silice[END_REF] where the use 𝑥 𝑖 𝑓 , an hysteretic relaxation distance, is significantly more explicit than the notion of reduced relaxation time. To illustrate how each branch contributes Fig. 16 displays the branch stress 𝐹 𝑖 and the instantaneous modulus first shown in Fig. 6. Two Jenkins cells are shown in blue. 𝐹 1 shown as a solid line does not slide and its instant modulus shows a nearly constant value with the small changes due to coupling with the hyperelastic behavior. 𝐹 2 shown with dashed lines has notable sliding and this induces a major drop of instant modulus during the period. The non-linear viscous version of the same load shown as a dot dashed line, follows [START_REF] Höfer | Modelling of frequency-and amplitude-dependent material properties of filler-reinforced rubber[END_REF]. The drop for higher velocities is still very much visible but the transition is much smoother.

Two Maxwell cells are shown in red. The solid curve 𝐹 3 has a relaxation frequency at 1 Hz, so the modulus is relatively high and nearly constant, especially for 15% static strain where the hyperelastic modulus varies much less than at 0% (see Fig. 6 left). The dashed 𝐹 4 has a much higher relaxation frequency at 25 Hz so the apparent modulus at 1 Hz is quite small. In conclusion, within all the possible forms coupling all three behavior of interest, the one retained in the end is a non linear relaxation of load fractions

Ḟ 𝑖 + ( 𝜔 𝑖 + ‖ ẋ‖ 𝑥 𝑖 𝑓 ) 𝐹 𝑖 = 𝑔 𝑖 𝑔 0 Ḟ 0 (𝑥) (41)
with controlling parameters being load fractions 𝑔 𝑖 , viscous frequencies 𝜔 𝑖 and hysteretic relaxation distances 𝑥 𝑖 𝑓 . Non-linear coupling of behavior occurs first because the hyperelastic force 𝐹 0 (𝑥) is a non-linear function of position, and second since the relaxation frequency in (41) has a non-linear dependence on velocity to provide an appropriate transition between viscoelasticity and rate-independent hysteresis. The later transition corresponds to known dependencies discussed in [START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF]: strain softening, anelastic deformation, stress activation, . . . Section 6.2 will now detail how test data was considered to support the choice of model form (41).

Combined form

The first coupling of interest is with hyperelasticity. When adding viscoelastic branches, it is important to test whether the strain or stress relaxation forms are more relevant. Using the strain relaxation (10), viscous loads are not proportional to hyperelastic stiffness, which changes with pre-strain. In Fig. 17 right it appears that the modulus increase due to viscoelasticity varies by nearly a factor 2 between -30% and 30% static strain. In Fig. 17 left, however, the ratio with the hyperelastic modulus still varies but only by less than 15% from the mean. The data thus indicate that stress relaxation [START_REF] Cappelli | Multi-scale identification of the elastic properties variability for composite materials through a hybrid optimisation strategy[END_REF] gives a more accurate representation, but the load fractions might need to be considered as slightly non-linear for even better accuracy.

The second coupling of interest is linked to the transition between rate-independent hysteretic and viscoelastic contributions. Using the instant modulus, Fig. 18 compares test and simulation results while focusing on amplitude variation at 1 Hz. The reproduction of trajectories is quite good and the best of all attempts: strain, strain rate, stress relaxation with various forms of non-linear load relaxation.

The instant modulus is clearly lower at nearly all instants and pre-strain levels for larger amplitude cycles. Fig. 1 illustrated this using the classical first harmonic extraction, but this figure indicates that the effect is present at all instants. When removing the dependency of poles on velocity, 𝑥 𝑖 𝑓 → ∞, the instant stiffness trajectories have the same mean. This confirms that the Payne effect is linked to rate-independent hysteresis, which here is possibly not distinguishable from non-linear effects on very low frequency viscoelasticity.

Fig. 19 focuses on frequency variations. The modulus cycles only shift up while keeping a nearly identical shape. This is consistent with the fact that at higher frequencies the behavior is essentially viscoelastic.

The results shown here use a piecewise cubic polynomial for hyperelasticity, three non-linear relaxation cells at low frequencies and three linear viscoelastic cells at higher frequencies. Refining the quality of instant modulus estimation and identification strategies is the next step. But the objective of this study was to propose objectives that would allow model discrimination of coupled non-linear models. Working with low order models showed that only model [START_REF] Österlöf | Modelling the Viscoplastic Properties of Carbon Black Filled Rubber[END_REF] had all the right trends in particular for the very discriminating instant modulus trajectory.

Conclusion

For bushing behavior, hyperelasticity is the underlying very long term effect and lowest stiffness value. Additional stiffness is first contributed by rate-independent hysteresis characterized by a hysteretic modulus converging to the hyperelastic behavior for large distances to the turning point. Close to the turning point and above moderate frequencies, viscosity cannot be ignored, as a range of combined effects exists. At higher frequencies, the modulus still increases and the linear time invariant system assumption or description using a complex modulus applies.

For non-parametric representations directly obtained from test using signal processing, hysteretic modulus obtained from triangular testing and instant complex modulus estimation were introduced and shown to provide relevant characterization of the non-linear behavior. In particular, the ordering hyperelasticity, hysteresis, viscoelasticity can be directly seen in measurements as illustrated in Fig. 1. The hysteretic modulus, while similar to approaches in [START_REF] Segalman | A four-parameter Iwan model for lap-type joints[END_REF], was measured and used to show that force rate integration with a hysteretic gain dependent on distance to turning point was most relevant. The novel instant modulus approach, where the system is considered time varying rather than non-linear, also gives a particularly detailed view of non-linear interactions present in sine tests used to characterize the Payne effect and was useful in discriminating between pertinent models.

For models, selection of order or number of internal states can be used to control discretization accuracy. For a given load fraction, added to the hyperelastic behavior, the evolution of relaxation frequencies with load gain characterizes viscoelasticity and while the evolution of turning point distances with load gain characterize hysteresis. A unified model considering both parameters can be obtained using non-linear variations of relaxation frequencies. Assuming these frequencies to be an affine function of the velocity modulus [START_REF] Österlöf | Modelling the Viscoplastic Properties of Carbon Black Filled Rubber[END_REF], as in [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF], gave good results but detailing relations with the interpretations of [START_REF] Heymans | Fractional calculus description of non-linear viscoelastic behaviour of polymers[END_REF] in terms of structural temperature, relaxation and annealing, . . . seems an important orientation of future work. While the Payne effect was clearly reproduced, improving identification procedures and analyzing other behavior such as the initial softening known as the Mullins effect are ongoing. Furthermore, many studies modify viscoelastic models to account for non-linear effects: first cycles softening [START_REF] Cantournet | Mullins effect and cyclic stress softening of filled elastomers by internal sliding and friction thermodynamics model[END_REF][START_REF] Vincent | Du Modèle Matériau à la Mécanique Des Systèmes : Étude Dynamique d'une Liaison Souple En Silicone Chargée de Silice[END_REF], plastic [START_REF] Bhuiyan | A rheology model of high damping rubber bearings for seismic analysis: Identification of nonlinear viscosity[END_REF][START_REF] Österlöf | Modelling the Viscoplastic Properties of Carbon Black Filled Rubber[END_REF][START_REF] Bardenhagen | Three-dimensional finite deformation, viscoplastic constitutive models for polymeric materials[END_REF], damage [START_REF] Diani | Observation and modeling of the anisotropic visco-hyperelastic behavior of a rubberlike material[END_REF] and amplitude effects [START_REF] Rendek | Amplitude dependence of filler-reinforced rubber: Experiments, constitutive modelling and FEM -implementation[END_REF], and clarifying the relation of the choices made here and those studies is clearly interesting.

Low parameter numbers such as fractional derivatives or continuous STS models were analyzed as intermediates allowing simple identification procedures with order selection being delayed to a later stage closer to use in simulations. Since such models, give extrapolations of relaxation frequencies and turning point distances it would seem quite interesting to analyze whether these extrapolations are relevant to estimate parameters in the intermediate range of low frequencies close to the turning point where the asymptotic viscoelastic and rate-independent hysteretic behavior cannot be separated.

This paper focused on illustrations rather than complete results. The same test campaign showed the proposed models to be relevant for shear tests, where hyperelasticity is less pronounced but viscoelastic damping higher, and part tests where geometric effects linked to large deformation play a major role. Using these tests to estimate the parameters of 3D material models is of clear interest, since this would allow prediction if part behavior in other geometric configurations or oblique loading conditions that were not tested. Finally, analyzing energy dissipation in time, seems important since hysteretic and viscous dissipated power are maximum at maximum displacement and velocity, respectively which has a visible impact on multibody simulations [START_REF] Penas | Dissipation in hysteretic rubber mount models[END_REF].

On the testing side, minimizing the viscoelastic effects at the velocity change in low speed triangular tests and better measurement sampling choices would improve hysteretic modulus estimation. For sine tests, obtaining an analytic signal model using transient Kalman filtering seems a promising approach that should be combined with the base new result that instant frequency response estimation gives a characterization allowing direct interpretations that seem much more relevant than the analysis of harmonics. The multi-step relaxation tests were considered the least useful even though very slow dependencies probably related to the Mullins effect were seen.

Fig. 1 .

 1 Fig. 1. Unified non-parametric representation of hyperelastic (solid red), hysteretic relaxation (green, dashed), viscoelastic complex moduli (blue to yellow map).

Fig. 2 .

 2 Fig. 2. Traction/compression specimen.

Fig. 3 .

 3 Fig. 3. Deformation profiles for each test: multi-step relaxation, triangular and sine testing.

Fig. 4 .

 4 Fig. 4. Multi-step test and hyperelastic behavior.

Fig. 5 .

 5 Fig. 5. Triangular traction tests. Left: stress. Right: instantaneous modulus.

Fig. 6 .

 6 Fig.6. Sine tests at 1 Hz and 1% (solid blue) or 5% strain amplitude (dashed red). Left: instantaneous modulus. Right: ratio of instantaneous and hyperelastic moduli.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Generalized Maxwell model and associated frequency domain response.

Fig. 9 .

 9 Fig. 9. Numeric implementation of truncated Grünwald's approximation of a fractional derivative model.

Fig. 10 .

 10 Fig. 10. Illustration of Madelung rules. Left: theoretic model. Right: large strain traction test.

Fig. 11 .

 11 Fig. 11. Hysteretic force/displacement cycle (left) and hysteretic stiffness (right) for model verifying Masing's law and Madelung rules.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Scheme for Iwan model and respective response in terms of hysteretic relaxation and on force/displacement domain.

Fig. 14 .

 14 Fig. 14. Simulated triangular test. Left: force displacement cycle. Right: hysteretic relaxation modulus.

Fig. 15 .

 15 Fig. 15. Standard triboelastic solid rheologic representation.

Fig. 16 .

 16 Fig. 16. Branch stress and stiffness contributions for 1 Hz excitation, 10% strain amplitude and 0 and 15% static strain.

Fig. 17 .

 17 Fig. 17. Left: viscoelastic and hysteretic modulus ratio with hyperelastic modulus. Right: modulus difference.

Fig. 18 .

 18 Fig. 18. Sine test different strain amplitudes and pre-strains at 1 Hz. Left: test. Right: model.

Fig. 19 .

 19 Fig. 19. Sine tests with different frequencies, amplitudes and pre-strains. Left: test. Right: model.

Acknowledgments

We wish to thank Vibracoustic NE, specially Pierre Charrier and Wilfried Hervouet for the realization of tests, and ANRT, France for funding the first author's thesis.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.