Introduction

Structural dynamics are used to determine the stability, safety and lifetime of mechanisms under normal operational conditions and to a great extent rely on modal analysis. Traditional modal analysis is a very well known tool used for identifying modal parameters of the structure (resonance frequencies, damping ratio and mode shapes) that consists in measuring the input and response signals, Ewins [START_REF] Ewins | Modal Testing Theory, Practice and Application[END_REF], Maia [START_REF]Theoretical and Experimental Modal Analysis[END_REF]. However, in certain cases, it is not possible to measure the excitation force, for example, during flight tests of aircraft and helicopters, and vibrations of buildings, bridges and offshore platforms.For these reasons, Operational Modal Analysis (OMA) has been developed under the following assumption: only structural responses are known, whereas excitation forces are unknown. Hermans [START_REF] Hermans | Modal testing and analysis of structures under operational conditions: industrial applications[END_REF] and Pintelon [START_REF] Pintelon | Continuous-time operational modal analysis in the presence of harmonic disturbances[END_REF], assume that excitations are random. This assumption allows to build an ideal DSP for the excitation in order to use classical modal analysis tools. Another method, usually performed to solve industrial problems, identifies peak Power Spectral Densities (PSD) and Cross-Spectrums (CS) under operating conditions [5]. These peaks are then used to estimate natural frequencies and normal modes. However this technique does not allow identifying damping with good accuracy. The "Natural excitation" technique (NEXT) [START_REF] James | Sandia National Laboratories, The natural excitation technique (NExT) for modal parameter extraction from operating structures[END_REF] considers the excitation force as a stationary white noise. Under this assumption, the Spectra and Cross Spectrum correspond to an impulse response. Modal parameters are determined from temporal responses with a complex exponential function [START_REF] Brown | Parameter estimation techniques for modal analysis[END_REF], or from the Ibrahim method (ITD) [START_REF] Ibrahim | A method for direct identification of vibration parameters from the free response[END_REF] or from the Auto Regressive -Moving Average Method (ARMAX) [START_REF] Moore | ARMAX modal parameter identification in the presence of unmeasured excitation-I: theoretical background[END_REF]. Several methods are proposed for separating periodic components from structural Natural modes, Mohanty [START_REF] Mohanty | Operational modal analysis in the presence of harmonic excitations[END_REF]. Filtering technics are often used in addition to these methods [START_REF] Parks | Digital Filter Design[END_REF]. Jacobsen [START_REF] Jacobsen | Using enhanced frequency domain decomposition as a robust technique to harmonic excitation in operational modal analysis[END_REF] defines an indicator using the spectral description of the signal. A linear interpolation is performed on spectral responses in order to eliminate the periodic components identified. Nonetheless, this method is not suitable when the harmonic component is close to a low-damped Natural mode. When Transmissibility is used in OMA, FRF must be known at several points and excitation force levels determined [START_REF] Devriendt | The use of transmissibility measurements in output-only modal analysis[END_REF], [START_REF] Devriendt | An operational modal analysis approach based on parametrically identified multivariable transmissibilities[END_REF]. One of the major difficulties in OMA is the coexistence in one signal of structural modal responses with periodic excitations due to the operation of machines fixed on the structure. This situation is often observed in transport technologies and large industrial structures.

Antoni et al. [START_REF] Garibaldi | New separation techniques for output-only modal analysis[END_REF] propose two methods for harmonic detection. The first one is based on envelope statistics and is especially applied on pure tone components. The second one is based on correlation in conjunction with an indicator built with a tuned time delay (correlation time indicator). This method is justified by the properties of the correlation (and the PSD) of random noise and periodic signals. This second method is especially suitable for non frequency modulated harmonic components. Antoni [START_REF]Blind separation of vibration components : principles and demonstrations[END_REF], [START_REF]The spectral kurtosis: a useful tool for characterizing non-stationary signals[END_REF] and Vrabie [START_REF] Vrabie | Spectral Kurtosis: From Definition To Application[END_REF], have defined and studied properties of Spectral Kurtosis respectively for non stationary signals and stationary signals. Both methods are based on the properties of Welch's method and compute the Spectral Kurtosis in the frequency domain. Another definition of the spectral Kurtosis, defined by Brincker [START_REF] Brincker | An Indicator for Separation of Structural and Harmonic Modes in Output-Only Modal Testing[END_REF] and Jacobsen [START_REF] Jacobsen | Using enhanced frequency domain decomposition as a robust technique to harmonic excitation in operational modal analysis[END_REF], [START_REF] Jacobsen | Using EFDD as a Robust Technique to Deterministic Excitation in Operational Modal Analysis[END_REF], [START_REF] Jacobsen | Eliminating the Influence of Harmonic Components in Operational Modal Analysis[END_REF] is based on the classical approach of the Kurtosis in a specific bandwidth. The present work describes a method developed in order to separate structural modal responses from periodic excitation due to rotor dynamics without assumptions regarding periodic excitation In the continuity of previous cited works, the objective of this paper is to study and define the necessary conditions of robust method for detection with Spectral Kurtosis. This method is applied to harmonic components detection, especially for frequency and amplitude modulated components produced by rotor dynamics. Theoretical definitions of the SK used in cited references are presented in section 2. Section 3 presents the proposed method so called OSK. The several SK definitions and methods are discussed and compared with the proposed Optimized Spectral Kurtosis method (OSK) in section 4. The proposed OSK and the different definitions of SK are computed from experimental results and presented in section 5. The proposed method is different from the methods used in references [START_REF]Blind separation of vibration components : principles and demonstrations[END_REF]- [START_REF] Vrabie | Spectral Kurtosis: From Definition To Application[END_REF] because of the definition of the Spectral Kurtosis. In these references, the Kurtosis is computed from spectral data: the temporal data are first computed with the FFT, and secondly the Kurtosis is computed from the DFT results. In the proposed method, the temporal signal is first filtered in the time domain with a series of sharp Pass Band Filters and secondly the Kurtosis is computed from the filtered temporal signals for each central frequency of PBF. In order to differentiate the two methods, we have called this new method the "Optimized Spectral Kurtosis" because in this new method the spectral resolution step of the FFT process is dissociated from the equivalent bandwidth used to compute the Kurtosis. This difference allows a better efficiency in the detection process. However, the efficiency of this new method depends on several parameters as the number of samples or the filtering technique. These parameters are discussed and optimized. This is the second main interest of this paper.

Statistical tools and Spectral Kurtosis

The properties of the probability density (PD) of a discrete signal make it possible to assume the nature of the initial temporal signal. However the same probability density could be produced by an infinite number of signals. Certain deterministic signals, such as sinusoidal ones, produce singular probability densities. The first step of the method is based on the differences observed on the PD between (Gaussian) random signals and sinusoidal signals.

In order to highlight these differences, two descriptors were selected: Peak Factor (PF) (Eq.1) and Kurtosis.

{ } { } { } ( ) N RMS n Max u[ n ] Max u[ n ] PF U u[ n ] u[ n ] N = = = ∑ 2 1 1 (1)
With: u[n]: the temporal sequence of the discrete signal U RMS : the Root Mean Square level of the discrete signal N: number of points of the temporal sequence The first four statistical moments allow defining the mean µ (expectation), variance (or σ standard deviation), Skewness and Kurtosis. These descriptors about its mean and central moments are reduced: central standardized moments are divided by the standard deviation (for orders greater than 2).

[ ]

E X = 1 (2a) E ( X -)   =   2 2 (2b) X - E σ     =           3 3 (2c) X - E σ     =           4 4 (2d)
With: µ 1 =µ: the mean σ: the standard deviation:

σ 2 =µ 2
Kurtosis is defined by K=µ 4 and excess Kurtosis is defined by:

N K = K -3 .
In the three signals studied, the expected values for these descriptors are presented in table 1 E{Estimator}\ Signals

Uniform Random [-0,5 0,5] Normal Random (0,1) Sinusoidal [-0,5 0,5] Mean: µ 1 0 0 0 Peak Factor: PF ( ) 1 2 12 - 2 Variance: σ² = µ 2 1/12 1 1/8 Kurtosis: K 9/5 3 3

/2 TABLE 1: EXPECTED VALUES OF STATISTICAL DESCRIPTORS FOR THE THREE CLASSICAL SIGNALS

Brincker [START_REF] Brincker | An Indicator for Separation of Structural and Harmonic Modes in Output-Only Modal Testing[END_REF] and Jacobsen [START_REF] Jacobsen | Using EFDD as a Robust Technique to Deterministic Excitation in Operational Modal Analysis[END_REF], [START_REF] Jacobsen | Eliminating the Influence of Harmonic Components in Operational Modal Analysis[END_REF] define the Spectral Kurtosis with the same definition used previously for the Kurtosis based on the forth statistical moment (Eq.2d) applied in a spectral bandwidth. For signal s, the Spectral Kurtosis of the frequency f has been previously defined and used by Antoni [START_REF]The spectral kurtosis: a useful tool for characterizing non-stationary signals[END_REF] and Vrabie [START_REF] Vrabie | Spectral Kurtosis: From Definition To Application[END_REF] in a similar way with spectral cumulants (C) and moments (S) of a signal divided into M unoverlapped blocks:

