Quentin Leclère
email: quentin.leclere@insa-lyon.fr.

Hugo André

Jérôme Antoni

A multi-order probabilistic approach for Instantaneous Angular Speed tracking debriefing of the CMMNO'14 diagnosis contest

Keywords: Instantaneous Angular Speed, wind turbines, order tracking, demodulation, virtual tacho

The aim of this work is to propose a novel approach for the estimation of the Instantaneous Angular Speed (IAS) of rotating machines from vibration measurements. This work is originated from the organisation, by the authors of this paper, of a contest during the conference CMMNO 2014, that was held in Lyon, december 2014. One purpose of the contest was to extract the IAS of a wind turbine from a gearbox accelerometer signal. The analysis of contestant contributions led to the observation that the main source of error in this exercise was the wrong association of one selected and tracked harmonic component with one mechanical periodic phenomenon, this association being assumed as an a priori hypothesis by all the methods used by the contestants. The approach proposed in this work does not need this kind of a priori assumption. A majority (but not necessarily all) periodical mechanical events are considered from a preliminary analysis of the kinematics of the machine (harmonics of shaft rotation speeds, meshing frequencies ...). The IAS is then determined from probability density functions that are constructed from instantaneous spectra of the signal. The efficiency and robustness of the proposed approach are illustrated in the frame of the CMMNO 2014 contest case.

Nomenclature Ω: IAS (Instantaneous Angular Speed) (in Hz) ω: argument of the probability density function (pdf) of the IAS (in Hz) A(f ): whitened instantaneous spectrum (dimensionless), depending on the frequency f (in Hz) Ω min : minimum value of the a priori IAS uniform law (in Hz) Ω max : maximum value of the a priori IAS uniform law (in Hz) H i (i ∈ [1...n]): orders of the mechanism, with respect to the IAS (dimensionless) [Ω|H i ]: pdf of the IAS knowing H i , function of ω (in Hz) [Ω j |Ω j+k ]: pdf of the IAS at time j knowing the IAS at time j + k, function of ω (in Hz) [Ω j ] j+k : pdf of the IAS at time j obtained from [Ω j+k ], function of ω (in Hz) γ: standard acceleration for IAS (smoothing parameter), in (Hz 2 ) [Ω j ] s : smoothed pdf of the IAS at time j, function of ω (in Hz)

Introduction

