Emmanuel Piat
email: emmanuel.piat@ens2m.fr

Joël Abadie
email: joel.abadie@femto-st.fr

Stéphane Oster

Nanoforce estimation based on Kalman filtering and applied to a force sensor using diamagnetic levitation

Keywords: micro and nano force sensor, magnetic spring, Kalman filtering, unknown input estimation

Nanoforce estimation based on Kalman filtering and applied to a force sensor using diamagnetic levitation

Introduction

The design of micro and nanoforce sensors is constrained by the fact that only force effects can be directly observed. Because of this, a transducer is necessary to convert the force into a measurable effect. The force is the unknown input to reconstruct and the effect is the measured output signal. Most of the time, the measured force effect is related to a displacement x and the usual scalar expression used to calculate the component F of the applied force F in one-direction x of space simply consists in the equation:

F = K x K > 0 (1)
in which K is the mechanical stiffness of the transducer along x (by convention x is set to zero when there is no displacement). This steady-state equation supposes that the transient response of the transducer can be neglected. This is usually considered to be the case for classical designs using monolithic elastic microstructures like microcantilevers which have a high resonant frequency [START_REF] Sepaniak | Microcantilever transducers: a new approach in sensor technology[END_REF]: AFM based microforce sensors [START_REF] Rougeot | Forces analysis for micro-manipulation[END_REF] [3], piezoresistive microforce sensors [START_REF] Arai | Micro tri-axial force sensor for 3d bio-micromanipulation[END_REF], capacitive microforce sensors [START_REF] Beyeler | A six-axis mems force-torque sensor with micro-newton and nano-newtonmeter resolution[END_REF], piezoelectric microforce sensors [START_REF] Shen | Contact and force control in microassembly[END_REF], etc. When the transient dynamic of the transducer due to the evolution of the successive derivatives of x is not negligible, Eq. [START_REF] Sepaniak | Microcantilever transducers: a new approach in sensor technology[END_REF] can not be used.

The general framework of the force estimation (shown in Fig. 1) is then equivalent to the deconvolution problematic of a noisy output signal or the design of a unknown input observer. In the specific case treated in this paper, the unknown input is a nanoforce that is applied to a macroscopic seismic mass that levitates passively thanks to the diamagnetic levitation principle. This seismic mass acts as a transducer that converts the unknown input force into a displacement that is the measured output signal corrupted by noise. This kind of macroscopic transducer has a under-damped and long transient response, thus this dynamic behaviour must be taken into account during the estimation process contrary to what is done in Eq. [START_REF] Sepaniak | Microcantilever transducers: a new approach in sensor technology[END_REF]. The estimation computation is based on a discrete Kalman filter that use an uncertain a priori model to represent the unknown force to be estimated. This paper begins by a short description of the force sensor and its dynamic behaviour (state-space modelling). The calibration process is then briefly presented and followed by the development of the unknown input estimation under Gaussian assumptions usually used to derive a Kalman filter. To be realistic, some performances of the force estimation obtained are then characterised in a non Gaussian framework and an analysis of the trade-off between resolution and bandwidth is developed using a steady-state Kalman filter. Finally, some experimental results are presented. A glossary of all mathematic symbols used is also present at the end of the article. Figure 1: Force estimation using a deconvolution approach or an unknown input observer (UIO).

2. Passive micro and nanoforce sensor prototype based on diamagnetic levitation

Sensor description

Microforce sensors based on "heavy" rigid seismic mass are really uncommon because of the inertia of their mass that should be a great handicap for time-varying force measurement. A force sensor with a range measurement of several millinewtons and based on a mass moving inside a pneumatic linear bearing is described in [START_REF] Fujii | Methods for generating and measuring the micronewton level forces[END_REF]. The mass is 21.17 g and the force resolution is 0.5 micronewton. The air friction inside the bearing is assumed small enough to be neglected. A force sensor based on magnetic springs and upthrust buoyancy is presented in [START_REF] Cherry | Analysis of a passive microforce sensor based on magnetic springs and upthrust buoyancy[END_REF][START_REF] Cherry | Modelling and optimization of a floating triangular platform used for nano and microforces sensing[END_REF]. The seismic mass is a float. Its weight is 4 grams. Resolution is around 10 nN and measurement range is ±50 µN. This force sensor can measure the horizontal component of an external force and the vertical component of an external torque. The design presented in this paper is also based on magnetic springs like [START_REF] Cherry | Analysis of a passive microforce sensor based on magnetic springs and upthrust buoyancy[END_REF] but uses repulsive diamagnetic force [START_REF] Boukallel | Levitated micro-nano force sensor using diamagnetic levitation[END_REF][START_REF] Boukallel | Passive diamagnetic levitation: theoretical foundations and application to the design of a micro-nano force sensor[END_REF] instead of repulsive upthrust buoyancy to stabilise the seismic mass. It reaches the same force resolution than an Atomic Force Microscope (AFM) but with a larger range measurement. This design is based on a macroscopic mass (≈ 70 mg) that is levitating passively thanks to to the diamagnetic levitation principle. This mass is a rigid 10-cm long capillary tube made of glass on which are stuck two small magnets M 2 . The whole structure is called maglevtube (see Fig. 2). As it is shown in Fig. 3, the maglevtube levitates passively around a given equilibrium state because of repulsive diamagnetic effects (generated by four graphite diamagnetic plates) coupled with attractive magnetic effects (generated by magnets M 1 and M 1 ). The maglevtube has a microscopic tip on which is applied the unknown external force F . The sensor is currently designed to only measure forces applied along the longitudinal axis x of the tube. Thus, the unknown force F is assumed to be collinear to x and has the following components in the global reference frame R 0 given in Fig. 3:

F F x 0 0 T (2)

F M 2

Graphite

M 1 M 1 z x y R o

Force sensing principle

The forces involved in the maglevtube dynamic along axis x are the external force F , the viscous friction due to the air and the magnetic force. The maglevtube weight and the diamagnetic force are not involved because they act along axis z and y respectively. When the force F x is applied to the tube tip, the displacement obtained corresponds to an under-damped behaviour because the viscous friction due to the air is very small. The simulated displacement computed with Matlab-Simulink and obtained with a step force (which amplitude F x is set to 1 µN) is given in Fig. 4 (step response). This simulation of the prototype presented in Section 2.3 is done with a complete computation of the internal magnetic and diamagnetic forces at each time step of the Simulink solver with an non-analytical physical model. Thus, the complete behaviour of the six DOF of the maglevtube can be plotted if necessary. In the case studied here, there is only a dynamic along the axis x due to the assumption done in Eq. ( 2). The settling time at 5% along x axis is typically 20 seconds. Overshoot is 97%. The nonlinear steady-state response of the maglevtube is given in Fig. 5. The slope of this curve corresponds to the magnetic stiffness K x m of the sensor that is equivalent to an invisible magnetic spring with a very small damping. One can notice that the linearity of the stiffness is good with displacements between ±1.5 millimetres. For such range of displacements, the maximum relative error between the linearised force and the nonlinear magnetic force is 0.63% in this simulation [START_REF] Abadie | Modeling and experimentation of a passive low frequency nanoforce sensor based on diamagnetic levitation[END_REF]. As a consequence, the behaviour of the maglevtube will be assumed linear for displacements between ±1.5 millimetres. Knowing the magnetic stiffness K x m , the force measurement is given by Eq. (1) in steady-state:

F x = K x m x K x m > 0. ( 3
)
Because the stiffness is equal to 0.0289 N/m in this simulation, the corresponding measured force range associated to a ±1.5 millimetres range displacement is ±43 µN. Eq. ( 3) can not be used for time-varying force measurement, especially if the seismic mass inertia is high. In this case, the dynamic response of the seismic mass must be taken into account in order to correctly estimate the unknown force applied on the mass. This specific problem is addressed in this paper. It is an extended analysis of the study briefly presented in [START_REF] Piat | Nanoforce estimation with kalman filtering applied to a force sensor based on diamagnetic levitation[END_REF] and it completes [START_REF] Abadie | Modeling and experimentation of a passive low frequency nanoforce sensor based on diamagnetic levitation[END_REF] in which only steady-state measurements are developed using Eq. (3).

Experimental prototype

Typical K x m stiffnesses obtained with the prototype shown in Fig. 6 are between 0.005 N/m and 0.03 N/m (same order of magnitude than for very flexible AFM cantilevers). The stiffness and the horizontal attitude of the maglevtube can be adjusted by changing the distances between magnets M 1 and M 1 . Lower is the stiffness and better is the sensitivity of the sensor. There is nevertheless a limitation on the lower value that can be reached for K x m because if magnets M 1 and M 1 are too far away from each other, the magnetic force along the vertical axis is not sufficient to compensate the maglevtube weight. Typical mass m for the maglevtube is around 70 mg. The sensor used to measure the maglevtube displacement is a confocal chromatic sensor (manufactured by STIL SA) that is aimed at a glass deflector stuck at the rear of the maglevtube (see Fig. 2). The force measurement range and the resolution obtained depends on the measurement range of the confocal head used. For instance, a CL4 (resp. CL2) head haves a ±1.5 mm (resp. ±150 µm) measurement range and covers the linear domain of the sensor. The associated typical standard deviation of a CL4 (resp. CL2) head is 75 nm (resp. 12 nm). Thus, without any signal processing and in steady-state with K x m set to 0.01 N/m, the minimal theoretical standard deviation that can be -0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 expected for the force is 0.75 nN with a CL4 head and 0.12 nN with a CL2 head. In practice, such small values can not be reached because of the seismic mass sensitivity to seismic disturbances (subsonic air disturbances can be avoided by enclosing the sensor with a chamber). Stochastic low frequency seismic vibrations of magnets M 1 and M 1 generate unwanted magnetic return forces that are applied on magnets M 2 . The result is a small stochastic oscillating behaviour for the maglevtube. All the forces induced by these environmental seismic and acoustic noises (see Fig. 1) along x axis are included in the force F x measured by the sensor. Thus care must be taken to reduce as much as possible this environmental noise. With a massive concrete ground to minimise seismic vibrations, a closed chamber and a passive anti-vibration table, the minimal standard deviation currently reached for the maglevtube displacement is around 30 nm (measured with a CL2 head).

Sensor modelling and calibration

One DOF linear modelling of the maglevtube dynamic

Let G the centre of gravity of the maglevtube and x its position in the frame R 0 attached to its base (cf. Fig. 3). The frame R 0 is supposed to be a Galilean reference frame, which implies that all ground vibrations are filtered by an antivibration table (no acceleration). This hypothesis will be later discussed in Section 6.1. If an external force F x is applied under the assumption done in Eq. (2) to the maglevtube tip, the dynamic of G along x can be modelled by [START_REF] Abadie | Modeling and experimentation of a passive low frequency nanoforce sensor based on diamagnetic levitation[END_REF]:

m ẍ = F x + F x mag + F x visc (4)
in which F x visc is the viscous friction force (mainly due to the air friction against the rear deflector) and F x mag is the return magnetic force along x (magnetic spring). The diamagnetic force is completely negligible along x. The force F x represents the sum of all others external forces applied on the tip along x. It includes all the environmental disturbance forces along x. Coordinate x is set to zero when the maglevtube is in steady-state with F x set to zero. If the displacement of the tube along axis x remains in the linear domain given in Section 2.2 and if the speed is small, Eq. ( 4) becomes:

m ẍ = F x -K x m x -K x v ẋ (5)
in which K x m is the magnetic stiffness and K x v the viscous damping coefficient. A possible state equation associated to Eq. ( 5) with F x (t) as input, x(t) as output and X(t) = x ẋ T as state is:

Ẋ(t) = A X(t) + B F x (t) (6)
x(t) = C X(t) (7)
A = 0 1 - K x m m - K x v m B = 0 1 m C = 1 0 (8)

Calibration

Calibration is usually a complex problem for micro and nanoforce sensors based on elastic microstructures because of the lack of standard forces at this scale. Stiffness absolute uncertainty is most of the time not specified and is still an open question on which are working international metrology laboratories [START_REF] Pratt | Progress toward système international d'unités traceable force metrology for nanomechanics[END_REF]. Calibrating micro force sensors based on macroscopic seismic mass is easier and several dynamic calibration methods have been investigated. They are based on particular external force generation like impact force [START_REF] Fujii | Proposal for an impulse response evaluation method for force transducers[END_REF], step force [START_REF] Fujii | Proposal for a step response evaluation method for force transducers[END_REF] and oscillating force [START_REF] Kumme | Investigation of the comparison methods for the dynamic calibration of force transducers[END_REF][START_REF] Fujii | A methode for calibrating force transducers against oscillation force[END_REF]. Because the maglevtube mass m can be easily measured with a precision balance, a Zero Input Response (ZIR) is another possible way to identify the two others parameters K x m and K x v with an ARMA identification process. It requires an unknown excitation force F x with the following dynamic:

F x (t) = 0 t 0 ≤ t < t 1 ∀ F x (9a) F x (t) = 0 t ≥ t 1 . (9b)
The output x after t 1 is the ZIR of the maglevtube with unknown initial conditions (position and speed of the maglevtube at t 1 ). The temporal shape of F x (t) before t 1 necessary to excite the maglevtube does not matter (see Eq. (9a)) but must be equal to zero after t 1 (see Eq. (9b)).

For the identification process, the behaviour of the maglevtube is taken into account only for time t superior to t 1 and is therefore a ZIR. In the prototype, F x (t) is generated using two coils located near the rear diamagnetic plates (see Fig. 6). Nevertheless in practice, Eq. (9b) can not be completely verified because of the environmental disturbance forces that always acts on the maglevtube (these forces can be considered as a noise on the input. Its amplitude is between ±1 nN (see Fig. 22) and the dynamic it induces needs to be negligible. This condition is achieved with an excitation of the maglevtube enough high). The acquisition of x starts at t 1 with a sampling time set to 0.01 seconds. Knowing the mass m, Fig. 7 shows the matching between both experimental ZIR and the continuous linear model after the parametric identification of K x m and K x v in Eq. ( 5) when F x is set to zero (done with Matlab ARMA identification). The frequency response of the identified linear model is shown in Fig. 8. Resonant frequency in displacement is around 3 Hertz. The low -3dB cutoff frequency (around 4 Hertz) will limit the resulting bandwidth during force measurements even if the force bandwidth can be reasonably extended to 20 Hertz because of the behaviour of the Kalman filter (see Section 5).

Unknown input force estimation

Despite the fact that the model-based deconvolution framework of a noisy output is an ill-posed problem with no exact solution, numerous approaches have been developed in the past such as for instance the Wiener deconvolution filter or deconvolution methods based on regularization [START_REF] Demoment | Image reconstrution and restoration: overview of common estimation structures and problems[END_REF]. In the specific context of micro or nanoforce measurement in microrobotics, the deconvolution problematic has been little addressed but some alternative approaches using unknown input observers have been recently published [START_REF] Rakotondrabe | Force estimation in a piezoelectric cantilever using the inverse-dynamics-based uio technique[END_REF]. These approaches generally requires to set several parameters and the level of noise in the input reconstructed can not be adjusted. The method proposed here only requires one parameter in order to adjust an intuitive trade-off (like in regularization methods) between the wished resolution of the force estimation (this point sets the level of noise in the force estimation) and the response time of the estimation.

A priory force modelling

The noisy measurement m x k of the maglevtube displacement along axis x with the confocal chromatic sensor is done with a sampling period T s at each sampling time t k = k T s . The estimation of the unknown force with the set {m x k } k≥0 is done at each sampling time t k . It is based on an a priori discrete-time stochastic model for the force evolution that will be processed inside a recursive discrete Kalman filter. In this approach, the maglevtube model is deterministic. The stochastic part only concerns the modelling of the uncertainty associated to the unknown input force and the modelling of the noise generated by the confocal chromatic sensor.

The working out of the a priori uncertain model used to represent the force evolution is based on a stochastic Wiener process:

Ḟ (t) = ω(t) (10)
F (t) in Eq. ( 10) is a model for the real force F x (t) and ω(t) is a zero-mean infinite-variance white Gaussian stochastic process representing the fact that the evolution of the force derivative is not known accurately. The autocorrelation function φ ω,ω of this process is characterised by its power spectral density (PSD) W Ḟ :

φ ω,ω (τ ) = W Ḟ δ(τ ) ∀τ ∈ IR (11)
The scalar term W Ḟ is a parameter to set by the enduser that will influence in a given way the dynamic of the unknown force estimation (see Section 5). Because the process is white, all random variables ω(t) for any t are independent and it is the same for the force increment between any time interval |t i -t j | [START_REF] Maybeck | Stochastic models, estimation and control[END_REF]. Thus knowing the force increment between two instants gives absolutely no information on the increment between two other instants. As a consequence, this a priori model on the force evolution is very little informative.

To estimate the force with a state estimator it is necessary to merge the uncertain modelling of the input force with the deterministic model of the maglevtube in order to obtain the following extended state:

X e (t) = x ẋ F T (12)
Its associated state-space model is obtained with Eqs. ( 10) and [START_REF] Shen | Contact and force control in microassembly[END_REF] in which the unknown input force F x (t) is replaced with the modelled random variable F (t):

Ẋe (t) = A X e (t) + Mω(t) (13)
x(t) = CX e (t) [START_REF] Pratt | Progress toward système international d'unités traceable force metrology for nanomechanics[END_REF] with

A = A11 A12 B11 A21 A22 B21 0 0 0 M = 0 0 1 C = 1 0 0 (15)
The state-space equation ( 13) is driven by the non-physical input ω(t) and thus by the parameter W Ḟ to set (see Fig. 9).

To estimate with a discrete Kalman filter the extended state X e and thus the external force F (third component of X e ) at each sampling time t k , a discretization of Eqs. ( 13)-( 14) is necessary. The use of a zero-order hold (ZOH) on ω(t) gives:

X e k+1 = F X e k + Ω k (16) x k = C X e k (17)
with:

X e k = x k ẋk F k T Ω k = ω x k ω ẋ k ω F k T (18) ω(t) Uncertainty modelling of the input force B F (t) + Ẋ(t) C X(t) x(t)
A State-space model of the maglevtube (Eqs.( 6)-( 7) with F x (t) replaced with F (t))

W Ḟ M ω(t)
Associated extented state-space model (Eqs.( 13)-( 14)) and

+ Ẋe (t) C X e (t) x(t) A W Ḟ
F = e A Ts =   F 11 F 12 F 13 F 21 F 22 F 23 0 0 1   (19)
The zero-mean white Gaussian process noise Ω k characterises uncertainties on x k , ẋk and F k due to the stochastic force model used and to the discretization of the maglevtube dynamic. Its 3 × 3 covariance matrix Q is:

Q = E Ω k Ω T k = Ts 0 e At MW Ḟ M T e A T t dt (20) = W Ḟ Ts 0 e At MM T e A T t dt (21) = W Ḟ η(T s ) (22)
It appears with Eq. ( 22) that the matrix Q is a tuning parameter derived from the scalar W Ḟ and the sampling period T s that are both chosen by the end-user. Thus the matrix Q is not associated to a physical process noise. This process noise Ω k related to X e k dynamic is only a mathematical consequence of the merging of the uncertain modelling of the input force with the deterministic model of the maglevtube followed by a discretization.

For a given sampling period T s in Eq. ( 22), Q is proportional to W Ḟ and it can be easily shown with Eq. ( 21) thanks to the specific shape of e At (see Eq. ( 19)) and M that the variance of ω

F k (third component of Ω k ) is equal to: σ 2 (ω F k ) = Q 33 = T s W Ḟ (23)
The evolution of F k (third component of X e k ) is obtained from Eqs. ( 16) and ( 18):

F k+1 = F k + ω F k k ≥ 0 (24)
Thus, if W Ḟ is set to a given value, longer is the sampling period T s in Eq. ( 23) and bigger is the variance of ω F k . That means in Eq. ( 24) that the uncertainty in the force evolution is larger. This is coherent with the fact that the longer sampling period will gives less information on the force evolution.

The statistic properties of the random process ω F k are:

E ω F k = µ(ω F k ) = 0 ∀k (25) E (ω F k -µ(ω F k )) 2 = σ 2 (ω F k ) = T s W Ḟ (26) E ω F i ω F j = σ 2 (ω F i ) δ ij = T s W Ḟ δ ij (27)
Eqs. ( 24) to ( 27) fully characterise the a priori discretetime Gaussian stochastic model on the force evolution that will be used inside the Kalman filter. The uncertainty model represented by Eq. ( 24) corresponds to a discretetime Gaussian random walk that is usually written as follow (index shift on ω F k ):

F k = F k-1 + ω F k k ≥ 1 (28)
It comes from Eq. ( 28) that:

F k = k i=1 ω F i + F 0 (29)
Because successive random variables ω F i form an a priori discrete zero-mean white Gaussian process (the white property is induced by Eq. ( 27)), F k in Eq. ( 29) is Gaussian if the uncertainty on F 0 is assumed Gaussian or if F 0 is supposed equal to some fixed value. Its a priori variance at each step k can be calculated thanks to Eq. ( 27):