( ) ( ) ( ) ( ) ( ) 4 4 2 2 2 2 2 = - ≜ S S S S S C S f f SK f S S f f (3)
Vrabie [START_REF] Vrabie | Spectral Kurtosis: From Definition To Application[END_REF] defines the unbiased estimation of the Spectral Kurtosis:

( ) ( ) ( ) ( ) 4 1 2 2 1 1 2 1 = =   +     - =   -           ∑ ∑ M i i S M i i S f M M SK f M S f (4)
This estimation allows defining the efficiency of detection for a signal composed of a random process and a harmonic component.

Description of the method

The main idea is based on a spectral observation of statistical descriptors. Previous works [START_REF]Blind separation of vibration components : principles and demonstrations[END_REF]- [START_REF] Jacobsen | Eliminating the Influence of Harmonic Components in Operational Modal Analysis[END_REF], have been performed in defining and using the Spectral Kurtosis. Detections are defined from numerical differences of statistical descriptors between random signals and sinusoidal signals. Also, although detections are not certain they are very likely.

Main assumptions

The excitations highlighting normal vibration modes are assumed to be mesokurtic or leptokurtik (positive or null excess Kurtosis).

In numerical applications, the method is still pertinent provided that excess Kurtosis does not exceed -1.

Global description of the method

If the signal is composed of a sinus and a random part, Kurtosis is assumed to be between 1.5 (sinus) and 3 (Normal Random). Under the same assumptions, the PF moves within a wider range (between 2 to 5 or more).

A numerical study of sensitivity of these descriptors highlights better detection with Kurtosis than with PF. Despite a narrower range of evolution, Kurtosis proves to be more robust and reliable which is why this study was performed with Kurtosis. The method is composed of two main steps: firstly, the harmonic components have to be detected (Section 4) and, secondly, they have to be eliminated.

The detection of harmonic components is based on combining two kinds of information: statistical and spectral. The OSK consists in establishing a spectral description of Kurtosis. First, the signal is filtered with a narrow Pass Band Filter (PBF) centered on the frequency studied. Secondly, the Kurtosis of the filtered signal is computed and the numerical result defines the assumed nature of the signal (included in the spectral bandwidth). The complete spectral description of Kurtosis is obtained by a series of PBF along the entire frequency bandwidth studied. The optimization of the harmonic component detection with OSK had to satisfy several questions given in section 4: 1 -What is the optimal number of samples? 2 -What kind of filter should be chosen? 3 -How can the PBF be translated along the entire bandwidth? 4 -What bandwidth is best? 5 -What are the limits of the method?

Detailed description and critical analysis of the OSK

First, the OSK was performed with a synthetic random signal added with a sinusoid. The random signal is a stationary Gaussian random white noise with mean µ=0. The main objective of this initial approach is to highlight the efficiency and limits of application of the method in terms of optimum filter bandwidth and Signal to Noise Ratio (SNR). Secondly, the OSK was performed with a synthetic random signal r(t) used as an excitation source applied on a system h(t) composed of two Natural modes in the frequency bandwidth studied. This system is considered as equivalent to a mechanical structure. The output signal s(t) is disturbed by a periodic signal composed of 3 harmonics. This disturbing periodic signal is considered equivalent to a rotating machine in operation and fixed on the structure (figure 1).

FIGURE 1: BLOCK DIAGRAM OF A TRANSFER FUNCTION WITH A POLLUTED OUTPUT p s [ n ] s[ n ] p[ n ] r [ n ] * h[ n ] p[ n ]

= + = + (

Field of the numerical study

The OSK was applied to a signal with a number of samples between 2 14 (16 284 samples) and 2 22 (4 194 304 samples).

An important interest of the present work consists in distinguishing the frequency resolution of the DFT process and the frequency bandwidth in which the OSK is applied. The relation between this two frequency bandwidths is studied and optimized in order to accurate the detection process. In previous studies [START_REF]The spectral kurtosis: a useful tool for characterizing non-stationary signals[END_REF], [START_REF] Vrabie | Spectral Kurtosis: From Definition To Application[END_REF] the Welch's method is used for spectral computations. In Welch's method, the two kinds of bandwidth (spectral resolution and filter bandwidth) are coupled (and equal). In that case, the Welch's Method allows a rapid computation with the FFT. In the present study signal had to be filtered in the time domain with an ARMA sequence before the Kurtosis computation. The Cauer Filter [START_REF] Parks | Digital Filter Design[END_REF] (used for the ARMA filter synthesis) has been chosen in order to keep a sharp bandwidth during the optimization process of the method. The sixth order results from a compromise between numerical costs and the efficiency of the filter. However if the PSD and the OSK have the same spectral resolution (equivalent bandwidth) and when the optimal filter is chosen the Welch's Method could be performed and thus allows an important reduction of the time computation. The PBF was applied with a constant frequency bandwidth. The OSK had already been performed with n th octave analysis (with frequency translation and down sampling) but, the detection obtained was less effective than with a constant spectral bandwidth. In order to reduce time computation the series of filters (with fixed bandwidth) are computed before the computation of the OSK.

The study was performed with a frequency bandwidth between ( )

0 0 2 , 2 2
∆ ∆   - +     = k f k f f f A SNR PSD r (6)
With sinusoidal amplitude A. The Power Spectral Density of the random signal r is defined for the k∆f Bandwidth.The Transfer Function used in the numerical study is shown in figure 2.

FIGURE 2: TEST SPECTRAL FUNCTION: A 2 DOF SYSTEM EXCITED WITH A RANDOM NOISE AND A PERIODIC SIGNAL COMPOSED WITH 3 FREQUENCIES REPRESENTED BY THE 3 ARROWS r[n] h[n] s p [n] s[n] p[n]

The harmonic components of the signal are 123Hz, 246 Hz and 492 Hz. The Sampling Frequency was 1600 Hz and the spectral resolution 6.1mHz with a signal composed of 2 18 samples (262 144 samples).

OSK in operation

Figure 3 represents the evolution of Kurtosis along the spectral domain studied. The frequency bandwidth of the filter is set at 16∆f which is less than 0.1Hz. In this case, the detection of the three harmonic components is clearly highlighted.

FIGURE 3: OSK FOR SEVERAL NUMBERS OF SAMPLES

The method was applied to a signal described with 2 18 samples for several values of SNR and several spectral bandwidths of the PBF (figure 3). Detection with OSK is represented on the map in figure 4. The efficiency qualification is based on the difference (in terms of standard deviation σ SK ) between the mean of the OSK µ SK and the measured OSK: "Perfect" is defined for OSK< µ SK -5 σ SK , "Very good", for µ SK -3 σ SK >OSK> µ SK -5 σ SK , "Good" for µ SK -2 σ SK >OSK> µ SK -3 σ SK, "Intermediate" for µ SK -σ SK >OSK> µ SK -2 σ SK , "Bad" for µ SK -σ SK /2>OSK<>µ SK -σ SK , "Very bad" for µ SK -σ SK /4 >OSK> µ SK -σ SK /2 "Null" for OSK in the neighborhood of the mean or greater.

Perfect

Very Good Good Intermediate Bad

FIGURE 4: EFFICIENCY OF HARMONIC COMPONENT DETECTION FOR 218 SAMPLES AGAINST FILTER FREQUENCY BANDWIDTH AND SNR. THE BEST OSK IS OBTAINED FOR 20 F BANDWIDTH (AND SNR GREATER THAN 3)

On figure 4, larger and darker the circle, better the expected detection. Figure 4 is built with a synthetic random signal added to a sinus and highlights: A field in which detection is impossible for a filter bandwidth lower than or equal to 4∆f, the best detection is obtained with a filter bandwidth close to 20 ∆f, certain detection of SNR greater than 3 and with a filter bandwidth equal to 20∆f. Figure 5 shows the same kind of chart as figure 4. The limit of detection with OSK is shown for several numbers of samples.

The larger the number of samples, the better the detection. However, accurate detection could only be obtained with a large number of samples. A reasonable compromise was found based on several numerical studies: the number of samples had to be chosen between 2 16 and 2 20 .

FIGURE 5: CHART OF THE DETECTION LIMITS FOR SEVERAL NUMBERS OF SAMPLES WHITE ZONES SHOW VERY GOOD DETECTIONS DARK ZONES SHOW POOR DETECTION

This chart has not been completely computed. Only few points with different number of samples have been computed in addition of figure 5. "Perfect", "Very good" and "Good" detections are considered with the same color in this chart. The more the number of samples increased, the better the detection performed.

The optimal bandwidth of the filter is close to 20∆f and do not depends on the number of samples. The detection could not be obtained for Bandwidth smaller than 4∆f. This dark zone does not depend on the number of samples. Over 20∆f, the larger the bandwidth, the worse the detection.

Experimental test bench

Description of the experiment

In order to improve the OSK in a real experiment, a test bench was set up. The experiment is composed of two floors insulated with flexible blades (Figure 6). The mechanical system could be described as a two DOF system in the bandwidth 0-100 Hz. The natural frequencies of the system are 9.7 Hz and 36.7 Hz. The structure is first excited with an electro-dynamic shaker and a random noise, then with an unbalanced electric engine rotating at fb=22Hz amplitude and frequency modulated in a short range. The both excitations are realized in the same time. The aim of this experimental study is to test the efficiency of the OSK. The expected results consists in detecting the harmonic component around fb=22Hz and in not detecting them at 9.7Hz and 36.7Hz. These two components had to be identified as natural frequencies and not as periodic excitations.

Experimental results

In agreement with previous numerical studies, the method is performed with 2 18 samples and Cauer filters with a bandwidth of 20∆f. For a sampling frequency of 2000 Hz, the filter bandwidth is around 0.15 Hz. Detection is defined for an OSK lower than 2.

FIGURE 7: DETECTION OF HARMONIC COMPONENTS USING OSK 1 FREQUENCY EXPECTED BUT 4 FREQUENCIES DETECTED

The results shown in figure 7 highlight the following expected observations: Frequency fb=22 Hz is well detected as a harmonic component and not as a structural response. The 2 structural Eigen modes (9.7 Hz and 36.7 Hz) are not identified as harmonic components. Nevertheless, other components, not initially foreseen, are identified as harmonic components: Electric network frequency fr=50 Hz, Two frequencies: fr-fb=28 Hz and fr+fb=72Hz. These components are very small but effectively included in the signal. The first is due to the electric power frequency used for the electric engine. The two others are due to magnitude modulation between the rotation frequency of the engine and its electric power frequency. These components are 100 to 100000 times smaller than the main excitation at fb=22 Hz. The unbiased Spectral Kurtosis SK* in figure 8 and the OSK figure 7 have been computed with the same record. SK* in figure 8 has been computed with 8, 32 and 128 unoverlapped blocs and with corresponding Hanning windows (respectively 32768, 8192 and 2048 samples). The choice of a large number of unoverlapped blocs lead, on the one hand, to an accurate detection but, on the other hand, to a worse spectral resolution for the DSP and the SK*. Due to little modulations of harmonic components the 2 main frequencies (22Hz and 50 Hz) are hardly detected and the two others (28 Hz and 72 Hz) are not detected : 2 minima are detected with the SK* in figure 8 and 4 minima are detected with OSK in figure 7. [START_REF] Vrabie | Spectral Kurtosis: From Definition To Application[END_REF] AND ANTONI [START_REF]The spectral kurtosis: a useful tool for characterizing non-stationary signals[END_REF] FOR N SEVERAL NUMBERS OF UNOVERLAPPED BLOCS.

FIGURE 8: DETECTION OF HARMONIC COMPONENTS WITH UNBIASED SPECTRAL KURTOSIS SK* USED BY VRABIE

In this experimental study, the lack of efficiency for the SK* is due to the use of the same bandwidth between the DFT (and PSD) and the SK. Separating the two kind of bandwidth needs to compute the SK in the time domain as it is proposed in the present paper for the OSK. The detection capability is more accurate with the OSK figure 7 than with the SK* in figure 8. This is the main interest of the proposed method. Antoni et al. [START_REF] Garibaldi | New separation techniques for output-only modal analysis[END_REF] propose two others detection methods for harmonic detection. The first one is based on envelope statistics and is especially applied on pure tone components. However this method could not be applied for amplitude modulated components. The second one is based on correlation in conjunction with an indicator built with a tuned time delay (correlation time indicator). This method is justified by the properties of the correlation (and the PSD) of random noise and periodic signals. However, this method is not appropriated in the case of frequency modulated signals (even with small variations of the instantaneous frequency). Due to small modulations of the 2 harmonic components (22Hz and 50Hz) this 2 previous methods are inappropriate for the experimental record presented in section 5.

The OSK shows that it is capable of performing beyond the initial detection target defined in this study and proves to be very efficient for harmonic component detection.

Conclusion

Harmonic components play a role in the dynamical response of many rotating mechanical systems. However certain OMA techniques are unable to detect these kinds of signal. The optimized OSK proposed here could detect harmonic components even when the latter were 1000 times smaller than the magnitudes observed around natural frequencies. The OSK does not require initial knowledge of excitation forces. What is more, the natural vibration mode parameters could be obtained from only one measurement (one location) of the dynamical response. The results of the experimental tests and comparisons with other techniques confirm the method's efficiency and the pertinence of the choice of optimum parameters in the detection process. However, as all cited methods, the OSK requires a large number of samples in order to perform the detection and its numerical computation is still costly.

5 )

 5 With: r[n]: Input Signal h[n]:Impulse Response of the system s[n]: Unpolluted Output Signal p[n]: Disturbing Signal pollution (periodic or pseudo-periodic) s p [n]: Polluted Output Signal 4.1.