In many aspects, the Instantaneous Angular Speed (IAS) is to the Condition Monitoring of rotating machines what the heartbeat signal is to medicine. Just as in medicine, the IAS provides vital information for diagnosis and it allows resynchronization of subsequent analyses. The importance of the IAS and of its measurement has been addressed in a recent special issue of MSSP [START_REF] Rémond | Editorial for the special issue on instantaneous angular speed (ias) processing and angular applications[END_REF][START_REF] Rémond | Instantaneous angular speed (ias) processing and related angular applications[END_REF] and in a review paper [START_REF] Lin | A review and strategy for the diagnosis of speedvarying machinery[END_REF]. Its use for order tracking is also discussed in the review paper [START_REF] Brandt | Main principles and limitations of current order tracking methods[END_REF]. As indicated by the adjective instantaneous, the IAS is a measure of the rotation speed of a rotating component in a machine with an angular resolution that corresponds to at least one value per revolution. The IAS enters as an essential ingredient to several processing techniques dedicated to the analysis of machines operating under nonstationary regimes. It has been for long an essential quantity in rotor dynamics and in NVH applications concerned with the testing of rotating machines in transient regimes such as runups and shutdowns. In this context, the IAS makes possible the synchronization of Fourier analysis on the shaft rotation -a practice known as order (spectral) analysis -or the tracking of the evolution of the Fourier coefficients as a function of speed -a practice known as order tracking [START_REF] Wang | The combined use of order tracking techniques for enhanced fourier analysis of order components[END_REF][START_REF] Wang | The application of order tracking for vibration analysis of a varying speed rotor with a propagating transverse crack[END_REF][START_REF] Borghesani | The velocity synchronous discrete fourier transform for order tracking in the field of rotating machinery[END_REF]. Early successes of order tracking are exemplified by the Vold-Kalman filter [START_REF] Vold | High resolution order tracking at extreme slew rates, using kalman tracking filters[END_REF][START_REF] Vold | Theoretical foundations for high performance order tracking with the vold-kalman tracking filter[END_REF][START_REF] Feldbauer | Realization of a vold-kalman tracking filter -a least squares problem[END_REF][START_REF] Bai | Fault diagnosis of rotating machinery using an intelligent order tracking system[END_REF][START_REF] Pan | Further exploration of vold-kalman-filtering order tracking with shaft-speed information-i: Theoretical part, numerical implementation and parameter investigations[END_REF][START_REF] Pan | Adaptive vold-kalman filtering order tracking[END_REF][START_REF] Bonnardot | Use of the acceleration signal of a gearbox in order to perform angular resampling (with limited speed fluctuation)[END_REF], a model-based approach well suited to slow-speed variations and high SNRs. It is traditionally accepted that order analysis and order tracking apply to deterministic signals which exhibit (modulated) harmonic structures in the order domain. It is only recently that these ideas have been extended to broadband signals. A typical example is provided by the vibrations of bearings under nonstationary regimes, where order analysis has been successfully applied to the signal envelope or to time-frequency energy distributions [START_REF] Urbanek | Time-frequency approach to extraction of selected second-order cyclostationary vibration components for varying operational conditions[END_REF][START_REF] Borghesani | A new procedure for using envelope analysis for rolling element bearing diagnostics in variable operating conditions[END_REF][START_REF] Cheng | Envelope deformation in computed order tracking and error in order analysis[END_REF]. The IAS also turns out essential to resample machine signals in the angular domain in all applications where synchronization to the machine kinematics is essential to correctly identify and analyze transient events in the machine cycle [START_REF] Antoni | Effective vibration analysis of ic engines using cyclostationarity. part i-a methodology for condition monitoring[END_REF][START_REF] Villa | Angular resampling for vibration analysis in wind turbines under non-linear speed fluctuation[END_REF][START_REF] Cocconcelli | An algorithm to diagnose ball bearing faults in servomotors running arbitrary motion profiles[END_REF]. The advantage of synchronizing signals to the angular domain has been fully highlighted within the cyclostationary framework, where the notion of machine cycle is central [START_REF] Antoni | Cyclostationary modelling of rotating machine vibration signals[END_REF][START_REF] Antoni | Cyclostationarity by examples[END_REF]. The cyclostationary framework happens to perfectly model most signals produced by rotating machines operating under stationary regimes. The latter condition has recently been extended to nonstationary regimes, thus generalizing the cyclostationary class of signals to the cyclononstationary class. One important contribution of the cyclononstationary framework is to extend the principles of order analysis to any type of signal, including both deterministic and random types. In that context, the IAS constitutes a fundamental variable for their description [START_REF] Abboud | Angletime cyclostationarity for the analysis of rolling element bearing vibrations[END_REF][START_REF] Abboud | Deterministic-random separation in nonstationary regime[END_REF]. The IAS is not only useful as a reference signal; just as the heartbeat signal in medicine, it may serve as a diagnostic quantity per se. This has been early recognized for the diagnostic of electric motors and of internal combustion engines [START_REF] Yousef | Instantaneous angular speed monitoring of electric motors[END_REF][START_REF] Yang | Fault detection in a diesel engine by analysing the instantaneous angular speed[END_REF][START_REF] Desbazeille | Model-based diagnosis of large diesel engines based on angular speed variations of the crankshaft[END_REF] for instance and, more recently and perhaps more surprisingly, for the diagnosis of rolling element bearings [START_REF] Renaudin | Natural roller bearing fault detection by angular measurement of true instantaneous angular speed[END_REF][START_REF] Bourdon | Introducing angularly periodic disturbances in dynamic models of rotating systems under nonstationary conditions[END_REF]. Research is also ongoing to use the IAS for the monitoring of machining processes [START_REF] Lamraoui | Indicators for monitoring chatter in milling based on instantaneous angular speeds[END_REF], in which case it directly reflects the cutting forces. The measurement of the IAS is ideally achieved by a tachymeter mounted on the machine. There exists various technologies, such as angular accelerometers or coder-based devices. In turn, the signals are either directly sampled in the angular domain by using the tachymetric signal as an external clock to the ADC or aposteriori resampled by using digital processing algorithms. The latter strategy is sometimes referred to as computer order tracking. These aspects are discussed in Refs. [START_REF] Potter | A new order tracking method for rotating machinery[END_REF][START_REF] Fyfe | Analysis of computed order tracking[END_REF][START_REF] Bossley | Hybrid computed order tracking[END_REF][START_REF] Li | The measurement of instantaneous angular speed[END_REF][START_REF] Saavedra | Accurate assessment of computed order tracking[END_REF][START_REF] Wang | The combined use of order tracking techniques for enhanced fourier analysis of order components[END_REF]. Technological aspects of the measurement of the IAS and associated errors are further investigated in Refs. [START_REF] Gu | An investigation of the effects of measurement noise in the use of instantaneous angular speed for machine diagnosis[END_REF][START_REF] Yu | A data processing method for determining instantaneous angular speed and acceleration of crankshaft in an aircraft engine-propeller system using a magnetic encoder[END_REF][START_REF] Janssens | Zebra tape butt joint detection and correction algorithm for rotating shafts with torsional vibrations[END_REF][START_REF] El | Impact of the non-uniform angular sampling on mechanical signals[END_REF][START_REF] Bourdon | Introducing angularly periodic disturbances in dynamic models of rotating systems under nonstationary conditions[END_REF][START_REF] Rivola | Zebra tape identification for the instantaneous angular speed computation and angular resampling of motorbike valve train measurements[END_REF][START_REF] André | Precision of the {IAS} monitoring system based on the elapsed time method in the spectral domain[END_REF][START_REF] Leclère | An analysis of instantaneous angular speed measurement errors[END_REF]. The provision of a tachymeter with sufficient resolution is not always possible for technological or cost reasons. In these situations, endeavors have been made to recover the IAS from the machine signals (e.g. vibrations). This is a most difficult situation, yet of major interest since it does not require any specific instrumentation and potentially returns the exact IAS of interest (e.g. the IAS used for subsequent resampling of the machine signals). One of the first successful attempt towards this direction is probably the non-parametric approach based on complex envelope demodulation proposed in Ref. [START_REF] Bonnardot | Use of the acceleration signal of a gearbox in order to perform angular resampling (with limited speed fluctuation)[END_REF] and further refined in Ref. [START_REF] Combet | An automated methodology for performing time synchronous averaging of a gearbox signal without speed sensor[END_REF]. It applies to low speed fluctuations characterized by a narrow-band spectrum as compared to the carrier frequency. Following this impetus, several demodulation methods have been proposed during the last six years in order to extract the IAS from complex vibration signals in more general scenarios [START_REF] Zimroz | Measurement of instantaneous shaft speed by advanced vibration signal processing-application to wind turbine gearbox[END_REF][START_REF] Urbanek | Time-frequency approach to extraction of selected second-order cyclostationary vibration components for varying operational conditions[END_REF][START_REF] Zhao | A tacho-less order tracking technique for large speed variations[END_REF][START_REF] Michael | Single and multi-stage phase demodulation based order-tracking[END_REF]. Model-based (e.g. Kalman filter) ap-proaches demonstrate good accuracy [START_REF] Dion | Tracking and removing modulated sinusoidal components: A solution based on the kurtosis and the extended kalman filter[END_REF][START_REF] Rodopoulos | A parametric approach for the estimation of the instantaneous speed of rotating machinery[END_REF][START_REF] Cardona-Morales | Nonlinear model for condition monitoring of non-stationary vibration signals in ship driveline application[END_REF], yet a patent limitation is that all harmonics have to be modeled -even those which are not of interestso as to produce a Gaussian residual. An important body of literature has also addressed the problem by using tailored time-frequency decomposition such as the Gabor transform and wavelet transforms [START_REF] Konstantinos | Estimation of the instantaneous rotation speed using complex shifted morlet wavelets[END_REF][START_REF] Shao | Order tracking by discrete gabor expansion. Instrumentation and Measurement[END_REF][START_REF] Cheng | Application of the improved generalized demodulation time-frequency analysis method to multi-component signal decomposition[END_REF][START_REF] Feng | Time-frequency analysis of timevarying modulated signals based on improved energy separation by iterative generalized demodulation[END_REF][START_REF] Czarnecki | The instantaneous frequency rate spectrogram[END_REF][START_REF] Ioannis | Ifesis: Instantaneous frequencies estimation via subspace invariance properties of wavelet structures[END_REF]. The variety of approaches that have been proposed so far reflects the inherent difficulty of the problem. One difficulty actually stems from the multi-component nature of the signals where different families of harmonics may coexist, each being characterized by numerous components (several tens or hundreds of harmonics are commonplace in some applications). Another difficulty stems from the interaction between the orders and the structural resonances of the machine. One last difficulty is to operate under low SNRs as is typically the case with vibration signals (the noise being understood as any component in the signal which is not of direct interest for the analysis). All these challenges jeopardize the application of standard signal processing tools dedicated to instantaneous frequency estimation [START_REF] Boashash | Algorithms for instantaneous frequency estimation: a comparative study[END_REF][START_REF] Boashash | Estimating and interpreting the instantaneous frequency of a signal. i. fundamentals[END_REF][START_REF] Boashash | Estimating and interpreting the instantaneous frequency of a ii. algorithms and applications[END_REF][START_REF] Barnes | The calculation of instantaneous frequency and instantaneous bandwidth[END_REF] and ask for a specific approach. In the presence of multi-component signals, it is actually advantageous to track all harmonics jointly instead of individually as is proposed by most of the methods of the state-of-the-art. Since more information is considered at once, the estimation of the IAS is expected to gain in accuracy. At the same time, a gain in robustness is also envisioned in scenarios where some of the harmonics happen to fade or cross with other families. It is specifically the object of this paper to propose a multi-harmonic tracking method. The present work finds its origins in the organization, by its authors, of a contest during the The fourth International Conference on Condition Monitoring of Machinery in Non-Stationary Operations (CMMN0 2014). Contestants were asked to extract an IAS from the vibration signal of a wind turbine's gearbox. The first section of this paper is dedicated to a debriefing of this contest. The analysis of contributions, also discussed in the first section, led the authors of the present work to propose an innovative approach for blind IAS estimation, based on a probabilistic formalism: this is the object of the second section. The third section illustrates the efficiency of the method, in the frame of the CMMNO'14 contest signals.

