Denis Efimov
email: denis.efimov@inria.fr

Christopher Edwards
email: c.edwards@exeter.ac.uk

Ali Zolghadri
email: ali.zolghadri@ims-bordeaux.fr

A note on improvement of adaptive observer robustness

In this paper the problem of adaptive observer design in the presence of disturbances is studied, and an augmented adaptive observer is proposed. First, the H ∞ gain of a conventional adaptive observer is estimated, which characterizes the effect of disturbances on output errors. Next, it is shown that if the disturbance is "matched" in the plant equations, then it is possible to introduce additional sliding-mode feedback, dependent on the plant output, improving the accuracy of observation. Simulation results confirm the improvement.

INTRODUCTION

The problem of adaptive observer design for nonlinear systems is very challenging [START_REF] Besanon | Nonlinear observers and applications[END_REF]. In this situation the observer has to generate estimates of the vector of unknown parameters and unmeasured state components. Almost all solutions are based on special canonical representations of the uncertain nonlinear systems, that allows the observer to be designed. An important obstacle is the relative degree between the output signal and the vector of unknown parameters (i.e. the number of derivatives of the output required, before the direct dependence on the vector of unknown parameters is obtained). The observers designed in the case when the degree is one [START_REF] Fradkov | Markov. Adaptive observer-based synchronization for communications[END_REF] and for the higher relative degree case [START_REF] Fradkov | Adaptive observers for nonlinear nonpassifiable systems with application to signal transmission[END_REF], [START_REF] Xu | Residual generation for fault diagnosis in linear time-varying systems[END_REF], [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF], [START_REF] Efimov | Robust adaptive nonlinear partial observers for time-varying chaotic systems[END_REF] have completely different structures, and the dimension of the observers in the latter case is much higher.

The structure of the observer proposed in [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF] for the high relative degree case is quite sophisticated, and strongly related to a canonical form of the plant equations. A small modification to the plant structure may render the observer design impossible. To overcome this limitation, a deviation of equations of the real system from the canonical form can be modeled by external disturbances. These disturbances have to be bounded, but in some situations it is hard to assume that the disturbance is sufficiently small. Thus it is necessary to ensure robustness of the designed observers with respect to such a disturbance.

There exist a number of solutions aimed at improving the robustness in nonlinear systems by applying dynamic or static feedback. Some very promising solutions have been obtained in the area of sliding mode theory since sliding mode feedback is able to fully compensate for matched disturbances even granting the closed loop system finitetime stability [START_REF] Shtessel | Sliding Mode Control and Observation[END_REF]. Recently the sliding mode approach has been successfully applied to adaptive observer design for the case of relative degree one [START_REF] Yan | Fault estimation for single output nonlinear systems using an adaptive sliding mode observer[END_REF], but the application of this theory for adaptive observer design with a high relative degree is somewhat complicated due to the fixed observer structure. In this paper a method is presented for modifying the adaptive observer from [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF] using a sliding-mode feedback in the spirit of [START_REF] Yan | Fault estimation for single output nonlinear systems using an adaptive sliding mode observer[END_REF], that allows us to improve the level of observer robustness with respect to some matched disturbances.

The paper is structured as follows. The adaptive observer equations and the problem statement are given in Section 2. The H ∞ gain of the adaptive observer from [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF] is computed in Section 3. Then the sliding-mode feedback methodology is applied in Section 4. An example of a computer simulation (for a nonlinear pendulum) is presented in Section 5 to demonstrate the efficacy of the approach.

PROBLEM STATEMENT

Consider the following uncertain nonlinear system: ẋ = Ax + φ(y, u) + G(y, u)θ + Bv, y = Cx, (1) where x ∈ R n , y ∈ R p , u ∈ R m are the state, the output and the control respectively, θ ∈ R q is the vector of unknown parameters, v ∈ R s is the vector of external disturbances and v : R + → R s is a (Lebesgue) measurable function of time; the matrices A, B, C are known and have corresponding dimensions, the functions φ : R p+m → R n and G : R p+m → R n×q are also assumed to be known and ensure uniqueness and existence of the system (1) solutions at least locally.

The symbol |x| denotes the Euclidean norm of a vector x (for a matrix A the symbol |A| denotes the induced matrix norm), and for the (Lebesgue) measurable functions v : R + → R s , the norm is defined as ||v|| = ess sup t≥0 {|v(t)|}. For a matrix function A : R + → R s×q we denote ||A|| = |A(t)| . The identity matrix of dimension n × n is denoted as I n . The symbols λ min (A), λ max (A) represent the minimal and maximal eigenvalues of a symmetric matrix A ∈ R n×n .

In this work we will assume that all external and internal signals in the system (1) are bounded:

Assumption 1. ||x|| < +∞, ||v|| < +∞, ||u|| < +∞.
Although assumed to be bounded, the disturbance v may have a large amplitude, and therefore special attenuation techniques have to be applied to ensure reliable estimates for the states in system (1).

The objective of this work is to design an adaptive observer for (1) under the conditions of Assumption 1. The observer has to provide estimates of the vectors x and θ with an improved degree of robustness with respect to the external disturbance v. The proposed design procedure is completed in two steps. Firstly, an adaptive observer is designed such that the H ∞ gain between the disturbances and the output errors is calculated (and thus it can be minimized). Secondly, an additional sliding mode output injection is applied to further reduce the influence of the disturbance components which cannot be completely rejected by the first step. The proof of stability is given for each case. An analysis of the structure of the resulting observers suggests that the augmented adaptive observer behaves better than the conventional adaptive observer with H ∞ gain optimization. The simulation results confirm this.

An adaptive observer for the system (1) has been proposed in [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF], and takes the form:

ż = Az + φ(y, u) + G(y, u) θ + L(y -Cz) + Ω θ, (2) where Ω = (A -LC)Ω + G(y, u), θ = γΩ T C T (y -Cz). (3)
In ( 2)-( 3), z ∈ R n is the estimate of x, θ ∈ R q is the estimate of θ, and Ω ∈ R n×q is an auxiliary/filter variable, that helps overcome possible high relative degree obstructions in system (1). In (3) γ > 0 is a design parameter, and L is the observer gain that is chosen to ensure a Hurwitz property for the matrix A -LC. The analysis of the estimation abilities of the observer in (2), ( 3) is based on the errors δ = x -z + Ω θ and θ = θ -θ whose dynamics can be shown to have the form:

δ = (A -LC)δ + Bv, (4) θ = γΩ T C T (Cδ -CΩ θ).
(5)

From equation ( 4) we conclude that δ → 0 for v = 0 and the variable δ stays bounded for any bounded disturbance v. The Hurwitz property of the matrix A -LC and Assumption 1 imply boundedness of the variable Ω. If the signal G(y, u) is persistently exciting (PE) [START_REF] Anderson | Exponential stability of linear equations arising in adaptive identification[END_REF], Yuan and [START_REF] Yuan | Probing signals for model reference identification[END_REF], then due to the filtering property of the variable Ω, the variable CΩ is also PE. Moreover, it is possible to show [START_REF] Efimov | Robust adaptive nonlinear partial observers for time-varying chaotic systems[END_REF] that for any bounded signal Cδ, the variable θ has a bounded response, and if Cδ → 0, then θ → 0 also. Since as a consequence of Assumption 1 the variables x, δ, θ and Ω are bounded, the variable z is also bounded according to the definition of δ. In addition, it is easy to verify that if θ → 0, then z → x. Thus, we have shown that the variables z, Ω and θ stay bounded for any bounded disturbance v, and for the case when v = 0 the estimates θ and z converge to their ideal values θ and x respectively. Thus the proof is based on general stability arguments, and no strict Lyapunov function has been proposed.

To improve the accuracy of estimation in the observer (2)-( 3), one can design the matrix L optimizing the H ∞ gain of the system (C, A -LC, B) [START_REF] Henry | Design of fault diagnosis filters: a multi-objective approach[END_REF]. In this case the influence of the disturbance v on the output Cδ of the system in ( 4) is minimized, and implies a proportional amelioration of the error θ in ( 5). This will be proven formally in Section 3 below. Moreover, it will be shown that if v is a Lipschitz function of x, then by proper selection of L it is possible to attenuate the disturbance v. If the amplitude of the disturbance v is large, then such an optimization could be insufficient. This problem can be addressed by modifying equation ( 2): this approach is later demonstrated in Section 4 using a sliding mode technique.

CONVENTIONAL ADAPTIVE OBSERVER

First let us show that the system in ( 5) is input-to-state stable with respect to the input Cδ. Lemma 1. Let the variable Ω T C T be PE and bounded, i.e. 0 < ρ = ||CΩ|| < +∞ and there exist constants ϑ > 0 and > 0 such that t+ t Ω(τ ) T C T CΩ(τ )dτ ≥ ϑI q ∀t ≥ 0. Then a) there exists a continuous symmetric matrix function P : R + → R q×q such that ρ -2 I q ≤ 2γP (t) ≤ αI q for all t ≥ 0, where α = γη -1 e 2η and η = -0.5 -1 ln(1

- γϑ 1+γ 2 2 ρ 4 ); b) for all t ≥ 0 Ṗ (t)-γP (t)Ω(t) T C T CΩ(t)-γΩ(t) T C T CΩ(t)P (t)+I q = 0; c) for S(t, θ) = θT P (t) θ we have for all θ ∈ R q , δ ∈ R n and t ≥ 0 Ṡ ≤ -γα -1 S + 0.5ρ 2 α 2 |Cδ| 2 .
In addition, for all θ(0) ∈ R q and all t ≥ 0 the following estimate is satisfied:

| θ(t)| ≤ ρ √ α[e -0.5γα -1 t | θ(0)| + ρα||Cδ||].
The proof is omitted for brevity.

Remark 2. This lemma also provides an estimate on the fastest rate of decrease of the parametric estimation error θ(t). Indeed, the rate of decrease equals 0.5γα -1 = 0.5ηe -2η = g(η), where the function g(η) for η > 0 has the unique maximum value g max = 0.25 -1 e -1 at η opt = 0.5 -1 . However, as shown in the proof of Lemma 1, 0 < η ≤ η max = 0.5 -1 ln 2. Since η opt > η max and the function g(η) is strictly increasing for 0 < η ≤ η opt , then the maximal rate of error decrease that can be reached is g(η max ) = 0.125 -1 ln 2.

Note that in order to use these estimates we have to know the value of the coefficient ρ. The existence of such a ρ follows from Assumption 1 for the Hurwitz matrix A-LC, but to compute it we have to know ||x||, which is assumed to be unavailable. Theorem 3. Suppose Assumption 1 is satisfied, the variable Ω T C T is PE, i.e. there exist constants ϑ > 0 and > 0 such that

t+ t Ω(τ ) T C T CΩ(τ )dτ ≥ ϑI q ∀t ≥ 0,
and there exists a n×n matrix W = W T > 0 and constants r > 0, h > 0 such that (A -LC)

T W + W (A -LC) + 0.5rα 2 C T C +hW BB T W + γα -1 W ≤ 0, (6)
where α = γη -1 e 2η , η = -0.5 -1 ln(1 -γϑ 1+γ 2 2 ||CΩ|| 4 ) and γ > 0. Then in system (1), ( 2), (3) all solutions are bounded for all t ≥ 0, and

|x(t) -z(t)| | θ(t)| ≤(1 + ||Ω||){ α 2 α 1 e 0.5γα -1 t δ(0) θ(0) + α α 1 γh ||v||},
where α 1 = min{λ min (W ), 0.5rγ -1 ||CΩ|| -4 }, α 2 = max{λ max (W ), 0.5r||CΩ|| -2 γ -1 α}.

Proof. It follows from Assumption 1 the variable Ω is bounded for any Hurwitz matrix A -LC (the stability of this matrix follows from the Riccati equation for the matrix W ), and it is also a continuous matrix function of time. Define ρ = ||CΩ|| < +∞. Consider as a candidate Lyapunov function

V (δ) = δ T W δ, then V = δ T [(A -LC) T W + W (A -LC)]δ + 2δ T W Bv. Note that 2δ T W Bv ≤ hδ T W BB T W δ + h -1 v T v, then V ≤ δ T [(A -LC) T W + W (A -LC) +hW BB T W ]δ + h -1 |v| 2 .
Define U (t, δ, θ) = V (δ)+rρ -2 S(t, θ), then using the result of Lemma 1 we have

U ≤ δ T [(A -LC) T W + W (A -LC) + hW BB T W ]δ -γrρ -2 α -1 S + 0.5rα 2 |Cδ| 2 + h -1 |v| 2 ≤ δ T [(A -LC) T W + W (A -LC) + 0.5rα 2 C T C +hW BB T W ]δ -γrρ -2 α -1 S + h -1 |v| 2 ≤ δ T [(A -LC) T W + W (A -LC) + 0.5rα 2 C T C +hW BB T W + γα -1 W ]δ -γα -1 U + h -1 |v| 2 ≤ -γα -1 U + h -1 |v| 2 .
Note that according to the result of Lemma 1 the following inequalities hold:

λ min (W )|δ| 2 + 0.5rγ -1 ρ -4 | θ| 2 ≤ U (t, δ, θ) ≤ λ max (W )|δ| 2 + 0.5rρ -2 γ -1 α| θ| 2 .
Therefore, for all δ(0) ∈ R n and θ(0) ∈ R q , and for all t ≥ 0 we obtain:

δ(t) θ(t) ≤ α 2 α 1 e 0.5γα -1 t δ(0) θ(0) + α α 1 γh ||v||.
Since the variables x, δ, Ω and θ are bounded, then the variable z has the same property (z = x -δ + Ω θ), and in addition

|x(t) -z(t)| ≤ |δ(t)| + ||Ω|| | θ(t)| ≤ (1 + ||Ω||) δ(t) θ(t) ,
which gives the required estimate.

This theorem estimates the H ∞ gain of the observer from the input v to the estimation errors θ and x -z. A minor modification of the LMI for W in Theorem 3 is needed in order to ensure the H ∞ gain from v to a generic output Z(x -z), where Z is a matrix (this modification is omitted). Next, a minimization of the H ∞ gain can be performed, using the nonlinear optimization routines (the detailed presentation of this part is skipped for brevity).

The Riccati equation ( 6) depends on α which is a function of the PE parameters , ϑ and ρ, which can be unknown a priori. However, this equation can be rewritten as

(A -LC) T W + W (A -LC) + r 1 C T C
+hW BB T W + r 2 W ≤ 0, where r 1 = 0.5rα 2 and r 2 = γα -1 are the new tuning parameters. If this Riccati inequality is satisfied for some h > 0, r 1 > 0 and r 2 > 0, and the PE property is true for Ω T C T , then there exist γ > 0 and r > 0 such that (6) is satisfied. Remark 4. It is straightforward to verify that if the disturbance v contains a globally Lipschitz nonlinear function of the state, then augmenting the value of h, and introducing a minor modification in (2), it is possible to attenuate this type of uncertainty. Indeed, let v(t) = µ(x) + ṽ(t), where µ : R → R s is a function which for any x, x ∈ R n satisfies the inequality |µ(x) -µ(x )| ≤ L µ |x -x | for some L µ > 0 and ṽ : R + → R s is a new bounded exogenous input. Then introducing in (2) an additional term Bµ(z) and solving the LMI of Theorem 3 for a sufficiently large value of h > 0, complete compensation of the influence of µ(x) can be achieved.

SLIDING-MODE ADAPTIVE OBSERVER

In this case the observer equation ( 2) is modified as follows:

ż = Az + φ(y, u) + G(y, u) θ + L(y -Cz) + Ω θ +kB sign[F (y -Cz)], (7)
where k > 0 is a design constant, F ∈ R n×p is a matrix to be designed and sign[v] = v/|v| for a vector v ∈ R m . Theorem 5. Let assumption 1 be satisfied and let k = ||v||. Assume the variable Ω T C T is PE, i.e. there exist ϑ > 0 and > 0 such that t+ t Ω(τ ) T C T CΩ(τ )dτ ≥ ϑI q ∀t ≥ 0, and suppose there exists a matrix W = W T > 0, a matrix

F ∈ R n×p and constants r > 0, ω > 0 such that (A -LC) T W + W (A -LC) + γα -1 W +(0.5rα 2 + ω)C T C ≤ 0, (8)
W B = C T F T , where α = γη -1 e 2η , γ > 0, η = -0.5 -1 ln(1 - γϑ 1 + γ 2 2 ||CΩ|| 4 ).
Then in the system (1), ( 3), ( 7), all variables z, Ω and θ are bounded and

|x(t) -z(t)| | θ(t)| ≤ (1 + ||Ω||)( α 2 α 1 e -0.5γα -1 κt δ(0) θ(0) +2 α α 1 ωγκ |F | ||v||).
where κ = 0.5r

0.5r+ω||CΩ|| 2 , α 1 = min{λ min (W ), γ -1 (0.5r||CΩ|| -4 + ω)}, α 2 = max{λ max (W ), γ -1 (0.5r||CΩ|| -2 + ω)α}.
Proof. As in the previous theorem the variable Ω is bounded for any Hurwitz matrix A -LC (the stability of this matrix follows from the Lyapunov equation for the matrix W ) and let ρ = ||CΩ||. The error δ = x -z + Ω θ for the system (1), ( 3), ( 7) has the following dynamics: δ = (A -LC)δ + B{v -k sign[F (y -Cz)]}.

Consider three Lyapunov function components given by V (δ) = δ T W δ, Z( θ) = γ -1 θT θ and S(t, θ) = θT S(t) θ. These have the property that