∆f

  and 50∆f (with ∆f being the spectral resolution of the Fast Fourier Transform (FFT) defined by ∆f=1/τ and τ the duration of the signal sampled) The signals were studied with an SNR η between 1 and 50. The SNR is defined only for the frequency bandwidth containing the sinusoidal component. (For other spectral bandwidth, the SNR tends towards 0).

FIGURE 6 . 1 :FIGURE 6

 616 FIGURE 6.1: PHOTOGRAPH OF THE EXPERIMENTS


Franck Renaud
email: franck.renaud@supmeca.fr

Gael Chevallier

Jean-Luc Dion

R Émi Lemaire

DRAFT: VISCOELASTICITY MEASUREMENT AND IDENTIFICATION OF VISCOELASTIC PARAMETRIC MODELS

Generally speaking, the behaviour of viscoelastic material is more complicated than the behaviour proposed by classical models as Voigt, Maxwell or Zener. The stiffness of such materials is a frequency dependent complex function. Above 1000Hz, classical measurements techniques are unable to achieve accurate measurements of the stiffness. In this paper, a new Dynamical Mechanical Analysis (DMA) tester is presented. It allows the characterization of the shear stiffness of preloaded viscoelastic materials between 200 and 3500Hz and without using timetemperature equivalences. Then the Generalized Maxwell model is used to describe behaviours measured with the DMA tester. A new iterative identification method of the parameter of the Generalized Maxwell model is described. This identification method is based on the asymptotes of the model.

NOMENCLATURE

σ Stress. ε Strain. h(t) Relaxation function. δ X Displacement. N Force factor. δV Tension. F Force. Q Power of the driver. f Frequency of excitation. C e Electrical capacity. M Mass. K Stiffness. C Damping. ω Angular frequency. χ = log(ω) Logarithm of the angular frequency. ω p Pole. ω z Zero. S Stiffening. A Area.

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 viscoelastic-

Copyright c 2011 by ASME ity 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 [START_REF] Caputo | Linear models of dissipation in anelastic solids[END_REF], Lakes [START_REF] Lakes | Viscoelastic Solids[END_REF], Chevalier [START_REF] Chevalier | Mechanics of Viscoelastic Materials and Wave Dispersion[END_REF]. Linear viscoelasticity assumes the existence of a relation between stress, σ , of a material and its strain, ε, history. Let us call h(t), the relaxation function which is also an element of the complex stiffness tensor. Linear viscoelasticity is defined by the Equation [START_REF] Ferry | Viscoelastic Properties Of Polymers[END_REF]. The Fourier transform is denoted with a hat and H = ĥ(ω).

     σ (t) = t -∞ h(t -τ)ε(τ)dτ σ (t) = H ε(ω) (1)