CMMNO 2014 diagnosis contest

The fourth International Conference on Condition Monitoring of Machinery in Non-Stationary Operations (CMMN0 2014) was held in Lyon (France) on the 15th and 16th of December 2014. More than one hundred people attended this event, from university and industry, among which several international level experts in the field of rotating machine monitoring. Several months before the event, a diagnosis contest was launched: a signal, as well as some kinematic data were provided to participants. The aim of the contest was to estimate the IAS of the rotating machine, and to propose a diagnosis of potential faults.

The contestants had to send back the estimation of the instantaneous speed, as well as a diagnostic report, a few weeks before the start of the conference. 23 people registered to the contest, but only 8 contributions from 7 different countries have been received in due form by the Organization Committee. The results were presented during the afternoon plenary session of the 16th of December, and the first prize was awarded to Bob Randall, Michael Coats and Wade Smith (University of new south Wales, in Sidney). They provided the best estimation of the IAS, as compared to a reference that was directly measured with an angle encoder. The next section briefly resumes the case study addressed during the contest.

The case study: a wind turbine gearbox signal

Though gearbox defects appear to be 50% less frequent than electric system faults, the main cause of overall wind turbine downtime is related to gearbox faults, as is can be seen in Fig. 1 and according to WindStats database analysis focused on more than 30000 wind turbines [START_REF] Sheng | Report on wind turbine subsystem reliability -a survey of various databases[END_REF]. This fact simply shows that these occurrences are longer to be dealt with. However, bearing or even gearbox replacement never exceeds a few days. Most of the time is actually spent to order new components, manufacture and bring them on site along with available technicians. In case cranes are needed, extra time must be spent preparing the platform and bringing a crane on site, which can be as long as obtaining a spare bearing. Hence, improving condition monitoring system ability to detect gearbox fault as soon as possible allows limiting turbine downtime to the duration of the replacement itself.

Although vibration monitoring benefits from wide academic and industrial experiences, most of the works are focused on time based methods. Industrial equipments are therefore mainly dedicated to time based observations and these simple systems are widely used in the wind industry for economical reasons. Inherently relying on unsteady operating conditions, wind turbines nevertheless presents the particularity to produce strongly non stationary measurements. Rotating speed and operating torque are constantly fluctuating and the use of angle domain observations seems unavoidable. This shows the necessity to propose an approach based on time based acquisition system to tackle non stationary problems. In particular, this paper focuses on a new method to resample the signal in the angle domain without the use of complex hardware or even a tachymeter. Common to any wind turbine gearbox vibration signal are 3 kinds of periodic phenomena:

• Time dependent phenomena produced by subsystems such as the converter chopper frequency, the electric gear pump to help lubricant flow through the system, the heat exchanger fan, the natural frequencies of the structure (tower, blades, gearbox housing ...) • Angle dependent, or more precisely, kinematic-dependent phenomena pro- duced mainly by gears and shafts unbalance. These phenomena have a fundamental frequency in constant ratio with the high speed shaft.

• Other phenomena such as those produced by the slip frequency of the asynchronous generator, and the yaw gearbox system (used to rotate the nacelle in the wind direction). These are neither time stationary (they depend on rotor angular speed for the first one and on resisting torque for the second one) nor angle stationary (they are not proportionally linked with angular speed). The signal is polluted by different sources of noise which all share the characteristic of not being angle dependent. Angular domain analysis will therefore not only focus on the observation on gearbox-kinematic dependent phenomena, but it will also blur all phenomena that are not steady in the angle domain as soon as the rotating speed varies. In the CMMNO contest, the signals are recorded with IFM-Eletronics VSE002 acquisition systems, sampled at 20kHz and quantified with a 12bits resolution. An accelerometric sensor is located on the gearbox housing, in the vicinity of the epicyclic gear train. A fault is located on the inner ring of planet bearing C (see Fig 2) and introduces an angular periodic impulsive phenomenon at 0.191 high speed shaft orders in the signal.

Analysis of the contributions of the contestants