V = δ T [(A -LC) T W + W (A -LC)]δ + +2δ T C T F T {v -k sign[F (y -Cz)]} ≤ -δ T [γα -1 W + (0.5rα 2 + ω)C T C]δ +2[y -Cz] T F T {v -k sign[F (y -Cz)]} +2 θT Ω T C T F T {v -k sign[F (y -Cz)]} ≤ -δ T [γα -1 W + (0.5rα 2 + ω)C T C]δ +4| θT Ω T C T | |F | ||v|| Ż = 2 θT Ω T C T (Cδ -CΩ θ) ≤ |Cδ| 2 -| θT Ω T C T | 2 Ṡ ≤ -0.5 θT θ + 0.5ρ 2 α 2 |Cδ| 2 .
Introduce a new candidate Lyapunov function given by U (t, δ, θ) = V (δ) + ωZ( θ) + rρ -2 S(t, θ), then we obtain

U ≤ -γα -1 δ T W δ -0.5rρ -2 θT θ +4| θT Ω T C T | |F | ||v|| -ω| θT Ω T C T | 2 ≤ -γα -1 δ T W δ -0.5rρ -2 θT θ + 4ω -1 |F | 2 ||v|| 2 ,
where the inequality 4s|F | |v||-ωs 2 ≤ 4|F | 2 ||v|| 2 /ω, which is satisfied for all s ∈ R + , has been used in the last step of transformations. Since

γ -1 (0.5rρ -4 + ω)| θ| 2 ≤ ωZ( θ) + rρ -2 S(t, θ) ≤ γ -1 (0.5rρ -2 + ω)α| θ| 2 , we finally obtain U ≤ -γα -1 κU + 4ω -1 |F | 2 ||v|| 2 ,
and for all t ≥ 0

δ(t) θ(t) ≤ α 2 α 1 e -0.5γα -1 κt δ(0) θ(0) +2 α α 1 ωγκ |F | ||v||.
According to the proof above, the variables x, δ, Ω and θ are bounded, consequently the variable z has the same property (z = x -δ + Ω θ), and again

|x(t) -z(t)| ≤ (1 + ||Ω||) δ(t) θ(t) ,
which gives the required estimate. Remark 6. Comparing with the case of Theorem 3, in Theorem 5, the accuracy of estimation is improved, since in part, the disturbance v is compensated by sign(•). Indeed, the dynamics of the state estimation error e = x -z yields the following differential equation ė

= (A -[L + γΩΩ T C T ]C)e -G(y, u) θ +B[v -ksign(F Ce)],
which for k = 0 reduces to (2). Consider the Lyapunov function V (e) = e T W e, then

V ≤ e T {(A -LC) T W + W (A -LC)}e

+2e T C T F T (v -ksign(F Ce)) -2e T W [γΩΩ T C T Ce + G(y, u) θ] ≤ e T {(A -LC) T W + W (A -LC)}e -2|F Ce|(k -||v||) -2e T W [γΩΩ T C T Ce + G(y, u) θ].
Therefore, for the observer in (7), with k = ||v||, the term proportional to |F Ce| disappears and the influence of v is cancelled. The Lyapunov equation for the matrix W is also simpler than in the case of Theorem 5 (and an additional decision variable F appears).

In order to find a solution to the LMI in (8) we have to reformulate it, to avoid direct dependence on α. To this end it can be rewritten as follows:

(A -LC) T W + W (A -LC) + r 1 C T C + r 2 W ≤ 0,
where r 1 = (0.5rα 2 + ω) and r 2 = γα -1 are new tuning parameters.

Remark 7. The introduction of the "structural constraint" W B = C T F T in (8) imposes limitation on the case of triple (A, B, C) to which these results are applicable. However a constructive approach to solving such problems is given in [START_REF] Edwards | On the solvability of the constrained lyapunov problem. Automatic Control[END_REF]. Equivalent conditions can be found in Lemma 1 of [START_REF] Efimov | Adaptive tuning to bifurcation for time-varying nonlinear systems[END_REF].

EXAMPLE

Consider a nonlinear pendulum: The problem studied in this paper is that of improving the performance of the adaptive observer originally proposed in [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF]. Following [START_REF] Yan | Fault estimation for single output nonlinear systems using an adaptive sliding mode observer[END_REF], a scheme augmented by including a sliding mode term is proposed. First, the H ∞ gain with respect to external disturbances is evaluated to minimize the effect of disturbances on the output errors in a H ∞ sense. Next, a slidingmode modification of the adaptive observer equations from [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF] is proposed in order to improve disturbance attenuation. The conditions for stability of the new scheme have been established. Numerical simulations confirm the accuracy improvement obtained from the newly proposed observer with respect to that proposed in [START_REF] Zhang | Adaptive observer for multiple-input-multipleoutput (mimo) linear time varying systems[END_REF].

n = 2, A = 0 1 0 -β , B = 1 0 , C = [1 0] , φ(y, u) = 0 u -2 sin(y) , G(y, u) = 0 sin(y) , v (

  Fig. 1. The estimates of unknown parameter

This work was supported in part by the Government of Russian Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031).


Juan Wang
email: juan.wang@ensam.eu

Florent Ravelet
email: florent.ravelet@ensam.eu

Farid Bakir
email: farid.bakir@ensam.eu

INFLUENCE OF DESIGN PARAMETERS ON THE GLOBAL PERFORMANCES OF LOW-SPEED COUNTER-ROTATING AXIAL-FLOW FANS

The present work aims at experimentally investigating the effects of some parameters on the performances of a counterrotating stage, and on the instationary flow between the rotors. Three counter-rotating fans, which have the same design point, have been designed. These systems differ by the distribution of the loading and of the ratio of angular velocity between the front rotor and the rear rotor. All the configurations have been tested in a normalized test rig, where the ratio of angular velocities and the axial distance between the two rotors can be varied. The influence of these parameters are then addressed by analysing the experimental results of the static pressure rise and static efficiency, as well as of the wall pressure fluctuations registered by a microphone at the wall. The three systems achieve the design point with a high efficiency. The counter-rotating systems lead to at least a 10 percentage points gain in static efficiency at the design flow rate, compared to the typical peak efficiency of a traditional rotor-stator stage. Meanwhile, counter-rotating systems display good working stabilities at very low volume flow rates. In addition, at the design speed ratio, the overall performance decreases almost monotonically with the axial distance. Nevertheless, an optimum in axial distance can be found for higher speed ratios. Finally, the investigations of the wall pressure fluctuations show that the amplitudes of power spectral density corresponding to the blade passing frequency of the rear rotor are significantly higher than that of the front rotor. The interaction peaks are also stronger for an equal distribution of the work on the two rotors.

NOMENCLATURE

INTRODUCTION

Adoption of a counter-rotating machine has opened a way to design high performance and compact turbomachines in various industrial domains. It has been already applied in the areas of subsonic fans, pumps and turbines [START_REF] Cho | Numerical and experimental analyses for the aerodynamic design of high performance counter-rotating axial flow fans[END_REF][START_REF] Shigemitsu | Internal flow measurement with ldv at design point of contra-rotating axial flow pump[END_REF][START_REF] Moroz | Comparison of counter-rotating and traditional axial aircraft low-pressure turbines integral and detailed performances[END_REF]. A Counter-Rotating System (CRS) is composed of a Front Rotor (FR) and a Rear Rotor (RR) which rotates in the opposite direction, as shown in Fig. 1. Compared with traditional rotor-stator system, the Rear Rotor is used not only to recover the static head but also to supply energy to the fluid. Therefore, to achieve the same performance, the use of a CRS may lead to a reduction of the rotational speed and may generate better homogeneous flow downstream of the stage.

On the other hand, the mixing area in between the two rotors induces complicated interacting flow structures. The understanding of this highly unsteady flow in the mixing area is an open problem. Moreover, the design method of such machines is still not sophisticated, due to a lack of systematic studies on the influence of free parameters, such as the distribution of load (L), the axial distance (S), the ratio of the rotation rates (θ ) and so on... In the Dynfluid Laboratory, series of experiments focused on axial counter rotating fans have been performed [START_REF] Nouri | Experimental study of the instationary flow between two ducted counter-rotating rotors[END_REF]. Based on this research, three CRS (JW1, JW2 and JW3) have been designed to attain the same design point, while varying the distribution of load. The main concern in this paper is the influence of the distribution of load (L) on the overall performances of a CRS.

The conception of the three CRS and their design parameters are presented in Section 2. Then the experimental facility, the measurement method and the analysis of uncertainty as well as repeatability are presented in Section 3. The results are presented in Section 4. The overall performances for the nominal S and θ are first presented in § 4.1, then the influence of the axial distance and of the variation of θ are discussed for each CRS, in respectively § 4.2 and § 4.3. In addition, the levels and spectral densities of the wall pressure fluctuations in the mixing area between the front and rear rotors are compared in § 4.4.

DESIGN OF THE COUNTER-ROTATING FANS SYS-TEMS 2.1 Methodology