Viscoelastic behaviour of materials is difficult to measure over a large frequential bandwidth. For that, several methods of characterization exist. Oberst and Frankenfeld [START_REF] Oberst | Damping of the bending vibrations of thin laminated metal beams connected through adherent layer[END_REF] proposed to study the first mode of a sandwich beam made up of metal skins and of a viscoelastic heart. Their method makes it possible to know the damping induced by the viscoelastic heart at the frequency of the mode. Several authors use this kind of technique for the identification of viscoelastic parameters, Barbosa and Farage [START_REF] Barbosa | A finite element model for sandwich viscoelastic beams: Experimental and numerical assessment[END_REF] and Castello et al. [START_REF] Castello | Constitutive parameter estimation of a viscoelastic model with internal variables[END_REF]. It is also possible to deduce the mechanical properties from a viscoelastic material starting from the measurement of the various Eigen frequencies of a sample having a simple form, a beam for example, see Chevalier [START_REF] Chevalier | Essais dynamiques sur composites. caractrisation aux basses frquences[END_REF]. These methods are equivalent with the deduction of the mechanical properties starting from Frequency Response Function (FRF). This type of methods allows only the characterization of the frequencies of the Eigen modes and not on a large frequential bandwidth. Moreover these methods are valid only under the linear material assumption in amplitude of excitation. Chen [START_REF] Chen | Determining a prony series for a viscoelastic material from time varying strain data[END_REF] proposed to directly measure the functions of relieving and creep to deduce the coefficients from them from a series of Prony. However, it is difficult to impose a level of perfect displacement, thus this method, which gives access to the behaviour on a large frequential bandwidth, is accurate only at low frequencies. For the characterization of viscoelastic materials, the Dynamic Mechanical Analysis (DMA) is the most suited technique. It allows one to measure the stress and the deformation of materials and thus to calculate the dynamic rigidity of the samples according to the frequency. DMA testers work into quasi-static, i.e. at frequencies lower than the first Eigen mode of the test bench. Currently, the best DMA testers allow characterizations, without using time-frequency equivalences, until approximately 1000Hz. Thus there is a need of a new DMA tester in order to measure the viscoelastic behaviour at higher frequencies in quasi-static excitation.

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 [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 accurate. Many rheological models have been proposed for H, as the constant complex modulus, Maxwell, Voigt, Zener, etc (see Renaud et al. [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF]). In the present paper Generalized Maxwell model is used. The parameters of this model will be identified from the curve of H against the angular frequency, see Renaud et al. [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF].

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. The identification method presented here has been detailed in the article of Renaud et al. [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF]. This method uses the characteristics of the asymptotes of Pole-Zero Formula- Copyright c 2011 by ASME tions and allows one to identify parameters from both the modulus and the phase curves.

PIEZOELECTRIC VISCOANALYZER (DMA) CONCEPT Shear measurement principle

The DMA tester has been designed to characterize the dynamic shear stiffness of the samples under normal preload and in the frequential band between 1 to 5000Hz. The DMA tester works in quasi-static excitation, which means at any frequency below the first Eigen mode of the tester. The DMA tester can reach a maximum normal preload of 4MPa. Moreover this preload must be uniform on the whole sample surface.

The design of the DMA tester is shown on the Figure 2. The samples are gripped between the jaws and the holders of samples. Due to symmetry, this tester needs four samples for characterization. In contrary to the Oberst beam test and thanks to the symmetry principle, this DMA tester allows pure dynamical shear excitation without inflection. The normal preload is applied with a unique M6 bolt. The preload is limited to 5000N in order not to plasticize the jaws. The sample size, 20*30mm, is thus selected to ensure a pressure of 4MPa with a preload of 5000N.

The shearing of the samples is obtained by moving the supports in phase opposition. For that, two chains of actuation are placed between the holders. Thanks to the symmetry, there is no need to clamp the chains of actuation on a frame. These chains consist in a piezoelectric actuator, a connector and a force sensor.

Quasi static excitation

The piezoelectric actuators make it possible to reach high frequencies higher than 10 kHz. The frequential limit is set by the frequency of the first Eigen mode of the DMA tester. So this frequency must appear as high as possible. Consequently, the DMA tester is suspended by rubber bands to avoid Eigen modes coming from the frame. Actually, the tested viscoelastic materials have so low shear stiffness that the first Eigen mode is due to the mass of the holders and the actuation chains suspended by the tested samples. This Eigen mode generally appears around 4100Hz and limits the characterization to a frequency of 3500Hz.

The jaws are made in aluminium. Their design has been optimized so that the unique preloading bolt generates a uniform preload of the samples on their whole surface. This design also ensures the decoupling between the normal load and the shear load on the samples. The first Eigen mode of this jaw appears around 7500Hz.

The Figure 3 shows the capabilities of each piezoelectric actuator according to the Equation (2), with δ X the displacement, F the force, K the stiffness of the actuator, δV the tension and N the force factor. The continuous line shows the forcedisplacement behavior of actuators, for the maximum electrical

δ X = (N δV -F)/K (2)
When their ends are blocked, the displacement is nil and the maximum pushing force they can generates is 800N. If their ends are free, the force is nil and the maximum displacement is δ X = 15 m. Unfortunately, the actuators are similar to electric capacities, thus the more the frequency of excitation is raised and the more the actuators need current. Since, the power of the actuators driver is limited to Q = 30 Watts, the current is also limited, see Equation ( 4) and the real behaviors of actuators at f = 1000 Hz and f = 5000 Hz are shown in Figure 3.

Q = 2π fC e V 2 p (3)

Test post processing

Basically, the complex shear stiffness of the samples is the ratio of the force spectra to the displacement spectra. However, the quasi-static excitation is applied until 3500Hz and some inertia effect must be taken into account in the post processing. Thus a simplified model of the DMA tester has been set, see Figure 4. The jaws are modeled by two identical masses, M P , connected between them by a spring, K jaw . The masses M P represent the moving part of the jaw which undergoes the accelerations read by the accelerometers. Since, the DMA tester is quasi-symmetrical between the sides left and right, the central part of the jaws which constitutes the symmetry plane, does not move. Consequently, the masses of the screw M6 and the sensor of force do not move and one has M P1 = M P2 = M jaw /3 + M accelerometer . The spring K jaw is the same one for the two jaws. It was evaluated starting from a calculation by finite elements. The two chains of actuation are identical and are made up of a sensor of force, a connector and a piezoelectric actuator. These two chains of actuation work in phase and push with the same force, thus they are regarded as one equivalent chain of actuation generating an equivalent force which is the sum of the forces of both F piezo = F 0 + F 1 . However, because of the efforts of inertia, the force measured on the level of the sensors, is neither equal to the force generated by the actuators nor equal to the force undergone by a sample. In order to take into account these efforts of inertia, the chains of actuation equivalent is modeled by three different masses M S1 , the mass of the support of left, M R the mass of the connector and M S2 the mass of the support of right-hand side. F FS is the force measured by the force sensors.

It is worth noticing that because of the efforts of inertia of the connector, the force generated by the piezoelectric actuators is F piezo = M R ü3 + F FS . The samples are modeled by pure stiffnesses, K S . The equations of the dynamics of the system are given by the Equations ( 4).

(a) (b) (c) (d) (e) ( f )                M P1 ü2 = -K jaw (u 2 + u 5 ) + K S (u 3 -u 2 ) M S1 ü3 = -K S (2u 3 -u 2 -u 4 ) + F FS M P2 ü4 = -K jaw (u 4 + u 7 ) + K S (u 3 -u 4 ) M P1 ü5 = -K jaw (u 2 + u 5 ) + K S (u 6 -u 5 ) M S2 ü6 = -K S (2u 6 -u 5 -u 7 ) + F FS + M R ü3 M P2 ü7 = -K jaw (u 4 + u 7 ) + K S (u 6 -u 7 ) (4)
By calculating the sum of the relations (b)+(e), one obtains the stiffness, K S , of the samples, see Equation ( 5)

K S = (M S1 -M R ) ü3 + M S2 ü6 -2F FS u 2 + u 4 + u 5 + u 7 -2u 3 -2u 6 (5)
Equation ( 5) is not used during post processing because displacements are calculated by double integration of accelerometer measurements. This double integration is carried out by the formula ü = -ω 2 u. Thus, the Equation is replaced by the Equation [START_REF] Oberst | Damping of the bending vibrations of thin laminated metal beams connected through adherent layer[END_REF].

K S = ω 2 (M S1 -M R ) ü3 + M S2 ü6 -2(F 0 + F 1 ) 2 ü3 + 2 ü6 -ü2 -ü4 -ü5 -ü7 (6)

Results on industrial specimens

The Figure 5 shows the complex shear stiffness of a multilayer viscoelastic material for two normal preloads : 500 and 5000N. This material is made of layers of glue, rubber and steel. One can see the stiffening effect of viscoelasticity. This stiffening effect cannot be taken into account by the common modal damping model. On the contrary, the Generalized Maxwell model is able to describe such behaviour.

IDENTIFICATION OF THE COMPLEX SHEAR MODU-LUS

The Generalized Maxwell model is a rheological model described by the Figure 6 and the Equation [START_REF] Barbosa | A finite element model for sandwich viscoelastic beams: Experimental and numerical assessment[END_REF]. It is composed of N Maxwell cells, with K 0 , the static stiffness, K i the stiffness in the i th Maxwell cell and C i the damping in the i th Maxwell cell.

Z(ω) = K 0 + N ∑ i=1 jωK i C i K i + jωC i (7) Z(ω) = K 0 N ∏ i=1 1 + ( jω/ω z,i ) 1 + ( jω/ω p,i ) (8)
       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 (9)
The identification of Generalized Maxwell parameters using a pole zero formulation has been fully detailed in the article of Renaud [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF]. Let us summarize the main steps. First, it is worth noticing that in the logarithmic scale, the behaviour of a pole zero formulation is the sum of the behaviour of each pole zero couple. The modulus and the phase generated by one pole zero couple in logarithmic scale is given by the Equation [START_REF] Chen | Determining a prony series for a viscoelastic material from time varying strain data[END_REF].

         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 (10)
Since this behaviour is non-linear against frequency, we will approach them with asymptotes, see Figure 7.

Let us define the stiffening S(|Z|, χ a , χ b ) and the area under the phase curve A(ϕ, χ a , χ b ) by the Equation [START_REF] Soula | Measurements of isothermal complex moduli of viscoelastic materials over a large range of frequencies[END_REF]. The identification of poles and zeros will be led by equating the stiffening and the area of experimentally measured curves and the stiffening and the area of the model.

5

Copyright c 2011 by ASME

∀ χ a = log(ω a ) ∈ R and χ b = log(ω b ) ∈ R,                      S(|Z|, χ a , χ b ) = χ b χ a ∂ ∂ χ log(|Z(χ)|) dχ = log(|Z(χ b )|) -log(Z|(χ a )|) A(ϕ, χ a , χ b ) = χ b χ a ϕ(χ) dχ (11)
The stiffening and the area of the asymptotes are given by the Equation ( 12) and are easy to calculate. They are proportional to the gap between the pole and the zero: log(ω p,i )log(ω z,i ). Thus, by the knowledge of the stiffening and the area of experimentally measured curves, it is easy to obtain the gap between the pole and the zero.

S(|Z asy

(χ)| i , χ a , χ b ) = log(ω p,i ) -log(ω z,i ) A(φ asy i (χ), χ a , χ b ) = π 2 (log(ω p,i ) -log(ω z,i )) (12)
The identification begins by an initialisation step. It consists in an apportionment of the frequency domain so that each sub domain has the same stiffening and phase area. Then the Equation [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF] allows the calculation of one pole-zero couple in each sub domain. The initialisation step gives a first set of parameters which is optimized in a second step.

The difference between experimental curves and the model obtained after initialisation can be seen as local lacks and local excesses of stiffening and area. These local lacks and excesses are translated into difference between pole and zero, log(ω p,i )log(ω z,i ) by the Equation [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF]. It allows the adjustment of initial parameters by iteration.

For example, a local excess of stiffening appearing on the right of a pole-zero couple means that the pole has to tend towards the zero. On the contrary if this excess of stiffening appears on the left of a pole-zero couple, the zero has to tend towards the pole. Moreover, all pole-zero couples of the model are able to correct a local excess of stiffening, so all pole-zero couples of the model have to be adjusted. Actually, the behavior of the model at a frequency far from a pole-zero couple is quasi independent of this couple. Indeed, the contribution of a pole-zero couple, at a frequency far from it, is very low, see Figure 7. So the adjustment of a pole zero couple, because of a local difference, is weighted according to its local contribution This identification technique is efficient, see Figure 8. Unlike some others techniques based on non-linear optimization, it allows constraining the values of the zeros to be lower than the values of the poles to ensure the causality of the model, see Renaud [START_REF] Renaud | A new identification method of viscoelastic behavior: Application to the generalized maxwell model[END_REF]).