The contestants were asked to provide an estimation in the angular domain of the rapid shaft IAS (carrying gear # 7 in Fig. 2), with a resolution of 8 samples by rotation. Results were compared to the reference speed directly measured with an angle encoder. All IASs (estimations and reference) were resampled in the time domain to ease the comparison. They are displayed in Fig. 3. It is clear from this figure that the results provided by the contestants are strongly dispersed around the reference (±15%), even if the general shape is correctly obtained. Different techniques were employed by the contestants to extract the IAS from the vibration signal, which can be classified in two families:

• single harmonic tracking in the time-frequency map,

• narrow-band filtering around one selected harmonic followed by instantaneous phase demodulation. The winning team used a method belonging to the second family (see [START_REF] Michael | Single and multi-stage phase demodulation based order-tracking[END_REF] for details). However, the reason of their success is less related to the method they used than to their ability to associate the selected harmonic to the correct periodic event of the kinematics. Indeed, all contestant have used more or less sophisticated techniques allowing the fine tracking of one given har-monic in the signal, but only the winning team has correctly associated this harmonic to the right phenomenon. This is the reason why the general shape of the IAS obtained by the contestants is correct (efficient tracking of the instantaneous frequency of the selected harmonic), yet with a scaling error due to the association to a wrong mechanical component. The difficulty of this second task is explained in the present case by the relative complexity of the provided signal : about 50 harmonics are clearly visible in the time frequency map (see Fig. 6 in section 4). Among these numerous harmonics, some happen to be related to mechanical elements that are not included in the main kinematic data provided for the contest. For instance, the strong harmonic starting at about 650Hz at t = 0 in the time-frequency map (Fig. 6), which is one of the strongest components is the signal, is not related to the main meshing frequencies of the gearbox. Moreover, it is important to note that it is not strictly a multiple of the rotation frequency, which is probably due to imperfect coupling between the main mechanism and the mechanical system responsible of these few asynchronous components. Unfortunately, most of the contestants choose this harmonic to estimate the IAS, probably because it is one of the most energetic. The dispersion of the results is due to the following step, which consisted in associating this harmonic to one of the meshing orders that where available through a basic analysis of the kinematics. This step was quite hazardous, when starting with a component that was not in the expected order chart. Finally, the extracted IAS profiles resulted from the association of this particular harmonic to various (but most of the time wrong) meshing orders of the main mechanism, which explains the different scaling errors in the results. It is important to note that this specific difficulty was not introduced purposely in the contest. Indeed, this particular case has been chosen because the vibration signal, as well as the angle encoder signal, were acquired in relatively good conditions, without too much background noise: this case was considered to be favorable by the organization committee. It is also realistic that signals obtained from industrial cases are most of the time complex, with some characteristics that are not always corresponding to a priori expectation, leading to potential misinterpretation. The bottom line of the CMMNO'14 contest is that the most delicate step of blind IAS estimation was to correctly associate one arbitrarily chosen harmonic to the right order of the kinematics. In the following, a robust probabilistic approach is proposed, which prevents the necessity of admitting this association as an hypothesis.

The multi-order probabilistic approach

The originality of the proposed approach is to use the whole spectrum harmonics, as well as all known potential orders, at each time step, to build a probability function of the IAS. The basic idea is to use, besides a response sig-nal, all the information that is available about the kinematics of the machine. The knowledge (eventually partial) of the kinematics leads to a list of potential orders that are supposed to strongly contribute to the vibration response, typically meshing frequencies and their harmonics. Each of these orders is used to transform a short time instantaneous spectrum into a probability density function of the fundamental rotation frequency. These pdfs are then combined to obtain a global pdf, allowing the extraction of the most probable IAS.

Expected shape of the short time instantaneous response spectrum

The response spectrum results from the combination of the structural response and the excitation spectrum. In the case of strongly periodical excitations, the effects of the structure and of the excitation are generally well separated: the former governs the general envelop of the spectrum, while the latter generates sharp peaks at integer multiples of the frequencies of periodic events constituting the excitation. The response of the structure is without interest considering the issue of identifying the IAS, and can actually be disturbing for the approach proposed in this work. A preliminary work therefore consists in suppressing, or at least attenuating the effect of the structure, by spectral whitening the signal. Several approaches can be used for spectral whitening -one of which being described in Sec. 4-yet this preliminary step is considered as out of the scope of this paper, and will not be developed further in this theoretical part.

Considering the instantaneous response spectrum as a probability density function of the IAS

Considering the knowledge of one given excitation order taken from the kinematics, noted H i (in events by revolution of the main shaft), the instantaneous spectrum of the response provides some information to build a pdf of the IAS Ω. Indeed, if the spectrum has a strong (resp. weak) amplitude at the frequency f , then there is a high (resp. low) probability that the shaft frequency is equal to f /H i . Moreover, the range of the IAS has to be defined. This range is an a priori information that will be formalized by a prior pdf. In this work, this function will be considered as constant between a left bound Ω min and a right bound Ω max . Finally, a conditional pdf can be defined as follows (case 1)

[Ω|H i ] = 1 ξ i A(H i ω) for Ω min < ω < Ω max (1) [Ω|H i ] = 0 for ω < Ω min | ω > Ω max
where A(f ) is a whitened version of the instantaneous spectrum of the response signal and ξ i a normalization factor:

ξ i = Ωmax Ω min A(H i ω)dω. ( 2
)
The necessity of the whitening operation is here obvious: if part of the spectrum is amplified by a strong modal behavior, it has to be compensated to avoid giving arbitrarily too much probability to the corresponding angular speed range.

Equation ( 1) is fundamental in this work because it establishes how a simple response spectrum, truncated between [Ω min ; Ω max ] and scaled by a factor H i , is considered as a probability density function of the IAS, noted Ω, expressed in Hz. Of course, the use of order H i only will conduct to a strongly multi-modal pdf (as many modes than there are harmonic components in the spectrum), that cannot be used alone. Several orders will lead to several pdfs, that have to be combined to obtain a meaningful result, i.e. an unimodal pdf (cf. Fig. 4). The definition expressed by Eq. ( 1) cannot be used in all cases, for the whole range of ω ∈ [Ω min , Ω max ]. The spectrum of the response is indeed known up to one given frequency noted f max . Then, for values of ω above f max /H i , the spectrum of the response doesn't bring any information. This situation appears if the considered harmonic H i is greater than f max /ω. In this case, the conditional pdf is made uniform above f max /H i (case 2):

[Ω|H i ] = 1 ξ i A(H i ω) for Ω min < ω < f max /H i (3) [Ω|H i ] = 1 (Ω max -Ω min ) for f max /H i < ω < Ω max [Ω|H i ] = 0 for ω < Ω min | ω > Ω max
where A(f ) is still the whitened instantaneous spectrum and ξ i is such that

ξ i = (Ω max -Ω min ) (f max /H i -Ω min ) fmax/H i Ω min A(H i ω)dω. ( 4
)
An illustration of the construction of [Ω|H i ] from A(f ) and H i in the two aforementioned cases is provided in Fig. 4. For orders H 1 and H 2 , f max /H i is higher than Ω max , Eq. ( 1) can be used (case 1). For order H 3 , f max /H i is lower than Ω max , Eq. ( 3) has to be used (case 2: a constant probability is considered between f max /H i and Ω max ). Finally, the probability density functions are multiplied to combine the information from all considered orders (Fig. 4, Subfig. (4)):

[Ω] ∝ n i=1

[Ω|H i ], (5) ω ω Ω ω
): case 1 (f max /H i > Ω max ). ( 2
): case 2 (f max /H i < Ω max ). (4): combined pdf. 3
a formula which reflects the assumption that all random variables Ω|H i are assumed independent.