The objective is to design three CRS which can achieve the same design point for various distributions of the total load between the front and rear rotors. The design point is presented in Tab. 1. In this table, the total pressure rise is the difference of total pressure between the inlet and the outlet of the CRS. One of the design constraint is to have a pure axial flow downstream of the CRS. Therefore, the corresponding expected static pressure rise is calculated as ∆P sC = ∆P tC -1 2 ρ a ( Q vC πD 2 /4 ) 2 ≈ 373 Pa. With this in mind, firstly, the front rotor was designed to achieve a part of the total pressure rise at design flow rate, by the in-house code MFT [START_REF] Noguera | Design and analysis of axial pumps[END_REF]. The details of the conception method can be found in Ref. [START_REF] Nouri | Experimental study of the instationary flow between two ducted counter-rotating rotors[END_REF]. The outlines could be depicted as follows: the geometrical parameters of the Front Rotor are designed by MFT with the inverse method. Then the axial and tangential velocities at the outlet of the Front Rotor are analysed and taken as the input conditions for the conception of the Rear Rotor. Therefore, the Rear Rotor is adapted to the outflow of the Front Rotor, and is designed such that the absolute tangential velocity at the outlet of the system vanishes.

Characteristics of the three different CRS

The main parameters can be found in Tab. 2 and Fig. 2. The distribution of load is defined as the ratio of the total pressure rise due to the Rear Rotor to that of the Counter-Rotating System at the design flow rate: L = ∆P t, RR ∆P t . All the three CRS have different L at the same design point. The Front Rotors of the three systems are designed with the same blade loading repartition, with a "Constant Vortex" Design (see Refs. [START_REF] Vad | Forward blade sweep applied to low-speed axial fan rotors of controlled vortex design: An overview[END_REF][START_REF] Sarraf | Experimental study of blade thickness effects on the overall and local performances of a controlled vortex designed axial-flow fan[END_REF]). The peculiar features of each system are the following: JW1. The Front Rotor of JW1 (JW1 FR ) is designed to have large stagger angles, in order to obtain a steep curve of static pressure rise ∆P s as a function of the volume flow rate Q v (see Fig. 2). Aside from this, the other parameters (rotation rates and ratio θ ) are very similar to those of the configuration that was studied in Ref. [START_REF] Nouri | Experimental study of the instationary flow between two ducted counter-rotating rotors[END_REF].

JW2. Among the three CRS, JW2 has the highest L = 52%, that is to say, in this Counter-Rotating System, the Rear Rotor transfers more energy to the fluid than the Front Rotor. Consequently, the Front Rotor JW2 FR possess the lowest ∆P s among the three Front Rotors. It is furthermore designed with low stagger angles and has a slowly decaying characteristics (see Fig. 2). The Rear Rotor rotates 1.44 times faster than the Front Rotor in that configuration. JW3. This is an extreme case where JW3 FR leads to the highest and steepest characteristics, as shown in Fig. 2. As a result, JW3 has the lowest L = 23% among the three configurations, and the Rear Rotor rotates much more slowly than the Front Fotor. Finally, please note that since the Rear Rotors are designed to rectify the outflow of the Front Rotors toward the axial direction, the shape of the Rear Rotors that are obtained is not usual, with non-monotonic stagger angle and blade camber profiles (see Fig. 3 for an example on JW1 RR ).

N FR /N RR θ C L Z FR /Z RR γ FR (

EXPERIMENTAL FACILITIES 3.1 Test rig

The experimental investigation are performed in a test rig, built according to the ISO-5801 standard. First, the air comes into the test pipe of diameter D = 380mm trough a bell mouth, then passes through the driving motor of the Front Rotor, and is homogenized by a honeycomb. Next, energy is transferred to the fluid by the two rotors. They are separated by a series of blocks for the purpose of varying the axial distance S between them. Then, the outflow passes the driving motor of the Rear Rotor and an anti-gyration device to remove the rotational component of the flow before the measurement of the static pressure by 4 pressure taps. After that, the fluid goes through an ISO-5167 orifice plate in order to measure the volume flow rate (Q v ). Finally, the fluid is regulated by an axial blower and an iris damper before being discharged into the ambient atmosphere.

Experimental method and Estimation of the uncertainties

The uncertainty are first estimated according to the repeatition of measurements at the design point. Specifically, ten measurements have been performed at the same rotation rates and for a fixed diameters of the iris damper.

Density of air ρ a . The static pressure rise and the power consumption vary with the density of air ρ a . The actual density is evaluated according to the the ISO-5801 standard, by measuring the atmospheric pressure P a , the dry temperature T ad and the wet temperature T aw , from which the partial pressure of the water vapour P av is computed:

ρ a = P a -0.378P av 287T ad (1)
According to the repetition of measurements and simple propagation of uncertainty rules, the relative uncertainty of ρ a is ±0.3%. Please note that all the presented results are rescaled to a reference density ρ a = 1.21 kg.m -3 .

Volume flow rate Q v . As presented previously, the volume flow rate is measured by an ISO-5167 orifice plate at more than 15D downstream of the CRS. According to the standard:

Q v = αεπd 2 4 2 ∆P q ρ a (2)
with ∆P q the pressure drop measured by a differential pressure transducer. According to the repetition of measurements and to the accuracy of the transducers, the uncertainty of Q v is ±0.5% at the design flow rate.

Correction for the losses. In order to eliminate the influence of the experimental facilities such as the honeycomb, the driving motor housing and the anti-gyration device, the static pressure drop was measured with both the Front and Rear Rotors removed. Meanwhile, the axial blower at the outlet of the test rig was used to create a flow. Then the correction Corr is modelled as a function of the orifice plate pressure drop ∆P q . At the design point, the uncertainty of Corr is about ±0.5%, which means ±0.6 Pa.

Static pressure rise ∆P s . The static pressure rise ∆P s is the difference between the static pressure downstream of the CRS and the inlet total pressure:

∆P s = ∆P v +Corr - 1 2 ρ a ( Q v πD 2 /4 ) 2 , (3)
with ∆P v being measured by the average of 4 pressure taps downstream of the anti-gyration device, with an uncertainty of ±1 Pa at the design point. The uncertainty of ∆P s is then ±4 Pa close to the design flow rate, that is a relative uncertainty of ±1%.

Power consumption P wt . The total power consumed by the CRS is defined as: P w,t = P w, FR + P w, RR = τ FR ω FR + τ RR ω RR [START_REF] Nouri | Experimental study of the instationary flow between two ducted counter-rotating rotors[END_REF] with τ the torque of motor, which is measured by the servocontrollers of the two brush-less AC motors and has been calibrated against a rotating torquemeter. It is corrected by the value measured when the rotor is removed from the shaft. The uncertainties of the total power consumption is ±4.5 W that is 0.8% of the power at the design point.

Static efficiency η s . The static efficiency is defined as:

η s = ∆P s × Q v P w,t (5)
The relative uncertainty of η s is ±2.3% (±1.5 percentagepoints at the nominal flow rate). . 3. Overall, the isolated Front Rotors of each system achieves the predicted values quite well. When coupled to their Rear Rotors to form the Counter-Rotating Systems, the ∆P s of the three are close to the design point, with at maximum a discrepancy of -11.0% for JW3.

RESULTS AND DISCUSSION

Furthermore, it can be observed from Fig. 4(a) that the characteristics can be divided into 4 regions where the slopes are different:

Region I, Q v = [0, 0.38].
In this region, the volume flow rates are very low. For each configuration, the characteristic curves of the Front Rotors alone and of the Counter-Rotating Systems show similar trend (i.e., flat curves for JW2 and JW2 FR , significantly negative slopes for JW3 and JW3 FR ). Pa.m -3 .s for JW1, JW2 and JW3). Please note that these values could change rather strongly according to the chosen points. They nevertheless illustrate the trend.

Region III, Q v = [0.6, 1]. In this region, of moderate partial flow rates, the curves of ∆P s have the smallest slopes. The values are roughly -110, 9 and -139 Pa.m -3 .s for JW1 FR , JW2 FR and JW3 FR , respectively. Please note the almost zero or even slightly positive slope of JW2 FR , that would lead to poor working stability for this fan if it were used alone in an air-loop. For the three CRS, the slopes are increased to -458, -435 and -208

Pa.m -3 .s for JW1, JW2 and JW3.

Region IV, Q v = [1, 1.3].
In this region of overflow, curves of ∆P s have the biggest slopes. The values are roughly -282, -133 and -288 Pa.m -3 .s for JW1 FR , JW2 FR and JW3 FR , respectively. Whereas the slopes are increased to -951, -836 and -468 Pa.m -3 .s for JW1, JW2 and JW3.

In short, the characteristics of the three CRS are steeper than that of the Front Rotors alone and are always negative on a wide range of partial and over flow rates. This is particularly impressive for the Counter-Rotating System JW2 for which the curve of JW2 FR working alone is quite flat. The presence of the Rear Rotors thus contributes to maintain the system stability even at extremely low flow rates in this type of machines. One can finally notice that the best system in terms of static pressure rise is JW1, with an intermediate repartition of the work on the two rotors (60% for JW1 FR , 40% for JW1 RR ). The values of the static efficiencies η s at the design flow rate are reported in Tab. 4. The three CRS are very efficient, the typical peak efficiency of a traditional rotor-stator stage being of the order of 55% according to Ref. [START_REF] Bakir | Efficient stator designed for automotive engine cooling fan systems[END_REF], and up to 60% for exceptional stages. For a single rotor stage, the typical static efficiency is 50%, up to 55% (see also Refs. [START_REF] Sarraf | Experimental study of blade thickness effects on the overall and local performances of a controlled vortex designed axial-flow fan[END_REF][START_REF] Bakir | Efficient stator designed for automotive engine cooling fan systems[END_REF][START_REF] Yen | Exit flow field and performance of axial flow fans[END_REF]). In view of this, the Counter-Rotating System is a promising solution for the designers who seek for high static efficiency turbomachines.