In general cases, several viscoelastic materials can appear in a same FE model. For some reasons concerning calculation 6

Copyright c 2011 by ASME time, that are not developed here, it is better to have the same relaxation times for all materials. A relaxation time is equal to the inverse of a pole. So, in the case of 3 pole-zero models with 5 pole-zero couple, it means it is better to have 5 poles common to the 3 models. The inclusion of this constraint requires to do 2 more steps. The third one consists in calculating the mean value of the poles and recalculating the zeros. Then, the fourth step is an optimisation step analogous to the second one. The Figure 9 gives the final result of the identification process in the case of 3 test functions. The identification is still good.

CONCLUSIONS

A new Dynamical Mechanical Analysis (DMA) tester has been developed for the characterization of the shear stiffness of viscoelastic materials which are known to be frequency dependent. It allows direct characterization under preload, for frequencies starting from 200Hz to 3500Hz without using the time-temperature equivalences.

The Generalized Maxwell model has been chosen to describe the experimentally measured behavior of viscoelastic material. This rheological model is equivalent to a pole-zero formulation. An identification method has been developed. It allows to compute the poles and the zeros from the behavior of materials measured with the DMA tester. This method gives good results and allows ensuring the causality of the model by constraining the zeros to be lower than the poles. In the case of several materials to identify, this identification method allows computing same poles for all material. Since a pole is the inverse of a relaxation time, each material is thus modeled with the same relaxation times. This characteristic is of great importance for computation time reduction in FE models.

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).

2

 2

FIGURE 2 .

 2 FIGURE 2. Principle and picture of the shear DMA tester.

FIGURE 3 .

 3 FIGURE 3. Real working characteristics of the piezoelectric actuators combined with a 30W driver.

3

  Copyright c 2011 by ASME alimentation: δV = 100 Volts.

FIGURE 4 .

 4 FIGURE 4. Model of the DMA tester for post processing.

FIGURE 5 .

 5 FIGURE 5. Modulus and phase of the shear stiffness of a multilayer viscoelastic material with layers of glue, rubber and steel. The blue curve has been obtained with a normal load of 5000N and the red one with a normal load of 500N.

FIGURE 6 .

 6 FIGURE 6. Rheological model of Generalized Maxwell.

FIGURE 7 .

 7 FIGURE 7. Behaviors of the modulus and the phase of one pole-zero couple and its asymptotes.

FIGURE 8 .

 8 FIGURE 8. Given a test stiffness function, the identification is led in two steps: initialization and optimization.

FIGURE 9 .

 9 FIGURE 9. Identification of 3 viscoelastic materials in the same time, with the constraint of common poles


Marco Rosatello
email: marco.rosatello@gmail.com

Jean-Luc Dion

Franck Renaud

Luigi Garibaldi

The skateboard speed wobble

Keywords: wobble, skateboard, stability

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The skateboard is a recreational vehicle whose origin extend back to the early 20's, when it consisted of a rudimentary crate scooter made of a wooden board, metal wheels and a wooden box. In the 50's the first commercially produced skateboards appeared on the market. The very turning point in the skateboard history took place in the 70's, thanks to the use of polyurethane wheels, which allowed more controllability and traction properties than the previous metal wheels. They also allowed higher speeds which brought to the emergence of a new type of skateboard, called longboard, specially designed for downhill skateboarding. They are longer than normal skateboards, they have different shapes and very high speed can be reached. The fastest skateboard speed achieved in a standing position is 129.94 km/h (80.74 mph) and was achieved by Mischo Erban at Les boulements, Quebec, Canada, on 18 June 2012 on a road descending at an 18% grade. To achieve high speeds the skateboard has to be set up correctly and the rider has to adapt his body to control the vehicle. If one of the this two aspects fails, at a certain speed the skateboardrider system can become unstable, showing an oscillatory behavior that can lead to loss of control and very dangerous falls. This phenomenon is known as speed wobble.

Literature Review

The existing literature about skateboard is limited and the speed wobble has never been treated. In [START_REF] Hubbard | Lateral Dynamics and Stability of the Skateboard[END_REF] Hubbard studied the lateral dynamics stability of a 1-DOF and a 2-DOF linearized skateboard model, but the speed wobble is not treated. Again Hubbard, in [START_REF] Hubbard | Human Control of the Skateboard[END_REF], after simplifying hypotheses and reducing the 2-DOF model to a 1-DOF model, implemented an optimal control to the skateboard in order to simulate the human control. However, even if the perfect control has been experimentally verified to be similar to the actual human control for low speeds, it's not interesting for higher speed ranges. Another major work in the skateboard field is [START_REF] Kremnev | Nonlinear Dynamics and Stability of a Simplified Skateboard Model[END_REF], in which Kremnev and Kuleshov, on the basis of [START_REF] Hubbard | Lateral Dynamics and Stability of the Skateboard[END_REF], study the stability of a more accurate 1 DOF skateboard model, obtaining the equations of motion through the Gibbs-Appell method. They offer a mathematically accurate nonlinear stability analysis, demonstrating a rattleback behavior for asymmetric skateboards around the equilibrium position. However, the speed wobble is never mentioned. The objective of this paper is to study the skateboard stability in a speed-wobble-focused way, analyzing the most affecting parameters, and to verify or disprove common skaters beliefs about speed wobble. In order to do that, a parametric stability analysis will be carried out on a 1 DOF skateboard model, with the intention of determining the parameters that has more influence on stability at higher speeds. Furthermore, a human control model based on Biomechanics researches on the human postural control (see [START_REF] Winter | Biomechanics and Motor Control of Human Movement[END_REF]- [START_REF] Chagdes | Dynamic stability of a human standing on a balance board[END_REF]) will be implemented on a second model of skateboard.

THE SKATEBOARD

The skateboard typically consists of three major components: a board, a set of two trucks and four wheels. A skateboard model is shown in Fig. 1. Modern boards are generally from 50 to 110 cm long, 15 to 30 cm wide and 1 to 2 cm thick. The stiffness of the board depends on the particular application desired: when a great manoeuvrability is required, a more flexible board is used, while stability at high speeds calls for more stiffness. In this study the board stiffness won't be considered in the equations of motion. The skateboard wheels are mounted on the axles on ball bearings and their stiffness is measured with the durometer A scale. The range of commercial wheels hardness is between 75A and 97A. Concerning the trucks, they are the elements that allow the skateboard to turn. The most common model is called Randal truck and it is characterized by three parameters: the wheelbase a, the pivot angle λ and the bushings. The wheelbase a is simply the distance between the front and rear truck. The pivot angle λ is the angle between the baseplate and the kingpin. The front truck pivot angle is λ f , while the rear one λ r , as shown in Fig. 1. Normally the mounted trucks are the same, thus λ f = λ r , but sometimes different configurations are used. The bushings exist in different shapes and hardness and, together with the tightening of the kingpin bolt, they determine the overall torsional spring stiffness k γ , shown in Fig. 3. The latter applies a restoring torque T γ between the wheel-set and the board, proportional to the tilt angle γ. In [START_REF] Hubbard | Lateral Dynamics and Stability of the Skateboard[END_REF][START_REF] Hubbard | Human Control of the Skateboard[END_REF][START_REF] Kremnev | Nonlinear Dynamics and Stability of a Simplified Skateboard Model[END_REF][START_REF] Ispolov | [END_REF][START_REF] Osterling | On the Skateboard, Kinematics and Dynamics[END_REF], the effect of different trucks tightening is not taken into account. In this study we will consider the possibility that the front and rear truck torsional spring levels of stiffness k γ F and k γ R are different, which is a very common set up among skaters. We can write the overall restoring torque as