A priori of continuity of the IAS

The method described in the previous section allows the construction of a pdf of the IAS at one specific time step from a short time Fourier transform. In non-stationary conditions, this operation is realized independently for each consecutive time steps and the continuity of the most probable IAS value is thus not guaranteed, although it is an important assumption for any mechanical system because of the inertia of rotating shafts preventing too strong accelerations/decelerations. The general concept of introducing an a priori of continuity is depicted in Fig. 5. Let's note [Ω j ] the pdf of the IAS at time step i. The a priori of continuity of the IAS will be enforced as follows:

[Ω j |Ω j+k ] = N (Ω j+k , σ k ) ∝ exp -(ω -Ω j+k ) 2 2σ 2 k ( 6
)
where σ k = |γk∆ t |, ∆ t the time step of the frequency map and γ the standard acceleration of the IAS. A pdf of the rotation frequency at time step j can thus be expressed from the pdf at time j + k:

+ k (k ∈ [-K...K]) [Ω j ] j+k = Ωmax Ω min [Ω j |Ω j+k ][Ω j+k ]dω ∝ exp -ω 2 2σ 2 k * [Ω j+k ], (7)
where [Ω j ] j+k is the pdf of the IAS at time j, obtained by convolving the pdf [Ω j+k ] at time j + k with a centered Gaussian. Several pdfs are thus available for the IAS at time j, from several time steps j + k before and after j (k ∈ [-K...K]). The variance of the Gaussian law (6) increases with k, which increases the smoothing effect of the convolution. Note that this formulation is also correct for k = 0, for which the standard deviation of the Gaussian law is null. In this case, the Gaussian is a Dirac distribution, and Eq. 7 becomes simply

[Ω j ] j = [Ω j ].
Finally, a new version of the pdf of the IAS at time j, including the a priori of continuity, is expressed as follows

[Ω j ] s ∝ K k=-K [Ω j ] j+k (8)
where K stands for the considered time interval around j.

Application to the wind turbine's gearbox

The proposed approach is now applied on the vibration signal of the CMMNO'14 contest. The meshing frequencies of the wind turbine gearbox are given in table 1, with respect to the rotation of the rapid shaft. This table is deduced from the information provided for the contest, given in Fig. 2. The vibration signal was recorded for about 10 minutes, at a sampling rate of 5kHz. A timefrequency map of the signal is given in Fig. 6 (left). The parameters chosen for this spectrogram are • T = 500ms short time window length,

• zero padding 300% of the short time window length (to achieve a good frequency resolution), • overlap: 50% of the short time window length,

• uniform weighting window (to ensure a minimum blurring of harmonic components). Firstly, the time-frequency map is whitened, in order to remove the response of the structure. The strong fluctuations of the IAS are used here to separate the harmonics from the wideband noise. For each frequency line of the spectrogram, it is considered that the less energetic time steps are not excited by harmonics. An estimation of the wideband noise is thus obtained by averaging at each frequency the less energetic time steps, the number of which being a priori defined as a fixed proportion of the total time steps. This operation can be written as follows :

A(f, t) = S(f, t) W (f ) with W (f ) = S(f, t) × {S(f, t) < S n (f )} t ( 9
)
where S(f, t) and A(f, t) are the initial and whitened spectrograms, t stands for an averaging operator over t, {C} returns 1 (resp. 0) if the condition C is true (resp. false), and S n (f ) represents the value exceeded by S(f, t) during n% of the time at the frequency f (this notation is equivalent to the n-Percent Exceeded Level in environmental acoustics). For the present case of the wind turbine, n is chosen equal to 75%, yet it is noteworthy that this parameter is not very sensitive, and similar results are obtained by setting it between 50 and 90%. This averaging operation over the less energetic time steps (defined at each frequency) gives an estimation of the average wideband spectrum, and the whole time-frequency map is multiplied at each time step by the inverse of this spectrum. It can be seen in Fig. 6 (right) that the wideband component of the time-frequency map is efficiently flattened and it can be therefore assumed that the effect of the structure has been strongly attenuated. An interest of this particular whitening approach, which is efficient only in case of non-stationarity of the rotation speed, is that fixed frequency components are removed because they contribute at all time steps to the same frequencies. Thus, these components are present in the estimation of the wideband spectrum and are thus removed by the whitening operation. See for instance, in Fig. 6, the constant harmonics between 1200 and 1400 Hz in the original map (left) that are no more visible in the whitened map (right).