FR

Influcence of the axial distance S on the performances of CRS at

θ = θ C and Q v = 1 m 3 .s -1
The variations of ∆P s and η s with the axial distance S between the Front Rotor and the Rear Rotor at the design flow rate are plotted in Fig. 5 for the three CRS working on their nominal speed ratios θ = θ c . As a reference, the values for Front Rotors alone are plotted at S = 0 with closed symbols. The values decrease monotonically with S for the three CRS, contrary to the results observed by Pundhir & Sharma [START_REF] Pundhir | A study of aerodynamic performance of a contra-rotating axial compressor stage[END_REF] on a transonic case where an optimum in distance was found. In the present

Q v = 1m 3 .s -1 , θ = θ c
case, the static pressure rise declines rapidly as S < 40 mm, then decreases slightly as S ∈ [40, 100] mm before decreasing significantly as S ∈ [100, 300] mm. When the distance is small, on the one hand the Rear Rotor recovers more swirl energy downstream of the Front Rotor. But on the other hand, the power consumption of the Front Rotor slightly increases. Nevertheless, in our case the relative increase in static pressure is greater than the relative increase in power consumption, which leads firstly, to a monotonic variation of the efficiency and secondly, to a stronger effect of the distance on the static pressure rise than on the efficiency.

Inluence of speed ratio θ on the performances of CRS 4.3.1 JW1

As illustrated in Fig. 6, the nominal flow rate is shifted towards higher flow rates as θ increases. In this section, the maximum static efficiency η s, max is presented and analysed instead of

η s at Q v = 1 m 3 .s -1 .
For JW1, the performances are improved up to a maximum static efficiency of 68% when θ is increased from 0.9θ c to 1.25θ c . This feature has also been observed for the system used in Ref. [START_REF] Nouri | Experimental study of the instationary flow between two ducted counter-rotating rotors[END_REF].

Due to the limitation of the rotational speed of the motors at maximum 3000 rpm, it is difficult to get higher θ and at the same time keep N FR = 2300 rpm. One solution to get higher speed ratios is to decrease the N FR . Figure . 7 presents the η s, max as a function of θ /θ c at N FR = 2300 rpm and N FR = 1400 rpm. Overall, the two curves have the same trend but are shifted by about 2 percentage points. This discrepancy could be explained as an effect of the Reynold's number Re. Lower Re induces larger losses in the flow which deteriorates the efficiency.

The Reynold's number based on the inlet relative velocity W and on the chord length L chord at mid-span, Re = ρ a W L chord µ is of the order of Re ≈ 1.2 × 10 5 for JW1 at N FR = 2300 rpm, and Re ≈ 7.7 × 10 4 at N FR = 1400 rpm. Aside from the influence of Re, it can be seen in Fig. 7 that η s, max reaches a peak value for θ ≈ 1.15θ c and decreases for θ 1.2θ c . Figure. 8 gives the trend of the variation of the maximum static efficiency η s, max with the speed ratio for the three CRS. Favourably, all the CRS could maintain the η s, max beyond 60% even working at far from design condition and low rotational speed. This shows that a CRS does not only improve the efficiency to a high level, but is also robust at maintaining its high performance at various off-set conditions.

JW1, JW2 and JW3

Additionally, it is obvious that η s, max increases significantly as θ is increased to θ c , then continues to rise slightly as θ reaches 1.2θ c , for all the three CRS. Nevertheless, the maximum efficiency drops down moderately as θ > 1.2θ c for JW1 and JW2, but on the contrary, for JW3 it continues to increase untill θ = 1.6θ c . This could own to the increased contribution of the Rear Rotor for this system.

Influence of both the axial distance S and θ on the performances of JW1

In the previous analysis, the influence of S and θ on the performances have been investigated separately. The figure. 9 0.85 0.9 0.95 demonstrates that the identified trend of the maximum static efficiency increase as θ increases is similar for all the axial distance in [START_REF] Pundhir | A study of aerodynamic performance of a contra-rotating axial compressor stage[END_REF]200] mm.

Analysis of the wall pressure fluctuations

The wall pressure fluctuations are measured by a microphone situated at a distance Z p = 5 mm downstream of the Front Rotor (see Fig. 1). The power spectral density (PSD) of these fluctuations for the three systems working at their design conditions and with S = 10 mm are plotted in the Figs. 10 to 12. For each system, two flow rates have been studied: the design flow rate Q v = 1 m 3 .s -1 and a partial flow rate Q v = 0.6 m 3 .s -1 . One can notice the presence of several peaks in the spectra. These peaks can be divided into three types: those that correspond to

Q v = 1 m 3 .s -1 and Q v = 0.6 m 3 .s -1 .
Several common features can be noticed. Firstly, the amplitudes of the peaks corresponding to f bp f , RR and its harmonics are always significantly higher than that of f bp f , FR . The influence of the Rear Rotor propagates upstream (potential effect) and is stronger than that of the Front Rotor (usually attributed to the wakes of the blades). Then, one can notice that the two rotors are in strong non-linear interaction at the design volume flow rate. However, for Q v = 0.6 m 3 .s -1 , the peaks corresponding to the interactions are dramatically attenuated, and in contrast, the amplitudes corresponding to both f bp f , FR and f bp f , RR are increased.

In order to compare now the three systems, the amplitude of the peaks corresponding to f bp f , FR and f bp f , RR are reported in Tab. 5 for the design volume flow rate. JW2 has the highest amplitude for f bp f , RR and the lowest for f bp f , FR among the three CRS. This is consistent with the ratio of rotational speeds. Finally, one can notice a correlation between the levels and the rotation rates of the rotors. The total level of the pressure fluctuations is thus obviously the lowest for JW3, and is almost similar for JW1 and JW2 when the sums of the two rotation rates are of the same order of magnitude.

CONCLUSIONS

Three subsonic Counter-Rotating axial-flow fans Systems (JW1, JW2 and JW3) have been designed. The design objective is to reach the same design point with various distributions of the work and rotation rates between the rotors. Series of experiments have been conducted to explore the influence of design parameters on the overall performances of these CRS. It can be concluded that:

1. All three CRS can achieve the design point, with maximum 11% of discrepancy. Besides, the best static efficiency is very high -above 65%for the three systems. The design method is thus satisfying for a pre-study, but may be improved, especially by including an analysis of the total system and developing an equation of the radial equilibrium for the whole stage. 2. The performance map shows that over the major part of operating range, with the contribution of the Rear Rotor, the slopes of the CRS are always negative, which stands for good working stability of this kind of axial turbomachine. 3. Short axial distances S are favourable for maintaining high performance. The results show that the ∆P s and η s diminished monotonically with the axial distance increase, in other words no optimum in distance is found contrary to what was reported for highly compressible CRS [START_REF] Pundhir | A study of aerodynamic performance of a contra-rotating axial compressor stage[END_REF]. 4. The maximum efficiency η s, max could be steadily increased as the speed ratio varies in [1, 1.2]θ C . Furthermore, η s, max could be kept above 60% for θ ± [0.8, 1.6]θ C . 5. The power spectral densities of the wall pressure fluctuations between the rotors show that the dominant frequency corresponds to the blade passing frequency of the Rear Rotor and that the zone between the two counter-rotating rotors is a zone of high interactions between instationary flows. This feature should also be taken into account for better design and optimization.

  ∆P t Total pressure rise (Pa) ∆P s Static pressure rise (Pa) D Pipe diameter (mm) Q v Volume flow rate (m 3 .s -1 ) R tip Blade tip radius(mm) R hub Blade hub radius (mm) N Rotational speed of rotor (rpm) θ Rotational speed ratio N RR N FR L Distribution of load ∆P t, RR ∆P t, FR +∆P t, RR Z Number of blades P w Power consumption of the rotor (W) Z p Axial position of microphone S Axial distance between front rotor and rear rotor (mm) L chord Blade chord (mm) W Relative velocity (m.s -1 ) U Rotor velocity (m.s -1 ) C a Axial velocity (m.s -1 ) p ′ Wall pressure fluctuation (Pa) γ Stagger angle ρ a Density of air (kg.m -3 ) τ torque (N.m) µ Dynamic viscosity of the fluid (Pa.s)

FIGURE 1 .

 1 FIGURE 1. Sketch of the CRS. The bold arrow stands for the microphone (position Z p = 5 mm downstream the front rotor).

FIGURE 2 .

 2 FIGURE 2. Conception by MFT [5] of JW1 FR (•); JW2 FR ( ) and JW3 FR (♦). (a): radial profile of the stagger angle and (b): static pressure rise calculated by a semi-empirical model [5].