σ 2 (F k ) = k i=1 σ 2 (ω F i ) + σ 2 (F 0 ) (30) = k T s W Ḟ + σ 2 (F 0 ) (31) = t k W Ḟ + σ 2 (F 0 ) (32)
Eq. (32) shows that bigger is the parameter W Ḟ to set and bigger is the a priori uncertainty (variance) of the modelled unknown force F k at time t k . The uncertainty growth in time is also linear with t k (see Fig. 10). Let us now assuming that the a priori uncertainty in the future on the value of a rapidly varying unknown force is necessarily higher than for a slowly varying unknown force. In this case, the model given by Eq. (32) involves that the parameter W Ḟ should be set higher for a rapidly varying force than for a slowly varying force. This primary qualitative analysis on the choice of W Ḟ will be completed in Section 5.

Finally, it is necessary to take into account the fact that the measurement m x k of x k is obtained with a confo- cal chromatic sensor that adds to x k a discrete-time bandlimited white Gaussian noise v k with zero-mean and variance R:

σ 2 (F k ) σ 2 (F 0 ) t t 0 t 1 t 2 t 3 t 4 t k T s T s W F .
m x k = x k + v k (33)
Eq. ( 33) combined with Eqs. ( 16)-( 17) gives the state-space model with output m x k that will be used by the discrete Kalman filter:

X e k+1 = F X e k + Ω k (34) m x k = C X e k + v k (35)
with

E[Ω k ] = 0 3×1 E Ω k Ω T k = Q = W Ḟ η(T s ) (36) E[v k ] = 0 E[v 2 k ] = R (37)
4.2. Force estimation using a time-varying Kalman filter Fig. 11 shows the force estimation using a Kalman filter. The input of the Kalman filter is the noisy measurement m x k of the maglevtube displacement. The output is Fk . If the parameter W Ḟ is changed by the end-user during the force estimation process, a time-varying Kalman filter must be used and a numerical computation of Q must be done each time W Ḟ is changed (nevertheless, the term η(T s ) in Eq. ( 22) can be pre-calculated for a given sampling period T s [START_REF] Loan | Computing integrals involving the matrix exponential[END_REF]). The prediction-estimation stages of the Kalman filter are derived from Eqs. ( 16) and (35). The prediction stage is:

Xe k|k-1 = F Xe k-1 (38)
P k|k-1 = FP k-1 F T + Q (39)
with P k|k-1 the covariance of the prediction error conditioned by all previous measurements done before time t k :

P k|k-1 = E (X e k -Xe k|k-1 )(X e k -Xe k|k-1 ) T |M k-1 (40) M k-1 = {m x 1 , . . . , m x k-1 } (41)
In Eq. ( 38), the third component of Xe k|k-1 is the force prediction Fk|k-1 that remains equal to the last estimation:

Fk|k-1 = Fk-1 (42) The estimation stage that gives Fk is:

W F , f s , R, X 0 , P 0 , , Maglevtube F x (t) x(t) f s time-varying kalman filter F k x k v k m x k e .
K k = P k|k-1 C T CP k|k-1 C T + R -1 (43) Xe k = Xe k|k-1 + K k m x k -C Xe k|k-1 (44)
Fk = C F Xe k ( 45
)
P k = (I 3×3 -K k C) P k|k-1 (46)
with P k the conditional covariance of the estimation error defined by:

P k = E (X e k -Xe k )(X e k -Xe k ) T |M k (47)
M k = {m x 1 , . . . , m x k-1 , m x k } (48)
and C F the output matrix giving the estimated force:

C F = 0 0 1 (49)
The Kalman filter is a priori initialised for instance with the maglevtube in its equilibrium state when no force is applied to it:

Xe 0 = E X e 0 = 0 0 0 T (50)
The covariance matrix P 0 of the initial estimation error is taken equal to:

P 0 = E (X e 0 -Xe 0 )(X e 0 -Xe 0 ) T = σ 2 (x0) 0 0 0 σ 2 ( ẋ0) 0 0 0 σ 2 (F0) (51)
In Eq. ( 51), each variance on the main diagonal represents the a priori uncertainty on x 0 , ẋ0 and F 0 . The term σ 2 (F 0 ) is shown in Fig. 10. These values are chosen to be coherent with the initial conditions associated to the experiment made. In practice, they have little importance if the user starts the Kalman filter with no force applied on the maglevtube and waits a few seconds such that the Kalman gain K k converges to its steady-state K ∞ (solution to the discrete Riccati equation that depends on W Ḟ , f s and R) before applying an unknown varying external force. Each time W Ḟ is changed, K ∞ (W Ḟ , f s , R) evolves to a new value.

Analysis of simulated input forces

Estimation of a force generated by a Wiener Process

The study in simulation of the estimation behaviour is first conducted with an external force F x (t) generated by a Wiener process (like Eq. ( 10)) considered as a shaping filter instead of an uncertainty model:

Ḟ x (t) = ω r (t) (52)
This specific case corresponds to a brownian evolution of the input force F x (t). Its interest is to observe qualitatively what happens if the parameter W Ḟ is not equal to the PSD W of the zero-mean white Gaussian process ω r (t):

φ ωr,ωr (τ ) = W δ(τ ) ∀τ ∈ IR ( 53
)
The PSD W Ḟ can be set equal, inferior or superior to W by the end-user who does not know the real value of W . Sampling frequency f s is 100 Hertz. The variance R of the measurement noise is 1.44 × 10 -16 m 2 (typical value for a CL2 confocal head provided by STIL SA).

Fig. 12 shows F x (t) and Fk when W Ḟ is chosen equal to W and both set to 10 -14 N 2 /Hz. The maglevtube displacement that remains in the linear domain is also plotted for reference. In this case, the error estimation between the sampled input force and Fk has a minimal variance conditionally to the measurements done (classical framework of Kalman filtering). Fig. 12(a), a slight smoothing induced by the sampling period T s can be seen between Fk and the non-sampled force F (t). Fig. 13(a) shows F x (t) and Fk if W Ḟ is inferior to W (the last one is still equal to 10 -14 N 2 /Hz). This choice of W Ḟ leads to a covariance matrix Q that does not reflect the real covariance matrix associated to the sampled force F x (t). In this case, Q can be seen as a parameter that influence the estimation dynamic in some ways that will be studied in the next Sections. The estimation is then much more smoothed and delayed (and it is not optimal in the sense of minimal variance). In Fig. 13(b) W Ḟ is chosen superior to W and the estimation is more noisy. Thus a primary observation is that W Ḟ has an effect on the resulting force bandwidth (smoothing), on the delay and also on the noise present in the estimation.

Estimation of a step input force

Because the real force F x (t) does not follow generally the uncertain model given by Eq. ( 10), the behaviour of the deconvolution block must be characterised with criterions well adapted to SISO transfer functions. This characterisation will start by a temporal analysis of the effect of W Ḟ on the estimation of a canonical step input force F x (t). This analysis will result in the study of a tradeoff between the resolution and the bandwidth of the force sensor. The step amplitude is set to 100 nN in all simulations. To be independent of P 0 specified in Eq. ( 51), a steadystate Kalman filter is used substituting K ∞ (W Ḟ , f s , R) to K k and using only Eqs. (38) (44) (45). The associated third-order state-space model (with measurement m x k as input and estimation Fk as output) is obtained with Eq. (44) reported in Eq. (38) and using Eq. ( 45):

Xe k+1|k = A K Xe k|k-1 + B K m x k (54) Fk = C K Xe k|k-1 + D K m x k (55)
with

A K = F(I 3×3 -K ∞ C) B K = FK ∞ (56)
C K = C F (I 3×3 -K ∞ C) D K = C F K ∞ (57)
The variance R of the measurement noise remains equal to 1.44 × 10 -16 m 2 . The maglevtube parameters are m = 74 × 10 -6 kg, K x m = 0.02818 N/m (in linear domain), K x v = 1.8 × 10 -5 N.s/m (the damping coefficient ζ is 6.23 × 10 -3 ). Identified values used in the Kalman filter are K x m = 0.02812 N/m, K x v = 1.772 × 10 -5 N.s/m (ζ is 6.14 × 10 -3 ). The sampling frequency f s = 1/T s is deliberately set to a high value (f s = 1000 Hertz) compared to the -3 dB cutoff frequency of the maglevtube (around 4 Hertz) in order to see the effects of this choice. Fig. 14 shows the force Fk estimated for W Ḟ set to 10 -18 N 2 /Hz and 10 -15 N 2 /Hz. In the last case, the response time is around 0.1 second and thus much shorter than the Maglevtube settling time at 5% (20 seconds). Smaller is W Ḟ and smaller is the noise on Fk but longer is the estimation response time. As a consequence, smaller is the amplitude to estimate and smaller must be W Ḟ to have a good signal to noise ratio in Fk . But in this case, the force bandwidth of the sensor is also reduced.

Harmonic frequency response of the Kalman filter

The previous behaviour can be explained with the frequency response of the deconvolution block which is the steady-state Kalman filter:

H K (ω) = F (e jω ) m x (e jω ) (58) = C K (e jω I 3×3 -A K ) -1 B K + D K ( 59
)
As it is shown in Fig. 15, this frequency response "inverts" the frequency response of the maglevtube with its resonance peak (see Fig. 8). Fig. 15 shows that bigger is W Ḟ and bigger is the gain in the high frequencies just after the antiresonance. Thus bigger is the amplification of the high frequency components in the measurement noise v k included in the input m x k . Nevertheless, whatever is W Ḟ the gain in the highest frequencies below the Nyquist frequency f s /2 is always decreasing to limit the amplification of the highest frequencies in the input noise. To reduce the averall noise level, it is necessary to reduce W Ḟ . But in this case, the high frequency components of the displacement m x k have a very small amplitude because the maglevtube acts as a low pass filter (see Fig. 8). These high frequencies are then insufficiently amplified by the Kalman filter (see Fig. 15(a) with W Ḟ set to 10 -18 N 2 /Hz) to correctly reconstruct the high frequency components of the input F (t) at the step time (time t = 1 sec in Figs. 14(a) and 14(b)). As a consequence, the response time increases (and the bandwidth decreases). This behaviour is interesting because the maglevtube displacement becomes very small for high frequencies, especially if f s is set to a high value as it is deliberately the case here. As a consequence, the signal to noise ratio of the maglevtube displacement measurement is also bad in the high frequencies but the Kalman filter takes care of this information worsening by adjusting a trade-off (driven by W Ḟ ) between the resolution of the estimated force (ie the noise level in Fk ) and the response time of the estimation.

Force sensor resolution

The additive coloured Gaussian noise n k presents in the estimation Fk (see Fig. 14(b) for instance) is the consequence of the measurement noise v k (present in the input m x k ) going through the transfer function H K (z). Statistical properties of n k can be determined using the Kalman filter state-space model given by Eqs. ( 54)-( 55). As it is shown in Fig. 16, thanks to additivity of linear systems, only input v k is considered instead of the sum m x k = x k + v k , which gives the following stochastic state-space model:

Xe k+1|k = A K Xe k|k-1 + B K v k ( 60
)
n k = C K Xe k|k-1 + D K v k (61)
The noise v k is the input and the noise n k is the output. X e k|k-1 is a random state with mean m k and covariance matrix S k given by:

m k = E[ Xe k|k-1 ] = (A K ) k m 0 (62) S k+1 = E[( Xe k+1|k -m k+1 ) 2 ] = A K S k A K T + B K R B K T (63)
with R defined in Eq. ( 37). Taking into account the fact that v k is a zero-mean Gaussian noise and that Xe k|k-1 and v k are independent at time t k , the mean µ k and variance Σ k of the Gaussian noise n k in Eq. ( 61) are given by:

µ k = E[n k ] = C K m k + D K E[v k ] = C K m k (64) Σ k = E[(n k -µ k ) 2 ] = C K S k C K T + D K RD K T (65)
Because of the sum m x k = x k + v k on the input, any initial condition different from zero for the (deterministic) state X e k|k-1 will be associated to the input x k . Thus, by (a) Amplitude of H K (ω) (fs = 1000 Hz). additivity, it remains only a null initial condition for the (stochastic) state X e k|k-1 associated to the input v k :

X e 0|-1 = [0 0 0] T (66)
Because the initial condition for this stochastic state X e k|k-1

is always the previous one, it comes:

m 0 = 0 3×1 S 0 = 0 3×3 (67)
Notice that m 0 and S 0 only influence the transient evolution of n k at the beginning of the measurement. With Eq. ( 62) the mean m k is always null and with Eq. ( 64) the statistical properties of n k are given by:

µ k = 0 ∀k (68) Σ k = C K S k C K T + D K RD K T (69)
The 99% confidence interval associated to Eq. ( 69) is plotted with red dotted lines in Fig. 17. The sampling frequency f s is set to 100 Hertz in this simulation. The mean in this Fig. corresponds to the value of Fk when there is no measurement noise v k . This unknown value belongs to the confidence interval with a 0.99 probability.

To avoid any misinterpretation, it does not mean that the real force F x (t) belongs to this interval, but only that Fk should belong to it if they were no measurement noise. The associated standard deviation √ Σ k is an image of the resolution (or the SNR ratio) of the sensor. It tends to √ Σ ∞ in a few sampling times t k (see the first 0.05 sec in Fig. 17). This value √ Σ ∞ can be given to the end-user to adjust W Ḟ . Bigger is W Ḟ and bigger is √ Σ ∞ (see Fig. 18), thus lower is the resolution of the sensor. Variance Σ k has no connection to the force sensor accuracy. This last one only depends on the estimation of m, K x m and K x v during the calibration stage. In steady-state, the accuracy only depends on the calibration of K x m .

Force sensor bandwidth

The bandwidth of the force estimation depends on the harmonic transfer between the input force F x (t) and output estimation Fk . Because time description is continuous for the input and discrete for the output, this transfer is difficult to model. Thus, it will be determined with only discrete time description. Because of this, the first transfer from F x (t) to x(t) that represents the dynamic of the maglevtube is discretized with a zero-order hold with sampling period T s added on the input F x (t). This will leads to a slightly erroneous result when a pure sinusoidal input F x (t) is applied on the maglevtube.

The discretized linear state-space representation of the maglevtube is obtained from Eqs. ( 6) and ( 7): and the measurement m x k is given by Eqs. ( 33) and (71):

X k+1 = A d X k + B d F x k (70) x k = C d X k (71)
m x k = x k + v k (73) = C d X k + v k (74)
The global discrete state representing the maglevtube state and the steady-state Kalman filter state is:

X k = X k Xe k|k-1 (75)
To obtain the associated global state-space model, Eq. ( 74) is merged in Eqs. ( 54) and (55):

Xe k+1|k = A K Xe k|k-1 + B K C d X k + B K v k (76) Fk = C K Xe k|k-1 + D K C d X k + D K v k (77)
and Eqs. ( 75), ( 76), ( 77), (70) leads to the following global representation of the force sensor with F x k and v k as input and Fk as output:

X k+1 = A g X k + B g F x k v k (78) Fk = C g X k + D g F x k v k (79)
with

A g = A d 0 2×3 B K C d A K B g = B d 0 2×1 0 3×1 B K (80)
C g = D K C d C K D g = 0 D K (81)
Setting v k to zero in the global state-space model (to study only the impact of F x k input) gives:

X k+1 = A g X k + B d 0 3×1 F x k (82) Fk = C g X k (83)
and the global harmonic transfer H g (jω) between the force to measure and the force estimated

H g (ω) = F (e jω ) F x (e jω ) (84)
is given by: Sampling frequency f s is set to 1000 Hertz to minimise the effect of the ZOH in Eq. ( 70). The step force response (with v k set to zero) is also plotted to see the matching between response time and frequency response of the estimation. Higher is W Ḟ and larger is the force bandwidth of the sensor. As it is shown, the force bandwidth can be extended beyond the -3dB cutoff frequency of the maglevtube (around 4 Hz) and then the phase difference decrease. But in this case the resolution also decreases (see previous Section 5. [START_REF] Arai | Micro tri-axial force sensor for 3d bio-micromanipulation[END_REF]). This frequency response can be provided to the end-user to help him adjusting W Ḟ . The assumption at the end of Section 4.1 ruling that W Ḟ should be set higher for a rapidly varying force than for a slowly varying one is coherent with the way the bandwidth of the sensor is varying with W Ḟ .

H g (ω) = C g (e jω I 5×5 -A g ) -1 B d 0 3×1 (85)
Fig. 20 shows the global amplitude frequency response when the calibration is not perfectly accurate. The values chosen are those given in Section 5.2. The scale on Yaxis is linear contrary to Fig. 19(b). Force measurement around the resonant frequency of the maglevtube can lead to non negligible errors. Fig. 21 shows the five eigen values of A g for W Ḟ set to different values between 10 -15 N 2 /Hz and 10 -18 N 2 /Hz.

Because of the open-loop design used, the two eigen values very close to the unit circle (inside a small red circle) are the maglevtube poles that are compensated by the zeros associated to the Kalman filter. These two maglevtube poles represent its oscillating and under-damped dynamic with long response time. They do not change if W Ḟ is modified. The three other poles depend on the value of W Ḟ and also on the calibration accuracy. Higher is W Ḟ and closer are the poles to the unit circle right boundary, thus slower is the estimation. The two complex conjugate poles generate an oscillation with a slower dynamic than the real pole that is the fastest pole. The location of the poles does not lead to an instability of the Kalman filter for choice of W Ḟ even smaller than 10 -18 N 2 /Hz with double precision floating-point calculation.

Experimental results

Disturbance forces

The force Fk is an estimation of F x (t) that represents the sum of all external forces applied on the Maglevtube. Fig. 22 shows Fk when no external force is artificially applied on the Maglevtube. In this case Fk is an estimation of all the external disturbance forces applied on the maglevtube along x by the environmental noise (see Fig. 1). Because of the maglevtube -3dB cutoff frequency, the most predominant disturbance forces are due to the seismic low frequency vibrations of the sensor base (vibration of magnets M 2 that induces a vibration of the magnetic field) and to air vibration in infrasound frequencies that generates an acoustic force that is mainly applied on the rear deflector. Extreme care must be taken to reduce as much as possible this environmental noise. With the setup used (sensor enclosed in a chamber put on a antivibration table inside a small room shown in Fig. 23), this disturbance force is probably mainly generated by the antivibration table that insufficiently filters the seismic low frequency vibrations of the massive concrete ground floor. Maximum disturbance force amplitude is inferior to 1 nN and unbiased standard deviation is inferior to 0.3 nN. The presence of vibrations involves that the reference frame R 0 (attached to the antivibration table) is non Galilean and the fundamental principle of dynamic can not be applied in R 0 to obtain Eq. ( 5). Nevertheless, under the hypothesis that the displacement sensor and the magnets M 2 vibrates in phase (realistic hypothesis in the low frequency domain), it can be shown using a new Galilean reference frame that:

m ẍ = F x + F x seismic -K x m x -K x v ẋ (86) m ẍ = F x -K x m x -K x v ẋ (87)
with

F x = F x + F x seismic (88)
The force F x seismic gives the low frequency seismic excitation of the maglevtube that depends on the sensor base low frequency vibration. The force F x represents all the external forces applied on the maglevtube except the seismic excitation. Thus Eq. ( 87) is the same than Eq. ( 5). The demonstration of this point and the taking into account of the seismic vibrations (measurement, compensation) is outside the scope of this paper.

Pull-off force measurement

Fig. 25 shows the evolution of the force Fk during a pull-off force measurement. A planar material is pushed

Maglevtube tip

Micro-sphere Planar material against a micro-sphere (see Fig. 24) stuck at the maglevtube tip (loading stage) and then moved back (unloading stage) until the contact is broken between the material and the tip. After this contact loss, the seismic mass is oscillating because of its big inertia but the unknown external force applied on the maglevtube becomes known: it is equal to zero and thus it can be compared with the force estimated with Eq. ( 3) or Eq. ( 45). Eq. ( 3) gives a bad estimation because it is proportional to the displacement of the ZIR. Kalman estimation associated to Eq. ( 45) gives a better result with a shorter and smaller oscillating transient response. The residual oscillating behaviour can be the result of several causes:

• the identification inaccuracy for the viscous damping coefficient K x v ,

• the push is not perfectly done in the horizontal direction thus the maglevtube is not aligned with the x axis when the contact is broken. The orientation of the maglevtube around y axis at this instant creates a complex behaviour in the plane xoz that is not taken into account in the one DOF model developed in Section 3.1, thus the Kalman filter does not react correctly,

• the horizontal and vertical seismic excitations of the maglevtube that also create a complex behaviour in the plane xoz defined in Fig. 3 (see previous example).

Fig. 26 also shows a pull-off force measurement but with a non-contact attractive force also applied on the maglevtube sphere before the loading, during the loading and after the contact break. The sphere used is an electrically charged macroscopic sphere (one millimetre diameter, Silicon nitride Si 3 N 4 ). The planar material is a cover glass (insulating material), thus the most predominant non-contact force is probably an electrostatic force because atmosphere relative humidity was low (no meniscus was observed during the pull-off measurement). The oscillations after the contact break represent the evolution of the non-contact force applied on the sphere when the cover glass is moved back. As the sphere moves towards and backwards the cover glass (because of the maglevtube inertia) the non-contact force is oscillating. The mean of these oscillations has the same mirroring shape than the non-contact force before the loading. Because the speed of the planar material is opposite in both cases, this probably means that the distribution of charges in both material remained the same during all the force measurement.

Conclusion