Once the whitening step is completed, each column of the spectrogram is used as one instantaneous spectrum to build a pdf of the IAS at each time step, as described in section 3.2. These pdfs will be presented in the following as pdf maps. The pdf map obtained using the 6 fundamental meshing orders given in table 1 is displayed in Fig. 7 (left). The time evolution of the IAS is visible in this figure, but the map is not really clean due to a lot of secondary local maxima still visible outside the main frequency trace. The 10 first harmonics of all fundamental meshing orders are then considered, that is to say a total of 60 orders, to obtain the pdf map given in Fig. 7 (right). The map is much clearer, and the dynamic much wider, which illustrates the interest to take into account as many orders as possible in the approach. This could be received as a satisfying result, but at this stage the a priori of continuity of the rotation frequency has still not been taken into account. In other words, the pdfs obtained at each time step are obtained from the short time Fourier transform: one column of the spectrogram is used to build one column of the pdf map. The last step of the probabilist approach, described in section 3.3, is to enforce an a priori of continuity to the IAS. At this stage, all pdfs in a given time interval around the considered time step are combined to obtained a "smoothed" one, embedding the a priori of continuity. The key parameter for this operation is the standard deviation of the Gaussian law (Eq. ( 6)), depending on the a priori tolerance γ on the acceleration. In this work, γ is choosen equal to 0.2Hz/s. This choice can be made by inspecting the previously obtained pdf map, and roughly estimating the rate of change in the IAS. However, it is noteworthy that this parameter is not too sensitive and equivalent results are obtained in a range from 0.1 to 0.4Hz/s. The parameter K in Eq. ( 8) is chosen so as to take into account, for the estimation of the smoothed pdf at time t 0 , the pdfs in the time window [t 0 ± 10s]. Considering that the time resolution of the spectrogram is 250ms, K is equal to 40. The pdf map resulting from this operation is drawn in Fig. 8. It is significantly cleaned and the resulting pdfs are characterized by clear and well defined maximum values. The expected values of pdfs are extracted, at each time step, to obtain the evolution as a function of the time (Fig. 9, left). This result is compared, in the same figure, with the IAS directly measured with the angle encoder (the reference of the CMMNO'14 contest). The estimation error is also calculated (Fig. 9, right), and it can be seen that this error, most often comprised between ±0.1% never exceeds ±0.4%.

Conclusion

The organization of a contest, during an international conference, on a specific application case, is a very interesting and original way of establishing the state-of-the-art in a particular domain. It gives an objective picture of what people actually use to solve one given problem. The short delay left to contestants to bring a solution is also a key point, because it requires them to use their existing techniques and algorithms, with limited possibilities to develop other specific approaches. It thus allows to assess the robustness and generality of up-to-date techniques. In the present case, the debriefing of contributions led to a valuable analysis, emphasizing a common weakness of existing methods: the fact that the IAS is generally tracked using only one harmonic component. The method proposed in this work has been conceived from this observation, and the problem of estimating the IAS has been reconsidered in the frame of a probabilistic formalism, trying to take advantage of any available a priori information on the studied mechanism. This approach has been evaluated in the frame of the contest case and showed a significant improvement of results as compared to the contest contributions.

Fig. 1 .

 1 Fig. 1. Wind stat data base, 2008-2012 aggregated downtime per turbine system[60].

Fig. 2 .

 2 Fig. 2. CMMNO 2014 contest: Kinematic data of the gearbox.

Fig. 3 .

 3 Fig. 3. Results of the CMMNO 2014 contest. Left: estimated IASs compared to the reference (solid blue), right: relative error with respect to reference.

Fig. 4 .

 4 Fig.4. Construction of the probability density functions of the rotation frequency, from the instantaneous spectrum and 3 orders of the fundamental rotation frequency. (1),(2): case 1 (f max /H i > Ω max ). (3): case 2 (f max /H i < Ω max ). (4): combined pdf.

Fig. 5 .

 5 Fig. 5. Estimation of the smoothed pdf at time j from pdfs at times j + k (k ∈ [-K...K])

Fig. 6 .

 6 Fig. 6. Time-frequency map of the vibration signal before (left) / after (right) whitening.

Fig. 7 .

 7 Fig. 7. Pdf of the IAS at each time step, considering fundamental meshing orders (left) and considering fundamental meshing orders and their 10 first harmonics (right).

Fig. 8 .Fig. 9 .

 89 Fig. 8. Pdf of the IAS at each time step, considering fundamental meshing orders and their 10 first harmonics, with smoothing.

Acknowledgements

The authors want to acknowledge the CMMNO 2014 Organization Committee for its help and support. All contestants are also warmly acknowledged for their contributions. The company OROS is acknowledged for providing the first prize award. This work was performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).


Joseph Gril

Modelling mechano-sorption in wood through hygro-locks and other approaches

Keywords: Wood, viscoelasticity, moisture, mechano-sorptive effect, rhéological model

Mechanosorption refers to the phenomena observed when wood is subjected simultaneously to mechanical loading and sorption resulting from changing hygrothermal conditions. A model based on the concept of hygro-locks has been recently reformulated. The hygro-lock spring further extends under load when softened, but is blocked when hardened. This approach offers an alternative to the widespread mechano-sorptive dashpot where the absolute value of moisture change replaces time increment. It will be presented and analysed in comparison with other models, and suggestions for further work will be discussed.

Introduction

Mechanosorption is a phenomenon resulting from the combined action of mechanical load and moisture content (m.c.) change in wood. Creep, relaxation or damping trigerred by changing humidity, expansion ratio modified by the load, permanent set after restrained expansion, are the most well-known expressions of the mechanosorptive (m.s.) effect. There is no clearly convincing explanation for it. The 'breaking of hygrogen bonds' resulting from sorption and reformed at different position under stress [START_REF] Grossman | Mechano-sorptive behaviour, in: General constitutive relations for wood and woo-based products[END_REF] does not explain why drying wood creeps more than humid wood, where bonds are also continuously broken an reformed. The transcient increase of molecular mobility resulting from a modification of m.c. was described as a sort of physical de-aging by [START_REF] Hunt | Evidence of a physical ageing phenomenon in wood[END_REF], a temporary disorder of the molecular structure similar to that existing in polymers heated above their glassy transition temperature and suddendly cooled. Such a physical aging in the classical sense, i.e. temperature induced, was actually characterised in wet wood (Dlouha et al 2009). [START_REF] Mukudai | Modeling and simulation of viscoelastic behavior (tensile strain) of wood under moisture change[END_REF] proposed for an explanation based on cell-wall delamination at the S1/S2 interface to explain mechanosorption. The sophisticated process of load transfer between inner and outer cell-wall layers resulted in complex equations that were hardly usable in practice. However it raised the important issue of the relevant structural scale to consider, in this case the cell-wall layers even, with an underlying molecular process concentrated at the interface between large blocks. Their model was also consistent with the quasi disappearance of mechanosorption in the wet range.

The phenomenology of mechanosorption, on the other hand, is abundantly documented. Ranta-Maunus (1975) modelled the m.s. creep of plywood using viscous-like elements where time has Gril J. (2015) International Symposium on Wood Science and Technology (IAWPS 2015) and 60th Anniversary of the Japan Wood Research Society,Tokyo,JP,Isogai A,Iwata T. (eds) been replaced by m.c., distinguishing the desorption (-), adsorption (+) and 'adsorption above highest level reached since loading' (++). Hunt in the 80ties introduced the m.s. 'trajectories' where the longitudinal creep compliance is plotted against versus m.c., and the concept of 'creep limit' reached after repeated humidity cycles under load (Fig. 1). His results suggested that the creep limit is approached both in sorption and desorption, and that the apparent recovery observed during adsorption was easily accounted for by assuming that the moisture expansion is modified by the strain [START_REF] Hunt | Present knowledge of mechano-sorptive creep of wood, Creep in Timber Constructions (a state of the art report)[END_REF]Shelton 1987, Hunt 1992). This allowed a much simpler rheological modelling, that will be recalled in the following section. In that approach, the (++) effect is considered as sufficiently taken into account by the higher creep of wetter wood (Montero et al 2012). [START_REF] Hunt | Present knowledge of mechano-sorptive creep of wood, Creep in Timber Constructions (a state of the art report)[END_REF], showing the apparent creep/recovery at the creep limit (explained by a modified hygroscopic expansion ε h ).