T γ = k γ γ = k γ f γ f + k γ r γ r ( 1
)
where γ f andγ r are the board tilt angles at the points F and R shown in Fig. 1. Assuming that

k γ f γ f = k γ r γ r (2)
and by replacing (2) in (1),we can find

γ f = 1 + µ 2µ γ γ r = 1 + µ 2 γ (3)
where µ = k γ f /k γ r .

Steering system

The skateboard is a negative-four wheel steering vehicle (for more details see [START_REF] Jazar | Vehicle Dynamics. Theory and Application[END_REF]): the angular motion of both the front and rear axles is constrained to be about their respective pivot axes, thus causing steering angles δ f and δ r , shown in Fig. 2, whenever the axles are not parallel to the plane of the board. The relationship between the board tilt angle γ and the steering angles δ f and δ r has been demonstrated in [START_REF] Kremnev | Nonlinear Dynamics and Stability of a Simplified Skateboard Model[END_REF] to be

tan δ f = tan λ f sin γ tan δ r = tan λ r sin γ (4)
Once the trucks are mounted, the pivot angles λ f and λ r are fixed. Furthermore taking into account (3), (4) becomes

tan δ f = c f sin 1 + µ 2µ γ tan δ r = c r sin 1 + µ 2 γ ( 5
)
where c f = tan λ f and c r = tan λ r . The last step to fully define the skateboard steering system is to assume that its wheels roll without lateral sliding. Looking at Fig. 2, this condition is satisfied through two nonholonomic constraints defined in [START_REF] Kremnev | Nonlinear Dynamics and Stability of a Simplified Skateboard Model[END_REF] as:

ẋo sin(θ -δ f ) -ẏo cos(θ -δ f ) - a 2 θ cos δ f = 0 ẋo sin(θ + δ r ) -ẏo cos(θ + δ r ) + a 2 θ cos δ r = 0 (6)

EQUATIONS OF MOTION

The skateboard equations of motion are obtained for two different models with different assumptions. In Model 1, outlined in Fig. 3, the rider is assumed to remain fixed and perpendicular to the board, so that the angle between the z-axis and the rider COM is the same board tilt angle γ. In Model 2, outlined in Fig. 4, the rider is assumed to move freely about the x-axis, thus forming an angle φ between the rider COM and the vertical z-axis. The rider is assumed to control its position through a torque exerted at the ankles level about the axis passing from the center of the board. In both models, when γ = 0, the projection of the rider COM is assumed to lie on the board centerline FR at a distance p from the board center B,as shown in Fig. 2. If p > 0 the rider COM is shifted in the front of the vehicle, while if p < 0 in the back. The equations of motion are obtained through Lagrange formalism with Lagrange multipliers to account for the nonholonomic constraints [START_REF] Jazar | Vehicle Dynamics. Theory and Application[END_REF]. Choosing the vector of N Lagrangian generalized coordinates q = (q 1 , q 2 , ..., q N ) T , the Lagrange energy method allows to write the equations of motion from the kinetic energy T and potential energy U expressions for the system under investigation. Being N C the number of constraints and

N ∑ i=1 c ji (q) • qi = 0 , j = 1, ..., N C (7)
the resulting first-order differential equations in Pfaffian form, following [START_REF] Janschek | Mechatronics Systems Design. Methods, Models, Concepts[END_REF], the equations of motion with N C Lagrange multipliers λ can be written as

d dt ∂ L(q, q) ∂ qi - ∂ L(q, q) ∂ q i = Q i + N C ∑ j=1 λ j c ji (q) i = 1, .., N (8)
where L(q, q) = T (q, q) -U(q) is called Lagrangian,Q i contains the generalized forcing terms f i and the dissipative forces terms d i .

The second-order differential system of equations ( 8) and the first-order differential system of equations ( 7) give a total of N + N C equations for computation of the N + N C unknowns q 1 , q 2 , ..., q N and λ 1 , ..., λ N C . Here the Lagrange multipliers represent the constraint forces required to maintain the constraints in Eq.( 7).

The final first-order differential system of equations is:

  I 0 0 0 A C T 0 C 0     q q -λ   =   q Q( q, q,t) + B( q, q,t) D( q, q,t)   (9)
where

D( q, q,t) = C q (10)
on which the Baumgarte stabilization method have been applied.

FIGURE 3: FRONT VIEW OF MODEL 1

Model 1 -Rider Constrained

The generalized coordinate vector for Model 1 is:

q = [x O , y O , θ , γ],
where x O and y O are the coordinates of point O, shown in Fig. 2,3. Considering a skateboard with a mass m B , moments of inertia I B x , I B z and its COM in point B and a rider with a mass m R , moments of inertia I R x , I R z and its COM in point R, the expression of the kinetic energy is: [START_REF] Herman | Physics of the Human Body[END_REF] where I x = I B x + I R x and I z = I B z + I R z , since the rotational speeds γ and θ are the same for the board and the rider. The potential energy is:

T = 1 2 m B ( ẋ2 B + ẏ2 B + ż2 B ) + 1 2 m R ( ẋ2 R + ẏ2 R + ż2 R ) + 1 2 I x γ2 + 1 2 I z θ 2
U = 1 2 k γ γ 2 + m B gz B + m R gz R ( 12
)
Concerning the generalized forces, there are no active forces or torques acting on the skateboard, thus the generalized forces vector is

Q 1 = [0, 0, 0, 0] T .
In order to write Eq.( 11) and ( 12) as function of the generalized coordinates, looking at Figs. 23, the following geometric constraints were found:

x B = x O + h sin γ sin θ , x R = x O + (h + l) sin γ sin θ + p cos θ y B = y O -h sin γ cos θ , y R = y O -(h + l) sin γ cos θ + p sin θ z B = h cos γ , z R = (h + l) cos γ (13)
and deriving them in respect of time we obtain

ẋB = ẋO + h( γ cos γ sin θ + θ sin γ cos θ ) ẋR = ẋO + (h + l)( γ cos γ sin θ + θ sin γ cos θ ) -p θ sin θ ẏB = ẏO -h( γ cos γ cos θ + θ sin γ sin θ ) ẏR = ẏO -(h + l)( γ cos γ cos θ + θ sin γ sin θ ) + p θ cos θ żB = -h γ sin γ , żR = -(h + l) γ sin γ (14)
Then we replace ( 14) and ( 13) in ( 11) and ( 12) At this point the first-order differential system (9) was derived with the Symbolic Math Toolbox of MATLAB R . The resulting matrices are omitted for brevity.

λ f 45 • λ r 45 • h 0.07 m p 0 m I B x 0.025 kgm 2 I B x 12 kgm 2 m f 3 kg I f x 0.45 kgm 2 k γ 250 Nm rad -1 µ 1

Model 2

The same procedure is followed for Model 2. The outline of the model is shown in Fig. 4. Concerning the human control, the latter is exerted through an ankles torque T B , assumed to act around the board x-axis. The generalized coordinate vector is q

= [x O , y O , θ , γ, φ ].
Since in this model the rider is assumed to move independently of the board, some considerations have to be made on the inertia terms. In particular it's important to pay attention to the feet inertia. In fact, even if the skater body moves independently of the board, the feet are fixed to the board and their rotation is equal to the board tilt angle γ and not the skater COM angle φ . Therefore, the feet inertia around the x-axis I f x has to be considered with the board inertia and not with the rider inertia. In the same way also the mass of the feet has to be considered with the board so that:

m ′ B = m B + m f , m ′ R = m R -m f ( 15
)
While the new inertia terms becomes then:

I ′ B x = I B x + I f x , I ′ R x = I R x -I f x ( 16
)
The expression of the kinetic energy is:

T = 1 2 m B ( ẋ2 B + ẏ2 B + ż2 B ) + 1 2 m R ( ẋ2 R + ẏ2 R + ż2 R ) + 1 2 I ′ B x γ2 + 1 2 I ′ R x φ 2 + 1 2 I z θ 2 (17)
while the potential energy is:

U = 1 2 k γ γ 2 + 1 2 k φ (γ -φ ) 2 + m B gz B + m R gz R (18)
Concerning the generalized forces, the only active force is the ankle torque exerted by the skater, so that the vector of generalized forces becomes