Despite the inertia of their mass that should be a great handicap for time-varying force measurement, nanoforce sensors based on a rigid macroscopic seismic mass are not disqualified compared to classical designs based on microscopic elastic cantilevers like AFMs. The design proposed in this paper is based on diamagnetic levitation and makes possible the adjustment of the stiffness of the magnetic springs. Stiffness as low as 0.005 N/m can be reached. The linearity of the dynamic behaviour is good for millimetric displacements of the seismic mass, which makes possible the measurement of a large range of force with the same seismic mass. Calibration is quite easy to perform because the mass can be measured with a precision balance. Nevertheless, this design imposes some constraints like the limited bandwidth and the deconvolution of the mass displacement to take care of the behaviour imposed by the mass inertia. This deconvolution is done with a discrete Kalman filter that is using a discretized Wiener process to model the unknown input force uncertainty. This processing requires the adjustment of a single parameter W Ḟ which directly adjusts a trade-off between the resolution (standard deviation) of Fk and the bandwidth of the force sensor. This parameter can be modified at any time by the end-user in accordance with its own knowledge on the force to measure. The resolution and frequency response of the sensor can also be plotted to make this choice easier. This deconvolution method is computationally cheap (third-order IIR filter if a steady-state Kalman filter is chosen) and can be implemented in cheap DSP or microcontrollers. Response time shorter than 0.1 second can be reached with a correct S/N ratio for forces in the micronewton scale despite the very long settling time of the transducer (20 seconds) and its under-damping. Compared to a simple low-pass filter added on the displacement measurement that will inevitably restrict the bandwidth, the force bandwidth can be extended reasonably four times higher than the displacement bandwidth. The main drawback of this open-loop design is its sensitivity to unwanted external disturbant forces like whose generated by the vibrations of the magnetic field and the vibration of the air in infrasound frequencies. Extreme care must be taken to reduce as much as possible this environmental noise thanks to an antivibration table and a closed chamber. As an outlook, a closed-loop design should be a alternative way to take into account these disturbances in order to cancel them.

Acknowledgment

This work was supported by the french National Research Agency under STILµFORCE ANR-07-ROBO-0005 contract.

Glossary of symbols

Order of appearance of mathematic symbols is the same than in the paper.

K x

m magnetic stiffness of the force sensor

K x v viscous

Figure 2 :

 2 Figure 2: Macroscopic seismic mass sensitive to external forces.

Figure 3 :

 3 Figure 3: Levitating seismic mass in the force sensing device.

Figure 4 :

 4 Figure 4: Simulated step response of the maglevtube with an external input force set to one micro newton.

Figure 5 :

 5 Figure 5: Simulated steady-state characteristic of the external force F x versus the displacement x of the maglevtube.

Figure 6 :

 6 Figure 6: Force sensor prototype.

Figure 7 :

 7 Figure 7: Measured and reconstructed zero input response (ZIR) of the maglevtube displacement.

Figure 8 :

 8 Figure 8: Simulated frequency response along x of the identified linear model.

Figure 9 :

 9 Figure 9: Blocks representation of the uncertain modelling of the input force and the deterministic model of the maglevtube.

Figure 10 :

 10 Figure 10: W Ḟ effect on the a priori variance σ 2 (F k ).

Figure 11 :

 11 Figure 11: Force estimation with a time-varying Kalman filter.

  real force F x (t) estimated force Fk (a) Real and estimated force.

  real displacement x(t) measured displ. m x (t) (b) Associated measured displacement.

Figure 12 :

 12 Figure 12: Force estimation if W = W Ḟ = 10 -14 N 2 /Hz.

Figure 13 :

 13 Figure 13: Real and estimated force if W Ḟ is different from W equal to 10 -14 N 2 /Hz.

  real force F x (t) Fk with W Ḟ = 10 -18 (a) Real and estimated force with W Ḟ = 10 -18 N 2 /Hz. real force F x (t) Fk with W Ḟ = 10 -15 (b) Real and estimated force with W Ḟ = 10 -15 N 2 /Hz. Real and Measured displacement with fs = 1000 Hz.

Figure 14 :

 14 Figure14:Step force estimation with F x set to 100 nN.

W

  Ḟ = 10 -18 W Ḟ = 10 -15 (b) Phase of H K (ω) (fs = 1000 Hz).

Figure 15 :

 15 Figure 15: Impact of W Ḟ on the transfer function F (e jω )/m x (e jω ).

Figure 16 :

 16 Figure 16: Output noise n k in estimation Fk due to the measurement noise v k .

Figure 17 :

 17 Figure 17: 99% confidence interval (±3 √ Σ k ) on the force estimation with W Ḟ = 10 -15 N 2 /Hz and fs = 100 Hz.

Figure 18 :

 18 Figure 18: Standard deviation √ Σ∞ versus W Ḟ for fs = 100 Hertz and R = 1.44 × 10 -16 m 2 . Small irregularities are mainly due to numerical errors during the solving of discrete Riccati to obtain K∞(W Ḟ , fs, R).

e

  withA d = e ATs B d = Ts 0 At dt B C d = C (72)

Fig 19

 19 Fig 19 shows the global frequency response of the sensor for different values of W Ḟ when the calibration is perfectly accurate (matrices A d and B d in Eq. (70) are replaced with the identified matrices given by the calibration stage).

W

  Ḟ = 10 -18 W Ḟ = 10 -15 (a)Step force response with noise v k set to zero.

W

  Ḟ = 10 -18 W Ḟ = 10 -15(b) Amplitude of Hg(e jω ).

W

  Ḟ = 10 -18 W Ḟ = 10 -15 (c) Phase of Hg(e jω ). 19: Impact of W Ḟ on the global frequency response Hg(e jω ) and the response time of the sensor (fs = 1000 Hertz -accurate calibration).

  scale) W Ḟ = 10 -18 W Ḟ = 10 -15

Figure 20 :Figure 21 :Figure 22 :

 202122 Figure 20: Global amplitude frequency response Hg(e jω ) with used calibration.

Figure 23 :

 23 Figure 23: force sensor enclosed in a chamber on a antivibration table.

Figure 24 :

 24 Figure 24: Micro-sphere stuck at the maglevtube tip.

Figure 25 :

 25 Figure 25: Experimental pull-off force measurement.

Figure 26 :

 26 Figure 26: Experimental pull-off force measurement with noncontact attractive force also applied on the Maglevtube.

  damping coefficient of the maglevtube (t) real external force applied to the maglevtube tip F (t) modelled external force applied to the maglevtube tipF k discretemodelled external force applied to the maglevtube tip ω(t) zero-mean white Gaussian stochastic process used to model uncertainty on Ḟ (t) W Ḟ power spectral density of stochastic process ω(t) X(t) continuous state of the maglevtube X e (t) extended continuous state that is including F (t) x k maglevtube discrete displacement if the external force (real or modelled) is discretized with a zero-order hold X k discrete state of the maglevtube e k by the Kalman filter M k set of measurements done between times t 1 and t k included P k|k-1 covariance matrix of prediction error X e k -Xe k|k-1 conditioned by measurements done before t k P k covariance matrix of estimation error X e k -Xe k conditioned by measurements done until t k K k time-varying Kalman gain K ∞ steady-state Kalman gain that depends on W Ḟ , f s and R (t) zero-mean white Gaussian stochastic process used to generate simulated derivative Ḟ x (t) W power spectral density of ω r (t) H K SISO transfer function (z-domain) of the steady-state Kalman filter with m x k as input and Fk as output m k mean of Xe k|k-1 when the steady-state Kalman filter is only excitated by noise v k S k covariance matrix of Xe k|k-1 when the steady-state Kalman filter is only excited by noise v k µ k mean of output noise n k present in Fk Σ k covariance of output noise n k present in Fk F x k F x (t) applied to a zero-order hold and sampled with frequency f s X k global state of the force sensor with F x k and v k as input and Fk as output H g SISO transfer function (z-domain) of the force sensor with F x k as input and Fk as output System, control, output and feed-forward matrices associated to state and output dynamics: K , C K , D K modelling of Xe k|k-1 and Fk dynamic A d , B d , C d modelling of X k and x k dynamic A g , B g , C g , D g modelling of X k and Fk dynamic

	m x(t) maglevtube displacement maglevtube mass A, B, C modelling of X(t) and x(t) dynamic A, M, C modelling of X e (t) and x(t) dynamic F, C modelling of X e k and m x k dynamic F x X e k extended discrete state including F k Ω k process noise associated to X e k dynamic Q covariance matrix of process noise Ω k T s sampling period used in the force sensor f s sampling frequency used in the force sensor v k scalar zero-mean white Gaussian measurement noise m x k noisy measurement of maglevtube displacement Xe k|k-1 prediction of X e k by the Kalman filter Xe Fk estimated force by the Kalman filter Xe 0 initial a priori estimation of X e k P 0 initial a priori covariance matrix of estimation error X e 0 -Xe R variance of noise v k C F output matrix giving Fk from Xe

k estimation of X 0 ω r k A K , B


Guillaume Olivier

Gaelle Penelet

Joël Poignand

Pierrick Gilbert

Lotton

C Olivier

G Poignand

J Gilbert

P Lotton

Weakly nonlinear propagation in thermoacoustic engines: a numerical study of higher harmonics generation up to the appearance of shock waves Come

à 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

Thermoacoustic engines (TAEs) are thermodynamic systems that make use of a temperature gradient along a porous material to generate self-sustained acoustic oscillations. This auto-oscillation occurs beyond a threshold temperature difference for which the thermoacoustic amplification taking place in the thermoacoustic core compensates exactly all losses in the system. This onset condition is well described by the linear theory of thermoacoustics [START_REF] Rott | Thermoacoustics[END_REF]. However, above the onset threshold, the amplitude of the acoustic oscillation grows to reach high levels. The saturation of the growth of amplitude is then controlled by nonlinear phenomena which are the cause for both acoustic and thermal power dissipation. Among the many effects usually described are the acoustically enhanced thermal pumping [START_REF] Rott | Thermoacoustics[END_REF][START_REF] Yazaki | Large heat transport due to spontaneous gas oscillation induced in a tube with steep temperature gradients[END_REF], acoustically induced DC flows [START_REF] Desjouy | Measurement of acoustic streaming in a closed-loop traveling wave resonator using laser doppler velocimetry[END_REF][START_REF] Thompson | Simultaneous measurement of acoustic and streaming velocities in a standing wave using laser doppler anemometry[END_REF][START_REF] Moreau | Measurements of inner and outer streaming vortices in a standing waveguide using laser doppler velocimetry[END_REF], complex aerodynamical effects due to geometrical singularities at interfaces between elements of different porosity or at geometric discontinuities [START_REF] Swift | Thermoacoustics: a unifying perspective for some engines and refrigerators[END_REF], or nonlinear propagation leading to higher harmonics generation [START_REF] Yuan | A simplified model for linear and nonlinear processes in thermoacoustic prime movers. part ii. nonlinear oscillations[END_REF][START_REF] Gusev | Asymptotic theory of nonlinear acoustic waves in a thermoacoustic primemover[END_REF].