At about the same period in the 80ties, the alternative 'hygro-locks' theory was proposed by the author [START_REF] Gril | Une modélisation du comportement hygro-rhéologique du bois à partir de sa microstructure[END_REF][START_REF] Gril | Hygro-mechanical threshold and memory effects in wood, 39th Annual Meeting of the Japan Wood Research Society[END_REF]. The initial aim of this model was to reproduce the features of the (-), (+) and (++) effects as introduced earlier by Ranta-Maunus and confirmed by Hunt's observations. A series of 'boxes' corresponding to increasing m.c. levels contributed to the strain through a spring, either blocked in its current configuration during drying or liberated during wetting. This approach was abandonned for two reasons: first, it was not able to account in a simple way to some loading cases such as m.s. relaxation; second, the success of the explanation given by Hunt for the apparent difference between (-) and (+) situation (pseudo creep/recovery) rendered other interpretations of little use.

A theory developped in the last 15 years by a group of Limoges University, France, offered a new opportunity to investigate the applicability of hygrolocks to wood rheology [START_REF] Dubois | Creep in wood under variable climate conditions: Numerical modeling and experimental validation[END_REF]. Their approach derived from the observation that when wood humidifies under tensile load, the deformation (corrected for the expansion) increases, whereas when it dries it remains blocked. The elongation can be easily explained by an increase of the elastic compliance with m.c.:

(1) ε=J e (h).σ

In the case of a m.c.decrease, this equation is not applicable, as it would predict a contraction. Instead, an incremental formulation must be used in this case:

(2) dε=J e (h).dσ

Based on this simple idea, they developped a sophisticated formulation of a 'hygrolocks spring' that includes the classical features of the +/-/++ cases, and can be combined with time-dependent dashpots (or even sliding frictional elements) to fit with the observed behaviour of wood [START_REF] Dubois | Modeling of the viscoelastic mechano-sorptive behavior in wood[END_REF].

Recently we were able to prove that their model is equivalent to a slightly modified formulation of the original hygrolocks theory [START_REF] Colmars | One-dimensional discrete formulation of a hygrolock model for wood hygromechanics[END_REF].

In this paper we will introduce various formulations of mechanosorption, including the hygrolock model. Although not limited to 1D formulation, we have here mostly in mind the L direction where the hygroexpansion does not dominate over other sources of deformation, at load levels small enough to avoid permanent damage or maked non-linearity relative to stress.

The mechanosorptive dashpot

In his list of 'requirements' [START_REF] Grossman | Requirements for a model that exhibits mechano-sorptive behaviour[END_REF] insisted that the m.s. effect is time-independent and driven by the moisture content. This led many authors to the use of a 'm.s. dashpot' to describe the m.s.

effect:

(2) dε ms /dt = k.σ.|dh/dt| or: dε/|dh| = k.σ where ε ms is the m.s. strain, σ the stress, h the moisture content and k a constant that does not need to be given different values in adsorption or desorption; in the case of 2D or 3D formulations, k can be a tensor. The second equation highlights the formal similarity with the time-dependent dashpot: dε/dt = η -1 .σ where η is the coefficient of viscosity.

To account for a creep limit, the m.s. dashpot can be associated with a spring in parallel, leading to the 'm.s. Kelvin' element:

(3) dε ms /dt = µ -1 .(J ms σ-ε ms )|dh/dt| or: dε/|dh| = µ -1 .(J ms σ-ε ms )

where µ is a characteristic moisture range and J ms the m.s. compliance. Such description is commonly associated with a decomposition of the total strain ε:

(4) ε = ε h + ε T + ε el + ε ve + ε ms
where ε h is moisture expansion, ε T thermal expansion, ε el elastic strain, ε ve viscoelastic strain and ε ms m.s. strain. In the models actually used to fit real data and to be implemented in numerical codes, ε ve and ε ms can be decomposed in several terms (generalized Kelvin-Voigt), and the constants are given tensorial forms. To account for the reversible pseudo-creep/recovery, the expansion can be written as:

(5)

ε h = (α -β.ε).(h -h 0 )
where h 0 is the initial moisture constant, α the initial moisture expansion ratio and β a constant.

Alternative representations using stress decomposition have been also proposed (generalized Maxwell).

In the case of time-dependent viscoelasticty, parallel or series arrangement of rheological blocks are theoretically equivalent, and can be related to each other through Laplace-Carson transformation. This certainly does not apply to models containing m.s. elements, as they introduce a non-linearity through the coupling between sorption and mechanical actions. Therefore model developers should not hesitate to explore all possible combinations of parallel and series representations in their search for better fitting to the observed behaviour in the widest range of conditions.

A major difficulty with the m.s. dashpot as defined by ( 2) or ( 3) is the presence of the term |dh|, which is poorly justified theoretically and makes the numerical integration unstable relative to the time discretization: a quasi stability with humidity fluctuations does not tend to the response at constant humidity when the level of fluctuation vanishes.

The mechanosorptive activation

The time-independent feature of m.s was mostly based on early results such as Schniewind (1966) who compared two drying schedules, slow and rapid but in the same range of magnitude. This question is still a matter of debate, as it is practically very difficult to separate the time-dependency of water sorption and diffusion, from that of mechanical viscosity. As a radically different approach, mechanosorption can be viewed as the acceleration of viscoelastic processes occurring at constant humidity and temperature. The formalism of thermal activation can be used to express this hypothesis.

Let consider, for instance, the case of a generalized Kelvin model:

(6) ε = ε 1 + … + ε N ; σ = σ 1 = … = σ N ; dε i /dt = (τ i ) -1 .(J i σ i -ε i ) (i=1..N)
where ε i , σ i , J i , τ i are the strain, stress, compliance and characteristic time of the i-th link, respectively, and N is the number of links. When only the effect of temperature T is considered like in typical polymers, the principle of time-temperature equivalence can be expressed by Arrhenius equations:

(7) τ i (T) = τ i ∞ exp[W i /RT
] where W i is the activation energy for the link i, R the gas constant and T the absolute temperature. The physical meaning of W i is the energetic barrier to overcome for configurational changes, while RT represents the thermal energy that allows them. The softening effect of humidity can be modelled by a lower activation energy:

(8) τ i (T, h)= τ i ∞ exp[(W i -a i h)/RT] (a i = cte)
leading to an equivalence between temperature T, moisture content h and time. To further account for the m.s. activation, one can subtract an additional term from the activation energy:

(9)

τ i (T, h, dh/dt) = τ i ∞ exp[(W i -a i h -b i |dh/dt|)/RT] (b i = cte)
In other words, the temperature controls the capacity of molecules to overcome energetic barriers, while the humidity (and its variations) controls the height of the barriers.

When the characteristic times τ i are regularly spaced in a log scale for reference hygrothermal condition, the generalized Kelvin-Voigt model can be used to discretize the compliance spectrum, and equations ( 8) or ( 9 Hanhijärvi (1995) proposed to add a m.s. dashpot in series to each time-dependent dashpot (Fig. 2). Gril (1987) gave another alternative to equation ( 9), with the drop of τ i (through that of W) operating for a fixed duration after every change of moisture content by a given amount (Fig. 3).

Fig. 3: Effects of a moisture content cycle on the activation energy: (i) direct effect, the energy barrier is higher at low moisture in order to predict a slower creep; (ii) creep activation due non-equilibrium state during drying periods. From [START_REF] Gril | A uniaxial model for hygro-thermo-mechanical loadings of wood[END_REF] The hygrolock formulations

The hygrolock theory is based on the 'hygrolock spring' (Fig. 4).

This body can be modelled by a combination of n Hooke springs, each spring k=1...n standing for a bounded moisture domain, a "box", of increasing level through the hygroscopic domain from ovendry to fiber saturation point (Fig. 5). The stress S and strain E applied on the body are the respective sums of the stresses σ k and strains ε k of all boxes, linked to each other locally by Hooke's law:

(9) S = Σσ k ; E = Σε k ; ε k = J k σ k (k = 1..n)
The boxes are neither disposed in series nor (10)

E t ( ) = J max ! t , t ( ) . ∂S ∂t d ! t 0 - t ∫
where J max (t',t) means the highest compliance of the material between times t' and t. Fig. 6 illustrates the way the m.s. springs Figure 6: Simulation test of the hygrolock spring. Stress and moisture content cycling is imposed, strain is calculated. Numbers refer to significant events: (1) HygroLock spring is wetted under load, compliance and strain increase.

(2) Spring is dried, compliance decreases but strain is blocked; (3) Higher moistening leads to strain increase. (4) Spring is unloaded at low humidity, remaining strain has been blocked during the previous drying. ( 5) New moistening phases allow to recover strain completely. (6) Now spring is loaded wet. (7) Unloading at dry state ; again strain has been blocked in the upper boxes. (8) Small moistening allows recovering partly residual strain. (9) A new moistening to the condition when loaded allows complete recovery. From Colmars et al 2014.

Hygrolock springs can be combined with dashpots to build generalized Kelvin or Maxwell models (Fig. 7). The only data required to determine the parameter values (the J k for each hygrolock spring) are creep or relaxation functions at a range of fixed humidity levels. In practice, viscoelastic data at lowest and highest humidity levels are sufficient, together with appropriate interpolation. Hence this type a model is essentially a reformulation of the (visco)elastic relationships, in order to take into account features of the m.s. effect.

τ 1 τ 2 τ j τ m J 1 J 2 J j J m J 0 ε 1 ε 2 ε j ε m ε 0
A simple example of such simulation is shown in Fig. 8. Fig. 8: Simulation of creep-recovery at changing humidity using a viscoelastic hygrolock model. From [START_REF] Dubois | Modeling of the viscoelastic mechano-sorptive behavior in wood[END_REF] More complex hygromechanical loadings have been considered by [START_REF] Saifouni | Modélisation des effets rhéologiques dans les matériaux: application au comportement mécanosorptif du bois[END_REF], who managed to obtain good agreement with experimental data (Fig. 9).

Conclusion

The challenge of a robust description of the hygro-thermo-mechanical behaviour of wood, applicable to a wide range of environmental and loading conditions, still exists after more than 50 years. The practical need for numerical implementation is important but should not prevent us from exploring alternative models that ultimately could better take into account the complexity of the material.

Hygrolock models offer a consistent formulation for the dependency of the viscoelastic behaviour on humidity. Without requiring any additional parameter, the produce naturally some features of the m.s. effect. They need to be developped further, especially with regards to 3D hygromechanical loading over a wide range of time scales and temperature and humidity range. It is also important, in order to stimulate the interest of the wood mechanics community, to compare the various types of approach through parallel simulations and comparative experimental data fitting.

Fig. 1 :

 1 Fig. 1: Construction of mechano-sorptive trajectories: (a) Humidity cycle under constant load; h is moisture content, ε strain and σ stress (b) Corresponding trajectory, relative compliance J versus h; (c) Representation of Ranta-Maunus model; (d) effect of repeated cycles under load according to[START_REF] Hunt | Present knowledge of mechano-sorptive creep of wood, Creep in Timber Constructions (a state of the art report)[END_REF], showing the apparent creep/recovery at the creep limit (explained by a modified hygroscopic expansion ε h ).

  ) result in the dependency of the effective compliance on humidity when observed within a given experimental window: the J i do not need to depend on h (or dh/dt), just due to τ i variations the range of links concerned is modified. This type of formulation offers a sounder theoretical basis than the mechanosorptive dashpot, however it does not contain any mechanism leading in a m.s. response proportional to the amount of humidity change, which is what is commonly observed at least approximately.

Fig. 2 :

 2 Fig. 2: Concept for describing the relation of viscoelastic (ve) and mechano-soprtive (ms) creep. From Hanhijärvi and Hunt (1998).

FigFig. 5 :

 5 Fig.4: symbol for hygrolock spring as proposed by[START_REF] Husson | Elastic response in wood under moisture content variations analytic development[END_REF]

Fig. 7 :

 7 Fig. 7: Hygrolock generalized Kelvin model. Hygrolock spring of compliance distribution J 0 is associated in series with m hygrolock Kelvin-Voigt bodies of compliances distribution J j and characteristic time τ j . From Colmars et al 2014.

Fig. 9

 9 Fig. 9 Simultaneous variation of stress and humidity. Left: Experimental data. Right: comparison of viscoelastic strain calculated from the data (black) and simulated with a hygrolock model (red, dotted line). From Saifouni (2014).