FIGURE 3 .

 3 FIGURE 3. Radial profiles of the stagger angle (•) and of the aerodynamic camber (♦, see Ref.[START_REF] Sarraf | Experimental study of blade thickness effects on the overall and local performances of a controlled vortex designed axial-flow fan[END_REF] for a definition) for the rotor JW1 RR .

4 . 1

 41 Overall performance of the three CRS working on design parametersThe static pressure rise and static efficiency as a function of volume flow rate, are presented in Fig.4, for the three systems. Closed symbols stand for the Front Rotors working alone and open symbols stand for CRS at their θ C and with S = 10 mm. The experimental data and the expected values at the nominal flow rate (Q v = 1 m 3 .s -1 ) are given in Table

FIGURE 4 .TABLE 3 .

 43 FIGURE 4. Overall performance of four systems: JW1(•), JW2 ( ) and JW3 (♦); JW1 FR (•), JW2 FR ( ) and JW3 FR ( ), design point,( ✡) (a): Static pressure rise ∆P s vs. volumetric flow rate Q v . (b): Static efficiency η s vs. volumetric flow rate Qv, at S = 10 mm.

FIGURE 5 .

 5 FIGURE 5. η s of JW1 (•), JW2 ( ) and JW3(♦) at S = 10, 20, 40, 100, 200, 250, 300 mm and Rront Rotor alone(S = 0, filled symbols) at design point Q v = 1m 3 .s -1 , θ = θ c

FIGURE 6 .

 6 FIGURE 6. Static efficiency of stage JW1 at Qv = 1 m 3 .s -1 , S = 10 mm, θ /θ c = [0.9 (•), 1, 0( ), 1.15(▽), 1.25(×)], θ c = 0.96.

FIGURE 7 .FIGURE 8 .

 78 FIGURE 7. Maximum static efficiency vs. θ . Stage JW1, at N FR = 1400 rpm (•), N FR = 2300 rpm ( ), S = 10 mm.

FIGURE 9 .

 9 FIGURE 9. Maximum Static efficiency of the stage JW1 at N FR = 2300 rpm, vs. θ /θ c for S = 10 mm (•), S = 20 mm ( ), S = 40 mm (♦), S = 50 mm (▽), S = 60 mm (×),S = 100 mm (+) and S = 200 mm (✡).

FIGURE 10 .

 10 FIGURE 10. PSD of the wall pressure fluctuations measured at Z p = 5 mm, for JW1 at N = 2300 -2200 rpm and S = 10mm, (a)Q v = 1 m 3 .s -1 ; (b)Q v = 0.6 m 3 .s -1 . (•): m f bp f , FR , (▽): n f bp f , RR , and ( * ): m f bp f , FR + n f bp f , RR with m = 0 and n = 0.

  FIGURE 11. PSD of the wall pressure fluctuations measured at Z p = 5 mm, for JW2 at N = 1800 -2600 rpm and S = 10mm, (a)Q v = 1 m 3 .s -1 ; (b)Q v = 0.6 m 3 .s -1 . (•): m f bp f , FR , (▽): n f bp f , RR , and ( * ): m f bp f , FR + n f bp f , RR with m = 0 and n = 0.

  FIGURE 12. PSD of the wall pressure fluctuations measured at Z p = 5 mm, for JW3 at N = 2600 -1100 rpm and S = 10mm, (a)Q v = 1 m 3 .s -1 ; (b)Q v = 0.6 m 3 .s -1 . (•): m f bp f , FR , (▽): n f bp f , RR , and ( * ): m f bp f , FR + n f bp f , RR with m = 0 and n = 0.

TABLE 1 .

 1 Design point for air at ρ a = 1.21 kg.m -3 .

TABLE 2 .

 2 Design parameters of the three CRS.

	rpm)	%
	JW1 2300/2200 0.96 41 10/7	large
	JW2 1800/2600 1.44 52 13/7	small
	JW3 2600/1100 0.42 23 10/7	large

TABLE 4 .

 4 Comparison of η s at design point Qv = 1 m 3 .s -1

  N CFR /N CRR PSD FR PSD RR Std(p ′ )

		(rpm)	(dB/Hz) (dB/Hz) (dB)
	JW1	2300/2200	17.4	36.5	40.9
	JW2	1800/2600	12.6	38.5	42.9
	JW3	2600/1100	20.2	29.8	35.1

TABLE 5 .

 5 The amplitude of the dominate frequency corresponded to f bp f , FR and f bp f , RR , for JW1, JW2 and JW3, at Q v = 1 m 3 .s -1 . Std(p ′ ) represents the power of the total signal.

Copyright c 2014 by ASME


Ewen Marechal
email: ewen.marechal@ensam.eusofiane.khelladi@ensam.eu

Sofiane Khelladi

Florent Ravelet
email: florent.ravelet@ensam.eu

Farid Bakir
email: farid.bakir@ensam.eu

Olivier Delepierre-Massue
email: olivier.delepierre-massue@snecma.fr

TOWARDS NUMERICAL SIMULATION OF SNOW SHOWERS IN JET-ENGINE FUEL SYSTEMS

Keywords: Finite volume, Simple algorithm, Moving Least Squares, Particle, Two-way coupling

Aircraft fuel systems are subject to icing at low temperatures. If the flow rate is increased, sudden releases of large quantities of ice may occur, called "snow showers". They threaten the safety of flights and have been the subject of several investigations over past years. Jet engines fuel system components may sensitive to clogging. When a snow shower happens, ice particles settle in seconds, forming a porous layer. Modelling such events involves transient hydraulics and solid dynamics. We propose to investigate numerically the dynamics of transient particle clogging. Equations of motion for the incompressible fluid phase are discretized in a high-order finite-volume context and solved using a pressure-based algorithm. The discrete phase is modelled in a Lagrangian frame. Contacts between solids are handled by a dedicated algorithm. Solid volume fraction is calculated in regions occupied by particles. Finally, two-way coupling is achieved by source terms for momentum exchange, viscous and inertial loss. 2D simulation of the clogging of an ideal filter is performed.

INTRODUCTION

Ice in flight can be disastrous. A thin layer on the airfoil degrades its ability to lift and increases its drag, leading the airplane to stall even at high speeds. Ice collected on wing can also broke off and be sucked into the engines, damaging the fan blades. An airplane can roll or pitch uncontrollably if icing occurs on the horizontal or vertical stabilizer. Obstruction of pitot tube will cause the airspeed indicator to give wrong information. These events are related to weather (moisture, snow, freezing rain, and drizzle) and caused hundreds of casualties over decades.

However, fuels used in jet aircraft contain significant quantities of water. Investigation of the Boeing 777-236ER, G-YMMM crash at Heathrow in 2008 revealed that this water can create ice structures in aircraft fuel ducts at low temperatures and moderate flow rates. If the latter is increased, shear stress may tear off the ice resulting in high concentrations of particles in the fuel. Such sudden releases of large amount of ice are often referred as a "snow showers" [START_REF] Sleight | Report on the accident to Boeing 777-236ER, G-YMMM[END_REF]. 1

Corresponding author

SimHydro 2014:Modelling of rapid transitory flows, [START_REF] Marechal | Dispositif de simulation d'une introduction de paquets de glace dans un moteur[END_REF][START_REF] Clift | Bubbles, Drops and Particles[END_REF][START_REF] Lefebvre | Modélisation numérique d'écoulements fluide/particles[END_REF] June 2014, Sophia Antipolis -Marechal E., Khelladi S., Ravelet F., Delepierre-Massue O,. Bakir F -Towards numerical simulation of snow showersin jet-engine fuel systems.

Jet engines have complex fuel systems, involving several hydraulic components sensitive to clogging e.g. filters, valves or heat exchanger inlet screens. When a snow shower occurs, ice particles settle in seconds to form a porous layer which is likely to occasion fuel flow restrictions. Nowadays, aviation safety authorities require that fuel systems must be designed so that any accumulation of water which may separate from the fuel will not cause engine malfunctioning [START_REF]Certification Specification for Engines CS-E, Amendment 3[END_REF]. Hence, understanding of such flow has become a rising matter of interest among jet engines manufacturers over last years.

In this paper, we propose a numerical approach to investigate dynamics of transient clogging of hydraulic components. Simulation strategies for particulate flows can be roughly classified in two categories. The first one is the Eulerian-Eulerian approach. The fluid phase and the dispersed phase are considered as interpenetrating continuum and updated using Mass, momentum and energy conservation laws. Additional constitutive equations may be required to describe interactions between particles. This approach is well suited for very small particles moving with the fluid, but is not adapted to packed beds involved in clogging.

Therefore we chose the second approach known as Eulerian-Lagrangian. Each particle, or group of similar particles (called parcel) is considered with its own position and velocity which are updated with the equation of motion. The effect of particles on the flow field is accounted through momentum and energy sources. Collisions between particles can be treated deterministically. Therefore it is possible to accurately describe how they stack on each other [START_REF][END_REF].

MODELING SNOW SHOWER 2.1 A multiphysic, multiscale phenomenon