Q 2 = [0, 0, 0, -T B , T B ] T .
Concerning the geometric constraints, for the point B they remain the same as for Model 1, while for point R they become:

x R = x O + (h sin γ + l sin φ ) sin θ + p cos θ y R = y O -(h sin γ + l sin φ ) cos θ + p sin θ z R = h cos γ + l cos φ (19)
while their derivatives become:

ẋR = ẋO + (h γ cos γ + l φ cos φ -p θ ) sin θ + θ (h sin γ + l sin φ ) cos θ ẏR = ẏO -(h γ cos γ + l φ cos φ -p θ ) cos θ + θ (h sin γ + l sin φ ) sin θ żR = -h γ sin γ -l φ sin φ (20)

Linearized Equations

To linearize the equations of motion several assumptions are made. The longitudinal speed ẋO is assumed to have a constant value u so that ẍO = 0. From [START_REF] Kremnev | Nonlinear Dynamics and Stability of a Simplified Skateboard Model[END_REF]

is demonstrated that θ = - u sin(δ f + δ r ) a cos δ f cos δ r (21)
and that, solving the two nonholonomic constraints ( 6) is possible to find

ẏO = - a θ 2 sin(δ f + δ r ) [cos δ f sin(θ + δ r ) + cos δ r sin(θ -δ f )] (22)
Furthermore, small angles assumption is made for the angles γ,θ and φ , so that cos

γ = cos θ = cos φ ≃ 1, sin γ ≃ γ, sin θ ≃ θ , sin φ ≃ φ .
Taking into account the assumption made and replacing (21),( 22) and their time derivatives in system (9), we obtain the linearized equation for Model 1:

I x + m R L 2 + m B h 2 γ + m R L ε + η p a + m B hε u γ + k γ -m R gL -m B gh + ηu 2 a (m R L + m B h) γ = 0 ( 23
)
where L = h + l,while

ε = µ + 1 4µ (c f -c r µ) , η = µ + 1 2µ (c f + c r µ) (24)
are the coefficients that take into account the eventuality of different trucks pivot angles and degrees of stiffness. They act as damping parameters but they don't involve any viscoelastic behavior of the trucks bushings.

The linearized system of equations for Model 2 is

I ′ B x + mh 2 γ + m ′ R hl φ + mε + m ′ B η p a hu γ + k γ -mgh + mh ηu 2 a γ = -T B I ′ R x + m r l 2 φ + m ′ R hl γ + m ′ R l ε + η p a u γ -m ′ R glφ + m ′ R lηu 2 a γ = T B ( 25
)
where m = m R + m B .

Human Control of the Skateboard

In order to design a human control for the skateboard, researches were carried out on the Biomechanics topic of human postural control, i.e. the ability of maintaining the balance. The latter is based on the pieces of information coming from three systems: visual, vestibular and somatosensory (see [START_REF] Mergner | Multisensory Control of Posture[END_REF]). The visual system provides information from the eyes in the absolute frame of reference. The vestibular system provides information about head position and acceleration in the absolute frame of reference through the inner ear vestibulum. The somatosensory system provides information about the position of the center of pressure under the feet. In [START_REF] Peterka | Sensorimotor Integration in Human postural Control[END_REF] David Peterka studied the postural control on subjects standing on an oscillating platform. He experimentally demonstrated that the postural control can be modeled through a sensory reweighting system. This means that the relevance of every sensory system changes with the external perturbation, accounting in this way for the nonlinearity of the human control. For low disturbance amplitudes the balance task relies more on the somatosensory system, while for higher disturbances the relevance of the vestibular system becomes greater. The control of the human balance consists of a torque exerted at the ankles level, which is made of two contributions: an Active Torque, based on the neuromuscular control, and a Passive Torque, result of the leg muscles intrinsic features. The control model designed for the skateboard is shown in Fig. 5. The core of the model is the skateboard-rider mechanical system, which contains the equations of motion (25) and takes as input the ankles torque T B and gives as output the angle φ and γ. The outputs are elaborated by the sensory systems: the visual and vestibular systems act on the same angle φ , which is the skater COM angle in the absolute frame of reference, while the somatosensory system acts on the skater angle relative to the board φγ. The red boxes contains the weighting factor W v and W s , representing the sensory reweighting system, such that W v +W s = 1. The error from the equilibrium position elaborated by the sensory reweighting system enters then in a PD controller, with K p as the proportional term and K d as the derivative term, and determines the active torque T a . The latter is summed to the passive torque T p , which depends on the muscular ankles stiffness K and damping C, in order to give the total ankles torque T B = T a + T p . The total ankle torque can then be written as:

T B = K p (φ -W s γ) + K d ( φ -W s γ) + K(φ -γ) +C( φ -γ) (26)
The values of the PD controller terms were taken from Peterka's paper [START_REF] Peterka | Sensorimotor Integration in Human postural Control[END_REF] as K p = 1000 Nmrad -1 and K d = 100 Nmsrad -1 . Concerning the passive torque, the value of the ankle passive stiffness was determined by Peterka to be K = 94 Nmrad -1 . Peterka also determined the value of C. However, its entity suggested him that the passive damping C may contribute more to overall damping when the surface is fixed than when it is moving. Therefore, considering the mobility of the skateboard, we will use the value found by Maura Casadio in her study [START_REF] Casadio | Direct measurement of ankle stiffness during quiet standing: implications for control modelling and clinical application Gait and Posture[END_REF] on the ankles impedance. The value is then assumed to be c = 5.8 Nmsrad -1 .

STABILITY ANALYSIS

A parametric stability analysis with the poles location method was carried out to investigate the influence of different skateboard setups. The Laplace transformation is applied to the linearized equations ( 23) and ( 25) to obtain the system charac-teristic polynomial, whose solutions determine its poles, thus its stability. All the figures show the poles going from a blue color to a red color as the chosen parameter value increases. The data of the chosen skateboard -rider configuration are shown in Table 1 and they represent an average skateboard set-up and an average human body in a standing position and were taken from previous studies and anthropometry sources (see [START_REF] Kremnev | Nonlinear Dynamics and Stability of a Simplified Skateboard Model[END_REF][START_REF] Winter | Biomechanics and Motor Control of Human Movement[END_REF][START_REF] Herman | Physics of the Human Body[END_REF]).

Model 1 -Results

The stability of the skateboard-rider model is studied by varying different parameters, while maintaining fixed all the others. The first considered parameter is the longitudinal speed u.

Figure 7 shows the eigenvalues position for a symmetric skateboard-rider configuration (λ f = λ r , µ = 1 , p = 0 m): at low speeds the system has a couple of real symmetric eigenvalues, thus it's unstable. With higher speeds the system has a couple of purely imaginary eigenvalues and shows a marginally stable oscillatory behavior, with increasing frequency. The latter aspect is due to the fact that the centripetal restoring torque, represented by the γ-term of equation ( 23), increases with the square of the speed, increasing the overall stiffness of the system.Thus the speed has a stabilizing effect,as it happens in bikes and motorbikes (see [START_REF] Ringwood | Bicycle Wheel Wobble International Conference on Informatics in Control[END_REF]). The characteristic equation of the model is of the form A 0 s 2 + B 0 s + C 0 = 0 and the critical speed, above which the motion becomes stable, corresponds to the value of u for which its discriminant ∆ = B 2 0 -4A 0 C 0 = 0 In Fig. 6 are shown the changes in critical speed due to the spring torsional stiffness k γ , the wheelbase a and the pivot angles λ . The critical speed u cr becomes lower as k γ become higher, until u cr = 0. This happens when, C 0 = 0 for u = 0, thus when k γ cr = g(m R L + m B h). Over this value, the skateboard is If we consider that the front and rear trucks are different we can fix the value of one of the two and let the other vary. This situation is represented in Fig. 8: it shows that if λ f > λ r the system is stable while if λ f < λ r the system becomes unstable. However, using different pivot angles is not very widespread among skaters.

A simpler way of acting on the trucks is to adjust the kingpin bolt tightening level. The latter aspect was modeled through the coefficient µ = k γ f /k γ r . Fig. 9 shows how the stability is affected by this parameter: the skateboard-rider system is stable for µ < 1, while it becomes unstable for µ > 1.

As far as the position of the rider on the board p is concerned, the latter is shown to be a key parameter for the skateboard stability.

Figure 10 shows that, if the rider shifts his weight in the front of the vehicle, the system is stable, otherwise the system becomes unstable. Except for these three parameters which can lead to an asymmetric configuration, the other skateboard and rider parameters, such as the wheelbase a, the trucks overall stiffness k γ and the COM height on the board l have no effects on the system stability. In fact their effect is only to increase or decrease the oscillation frequency, as it happens for the speed in Fig. 7. Further investigation of these parameters will be offered in the human controlled skateboard of Model 2.

The results obtained from the stability study of the linearized system were then verified by solving the nonlinear full system (9) by imposing a small initial board tilt angle perturbation γ 0 and a an initial longitudinal speed u 0 . The results are shown in Fig. 11: Fig. 11a shows the system time response of a marginally stable skateboard-rider configuration, while Fig. 11b shows the system response for an unstable skateboard-rider configuration at the same speed u 0 = 18 km/h.