Thermoacoustic pumping is well described analytically for simple pore shapes [START_REF] Arnott | General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross sections[END_REF], or with hybrid analytic/experimental approaches for more complex materials [START_REF] Guedra | Theoretical prediction of the onset of thermoacoustic instability from the experimental transfer matrix of a thermoacoustic core[END_REF][START_REF] Bannwart | Measurements of the impedance matrix of a thermoacoustic core: Applications to the design of thermoacoustic engines[END_REF]. Acoustic streaming is accounted for in models by higher-order developments of the basic equations [START_REF] Gusev | Acoustic streaming in annular thermoacoustic prime-movers[END_REF][START_REF] Bailliet | Acoustic streaming in closed thermoacoustic devices[END_REF], and experimental solutions are implemented to limit its effects [START_REF] Backhaus | A thermoacoustic-Stirling heat engine: Detailed study[END_REF][START_REF] Tijani | A high performance thermoacoustic engine[END_REF]. The aerodynamical effects at geometrical singularities have been extensively studied numerically [START_REF] Marx | Numerical simulation of stackheat exchangers coupling in a thermoacoustic refrigerator[END_REF] and experimentally [START_REF] Wakeland | Measurements of the resistance of parallel-plate heat exchangers to oscillating flow at high amplitudes[END_REF][START_REF] Berson | Nonperiodicity of the flow within the gap of a thermoacoustic couple at high amplitudes (L)[END_REF] in thermoacoustic systems, but are often analytically disregarded because of the complexity of the involved phenomena. Some of these effects may be minimized empirically by the use of artifacts such as flow straighteners [START_REF] Backhaus | A thermoacoustic-Stirling heat engine: Detailed study[END_REF].

Nonlinear propagation due to the high acoustic amplitudes can be dealt with by a proper design of resonators [START_REF] Ilinskii | Nonlinear standing waves in an acoustical resonator[END_REF]. It is usually ignored in the description of thermoacoustic engines for several reasons. Firstly, when present, the magnitude of the losses due to nonlinear propagation is small compared to the one produced by other nonlinear phenomena. Also, most engines have complex geometries making their resonances inharmonic [START_REF] Gaitan | Finite amplitude standing waves in harmonic and anharmonic tubes[END_REF]. Even when built with a simpler shape of resonator such as a straight tube, the presence of the steep temperature gradient tends to make the frequencies of the higher modes non-integer multiples of the fundamental resonant frequency f 0 of the engine. Therefore, the formation of shock waves or a strong deformation of the wavefront is rarely observed in thermoacoustic engines, although the acoustic pressure oscillations may reach an amplitude of 10% of the static pressure [START_REF] Swift | Analysis and performance of a large thermoacoustic engine[END_REF][START_REF] Olson | A loaded thermoacoustic engine[END_REF]. Nevertheless, it has been reported by Biwa et al. [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF] that is it possible to obtain an experimental traveling shock wave in an annular thermoacoustic engine under specific conditions. Furthermore, the same team showed more recently that shock wave formation is also possible in a standing wave engine [START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]. This paper proposes a numerical investigation of the nonlinear behavior of simple standing or traveling wave thermoacoustic engines such as the ones presented schematically in Fig. 1. In Section 2, the analytic tools used in this study are briefly presented, for the description of the propagation in the thermoacoustic core with the linear equations of thermoacoustics [START_REF] Rott | Thermoacoustics[END_REF] and for the nonlinear propagation in the waveguides with generalized Burgers equation in the form introduced by Sugimoto [START_REF] Sugimoto | Nonlinear Waves Motion, ch. Generalized Burgers equations and fractional calculus[END_REF]. A numerical solving process of these equations of propagation is proposed, based on the one used by Menguy and Gilbert [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF] for the case of forced oscillations and adapted

l h l s l h l w T h T c x x = 0 (A) (B) x = L x = 0 or L x 1 x 2 R x s
Figure 1: Architectures of the thermoacoustic engines under consideration in this study: (A) Standing-wave engine made of a thermoacoustic core (TAC) in a resonant waveguide; (B) Traveling-wave engine made of a similar TAC in a looped tube. On top is the temperature distribution imposed in both engines, assuming linear variations. The TAC consists of two heat exchangers (striped), a stack (hatched) and the part of the resonator comprising a temperature gradient, referred to as the thermal buffer tube (TBT).

to the case of self-sustained oscillations of undetermined frequency. Section 3 is devoted to presenting numerical results for different engine configurations. We show that higher harmonics generation can be favored or inhibited for various parameters of the thermoacoustic core (geometry, temperature distribution T 0 (x), position in the resonator x s /L, etc.). Computed steady-state waveforms are presented that show the possibility of shock wave formation for both standing and traveling wave configurations under particular conditions. These results are discussed in relation with available experimental results [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF][START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]. A brief discussion is also given to interpret these results. The conclusion presented in Section 4 highlights the main results and proposes perspectives to this study.

Simplified approach for the description of TAEs

A basic treatment of the nonlinear propagation in thermoacoustic engines (TAEs) suggests a separation in two sub-problems: on one hand, the simplified description of the nonlinear propagation in the temperaturehomogeneous waveguide and, on the other hand, the description of the propagation in the idealized thermoacoustic core comprising an inhomogeneous temperature distribution. Analytic tools are available for both problems: nonlinear propagation in waveguides is expressed making use of generalized Burgers' equations [START_REF] Nayfeh | Nonlinear oscillations[END_REF][START_REF] Hamilton | Nonlinear acoustics[END_REF], and the propagation in the thermoacoustic core (TAC) is given by Rott's linear theory of thermoacoustics [START_REF] Rott | Thermoacoustics[END_REF]. The coupling of the subsystems and the research of a solution with an iterative harmonic balance method allows to compute the amplitudes for all harmonics of the periodic steadystate regime, for a given set of geometrical parameters, thermophysical properties and an assigned temperature distribution through the thermoacoustic core.

Nonlinear propagation in waveguides

The nonlinear propagation of waves in ducts has since long been described, with application to various fields [START_REF] Chester | Resonant oscillations in closed tubes[END_REF][START_REF] Sugimoto | Burgers equation with a fractional derivative; hereditary effects on nonlinear acoustic waves[END_REF][START_REF] Chester | Resonant oscillations of a gas in an open-ended tube[END_REF]. In the context of this study, the work of Menguy and Gilbert [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF] is worth considering. Their use of the formalism of generalized Burgers' equations to describe the lossy, weakly nonlinear acoustic propagation allowed them to emphasize the shock wave formation in wind instruments. This phenomenon, previously observed experimentally by Hirschberg et al. [START_REF] Hirschberg | Shock waves in trombones[END_REF] is now known to explain the brassiness of their tone [START_REF] Gilbert | A simulation tool for brassiness studies[END_REF]. The same approach is used in this study to describe the cascade process of higher harmonics generation, applied to the thermoacoustic auto-oscillation whose frequency is likely to drift with the variations of amplitude of the acoustic variables. In a fluid characterized by its shear dynamic and kinematic viscosities µ and ν, its bulk viscosity η, its thermal conductivity κ and its specific heat capacity c P (all five parameters being assumed constant in absence of a temperature gradient), the acoustic plane wave of angular frequency ω propagating in a resonator of radius R at the typical value c 0 of the sound speed in that fluid may be characterized by three numbers:

• M = u 0 /c 0 , the acoustic Mach number, which is the ratio of the peak amplitude of the particle velocity to the speed of sound;

• Re = u 0 / √ νω, an acoustic Reynolds number [START_REF] Merkli | Transition to turbulence in oscillating pipe flow[END_REF], which is the ratio of inertia to viscous forces and qualifies the type of flow in the boundary layer thickness δ ν = ν/ω;

• Sh = ν/ω/R, the shear number, which is the ratio of the viscous boundary layer thickness to the characteristic transverse dimension of the waveguide.

Thermoacoustic engines usually operate at frequencies from tens to a few hundreds of hertz, with waveguides of a few centimeters in radius, and high sound levels (145-175 dB SPL with a working fluid at atmospheric pressure, that is drive ratios -ratio of acoustic pressure amplitude |p| to the static pressure P 0 -from 0.5% up to 10%). Under these conditions, the three dimensionless numbers M , 1/Re and Sh are small compared to unity, indicating that the nonlinear effects of propagation are not local but cumulative [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF], with two consequences. Firstly, any acoustic wave may be considered as the superposition of two noninteracting counterpropagating waves, whose propagation is described by two independent equations. Secondly, this allows to separate two scales of time and space for the description of the propagation, one fast scale (resp. short) for the quasi-linear local acoustic propagation, and one slow scale (resp. long) for the cumulative effects. Furthermore, the comparison of the relative amplitudes of the dimensionless numbers allows to evaluate the importance of the different phenomena taking place during propagation. A diagram of predominance is given in Fig. 2, allowing to visualize the domain of validity of the chosen filled with air at ambient pressure. Each straight line represents the limit of predominance of one phenomenon over another: Merkli condition is the upper limit for laminar flow in boundary layers; above M 2 = Sh, second order effects take over boundary layer effects; M = Sh 2 and M = Re 2 Sh are two limits under which a nonlinear description of the propagation is not necessary. Parallel dashed lines represent the same conditions in a pore of the stack. Adapted from Fig. 1 of [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF]. approach as a function of frequency and drive ratio. Conditions M = Sh 2 and M = Re 2 Sh represent the limits under which nonlinear phenomena are small compared to the neglected linear effects (e.g. the bulk losses or the effects of the curvature of the waveguide), and a linear description of acoustic propagation would suffice. Condition M 2 = Sh stands as the limit above which boundary layer effects are to be ignored compared to second order effects, better described by another approach than the one used here. Merkli condition Re = u 0 / √ νω < 400 is the limit under which the flow in the boundary layers is laminar, a necessary condition for the expression of the boundary layer losses used thereafter. The cut-off frequency of the waveguide for plane waves is given as the vertical dotted line for information. The description given here is that of a one-dimensional system (only plane waves propagate and the amplification from the thermoacoustic core is 1-D) and therefore does not account for transverse modes, though the frequency of the higher harmonics considered can be above that cut-off frequency.

These limits border a domain, represented by the shaded area in Fig. 2, which comprises the condition of nominal operation of the two studied thermoacoustic engines. The acoustic propagation in the resonators may therefore be described by the two independent generalized Burgers' equations [START_REF] Sugimoto | Nonlinear Waves Motion, ch. Generalized Burgers equations and fractional calculus[END_REF]

∂ σ p ± ∓ p ± ∂ θ± p ± = ± St ǫ ∂ 2 θ± p ± ∓ T ǫ ∂ 1 2 θ± p ± ( 1
)
where σ is a long spatial scale and θ a delayed time scale [START_REF] Rudenko | Theoretical foundations of nonlinear acoustics[END_REF]. Subscript + (resp. -) stands for propagation in the positive (resp. negative )