Snow shower are tree-stage events. The first step is a slow (several hours) growth of ice crystal within the aircraft fuel pipes. Then, accumulated snow is unexpectedly released and carried by the fuel into the jetengine fuel system. At last, ice particles will choke fuel systems components in seconds.

Ice accretion onto sub cooled surfaces

Icing is influenced by several parameters, e.g. surface properties, temperature variations of fuel and fluid velocity. Laboratory observations show that accretion process involves the combination of growth of ice crystal from the water dissolved in the fuel (frost) and deposit of supercooled droplets (rime) [START_REF] Baena | Cold fuel test rig to investigate ice accretion on different pump inlet filter-mesh screens[END_REF][START_REF] Murray | Supercooling of water droplets in jet aviation fuel[END_REF]. Cristal growth is driven by thermodynamics: higher molecular diffusion favours adsorption of water molecules to the growing crystal while thermal diffusion removes latent heat released by water solidification [START_REF] Hobbs | Ice Physics, Oxford Classic Texts in the Physical Sciences[END_REF].

Ice shedding

The mechanical strength of an ice sample at a given temperature depends of the shape and size of the specimen, the way ice is formed, and how the load is applied. The fracture is initiated by the formation of a crack, depending on the crystallographic orientation, which then spread through the specimen. A generally fracture law has yet to be agreed upon, but correlations were proposed by several authors. Ice present weak tensile and shear strengths [START_REF] Butkovich | The ultimate strenght of ice[END_REF][START_REF] Jellinek | The influence of imperfections on the strength of ice[END_REF]. Investigations following the 2008 accident demonstrated that release of accretion occurs when the fuel velocities is increased to the range of 0.7 to 0.8 m/s. Hence, ice shedding is likely caused by shear stress resulting from a drag increase [START_REF] Reid | Engine fuel system tolerance to fuel born ice. Managing water and ice in aviation fuel under low temperature conditions[END_REF].

Clogging

Components with small holes are obviously threatened by ice particles. This includes filters, strainers, injectors and heat exchanger inlets. It is usual to provide by-pass means to such components so that if they were to be obstructed, fuel will continue to flow at an acceptable rate through the rest of the system. Accumulation of a significant quantity of ice will form a thin layer. Resulting pressure drop may cause bypass valves to open. Therefore it is critical to know in which extent ice particles will follow the flow through the by-pass and propagate in the fuel system. Given a carrier fluid, the dynamic is driven by the shape, size and density of ice particles. According to recent studies, ice particles size is in the range of to 1-10 μm and may be considered almost neutrally buoyant in fuel [START_REF] Lam | WAFCOLT -Water behavior in aviation fuel under low temperature conditions[END_REF].

Replicating snow showers

The European Aviation Safety Agency (EASA) issues certification specifications to ensure the tolerance of fuel systems to high concentration of ice particles. CS-E 560 states: «The fuel system must be designed so that any accumulation of likely quantities of water which may separate from the fuel will not cause Engine malfunctioning. » [START_REF]Certification Specification for Engines CS-E, Amendment 3[END_REF].

Experimental setup

In order to comply with international safety standards, an experimental setup was designed by the French jetengine manufacturer SNECMA. The equipment used to make such tests is described in more detail previously [START_REF] Marechal | Dispositif de simulation d'une introduction de paquets de glace dans un moteur[END_REF]. Briefly, pure water is introduced in cold kerosene fuel flow through a homemade injector. The system regulation ensure that a given quantity of water is injected at the desired volumetric flow rate (around 10 000 ppm). It was observed that injected water froze almost immediately when the fuel temperature is below -20°C.

Collected ice

80mL of water were introduced in a -45°C jet-engine fuel and ice was collected downstream by a 25 μm filter. The flow decreased by about 20% though the filter appeared to be completely obstructed. Figure 1 shows the layer of ice accumulated on the filter and a close-up on ice particles. Ice was collected as dry as possible from the deposit and left to melt. A volume fraction up to 50% of fuel was still observed. This implies that permeability of the ice layer is a combination of compactness from particle stacking and porosity of particles themselves.

Assumptions

Complete modelling of a snow shower is outside of this research scope. This paper focuses on the clogging of system components by ice, which is the most critical step of the phenomenon. Furthermore, we mean to keep the model close to what have been observed in previously described tests in the perspective of future validation of the method. Hence, several assumptions are made in this study:

1. Supercooling droplets are not considered. 2. Injected water is considered as already frozen. We do not consider any freezing or melting.

NUMERICAL METHOD

Numerical simulation of solid dynamics

Particle equation of motion

Let consider a rigid body of mass m and inertia matrix J submitted to total forces F and torques T . Its center of mass has for coordinate and orientation  

W X,

. Its motion is governed by the following equations:

         T W J F X dt d dt d m   , (1)
Modelling total forces F and torques T applied on a solid particle in an infinite medium has been studied

extensively. Various forces are involved: drag, buoyancy, Magnus, Suffman and Basset forces, virtual mass force. In this study we only consider drag D F , buoyancy B F and contact forces C F . Rotation of the particle is discarded. Let n denote the simulation time step of length t  . The discretized equation of motion for a particle is:

  2 1 1 2 1 t m t n n n n          C B D F F F X X X  , (2)

Drag in a multiparticle system

In the following, subscript M denotes the mixture, C the fluid continuum, and P the particle. Drag arises from the relative motion R V between the particle and the flow:

C P R X X V     , (3)
A standard relation for drag is force is:

R R P C D S C V V F D  2   , (4)
In a multiparticle system, the effect of others particles on the drag force is modelled as an increase of viscosity resulting from additional stress in the flow. We use the mixture viscosity model proposed by Clift & Al [START_REF] Clift | Bubbles, Drops and Particles[END_REF]:

  5 . 2 1       C M , (5)
Where  denotes the mixture viscosity, and  the fraction of volume occupied by particles. Two similarity hypotheses are introduced. First, the particle Reynolds number is defined as:

M C P P C P d   X X     Re , ( 6
)
Then it is assume that a similarity exists between a single-particle system and a multi-particle system. The relation for the drag coefficient is given by:

          10 Re 1 Re 24 75 . 0 P P D C , (7)
This equation holds for particle Reynolds number lesser than 1000. Above this value, the drag coefficient D C should be set to a constant value of 0, 44 [12]. It is also necessary to set a maximum value for D C .

Otherwise, particle velocity may exceed fluid velocity at the end of a step. Maximum drag coefficient value is related to the time step. It has been found that sedimentation speed is sensitive to this parameter.

Particle interactions

The mixture viscosity model only holds for moderate solid volume fractions. In dense flow and packed bed we need to take into account contact forces between particles. In this paper, we derive a discrete contact model from equation of motions. Let i and j denotes two particles of radius i r and j r . Let n i X and n j X be their centers of mass of mass. We define the directional unit vector between i and j as:

  1     n i n j n i n j X X X X E ij , ( 8
) Let ji ij F F  
be a force exerted by particle j on particle i over a time step. This force is directed along ij E .

Writing equation of motion (2) for each particle and taking into account ij F give the following system:

                                   ij j n j n j ij j j n j n j n j ij i n i n i ij i i n i n i n i m t t m t t m t t m t t F V X F F X X X F V X F F X X X 2 2 2 2 2 2 1 2 2 1   , ( 9
)
Where n i V ˆ and n j V ˆ are called "a priori" velocity of particles i and j. These velocities can be such that if no force is exerted between particles, there may be interpenetration.

    ij ij n i n j n i n j n i n j m t t F V V X X X X 2 ˆ2 1 1           , (10)
Where

  1 1 1      j i ij m m m
.Now we perform the dot product of (10) with the directional vector ij E :

      ij ij n i n j n i n j n i n j m t t F E V V E X X E X X ij ij ij 2 ˆ2 1 1              , (11)
Let denote n ij l the projected distance between particle centers of mass and n ij V ˆ the relative a priori velocity.

To guarantee that no interpenetration occurs at the end of a time step, we need to set a constraint

j i n ij r r l   1
. Introduction of this constraint in equations [START_REF] Marechal | Dispositif de simulation d'une introduction de paquets de glace dans un moteur[END_REF] give:

0 2 ˆ2         ij ij ij n ij j i n ij m t t r r l F E V , (12)
According to the velocity projection method introduced by Lefebvre in [START_REF] Lefebvre | Modélisation numérique d'écoulements fluide/particles[END_REF], finding the non-overlapping velocity field is equivalent to minimizing the following functional:

              j i ij n ij j i n ij ij n t r r l J E V V V V , 2   , (13)
Solution of this saddle-point problem is obtained iteratively by the Uzawa algorithm:

                             t r r l t ij n ij j i n ij k ij k ij j i ij k ij n i k i E V E V V , 0 max 1 1     , ( 14
)
The real is the non-dimensional artificial time parameter. Comparison of equations ( 12) and [START_REF] Khelladi | Toward a Higher-Order Unsteady Finite Volume Solver Based on Reproducing Kernel Particle Method[END_REF] shows the relation between contact force and Lagrange multipliers:

ij ij ij m t F 2     , (15)
       Ns j j n i j n n i n u x N x u 1 ˆ, (17)