Model 2 -Results

To analyze the system of equations (25), the ankles torque shown in Eq.( 26) was replaced in the system and the Laplace transform was applied, obtaining:

A 1 s 2 + B 1 s +C 1 A 2 s 2 + B 2 s +C 2 A 2 s 2 + B 3 s +C 3 A 3 s 2 + B 4 s +C 4 L γ φ = 0 0 (27)
The characteristic equation is then a quartic equation. The main task is to analyze the stability of the human control and in particular of the sensory reweighting system. In [START_REF] Peterka | Sensorimotor Integration in Human postural Control[END_REF], David Peterka demonstrated that that, for low perturbations of his motion board (0.5 • ), the sensory weight distribution is the following: 50% somatosensory and 50% visual+vestibular. Instead, for larger amplitude perturbations (8 • ), the sensory weight distribution is 15% somatosensory and 85% visual+vestibular. Thus, considering this aspect, we analyzed several weighting factors situations. Figure 12 shows two different situations of weighting factors distribution. Figure 12a shows that, for a somatosensory weighting factor W s = 0.5, the human control is stable at every speed. Instead, Fig. 12b shows that, for W s = 0.15, at a certain point, the real part of a couple of eigenvalue becomes positive and the system become unstable. The value of the speed for which the latter happens is called critical wobble speed u cr W . Figure 13 shows the dependency between the somatosensory weighting factor W s and the critical wobble speed. It can be seen that, under a certain value of W s , the speed wobble becomes possible and the speed at which it happens decreases with W s . In particular it can be said that, if the weight of the somatosensory system is larger than 30%, the speed wobble instability can't take place.

At this point, assuming a symmetric skateboard-rider configuration, the parameters whose influence on the system stability of Model 1 was uncertain will be analyzed in this context. when the pivot angle λ f = λ r or the feet inertia I f x increases. The latter parameter has been considered because the position of the feet on the board can be changed: in fact, if we consider a professional skater position on the board at high speed, the feet are positioned longitudinally on the board center line. In this way the feet inertia around the x-axis I f x is reduced. The effect of the different parameters on the critical wobble speed are important because they show that, the same parameters whose influence on Model 1 was not clear, are instead important from a human postural control point of view. The results obtained from the stability analysis of the human control of the skateboard were then verified on the nonlinear full system by imposing an initial perturbation of the board tilt angle γ 0 and an initial speed u 0 . Figure 15 shows the time response for two different configurations with different wobble critical speed.

It can be seen that, higher the critical wobble speed, higher the oscillation frequency of the wobble instability. The frequency range of the instability is between 4 and 9 Hz.

DISCUSSION

The stability analysis on the two skateboard-rider models show different results. Model 1 showed that the there is a critical speed above which the motion is marginally stable if the skateboard-rider configuration is symmetric, i.e. when λ f = λ r , µ = 1, p = 0. If the configuration is not symmetric, the sys-tem can become stable or unstable. By varying one parameter at a time has been discovered that the skateboard is stable over the critical speed when λ f > λ r , i.e. the front truck pivot angle is greater than the rear one; when µ < 1, i.e. the rear truck is more tighten then the front one; and when p > 0, i.e. the rider COM is shifted more in the front of the vehicle. In the opposite cases, instead,the skateboard becomes unstable. Concerning the wobble, we were expecting an instability arising with speed increase, but in the obtained results there is no trace of such behavior. Furthermore, considering that the skateboard parameters k γ , λ and µ can't be changed during the motion, the skateboard will always be either in a stable or unstable configuration.The only aspects that can be changed during the motion are the position p and the height l of the rider COM. The most important parameter between the two is the position p. Looking at professional long-boarders, they try to shift their weight in the front as much as possible and their COM as low as possible. In fact, as demonstrated by a few wobble videos, if the weight is shifted in the back even for a short moment, the wobble occurs and in most cases the rider loses the control of the vehicle. However, if the skateboardrider is in an unstable configurations, this doesn't mean that the rider can't control the vehicle, it will only be more difficult as the speed increases. For these reasons a human control was applied on the skateboard in Model 2. The results demonstrate that the speed wobble is due to the sensory reweighting system of the human postural control, and in particular to the loss of relevance of the somatosensory system. In fact, under a certain value of the somatosensory weigthing factor, it's possible to have a critical wobble speed, above which the human control becomes unstable. However, the evolution of the weighting factors can be very different from person to person and even within the same person, and it depends on many aspects, such as experience, training, concentration, intrinsic abilities to maintain the balance and so on. For example, a professional skater will be able to maintain the somatosensory weighting factor higher than a beginner skater for the same external perturbations. Therefore, the speed wobble effect is due in part to the skateboard setup, but in a major part to the rider skills in maintaining its balance. Furthermore, the results of the parametric stability analysis have pointed out that the human postural control is influenced by several skateboard parameters such as the overall torsional truck stiffness k γ , the wheelbase a and the trucks pivot angle λ . In order to increase the wobble critical speed, the skater should use longer wheelbase skateboards, tighten more the trucks and use trucks with low pivot angles. Moreover, also the feet inertia has been evaluated as a very influential parameter, and the rider should put them in a longitudinal position instead of perpendicularly. However, by doing this, the ability to turn is lowered and the available base area in which the COM projection has to fall in order to maintain the equilibrium is lower. Therefore the rider should also lower its COM on the board.

CONCLUSION

The equations of motion of two models of skateboard were obtained through Lagrange energy method and numerically solved in MATLAB R . A human postural control was implemented on the second model and the stability of the two systems evaluated through a parametric analysis in order to investigate the so called speed wobble instability. Model 1 showed that there are stable and unstable skateboard-rider configurations, but they don't show a clear dependency on the speed. Model 2 demonstrates that the speed wobble instability depends on the sensory reweighting system of the human postural control and that the value of the critical wobble speed varies according to skateboard and rider parameters. The results seem to confirm skaters beliefs and opinions about the speed wobble. Therefore, even if it was demonstrated that the speed wobble depends on the rider skills, the factors that influence the instability were discovered and theoretically demonstrated. In order to validate the theoretical results obtained in this paper, a series of experiments will be needed.

NOMENCLATURE a

  Skateboard wheelbase COM Center of Mass µ Ratio between front and rear truck torsional stiffness γ Board tilt h Board height from axles plane I B x Board moment of inertia about x-axis I R x Rider moment of inertia about x-axis I B z Board moment of inertia about z-axis I R z Rider moment of inertia about z-axis λ f Pivot angle front truck λ r Pivot angle rear truck k γ Overall skateboard truck torsional stiffness k γ f Front truck torsional stiffness k γ r Rear truck torsional stiffness l Rider COM height from board m B Skateboard mass m R Rider mass p Distance from the board center of the COM projection T B Total ankles torque W s Somatosensory weighting factor W v Visual+Vestibular weighting factor u Longitudinal skateboard speed

FIGURE 1 :

 1 FIGURE 1: LATERAL AND TOP VIEW OF A SKATEBOARD

FIGURE 2 :

 2 FIGURE 2: SKATEBOARD COORDINATE SYSTEM

FIGURE 4 :

 4 FIGURE 4: FRONT VIEW OF MODEL 2

FIGURE 5 :

 5 FIGURE 5: HUMAN CONTROL MODEL OF THE SKATEBOARD

FIGURE 7 :FIGURE 8 :FIGURE 9 :FIGURE 10 :

 78910 FIGURE 7: POLES LOCATION FOR SPEED 0 < u < 72km/h

FIGURE 11 :

 11 FIGURE 11: NONLINEAR SYSTEM RESPONSE (γ 0 = π/70,u 0 = 18Km/h)

FIGURE 12 :

 12 FIGURE 12: HUMAN CONTROL STABILITY FOR DIFFERENT SENSORY WEIGHTING FACTORS 0 < u < 72km/h

FIGURE 13 : 2 (FIGURE 14 :

 13214 FIGURE 13: WOBBLE CRITICAL SPEED DEPENDENCY ON SOMATOSENSORY WEIGHTING FACTOR

FIGURE 15 :

 15 FIGURE 15: NONLINEAR SYSTEM RESPONSE MODEL 2 γ 0 = π/70

TABLE 1 :

 1 AVERAGE SKATEBOARD-RIDER DATA

		Skateboard -Rider Parameters
	m B	3	kg	m R	75	kg
	a	0.6	m	l	0.7	m

:

  STIFFNESS k γ , PIVOT ANGLE λ AND WHEELBASE a

		8
	[km/h]	4 6
	u cr	2
		0	100	200	300	400	500	600	700	800	900	1000
						k γ	[Nm/rad]
	10
	[km/h]	5
	u cr
		20 0	25	30	35	40	45	50	55	60	65	70
							λ [m]
		8
	[km/h]	6
	u cr	4
		0.3 2	0.4		0.5	0.6	0.7	0.8	0.9		1
							a [m]
	FIGURE 6INFLUENCE ON CRITICAL SPEED