x direction. ǫ = M γ+1
2 is an evaluation of the importance of the nonlinear effects.

Stokes' number St = M 2 2Re 2 4 3 + η ν γ-1
Pr , where Pr = µcP κ , stands for the importance of bulk losses, that are proportional to the square of the frequency (as is emphasized by the double derivative in the time domain). This last term is negligible for small deformations of the wavefront, but cannot be ignored any more once the shock is formed. The term with the fractional derivative stands for the thermoviscous losses taking place in the boundary layer, whose relative importance is designated by

T = Sh 1 + γ-1 √ Pr
and proportional to the square root of the frequency. The 1 2 fractional derivative is defined as the convolution product

∂ 1 2 θ p = ∂ θ p * 1 √ πθ = 1 √ π θ -∞ 1 √ θ -θ ′ ∂ θ ′ p(θ ′ , σ)dθ ′ . ( 2
)

Propagation through the TAC

By thermoacoustic core (TAC), we designate the parts of the engine between x = x 1 and x = x 2 (see Fig. 1) comprising an inhomogeneous temperature distribution and consisting of an ambient heat exchanger (AHX), a porous medium designated as the stack, submitted to a temperature gradient by a hot heat exchanger (HHX) and a thermal buffer tube (TBT) which provides a thermal buffer between the HHX and the ambient temperature waveguide. On the diagram of predominance (Fig. 2), the limits of predominance are also represented in dashed lines for the propagation of an acoustic wave in a pore of the stack (or of one of the heat exchangers, of equivalent hydraulic radius). The nominal working conditions of the engine are below the validity domain of the approach presented above, represented by the dotted area. More particularly, the working conditions are below the M = Sh 2 condition, indicating that a nonlinear description of the propagation in the stack is not necessary as long as the predominant boundary layer effects are well accounted for. Moreover, in the engine core, the presence of a strong temperature gradient T 0 (x) does not allow the decoupling of the acoustic propagation into two counterpropagating waves and the temperature dependence of the thermoviscous properties of the fluid must be accounted for. Therefore, the TAC will be described with a different approach than the rest of the engine. In the frequency domain, it can be described by a scattering matrix [START_REF] Penelet | Dependence of sound amplification on temperature distribution in annular thermoacoustic engines[END_REF] obtained from the classical equations of thermoacoustics [START_REF] Rott | Thermoacoustics[END_REF], and expressing the counterpropagating components for each harmonic p