We refer the reader to [START_REF] Khelladi | Toward a Higher-Order Unsteady Finite Volume Solver Based on Reproducing Kernel Particle Method[END_REF] for a complete description of the computation of MLS shape functions and derivatives.  the fluid density (kg.m -3 ), and p the pressure (Pa). In order to achieve two-way coupling between particles and fluid, we replace the fluid viscosity by previously defined mixture viscosity and introduce a source term V S to account for viscous loss within the porous layer [15]:

Navier-Stokes equations

U S V   M  , (18)
Let P d denotes the mean particle diameter and the volumetric solid fraction. Parameter  is defined:

  2 3 2 150 1      P d , (19)
The continuity and momentum conservation equations integrated over a control volume  can be written:

0 0      d n U  , ( 20
)                                d d d p d d t M V S n U n n U U U     0 0 , (21)
Equations ( 20) and ( 21) are discretized using a novel high-order formulation based on MLS approximants introduced by Ramirez & Al solved. The solution procedure is based on the SIMPLE algorithm [START_REF] Patankar | A Calculation Procedure for Heat, Mass and Momentum Transfer in Three Dimensional Parabolic Flows[END_REF]. We refer the interested reader to [START_REF] Ramirez | A new higher-order finite volume method based on Moving Least Squares for the resolution of the incompressible Navier-Stokes equations on unstructured grids[END_REF] for a complete description of the method.

APPLICATION: CLOGGING OF A FILTER

By way of application we consider the simple 2D case of a filter clogged by a snow shower. An ideal filter (permeable to the flow but blocking any kind of particle) is located in the middle of rectangular domain representing a segment of pipe. Domain dimensions are 2.5 cm (the diameter) by 10 cm (the length).

As inlet condition, we set the fluid velocity to 0.1 m/s and up to 1000 particles are seed in order to have a concentration of water close to 10000 ppm. The carrier fluid has a density of 850 kg/m3 and a viscosity of 1.10-3 Pa.s. Ice particles have a density of 917 kg/m3 and porosity is set to 0.5. We set the time step to 1ms and perform up to 3000 iterations. Figure 3 shows the evolution of the velocity field U and solid fraction  over time. Particles seem to be by influenced by buoyancy, though this could be due to the maximum value imposed on drag coefficient. The deposit starts building on the bottom of the pipe, but contact forces push particles upward until equilibrium is reached. The velocity field is disturbed by particles, and a global decreased of the velocity is observed as particles settle to form a porous layer.

CONCLUSION

We proposed a method to simulate clogging of components by snow showers. A standard formulation is used to describe particle dynamics and retroaction of particles on the flow. Contacts forces between particles are introduced to describe the behavior of the packed bed formed by accumulating particles. A first simulation performed on the simple case of a 2D ideal filter in a pipe is in qualitatively consistent with experimental observations of snow shower dynamics. However, the method has to be extended in 3D and confronted quantitatively to available experimental data. The novel high order formulation used for computation of fluid dynamics should allow handling of large and complex geometries.

SimHydro 2014 :

 2014 Modelling of rapid transitory flows,11-13 June 2014, Sophia Antipolis -Marechal E., Khelladi S., Ravelet F., Delepierre-Massue O,. Bakir F -Towards numerical simulation of snow showersin jet-engine fuel systems.

Figure 1 :

 1 Figure 1: Left: clogging of a filter -Right: picture of particles

3 .

 3 Particles are considered already released and created at flow inlet. 4. Particle initial velocity is set equal the fuel one. 5. Particles radius is evenly distributed within the range of 150-250 μm. 6. Particles porosity is in the range of 0.1 to 0.5. 7. The carrier fluid is incompressible. 8. We restrict our study to low Reynolds flow.

SimHydro 2014 :

 2014 Modelling of rapid transitory flows,11-13 June 2014, Sophia Antipolis -Marechal E., Khelladi S., Ravelet F., Delepierre-Massue O,. Bakir F -Towards numerical simulation of snow showersin jet-engine fuel systems.

SimHydro 2014 :

 2014 Modelling of rapid transitory flows,11-13 June 2014, Sophia Antipolis -Marechal E., Khelladi S., Ravelet F., Delepierre-Massue O,. Bakir F -Towards numerical simulation of snow showersin jet-engine fuel systems.

Figure 2 :F

 2 Figure 2: Particles initial position, directional vector, and a priori velocities and particles final positions given

Fluid

  motion is governed by the well-known Navier-Stokes equations. Let

SimHydro 2014 :

 2014 Modelling of rapid transitory flows,11-13 June 2014, Sophia Antipolis -Marechal E., Khelladi S., Ravelet F., Delepierre-Massue O,. Bakir F -Towards numerical simulation of snow showersin jet-engine fuel systems.

Figure 3 :

 3 Figure 3: 2D simulation of transient clogging of a filter. Left: solid volume fraction field. Right: velocity magnitude field. Snapshots taken at time = 0.5, 1.0, 2.0 and 3.0 seconds.

Figure 4

 4 Figure4shows the network of contact forces between particles. Because diameters are randomly generated between of 150 and 250 μm, the stack is irregular; Therefore some particles close to the filter are under the full load of the drag cumulated by the stack of particles while other are not.

Figure 4 :

 4 Figure 4: 2D simulation of transient clogging of a filter. Network of contact forces between particles.

ACKNOWLEDGEMENTS

This research was carried out within the framework of a joint supervised doctoral thesis (grant #2011/0411). The first author is grateful to SNECMA for funding and would like to thank colleagues of DynFluid laboratory for their support.

Delepierre-Massue O,. Bakir F -Towards numerical simulation of snow showersin jet-engine fuel systems.

Numerical simulation of fluid dynamics

Flow simulation is performed using a homemade finite volume code. The geometry is meshed so that computational cells remain greater than particles. This is essential to get cell-averaged quantities as solid volume fraction, but can lead to rather coarse meshes.

Moving Least Squares approximation

In order to have a good accuracy of the solution, we use a high order formulation based on Moving Least Squares (MLS) approximants, a technique widely used in the meshless community. Let ) (x u be a function. Its MLS approximation ) ( ˆx u is computed from a set of Ns neighbouring points called stencil where the value j u of ) (x u is known. The number of neighbours depends of the order of accuracy and the stencil should be as compact as possible, centred around the node. The MLS approximation ) ( ˆx u is written:

High order approximate derivatives can be expressed in terms of the derivatives of the MLS shape function.

The n-th derivative of ) (x u can be expressed:


Pascal Audigane

Anne-Gaëlle Bader

Sylvie Gentier

Laurent Beccaletto

Gaël Bellenfant

The role of the underground for massive storage of energy: a preliminary glance of the French case

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

The role of the underground for massive storage of energy: a preliminary glance of the French case Audigane P., Bader A., Gentier S., Bellenfant G., Beccaletto L.

The question of storing energy in France has become of primary importance since the launch of a road map from the government which places in pole position this topic among seven major milestones to be challenged in the context of the development of innovative technology in the country. The European objective to reach 20% of renewables in the energy market, from which a large part would come from wind and solar power generation, raises several issues regarding the capacity of the grid to manage the various intermittent energy sources in line with the variability of the public demand and offer. These uncertainties are highly influenced by unpredictable weather and economic fluctuations. To facilitate the large-scale integration of variable renewable electricity sources in grids, massive energy storage is needed. In that case, electric energy storage techniques involving the use of underground are often under consideration as they offer a large storage capacity volume with a adapted potential of confining and the space required for the implantation. Among the panel of massive storage technologies, one can find (i) the Underground Pumped Hydro-Storage (UPHS) which are an adaptation of classical Pumped Hydro Storage system often connected with dam constructions, (ii) the compressed air storage (CAES) and (iii) the hydrogen storage from conversion of electricity into H2 and O2 by electrolysis.

UPHS concept is based on using the potential energy between two water reservoirs positioned at different heights. Favorable natural locations like mountainous areas or cliffs are spatially limited given the geography of the territory. This concept could be extended with the integration of one of these reservoirs in an underground cavities (specifically mined or reuse of preexisting mines) to increase opportunities on the national territory. Massive storage based on compression and relaxation of air (CAES) requires high volume and confining pressure around the storage that exists naturally in the underground and which increases with depth. However, the move to an interesting efficiency requires that the heat generated during compression can be stored and used during expansion. This storage can be also underground. H2 underground storage is part of the "Power to gas" concept which allows for converting electricity into a gas available for either electrical or gas grid. Each of these techniques requires the selection of appropriate geological formations which contains specific characteristics in agreement with several criteria under consideration when choosing electric energy storage methods for application (lifetime, life cycle, discharge rate, environmental impact, public acceptance …). We propose in this paper a preliminary review of the potential massive electric energy storage capacities in France of using specific geological formations (salt, basement) and the various physical phenomena linked to the couple geology/technology. Several approaches and methodologies developed formerly with other applications (geothermal, CO2 storage, heat storage …) will be used to investigate rock mechanical integrity, geochemical and thermal environmental impacts associated to these innovative technologies.