(n) ± of frequency nf 0 , (n ∈ [[1, ∞[[) of
the acoustic pressure variations at one extremity of the core (x = x 2 ) as a function of the same variables at the other end (x = x 1 ), with help of the reflection and transmission coefficients R ± and T ± of the core

p(n) + (x 2 ) p(n) -(x 1 ) = S (nf 0 , T 0 (x)) p(n) + (x 1 ) p(n) -(x 2 ) = T + R - R + T - p(n) + (x 1 ) p(n) -(x 2 ) . ( 3
)
Coupling this description of the acoustic propagation through the thermoacoustic core together with the ex-pression of the nonlinear propagation in the temperature homogeneous waveguides is sufficient to obtain the saturated waveforms of a thermoacoustic engine. However, as stated in Section 1, various phenomena of heat and mass transfer work together to dissipate energy in TAEs.

Ignoring the additional losses generated by these phenomena would yield numerical amplitudes obtained with the current linear description of the thermoacoustic core beyond the observed experimental saturation levels of the acoustic amplitudes. It is therefore necessary to include some nonlinear effects in the description of the TAC for the solution to converge to realistic amplitudes of saturation. Among the previously cited effects, sudden changes of the effective cross-section area at porous elements interfaces are the source of localized dissipation of acoustic energy that results in a pressure drop ∆p through the TAC. These so-called minor losses may be evaluated roughly by extrapolating to oscillatory flows [START_REF] Penelet | Experimental and theoretical study of processes leading to steadystate sound in annular thermoacoustic engines[END_REF] the DC behavior of such interfaces, yielding an average pressure drop for each harmonic

∆p (n) (t) ≃ - 4 3π ρ 0 2 (K 12 φ 1 + K 21 φ 2 )|u (n) | 2 (t), ( 4
)
where ρ 0 is the density of the fluid at the temperature of the interface, φ 1 and φ 2 the porosity of the elements on each side of an interface, and K 12 and K 21 are tabulated losses coefficients for a sudden section expansion or contraction [START_REF] Idelchik | Handbook of Hydraulic Resistance[END_REF]. The acoustic velocity u is evaluated from Euler's equation:

u(t) ≃ (p + (t) -p -(t))/(ρ 0 c 0 ).
Taking into account those losses proportional to u 2 yields the nonlinear scattering matrix for the thermoacoustic core

p(n) + (x 2 ) p(n) -(x 1 ) = S nl nf 0 , T 0 , |ũ (n) | 2 p(n) + (x 1 ) p(n) -(x 2 ) . (5)
One should note that such a description is not intended to be accurate, but should give an order of magnitude of the encountered losses in the TAC, allowing realistic levels of saturation. Neglecting the other nonlinear saturating effects is a choice to keep the description simple. It can however be justified by the following reasoning: aerodynamical edge effects, acoustic streaming and thermoacoustic pumping reduce the efficiency of thermal-to-acoustic conversion by causing heat transfer via convection and/or conduction. Thus, their effect would merely translate into a variation of the temperature of the hot heat exchanger. However, the description adopted in the frame of this paper concentrates on steady-state regime and therefore the imposed temperature of the hot heat exchanger can be considered as the resulting one, accounting for these effects.

Solving process

Though approximate solutions are available [START_REF] Bednarik | Propagation of quasiplane nonlinear waves in tubes and the approximate solutions of the generalized burgers equation[END_REF], equations (1) describing the nonlinear propagation of quasiplane waves in a thermoviscous fluid do not have any exact analytic solution. They are solved numerically in the frequency domain for the different sections of waveguide [0; x 1 ] and [x 2 ; L], coupled by the scattering matrix (5), and with boundary conditions depending on the geometry of the engine. The solution is sought in the form of a harmonic series

p = +∞ n=1 [a n (σ) sin nθ + b n (σ) cos nθ] , (6)
the introduction of which in the lossy Burgers' equation ( 1) and term-by-term equalization yields a set of equations for the coefficients a n and b n . This set is solved numerically following the steps stated here after, which consist in an extension to self-sustained oscillations of the work of Menguy and Gilbert [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF] for forced oscillations.

• (i) The initial harmonic spectrum is defined with a fundamental frequency close to the resonant frequency of the engine (e.g., at the resonance frequency of the empty resonator) with an amplitude of the fundamental of a few hundreds of Pa and 0 for all harmonics, at one extremity of the engine.

• (ii) Each harmonic of the signal is propagated through the whole engine back and forth to the starting point by applying a backwards finite differences scheme on equations (1) in the resonator parts, and thanks to the scattering matrix (5) through the TAC.

• (iii) After one complete cycle, the amplitudes of each harmonic of the final spectrum are compared elementwise to the amplitudes of the initial spectrum. In the same manner, the phase of the fundamental component is compared between the initial φ i and final states φ f .

• (iv) If the phase difference between initial and final states of the fundamental is not zero (i.e, if φ f -φ i = 0 mod [2nπ], with n depending on the working mode), the frequency of the fundamental is adjusted. Due to the nonlinear drift of the resonant frequency with the amplitude of the acoustic pressure (sometimes called nonlinear detuning [START_REF] Gusev | Asymptotic theory of nonlinear acoustic waves in a thermoacoustic primemover[END_REF]), this correction has to be repeated at each iteration of the solving process, contrarily to the case of forced oscillations where the frequency is imposed by a driving source.

• (v) If the difference e of amplitudes (in the sense of the absolute deviations) between initial and final states is smaller than a predefined small threshold ε, the obtained solution is considered to be the steady-state saturated spectrum and the search of a solution stops here. If e > ε, the steady-state has not been reached yet. The spectrum obtained after this iteration is used to initiate the next iteration and steps (ii) through (v) are repeated until a satisfying solution is found (see Ref. [START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF] for a more detailed description of the procedure in the case of forced oscillations).

The definition of the system is such here that it tends to the point where losses compensate the thermoacoustic amplification for all harmonics after a sufficient number of iterations, whatever the initial state provided. In practice, in order to reduce the computation time, the number of computed harmonics has to be limited to a reasonable number. However, since the energy content of each harmonic decreases with its rank due to volumetric losses proportional to ω 2 , the truncation error drops quickly when the number of harmonics is increased over a few tens.

Results and discussion

Two classical engines are modeled for this study. The configurations are chosen to describe actual engines for which experimental results showing shock wave formation are available [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF][START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]. The first one (A) is a standing wave engine closed at both ends [START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]. The second one (B) is a quasi-traveling wave engine with a looped resonator, using a core of similar construction, but with a stack with slightly different specifications [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF]. The dimensions of both engines are provided in Table 1. The length l w of the TBT is not defined in either reference, because no secondary heat exchanger is used to constrain the cold temperature in the resonators. In the quasi-traveling wave engine (B), the amplification of the TAC can be highly sensitive to the temperature distribution [START_REF] Penelet | Dependence of sound amplification on temperature distribution in annular thermoacoustic engines[END_REF]. Assuming a linear temperature distribution T 0 (x), l w is then defined to match the given experimental onset threshold of the experimental apparatus. In the standing wave engine however, the sensitivity to the temperature distribution appears to be much less important. The TBT is chosen to be short enough to allow the placement of the core close to the extremity of the engine. It is to be noted that, though the data provided in Ref. [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF][START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF] is admittedly sufficient for the description of the linear behavior of the engines (i.e., their onset threshold), it is not enough for a precise rendering of their operation above threshold. Some parameters are not provided in the description of the experimental configurations, either because they were not relevant in the frame of the presented work or because they were not available to the authors (e.g., the temperature distribution in the whole TAC); these parameters are nonetheless important to describe theoretically the steady-state operation above threshold. These above-mentioned experimental configurations are however used as reference for this study, still allowing for a qualitative comparison. The description of the nonlinear effects controlling the amplitude of saturation in the model presented in this paper provides only an estimate of their amplitude, and therefore does not allow for a quantitative comparison.

The computed periodic steady state waveform in a thermoacoustic engine described with the previous tools is the result of the balance of the phenomena taking place in the two sections of the engine. On one hand, the nonlinear propagation of high amplitude acoustic waves induces the cascade of energy from the fundamental frequency of the periodic oscillation regime to its higher harmonics. This phenomenon prevails in the temperature homogeneous resonant parts of the engine. With the generation of higher harmonics also comes additional thermal and viscous dissipation, increasing with frequency. These losses are the main source of saturation of the amplitude, together with the minor losses in the TAC. On the other hand, from the linear theory of thermoacoustics, a thermoacoustic core with a given set of parameters (geometry, temperature distribution T 0 (x), position in the resonator x s /L, etc.) can be seen either as a source or a sink of volume flow rate [START_REF] Swift | Thermoacoustics: a unifying perspective for some engines and refrigerators[END_REF]. Its amplification gain (resp. positive or negative) is dependent on the acoustic impedance within the stack [START_REF] Ceperley | Gain and efficiency of a short traveling wave heat engine[END_REF], and therefore is different for each harmonic of a wave. Indeed, the acoustic impedance is dictated by geometrical considerations, including the boundary conditions at the extremities of the engine.

To favor the wavefront deformation and possibly lead to the formation of a shock wave, balance has to be achieved for each harmonic between the possible amplification from the core and the increase of energy caused by the cascade from the lower harmonics on one hand, and the losses, the decrease of energy due to the cascade to higher harmonics and the possible attenuation from the core on the other hand.

Quasi-standing wave engine

Figure 3 shows the progressive deformation of the steadystate wavefront with the increase of the imposed temperature difference ∆T = T h -T h,o above the onset threshold for two different positions of the core in the resonator of the standing wave engine (configuration (A) of Fig 1). In configuration (A1) the core is placed close to the extremity (x s /L = 0.14) (which corresponds to the final experimental configuration of Ref. [START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]), whereas in configuration (A2) the core is closer to the center (x s /L = 0.22). At first glance, one can see that in both cases the wavefront distorts progressively as the imposed temperature difference increases. However for an amplitude of the pressure oscillations corresponding to a drive ratio (that is ratio of the peak-to-peak amplitude of pressure oscillations to the static pressure) of 10 %, configuration (A1) clearly exhibits formed shocks -the maximum slope shown by the vertical dashed line tends to infinity -when configuration (A2) does not -the maximum slope is finite -, as one can see for the cases in bold black on Figure 3. Configuration (A1) appears more suitable for shock formation. 1) for increasing temperature differences beyond threshold from 5 K to 50 K with a step of 5 K. (A1): x s /L = 0.14, (A2): x s /L = 0.22. The case highlighted in black is for the same drive ratio of 10 %. Dashed lines stand for maximum slope in each case, showing that a shock is formed in the upper case (infinite slope), and not in the lower one (finite slope). Figure 4.(a) shows the saturated waveforms for higher temperature differences in the configuration (A1) where the shock forms more easily (i.e. when the thermoacoustic core is closer to the extremity, x s /l = 0.14). Once the shock is formed, its amplitude continues to increase. The peak amplitude however grows more slowly as more energy is dissipated by the shock for higher ∆T . At high amplitudes, one should note the presence of irregularities in the waveform. The large sharp accidents (circled on Fig. 4.(a)) are not numerical artifacts. They remain present no matter how many harmonics are taken into account in the computation and are independent of the spatial discretization of the problem; they are more likely due to the multiple reflections caused by the presence of the steep temperature gradient and the porous elements interfaces in the TAC. One can also distinguish small "high frequency" oscillations (emphasized by the zoomed portion on Fig. 4.(a)). They are the result of the harmonic series truncation; the energy carried by the truncated part of the spectrum is rejected on the last computed harmonics. Figure 4.(b) illustrates the temporal evolution of the saturated shock wave at regularly spaced points in the engine. The spatio-temporal evolution of the shock in the quasistanding wave engine is the result of the two counterpropagating waves of equivalent amplitudes propagating back and forth, being reflected on its rigid terminations. The propagation of the shock in the engine is highlighted by the arrows. One acoustic period last for one round trip of the shock in the engine, showing that the engine works on its first mode -the half-wavelength mode. This results are to be compared with those obtained experimentally by Biwa et al. [START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]. The computed waveforms present a similar temporal evolution, and amplitude of the same order of magnitude to the measured ones; the model seems to have captured the nonlinear behavior of the wave propagation.

Quasi-traveling wave engine

Figure 5 shows in the same fashion the numerical results for the configuration of quasi-traveling wave engine (B). In this case however, it is worth mentioning that the engine is spatially translation invariant; that is, the position of the core is not relevant to describe the behavior of the engine. Fig. 5.(a) shows the progressive formation and growth of the shock with the increase of the imposed temperature difference above threshold ∆T . Fig. 5.(b) shows a formed traveling shock wave going around the looped engine as emphasized by the black arrow, operating one rotation around the engine per period because the engine is working on its first mode -the one-wavelength mode. Signals at positions 0 and 1 are identical. Since the wave is not purely traveling [START_REF] Hyodo | Stability analysis for the onset of thermoacoustic oscillations in a gas-filled looped tube[END_REF], a small counterpropagating wavefront can also be observed, emphasized by the gray line. These results are to be compared with those obtained experimentally by Biwa et al. for a similar configuration (B)

T h -Tc [K] α [%]
Figure 6: Amplification rate α of the three first modes of the engines versus the temperature gradient imposed for the three configurations (A1), (A2), and (B). Modes 4 to 6 are also given in gray. A positive amplification rate means that a wave at the frequency of the corresponding mode will be amplified by the TAC. α = 0 for mode 1 represents the onset threshold, for which T h -T c = ∆T o .

[23], where a good qualitative agreement is observed.

A linear interpretation

To attempt an interpretation of the results presented above, one can rely -at least partially and cautiouslyon a linear analysis of the amplification of the thermoacoustic engines. To this end, we will rely on the method presented in Ref. [START_REF] Guedra | On the use of a complex frequency for the description of thermoacoustic engines[END_REF]: the amplification rate α is given by the imaginary part of the complex frequency obtained from the linear transfer matrix of the thermoacoustic engine. Figure 6 shows the amplification rate α for the first modes of the three different configurations of engines. If α > 0, the mode is amplified by the thermoacoustic core; that is the core acts as a source of energy at that frequency. If α < 0, the thermoacoustic core behaves as an energy sink and the wave is attenuated. Though the frequencies of the higher modes do not correspond strictly to those of the harmonics of the working mode -mainly because of inharmonicity induced by the presence of the temperature gradient -, they are still close enough for a qualitative analysis of the behavior of the harmonics in the engine.

One can notice that the first modes are more amplified with the increase of the imposed temperature difference along the stack in both cases (A1) and (B), corresponding respectively to the standing wave engine with the core close to the extremity and the traveling wave engine. This means that the greater the temperature difference is, the less the energy cascaded from the fundamental frequency will be attenuated in the TAC, thus facilitating the wave-front deformation. Case (A2) however shows a different behavior: though the fundamental frequency of the working mode seems more easily amplified in this configuration, the second and third harmonics are always (linearly) attenuated. The energy cascaded from the fundamental frequency will then be dissipated in the thermoacoustic core, limiting its transit to the superior harmonics and consequently prevent the deformation of the wavefront. Finally, the shock wave generation in a thermoacoustic engine appears to be related to the attenuation or amplification in its core of the higher harmonics produced by the energy cascade from the fundamental frequency. In cases (A1) and (B), the thermoacoustic core acts as an amplifier for the higher harmonics, supporting the wavefront deformation, whereas in case (A2), it behaves as an energy sink, creating a barrier for energy cascade impeding the deformation. Because such an interpretation is based on observing the variations of an amplification rate defined in the frame of a linear theory, it has to be considered with caution. Nevertheless, this interpretation is also the one supported by Biwa et al. [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF][START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF] from the measurements of intensity at both side of the TAC. This work furthermore brings a numerical approach that successfully renders the nonlinear behavior of the studied engines, as well as an a priori reading of the linear model of the acoustic propagation to support this interpretation.

Conclusion

In this paper, a numerical simulation of the nonlinear propagation in thermoacoustic engines of basic geometry is presented. It uses the formalism of generalized Burgers' equations for nonlinear propagation in the temperature homogeneous waveguide sections and the linearized equations of thermoacoustics in the engine's core, solved numerically to obtain the periodic steady-state regime of auto-oscillation. This expressions proved to be sufficient for a qualitative description of the wave steepening in thermoacoustic engines of simple geometry, for both standing and quasi-traveling wave configurations. Numerical results show good qualitative agreement with experimental results obtained in [START_REF] Biwa | Observation of traveling thermoacoustic shock waves (L)[END_REF][START_REF] Biwa | Observation of thermoacoustic shock waves in a resonance tube (L)[END_REF]. Furthermore, a more accurate description of the thermoacoustic core should help to obtain a better quantitative match. Notably, using a nonlinear profile temperature, both in the stack and the TBT, will influence the amplification rate of the core, particularly in the quasitraveling wave engine. Besides, the description of minor losses is quite rough in this work, only aiming at getting an order of magnitude of the losses. A more accurate description is important for the prediction of the levels of saturation, and should be accompanied by more exhaustive evaluation of the thermophysical parameters in the experimental configuration, such as the presence of acoustically induced streaming or the spatial distribution of temperature in the different elements composing the core. A quantitative comparison would be possible only under these conditions.

3 M = Re 2 ShM = Sh 2 M 2 Figure 2 :

 32222 Figure2: Predominance diagram in the amplitude vs. frequency plane for a cylindrical tube of 4 cm in diameter filled with air at ambient pressure. Each straight line represents the limit of predominance of one phenomenon over another: Merkli condition is the upper limit for laminar flow in boundary layers; above M 2 = Sh, second order effects take over boundary layer effects; M = Sh 2 and M = Re 2 Sh are two limits under which a nonlinear description of the propagation is not necessary. Parallel dashed lines represent the same conditions in a pore of the stack. Adapted from Fig.1of[START_REF] Menguy | Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments[END_REF].

Figure 3 :

 3 Figure3: Computed deformation of the wavefront for two different positions x s /L of the core in the resonator of the standing wave engine (configuration (A) of Fig.1) for increasing temperature differences beyond threshold from 5 K to 50 K with a step of 5 K. (A1): x s /L = 0.14, (A2): x s /L = 0.22. The case highlighted in black is for the same drive ratio of 10 %. Dashed lines stand for maximum slope in each case, showing that a shock is formed in the upper case (infinite slope), and not in the lower one (finite slope).

Figure 4 :Figure 5 :

 45 Figure4: (a): Acoustic pressure waveform computed in the standing wave engine (A1) when the TAC is placed at position x s /L = 0.14, for increasing temperature differences beyond threshold from 5 to 100 K with a step of 5 K up to 50 K, then 10 K, shown at position x = 0. (b): Spatio-temporal evolution of a quasi-standing shock wave computed for ∆T = 100 K in the engine configuration (A1). One can clearly see a shock wave going back and forth in the engine as emphasized by the arrows.

Acknowledgments

The authors wish to acknowledge A. Hirschberg for his help and precious advice in the redaction of this paper.