David Marx
email: david.marx@univ-poitiers.fr

Yves Aurégan
email: yves.auregan@univ-lemans.fr

Effect of turbulent eddy viscosity on the unstable surface mode above an acoustic liner

Keywords: duct acoustics, acoustic liner, instability PACS: 43.28.Py, 43.20.Mv, 43.28.Ra

Lined ducts are used to reduce noise radiation from ducts in turbofan engines.

In certain conditions they may sustain hydrodynamic instabilities. A local linear stability analysis of the flow in a 2D lined channel is performed using a numerical integration of the governing equations. Several model equations are used, one of them taking into account turbulent eddy viscosity, and a realistic turbulent mean flow profile is used that vanishes at the wall. The stability analysis results are compared to published experimental results.

Both the model and the experiments show the existence of an unstable mode, and the importance of taking into account eddy viscosity in the model is shown. When this is done, quantities such as the growth rate and the velocity eigenfunctions are shown to agree correctly.

Introduction

Lined ducts are widely used to reduce noise radiation from ducts, for example in air conditioning systems and turbofan engines. In the presence of flow, the interaction of sound and flow above the liner can be complex.

One of the main feature of this interaction are nonlinear effects, which have an effect on peak absorption frequency and on sound transmission when the Mach number increases [START_REF] Guess | Calculation of perforated plate liner parameters from specified acoustic resistance and reactance[END_REF][START_REF] Motsinger | Design and performance of duct acoustic treatment, Aeroacoustics of Flight Vehicles[END_REF]. Another possible effect is the presence of an instability over the liner, which is the subject of the present paper.

In some circumstances, sound amplification by a liner has been observed [START_REF] Meyer | Experiments on the influence of flow on sound attenuation in absorbing ducts[END_REF][START_REF] Brandes | Sound Amplification in flow ducts lined with a periodic sequence of resonators[END_REF][START_REF] Aurégan | Experimental evidence of an instability along an impedance wall with flow[END_REF]. This has been the case mostly for plane waves, for low-resistance liners, and for a high enough mean flow velocity. This amplification has been ascribed to an underlying flow instability. In recent experiments [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF] the flow has been measured over a liner in such a situation when the liner amplifies sound at some frequency close to the resonance frequency of the liner.

The presence of a spatial instability has been confirmed, and quantites such as growth rate, spatial wavelength, and eigenfunctions have been provided. Some numerical evidence of instability over a liner also exist [START_REF] Li | Time-domain impedance boundary conditions for surfaces with subsonic mean flows[END_REF][START_REF] Burak | Validation of a Time and Frequency Domain Grazing Flow Acoustic Liner Model[END_REF]. Many analytical works have been performed on propagation modes in lined ducts with flow [START_REF] Tester | The propagation and attenuation of sound in lined ducts containing uniform or plug flow[END_REF][START_REF] Koch | Eigensolutions for liners in uniform mean flow ducts[END_REF][START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF][START_REF] Brambley | Classification of aeroacoustically relevant surface modes in cylindrical lined ducts[END_REF]. They often rely on a constant mean flow (or plug flow) assumption. As a result, there is a sheet of vorticity at the wall, and it becomes necessary to introduce a model for enforcing the liner impedance at the wall. A widely used model is that by Myers, which enforces the continuity of particle displacement across the vorticity sheet [START_REF] Ingard | Influence of fluid motion past a plane boundary on sound reflection, absorption, and transmission[END_REF][START_REF] Myers | On the acoustic boundary condition in the presence of flow[END_REF]. Us-ing this model, Tester [START_REF] Tester | The propagation and attenuation of sound in lined ducts containing uniform or plug flow[END_REF] showed that the response of a lined duct to an impulse contains some modes that are amplified during their propagation. This instability was found to be of the convective type. In his classification of modes in lined duct with flow, Rienstra [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF], besides acoustic modes, also predicts some unstable hydrodynamic surface modes. Though, the effect of the Myers boundary condition has been questioned [START_REF] Tester | The propagation and attenuation of sound in lined ducts containing uniform or plug flow[END_REF][START_REF] Rienstra | Spatial instability of boundary layer along impedance wall[END_REF][START_REF] Brambley | Fundamental problems with the model of uniform flow over acoustic linings[END_REF]. In particular, Brambley [START_REF] Brambley | Fundamental problems with the model of uniform flow over acoustic linings[END_REF] has shown that the Myers' boundary condition is ill-posed. The use of classical stability criteria such as Briggs-Bers' [START_REF] Briggs | Electron-Stream Interaction with Plasmas[END_REF][START_REF] Huerre | Hydrodynamic Instabilities in Open Flows, Hydrodynamics and Nonlinear Instabilities[END_REF] then becomes invalid. Such a criterion was used by Rienstra [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF] to prove the existence of spatially amplified waves. This renders the conclusions about the stability of flow in lined flow ducts fragile. Fortunately, it is possible to avoid the use of the Myers' boundary condition: one has to use a velocity profile that accounts for the boundary layer and vanishes to zero at the wall. It is then not possible to tackle the problem analytically and one has to rely on a numerical solver to calculate the eigenmodes [START_REF] Shankar | Acoustic refraction and attenuation in cylindrical and annular ducts[END_REF][START_REF] Agarwal | Acoustic wave propagation in a pipe with fully developed turbulent flow[END_REF][START_REF] Kousen | Eigenmodes of ducted flows with radially-dependent axial and swirl velocity components[END_REF][START_REF] Rienstra | Spatial instability of boundary layer along impedance wall[END_REF][START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF][START_REF] Boyer | Theoretical investigation of hydrodynamic surface mode in a lined duct with sheared flow and comparison with experiments[END_REF], which will be the case here.

The effect of dissipative mechanisms on acoustic propagation and the resulting acoustic damping has been the object of some researches. This effect is well known in the absence of flow, where the wavenumber becomes complex with an imaginary part given by Kirchoff's formula [START_REF] Davies | Measurement of plane wave acoustic fields in flow ducts[END_REF]. The presence of flow complicates the matter somewhat [START_REF] Peters | Damping and reflection reflection coefficient measurements for an open pipe at low Mach and low Helmholtz numbers[END_REF][START_REF] Howe | The damping of sound by wall turbulent shear layers[END_REF], and there are still different models for the wavenumbers with flow when dissipation is present [START_REF] Davies | Measurement of plane wave acoustic fields in flow ducts[END_REF][START_REF] Dokumaci | A note on transmission of sound in a wide pipe with mean flow and viscothermal attenuation[END_REF][START_REF] Allam | Investigation of damping and radiation using full plane wave decomposition in ducts[END_REF].

Few works have focused on the effect of viscosity over the modes above an acoustic liner, even is some studies have considered the effect of viscosity on the nature of the boundary condition at a liner [START_REF] Aurégan | Influence of grazing flow and dissipation effects on the acoustic boundary conditions at a lined wall[END_REF]. The effect of turbulent eddy visosity on the other hand has been taken into account in several studies dealing with the mechanics of an organized wave in a turbulent shear flow [START_REF] Reynolds | Stability of turbulent channel flow, with application to Malkus's theory[END_REF][START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3: Theoretical models and comparisons with experiments[END_REF][START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF], and the same kind of formalism will be used here.

In the present paper, a local 2D linear stability analysis is performed to investigate the presence of an unstable mode over a lined wall in a duct with flow, and the effect of viscous dissipation on the mode is taken into account, which is new to the author knowledge. The different parameters are taken to be identical to those in reference [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF], so that the models can be tested against experiments. A brief summary of the experimental results are recalled in section 2. The modelling of the turbulent flow is provided. Several models for the stability analysis are then presented in section 3. These are an inviscid compressible model (Pridmore-Brown), an inviscid incompressible model (Rayleigh), and a viscous incompressible model which takes into account both the molecular and the turbulent eddy viscosity. The numerical solution method is also presented. Section 4 presents the results for the inviscid compressible and incompressible models. These models predict absolute instability and do not match the experiments. Section 5 presents the results for the incompressible but viscous model. This model predicts spatial instability and provides a correct match with the experiments. Finally conclusions are given.

Summary of experimental results

In Marx et al [START_REF] Marx | Evidence of Hydrodynamic instability over a liner in a duct with flow[END_REF][START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF] an instability has been measured in a lined flow duct.

The setup and main results are now briefly recalled, and the modelling of the mean turbulent flow is presented.

A 2D cut of the experimental configuration is shown Fig. 1. The channel has a height 2 H=0.02m, with a rigid top wall and a partly-lined bottom wall. In following the tilde˜is used to indicate physical quantities, while no tilde is used for normalized quantities. Quantities used for normalizing are introduced in the next section. The duct width (0.08 m) gives an aspect ratio of 4 which, while not very large, we assume is large enough for a 2D model to be suitable. The liner on the bottom wall covers the whole span of the section

Y = -iQ tan Ãω + i 2 = -iQ tan Aω + i 2 (1)
with =0.3, Q=0.8, Ã=2.5 x 10 -4 s, A = Ãc 0 / H=8.5. For the stability analysis performed below, the expression for the impedance needs to be continued in the complex ω-plane, and this continuation is valid for the present model [START_REF] Rienstra | Impedance models in time domain including the extended Helmholtz resonator model[END_REF].

The flow velocity Ũ (ỹ) is along the x-direction and depends on ỹ. The velocity profile has been measured at several stations above the liner and its shape changes with streamwise location, with a displacement thickness changing from about 0.9mm to 1.5mm. While the profile just upstream of the liner (above the rigid wall) does not represent the flow at all streamwise stations above the liner, it will be used for the stability analysis. This profile is given in Fig. 2. The section-averaged mean Mach number is defined as:

M 0,avg = 1 2 1 -1 M 0 (y)dy , (2)
where M 0 (y) is the Mach number. Its value is: M 0,avg =0.27.

When performing the stability analysis, it is desirable to have an analytic form of the velocity profile. In Marx and Aurégan [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF], the profile was fitted by using polynomials: a quadratic polynomials in the near wall region, smoothly connected to a quadratic one in the central region. Such a fit ignores the physics of the flow in the near wall region. This region is well known to exhibit a viscous sublayer, a buffer layer, and a logarithmic layer.

Ignoring these regions seems irrelevant when using a turbulent eddy viscosity model as will be done below. Thus a profile physically acceptable in the near wall region will be used in the following even if it is not as good as the fit in [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF] for the central region of the channel. Note that in the literature the flow in the wall region has been measured mainly for flow over rigid or possibly rough surfaces, but not for the kind of wall corresponding to the acoustic liners. This is another reason for using the velocity profile just upstream of the acoustic material, since it is is measured over a rigid wall.

The mean flow used in the present work is obtained by integrating the following expression [START_REF] Reynolds | Stability of turbulent channel flow, with application to Malkus's theory[END_REF][START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF]:

dU + 0 dy = - Re τ y ν + T (y) (3)
The + superscript indicates a normalization in wall units. In such a normalization the velocity is normalized by the wall friction velocity, ũτ ; the total viscosity, νT (y), is normalized by the physical viscosity, ν; and the Reynolds number is given by Re τ = ũτ H/ν. The total viscosity is the sum of the physical viscosity and of the turbulent eddy viscosity, νt , that is:

νT = ν + νt or ν + T = 1 + ν + t
It is given by Cess' expression [START_REF] Reynolds | Stability of turbulent channel flow, with application to Malkus's theory[END_REF][START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3: Theoretical models and comparisons with experiments[END_REF][START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF]:

ν + T (y) = 1 2 + 1 2 1 + κ 2 Re 2 τ 9 (1 -y 2 ) 2 (1 + 2y 2 ) 2 1 -e (|y|-1)Reτ /A V D 2 (4)
κ=0.42 is the von Karman constant, and A V D =25.4 is a constant. This expression is valid for Re τ up to a few thousand, which is the case here. This model is based on a turbulent eddy viscosity using a Prandtl mixing layer hypothesis and a Van Driest correction. It includes the viscous sublayer, the buffer layer, the logarithmic layer and a wake correction at the center of the channel. These features have also been used by Agarwal and Bull [START_REF] Agarwal | Acoustic wave propagation in a pipe with fully developed turbulent flow[END_REF] in their study of acoustic modes in a cylindrical duct, while with a slightly different model. These authors did not include the effect of viscosity on the perturbation though. The previous expressions allow for a degree of freedom, which is Re τ , or likewise, ũτ . These are determined so that the experimental profile and the analytical profile correspond to the same sectionaveraged Mach number of 0.27, which provides: Re τ =2760, corresponding to ũτ =4.14 m.s -1 . The analytical velocity profile obtained by integrating Eq. ( 3) together with Eq. ( 4) using this value for Re τ is given in Fig. 2. It matches correctly the experimentally measured profile. The total eddy viscosity is given in Fig. 3. It is seen that the turbulence is generating viscosity more than 200 times that of the molecular viscosity.

For a high enough flow speed, the acoustic transmission by the liner of acoustic waves in a frequency range around the liner resonance frequency was abnormally large [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF]. Some phase-locked laser measurements of the flow were performed and structures convected by the flow and growing spatially were observed above the liner. An example of fluctuating axial velocity, ũ , and fluctuating transverse velocity, ṽ , at one phase is shown in Fig. 4. This was obtained for a section-averaged mean Mach number M 0,avg =0.27 and an incident acoustic wave with frequency fa =1060Hz at sound pressure level L i = 133 dB. The normalized fluctuating velocity could be put into the form of a travelling wave amplifying in the downstream direction. That is:

u (x, y, t) = û(y, α, ω a )e i(αx-ωat) v (x, y, t) = v(y, α, ω a )e i(αx-ωat) (5)
Here α = α r + iα i is the complex wavenumber (-α i is the growth rate). The normalized angular frequency ω a corresponds to the particular frequency of the exciting incident acoustic wave (ω a ∼ 0.196). The eigenfunctions û and v can also be written:

û(y) = |û(y)|e iφu(y) v(y) = |v(y)|e iφv(y) (6)
In Marx et al [START_REF] Marx | Evidence of Hydrodynamic instability over a liner in a duct with flow[END_REF][START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF], the wave number has been given to be α = 1.3 -0.3i

(corresponding to αr ∼130 m -1 and αi ∼ 30 m -1 ), and the eigenfunctions |û(y)| and |v(y)| have also been given. These will be compared to the result of the stability analysis in the following.

Model and stability analysis

Compressible inviscid model

A usual model for investigating the modes in duct acoustics are the compressible Euler, the mass conservation, and the entropy conservation equations, assuming isentropic movement [START_REF] Tester | The propagation and attenuation of sound in lined ducts containing uniform or plug flow[END_REF][START_REF] Koch | Eigensolutions for liners in uniform mean flow ducts[END_REF][START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF][START_REF] Brambley | Classification of aeroacoustically relevant surface modes in cylindrical lined ducts[END_REF]:

∂ ρ ∂ t + ∇ • (ρũ) = 0 ρ ∂ũ ∂ t + ρũ • ∇ũ = -∇p dp d t = c2 0 dρ d t (7)
where ũ = (ũ, ṽ) is the velocity vector, ρ is the density, p is the pressure.

According to the previous section, a 2D model is considered. In order to perform a linear stability analysis, these equations are linearized. For this we introduce the following splitting of the variables: ũ(x, ỹ, t) = Ũ0 (ỹ) + ũ (x, ỹ, t); p(x, ỹ, t) = p0 + p (x, ỹ, t) ; ρ(x, ỹ, t) = ρ0 + ρ (x, ỹ, t). The subscript 0 refers to the base flow and the prime represents a fluctuation upon this base flow. The mean base flow Ũ0 = ( Ũ0 (ỹ), 0) is in the x-direction. It is allowed to depend on ỹ so that a mean shear is possibly taken into account. The set of linearized equations for the perturbations ũ , ṽ , p can be reduced to:

∂ ũ ∂ t + Ũ0 ∂ ũ ∂ x + ṽ d Ũ0 dỹ = - 1 ρ0 ∂ p ∂ x ( 8
)
∂ṽ ∂ t + Ũ0 ∂ṽ ∂ x = - 1 ρ0 ∂ p ∂ ỹ (9) ∂ p ∂ t + Ũ0 ∂ p ∂ x + ρ0 c2 0 ∂ ũ ∂ x + ∂ṽ ∂ ỹ = 0 (10)
For every variable a modal solution of the form ũ (y) = û(ỹ, α, ω)e i( αx-ω t) is introduced, where the wavenumber, α = αr + iα i , and the angular frequency, ω = ωr + iω i , are both complex. All variables are then normalized using for the length scale the channel half-width H; for the velocity scale the speed of sound c0 ; for the density scale the mean unperturbed density ρ0 ; and for the pressure scale ρ0 c2 0 . This gives:

y = ỹ/ H ; x = x/ H ; α = α H ; ω = ω H/c 0 ; û = û/c 0 ; v = v/c 0 ; p = p/(ρ 0 c2 0 ).
The normalized modal solution is then u (y) = û(y, α, ω)e i(αx-ωt) and the following system is finally obtained for the modal amplitudes:

i(αM 0 -ω)û + v dM 0 dy = -iαp (11) i(αM 0 -ω)v = - dp dy (12)
i(αM 0ω)p + iαû + dv dy = 0 [START_REF] Ingard | Influence of fluid motion past a plane boundary on sound reflection, absorption, and transmission[END_REF] where M 0 (y) = Ũ0 (y)/c 0 is the Mach number. Note that it is possible to combine these equations to obtain the well known Pridmore-Brown equation for calculating compressible modes uppon a sheared mean flow:

d 2 p dy 2 - 2α (αM 0 -ω) dM 0 dy dp dy -α 2 -(αM 0 -ω) 2 p = 0 (14)

Incompressible inviscid model

Let us introduce the incompressible equivalent to the compressible model of the previous section. It is obtained by replacing ∂ ρ ∂ t + ∇•(ρũ) = 0 in Eq. ( 7) by ∇ • (ρũ) = 0. This gives for the modal amplitudes:

i(αM 0 -ω)û + vM 0 = -iαp (15) i(αM 0 -ω)v = - dp dy (16) iαû + dv dy = 0 (17)
The velocity used for the normalization is still the numerical value of the speed of sound c0 of the medium if it were compressible. This is the reason why the normalized velocity still appears as a Mach number, M 0 . At first sight this normalization may appear not very relevant for an incompressible study but this choice will allow a direct comparison with the inviscid compressible case (and of course it is correct in terms of numerical value).

It is also possible to combine these equations to obtain the Rayleigh equation:

d 2 p dy 2 - 2α (αM 0 -ω) dM 0 dy dp dy -α 2 p = 0 (18)
Note that the equation given here involves the pressure variable and is a bit different from the more frequently encountered form involving the streamfunction and a double derivative of the mean velocity with respect to y.

Incompressible viscous model

As found in [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF] and as shown below, it is important to account for the effect of viscosity on the unstable mode. When dealing with acoustic modes, viscous and thermal losses are both equally important. Taking into account losses requires one to use an energy equation. As will be argued below, the unstable mode is hydrodynamic in the sense that it is not much affected by compressibility. This is the reason why an incompressible model will be used for studying the effect of viscosity, and no energy equation is required. As in [START_REF] Reynolds | Stability of turbulent channel flow, with application to Malkus's theory[END_REF][START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3: Theoretical models and comparisons with experiments[END_REF][START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF], the stability analysis concerns some coherent modes superimposed on a turbulence whose effect is to add a kinematic (eddy) viscosity. As seen in Fig. 3, at the high Reynolds number investigated here this eddy viscosity is large compared to the molecular viscosity except in the close vicinity to the wall. Of course, the turbulent eddy viscosity depends on y. Hence, the 2D incompressible viscous linearized Navier-Stokes equations with a variable viscosity are introduced. For the physical variables, these are [START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF]:

∂ ũ ∂ t + Ũ0 ∂ ũ ∂ x + ṽ d Ũ0 dỹ = - 1 ρ0 ∂ p ∂ x + d νT dỹ ∂ ũ ∂ ỹ + d νT dỹ ∂ṽ ∂ x + νT ∂ 2 ũ ∂ x2 + ∂ 2 ũ ∂ ỹ2 (19) ∂ṽ ∂ t + Ũ0 ∂ṽ ∂ x = - 1 ρ0 ∂ p ∂ ỹ + 2 d νT dỹ ∂ṽ ∂ ỹ + νT ∂ 2 ṽ ∂ x2 + ∂ 2 ṽ ∂ ỹ2 (20) ∂ ũ ∂ x + ∂ṽ ∂ ỹ = 0(21)
where as before νT (ỹ) is the total viscosity, including molecular and turbulent contributions. By combining these equations, one could obtain the modified Orr-Sommerfeld equation initially proposed by Reynolds and Hussain [31, Eq 6.1] for calculating the modes of an organized wave in a turbulent medium.

A model is needed for this viscosity, and the expression by Cess, Eq. ( 4), will be used. It is recalled that Ũ0 is obtained by integrating Eq. ( 3), so that the mean velocity and the total viscosity model are consistent, which was not the case in [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF].

As before modal solutions are introduced and the previous equations are normalized. As in the previous section (see remark below Eqs. [START_REF] Rienstra | Spatial instability of boundary layer along impedance wall[END_REF][START_REF] Brambley | Fundamental problems with the model of uniform flow over acoustic linings[END_REF][START_REF] Briggs | Electron-Stream Interaction with Plasmas[END_REF]) the velocity for the normalization is the numerical value of the speed of sound c0 if the medium were compressible. By gathering the powers of α the following set of equations is obtained:

-iωû + vM 0 - 1 Re ν + T d 2 û dy 2 - 1 Re dν + T dy dû dy + α iM 0 û + ip - 1 Re dν + T dy iv + α 2 1 Re ν + T û = 0 (22) -iωv + dp dy - 1 Re ν + T d 2 v dy 2 - 2 Re dν + T dy dv dy + α (iM 0 v) + α 2 1 Re ν + T v = 0(23) dv dy + α (iû) = 0(24)
The Reynolds number is defined by Re = Hc 0 /ν. The total viscosity is normalized in the same way as in section 2.

Boundary conditions

To solve the equations governing the flow, boundary conditions are needed.

The liner at the bottom wall is characterized by its admittance Y (ω) (or its impedance Z(ω)), which provides a relation between the pressure and the wall normal velocity. Using the same normalization as before, as well as the following one:

Y = ρ0 c0 Ỹ
we have at the bottom wall at y = -1:

v(-1, α, ω) = -Y (ω)p(-1, α, ω) ∀α, ω (25)
This relation can be used when the mean flow velocity satisfies M 0 (-1) = 0.

If there is a slip velocity at the wall, that is if M 0 (-1) = 0 (as is the case when the plug flow assumption, M 0 =cst, is used), one needs to account for the vorticity sheet at the wall, and the boundary condition becomes the Myers' boundary condition:

ωv(-1, α, ω) = Y (ω)(αM 0 (-1) -ω)p(-1, α, ω) ∀α, ω (26)
When M 0 (-1) = 0 relations ( 25) and ( 26) are equivalent. Hence, we can use only the latter one. But according to the introduction, when M 0 (-1) = 0, the Myers boundary condition is not "well-posed", and usual criteria for stability are not applicable. We will nevertheless use this condition for purposes of validation.

At the upper rigid wall (at y = 1), the admittance is zero, we have:

v(1, α, ω) = 0 ∀α, ω (27)
In the viscous case, a non-slip velocity is enforced at both walls:

û(±1, α, ω) = 0 ∀α, ω (28)
3.5. Stability problem and solution method

Generalized eigenvalue problem

In the inviscid compressible case, we have to solve the system (11-13) with boundary conditions [START_REF] Howe | The damping of sound by wall turbulent shear layers[END_REF][START_REF] Dokumaci | A note on transmission of sound in a wide pipe with mean flow and viscothermal attenuation[END_REF]. In the viscous incompressible case, we have to solve the system [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF][START_REF] Boyer | Theoretical investigation of hydrodynamic surface mode in a lined duct with sheared flow and comparison with experiments[END_REF][START_REF] Davies | Measurement of plane wave acoustic fields in flow ducts[END_REF] with boundary conditions [START_REF] Howe | The damping of sound by wall turbulent shear layers[END_REF][START_REF] Dokumaci | A note on transmission of sound in a wide pipe with mean flow and viscothermal attenuation[END_REF][START_REF] Allam | Investigation of damping and radiation using full plane wave decomposition in ducts[END_REF]. In either case, the problem is an eigenvalue problem, which can be solved only for some couples (α,ω), that is, when α and ω satisfy the dispersion relation D(α, ω) = 0.

As the problem cannot be solved analytically, a numerical solution method is used. The problem is discretized in the y-direction using a spectral Chebyshev collocation method [START_REF] Trefethen | Spectral methods in Matlab[END_REF][START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF] and N+1 collocation points. The collocation points are of the Gauss-Lobato type [START_REF] Trefethen | Spectral methods in Matlab[END_REF][START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]: y j = cos(jπ/N ) for j = 0...N .

On the grid thus defined, we define the vectors of unknowns:

U = [û 0 û1 ... ûN ] t ; V = [v 0 v1 ... vN ] t ; P = [p 0 p1 ... pN ] t (29)
where [ ] t denotes the transpose. These vectors contain the value of the variables at the collocation points. Hence, for example, û0 is the value of û at collocation point y 0 = 1, and ûN is the value of û at collocation point

y N = -1.
The y-derivative of a vector is obtained by multiplying this vector by a suitable differentiation square matrix D of size N + 1 [START_REF] Trefethen | Spectral methods in Matlab[END_REF][START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]. This procedure yields a generalized eigenvalue problem that is now specified for the different cases under consideration.

In the inviscid compressible case, discretizing the system (11-13), one obtains the following generalized eigenvalue problem:

              -iωI M 0 -iωI D D -iωI      A +α      iM 0 iI iM 0 iI iM 0      B          •      U V P      = 0 (30)
A and B are square matrix of size 3N +3; M 0 is the diagonal matrix (of size N +1) with the values of the mean velocity at the collocation points on the diagonal; M 0 is also a diagonal matrix (of size N +1) containing the values of the derivative of the mean velocity at the collocation points; I is the identity matrix (of size N +1). Parts left blank are filled with zeros. One has still to take into account the boundary conditions [START_REF] Howe | The damping of sound by wall turbulent shear layers[END_REF][START_REF] Dokumaci | A note on transmission of sound in a wide pipe with mean flow and viscothermal attenuation[END_REF]. This is done by replacing the rows N+2 and 2N +2 of the matrices A and B by the discretized form of the boundary conditions.

In the inviscid incompressible case, discretizing the system [START_REF] Rienstra | Spatial instability of boundary layer along impedance wall[END_REF][START_REF] Brambley | Fundamental problems with the model of uniform flow over acoustic linings[END_REF][START_REF] Briggs | Electron-Stream Interaction with Plasmas[END_REF] gives:

               -iωI M 0 -iωI D D      A inc +α      iM 0 iI iM 0 iI      B inc           •      U V P      = 0 (31)
In the viscous incompressible case, the system given by Eqs. [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF][START_REF] Boyer | Theoretical investigation of hydrodynamic surface mode in a lined duct with sheared flow and comparison with experiments[END_REF][START_REF] Davies | Measurement of plane wave acoustic fields in flow ducts[END_REF] becomes quadratic in α. To get a system that depends only on the first power of α, the two following auxiliary variables are introduced [START_REF] Bridges | Differential eigenvalue problems in which the parameter appears nonlinearly[END_REF]: φu = αû and φv = αv

Using the same Chebyshev collocation method as before, the system can then be put in the form of an eigenvalue problem of the type:

(A v + αB v ) [U V P Φ u Φ v ] t = 0 ( 32
)
where Φ u and Φ v are the vectors containing the φu (y j ) and the φv (y j ) at the collocation points. The matrix A v and B v are now square matrices of size 5N +5 and are given by:

A v =            -iωI -1 Re (ND 2 + N D) M 0 -iωI -1 Re (ND 2 + 2N D) D D I I            (33)
B v =            iM 0 -1 Re iN iI 1 Re N iM 0 1 Re N iI -I -I            (34)
The matrix N and N are diagonal matrices of size N +1 that contain respectively on their diagonal the values of the viscosity ν + T (y j ) and of the first derivative of the viscosity dν + T /dy(y j ) at the collocation points. As for the inviscid compressible case, boundary conditions are enforced by replacing the suitable entries of the matrices.

Shooting methods

The Chebyshev solution method is a global method that allows obtaining the full spectrum by solving a matrix system in one go. It is very useful but may sometimes be unpractical when the interest lies on a particular eigenvalue and on the tracking of this eigenvalue during the variation of a given parameter. In the later case it is convenient to obtain α(ω) by a shooting method. The shooting method has been used in the compressible and incompressible inviscid cases only, using the Pridmore-Brown equation, Eq. ( 14), and the Rayleigh equation, Eq. ( 18), together with boundary conditions [START_REF] Howe | The damping of sound by wall turbulent shear layers[END_REF][START_REF] Dokumaci | A note on transmission of sound in a wide pipe with mean flow and viscothermal attenuation[END_REF][START_REF] Allam | Investigation of damping and radiation using full plane wave decomposition in ducts[END_REF] expressed in terms of pressure only. For a fixed value of ω the shooting method determines the branch α(ω) in the vicinity of a guessed value of α (which is provided by the Chebyshev method). This method will complement the Chebyshev method occasionally.

Stability analysis

The generalized eigenvalue problems (30), ( 31) and ( 32) have been derived above. These are of the form:

A ω f = αB ω f (35)
where A ω and B ω are matrices whose dependence on ω has been stressed. Solving these problems using standard libraries (in the present case, the function eig of Matlab is used) provides both the eigenvalues and the eigenfunctions. That is, fixing ω, it is possible to obtain the branches α (n) (ω)

satisfying the dispersion relation. Knowing these, it is possible to apply the classical methods of stability such as the Briggs-Bers criterion to classify the modes as stable, convectively unstable, or absolutely unstable [START_REF] Briggs | Electron-Stream Interaction with Plasmas[END_REF][START_REF] Huerre | Hydrodynamic Instabilities in Open Flows, Hydrodynamics and Nonlinear Instabilities[END_REF][START_REF] Huerre | Absolute and convective instabilities in free shear layers[END_REF].

The reader is referred to these references for a detailed discussion of the method; only the essential material is now outlined. First, let us define in the ω-plane a contour L ω = {ω; ω i = cst = ω L,i }. Such a contour is a line parallel to the real axis. For a well-posed problem and a L ω contour having a sufficiently large and positive value of ω L,i , the branches α (n) (ω) for ω ∈ L ω are well separated: some are lying entirely in the half-plane α i > 0, the others are lying entirely in the half-plane α i < 0. The method consists in tracking the behaviour of the branches α (n) (ω) when decreasing the value of ω L,i from its large initial value towards 0. There are two possibilities:

• (C1) If in the process of decreasing ω L,i , and before ω L,i has reached 0, two opposite spatial branches merge (or pinch) there is an absolute instability consisting of an exponential growth in time in the whole space. By opposite spatial branches, it is meant that one branch is lying entirely in the half-space α i > 0 while the other is lying entirely in the half-space α i < 0 for the initial value ω L,i .

• (C2) If such a pinch does not occur there is no absolute instability and the system is at most convectively unstable. It is valid to perform a spatial analysis whereby ω is real (ω then belongs to the L ω contour having ω L,i =0) and the wavenumber is complex. One still needs to identify if there are convectively unstable modes. Any mode having α i =0 is a candidate and stability is ultimately determined by the direction of propagation of that mode. To obtain this direction, one has to vary the imaginary part of ω from 0 to a large value and track by continuity the corresponding branch α(ω). If, for a very large value of ω i , the mode α(ω) ends up in the half-plane α i >0, this mode is propagating toward x > 0 (it contributes to the impulse response only in the region x > 0). If it ends up in the half-plane α i <0 it is propagating in towards x <0. For example suppose a mode has α i <0 for real ω and ends up with α i > 0 after increasing the imaginary part of ω to a large value:

then this mode would be a convective instability growing towards x >0.

Both criteria (C1) and (C2) will be used below to determine the kind of stability predicted by the different models.

Validation of the numerical solver

Since the focus will be on unstable modes, it is proposed to validate the method against the solution Rienstra called an unstable hydrodynamic mode [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF], which is one of the different surfaces modes for which he obtained an asymptotic analytical expression. In the limit ω → ∞, surface waves are solution to the following equation [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF]:

(1 -M 0 σ) 4 = (1 -M 2 0 ) 3 Y 2 (1 -σ 2 ) ( 36
)
where σ is the reduced wave-number defined by: α = (σ

-M 0 )ω/(1 -M 0 ) 2 .
One has to keep only solutions that verifies (γ) > 0, where γ is the reduced transverse wavenumber defined by

γ = Y (1 -M 0 σ) 2 /(1 -M 2 0 ) 3/2
, and designates the imaginary part. Initially, results were given by Rienstra in cylindrical coordinates (and for the opposite time convention) but it can be checked that they remain valid in 2D cartesian coordinates. Note also that these results are for a plug flow and the Myers' boundary condition. This is not a problem since we are not interested here in the validity of the boundary condition but only in the validation of the numerical solver.

In Fig. 5, the different modes are given for the following conditions: M 0 =0.5, ω=2, Y =2.5-2.5i. N =100 modes are used. Two modes are surface modes as calculated by [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF], the one labelled "HI" represents the hydrodynamic instability (according to Rienstra). The precise value obtained for the corresponding wavenumber is given in table 1: the value obtained by the present Chebyshev integration using N +1=101 collocation points is in perfect agreement with the asymptotic analysis (which appears to be valid for a value of ω as low as 2). The agreement is equally good with the shooting method.

With a more complicated velocity profile such as the Cess profile that will be Method α Asymptotic analysis 7.62064475 -2.04970330i (Eq. ( 36)) Numerical integration 7.62064475 -2.04970330i

Chebyshev using Eq. ( 30), N =100 Numerical integration 7.62064475 -2.04970330i

Shooting method using Eq. ( 14)

Table 1: Wave-number for the HI mode, as calculated using different methods.

used below, the Chebyshev method requires more points, typically N =200. The numerical value for the Hydrodynamic Instability (HI) mode is given in table 1.

Results for the inviscid models

Spatial analysis

Many studies of modes in ducts have involved a spatial stability analysis [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF][START_REF] Brambley | Classification of aeroacoustically relevant surface modes in cylindrical lined ducts[END_REF][START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF] whereby the angular frequency is real while the wavenumber is complex. This approach corresponds to point (C2) in section 3.5.3; it is a practical one in the sense that it corresponds to the physical situation whereby one excites the flow at some real frequency. The spatial analysis provides the spatial modes at this frequency. This approach is valid as long as there is no absolute instability, which we are now going to suppose for a moment.

The spatial spectrum is calculated by solving the compressible inviscid problem, Eq. ( 30). The conditions are identical to the experiment: the flow is approximated by the Cess profile as in Fig. 2, and the angular frequency is the real frequency ω = ω a =0.196 of the incident wave. The wave number spectrum is complex; it is shown in Fig. 6(a). Since a convective instability was observed in the experiments, it is appealing to attempt to identify such an instability by applying the criterion (C2) given in section 3.5.3. The imaginary part of the angular frequency is thus increased from 0 to a very large value. Here ω is varied between ω a and ω a + 0.4i; the later value is not extremely large but is sufficient to observe the main move of the spectrum (a further increase does not change the sign of α i ). The α-spectra thus obtained are shown in Fig. 6 on Fig. 6(b). Hence, if the spatial stability analysis were valid this mode would be a convective instability growing towards x> 0 with a growth rate -α i ∼1.1. This reasoning was done in reference [START_REF] Marx | Comparison of experiments with stability analysis predictions in a lined flow duct[END_REF] assuming that a spatial analysis was valid, that is, assuming that there is no absolute instability. Such a kind of reasoning was also done in [START_REF] Rienstra | A classification of duct modes based on surface waves[END_REF] with the ill-posed Myers boundary condition. Unfortunately, it is shown below that in our specific case there is an absolute instability, meaning that the spatial analysis is indeed not valid.

Spatio-temporal analysis

In the previous section the criterion (C2) was applied directly without checking condition (C1) (see section 3.5.3). A safer way to investigate sta-bility is to start with L ω contours having a high value of ω L,i , as indicated in condition (C1). Then the value of ω L,i is decreased and the branches α are monitored. The result is shown in Fig. 7 where only the α-branch that is suspected to be an instability has been tracked. This figure features a typical branch swap as ω L,i is decreased from 0.11 to 0.10, indicating the presence of a saddle point of the mapping α(ω). Moreover the branches that swap are opposite in the sense given in C1. This reveals the presence of a pinch point for a positive value of ω L,i ∼0.107. Hence, there is indeed absolute instability, and the spatial analysis of the previous section was not valid.

Note that the sole trajectory, when the imaginary part of ω is increased, of ω L,i =0.11; ω L,i =0.12 ; ω L,i =0.18; ω L,i =0.41. The curve is the trajectory of the branch for ω r =ω a and ω i increasing from 0 to positive values (this is the same as in Fig. 6(b)).

the branch starting at α(ω a ) ∼ 0.79 -1.1i when ω = ω a (square symbol), does not allow identifying such a branch swap in Fig. 6(b). This trajectory is again represented in Fig. 7, and combining this with the mappings α(L ω )

shows that this is indeed the case.

It was shown in [START_REF] Rienstra | Mean flow boundary layer effects of hydrodynamic instability of impedance wall[END_REF][START_REF] Rienstra | Boundary-layer thickness effects on the hydrodynamic instability along an impedance wall[END_REF] using a piecewise linear mean flow profile and an inviscid model that absolute instability is unlikely in many practical aeronautical applications. A similar kind of study with a piecewise linear flow and an incompressible inviscid model was performed in a channel for conditions similar to the present study [START_REF] Marx | A piecewise linear mean flow model for studying stability in a lined channel[END_REF]. It concluded that there was an absolute instability of the flow in our configuration. This conclusion remains valid with the present model despite its using a more realistic velocity profile.

Incompressible case

Below, viscous processes are to be taken into account together with an incompressible model. In the present section, the use of an incompressible model is justified by comparing the spectra provided by compressible and incompressible inviscid models. These are respectively given by the solution of Eq. ( 30) and Eq. [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3: Theoretical models and comparisons with experiments[END_REF]. The spectra for both models obtained for ω = ω a are provided in Fig. 8(a). It is seen that the general agreement is acceptable. In particular the branch that is involved in instability, marked by an arrow, is well described by an incompressible model. Figure 8(b) compares the branches implied in the pinching. The compressible and incompressible models compare very well. The incompressible model predicts absolute instability, and the branch obtained for the contour ω L,i =0.12 (that is, before pinching) are almost indiscernible. This is important since this branch (also marked by an arrow) is basically responsible for instability. Hence, as far as the unstable mode is concerned, the incompressible model seems to be a satisfying approximation. ω L,i =0.11 compressible; ω L,i =0.12 compressible; ω L,i =0.10 incompressible; ω L,i =0.11 incompressible; ω L,i =0.12 incompressible.

Results for the viscous incompressible model

In the previous section, the inviscid compressible or incompressible models applied to our experimental configuration predicts absolute instability.

This conclusion does not match the experiments, the later rather displaying spatial instability. In the present section, the effect of viscosity is taken back into account, and this will be seen to change the type of instability, from absolute to convective. The model to be used has been presented in section 3.3.

Spatially unstable mode

The spectrum obtained at ω=ω a by solving Eq. ( 32) is given in Fig. 9(a).

Some differences are to be noted compared to the spectra obtained at the same frequency with the inviscid, incompressible or compressible, spectra displayed in Fig. 8(a). First in the viscous approximation, there is no continuous branch along the α r -axis, since the system is no more singular for α(ω)U (y) = ω. Second there are additional branches due to viscosity located along α r ∼0.8. The evanescent modes along α r =0 are present in all cases. Finally, there is a mode with α num ∼0.94-0.22i indicated by the square symbol, which as before, belongs to an unstable branch. The branch is un- stable because it lies entirely in the upper α-plane for large values of ω i , and moves into the lower α-plane when ω i decreases to 0. This is shown in There is a small discrepency here with the experiment since experimentally one would expect the growth rate to be maximal at ω = ω a . Indeed, in the measurements the acoustic transmission was maximal at this frequency [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF], and this should also correspond to a maximal growth rate of the instability. This discrepency is probably due to the idealized modelling as well as to the sensitivity of the growth rate to both flow velocity and frequency (see next subsection).

Here a turbulence model has been used to account for the turbulent flow.

Calculations done with a uniform molecular viscosity with a value increased by 30 times gave comparable results. Interestingly, Reynolds and Hussain [START_REF] Reynolds | The mechanics of an organized wave in turbulent shear flow. Part 3: Theoretical models and comparisons with experiments[END_REF] also found that their spatial modes calculated either with an eddy viscosity model or with a uniform viscosity multiplied by 30 were very similar.

Hence, for preliminary studies about the effect of viscosity on incompressible modes in a channel, one could possibly use a molecular viscosity multiplied by 30.

In conclusion the instability is a spatial one, growing in the x-direction with a growth rate -α i =0.22. The computed value α num ∼0.94-0.22i has to be compared to the experimental value α exp =1.3-0.3i [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF]. This is less than a 30% error, which is acceptable given the approximations that have been made.

The numerical eigenfunction corresponding to this unstable mode are also provided by the solution of Eq. [START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF]. They are compared to the eigenfunctions measured experimentally [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF] in Fig. 10(a-d). The amplitudes are normalized with respect to the maximum of the transverse velocity and the phase is normalized with respect to the phase of the axial velocity at the channel center.

There is a good agreement between the numerical and experimental eigenfunctions, with an error more pronounced on the axial velocity magnitude.

Effect of mean flow velocity

The only parameter that could be changed during the experiments reported in reference [START_REF] Marx | PIV and LDV evidence of hydrodynamic instability over a liner in a duct with flow[END_REF] was the flow Mach number. While the laser velocity measurements were performed only at M 0,avg =0.27, acoustic transmission were measured over a range of Mach number below this value. Indeed, the instability was first uncovered from high acoustic transmission rates by the liner [6, Fig. 5]. Sound transmission by the liner displays a peak for Mach numbers larger than about 0.2. Above this value, the frequency for maximal transmission increases with Mach number. If we suppose that the sound transmission is proportional to the unstable wave amplification rate, this means that: an unstable wave should exist only above some threshold value of the Mach number; and past this value the frequency for maximal amplification should increase. To test these statements, the effect of velocity is now investigated. Notice that the model is incompressible, but the mean sectionaveraged velocity is normalized by the speed of sound in the compressible fluid and is thus presented as a section-averaged Mach number, M 0,avg . In the mean velocity range investigated no pinching point has been observed and the instability is spatial whenever it exists, which legitimates a spatial A wider picture is obtained by varying the (real) frequency for several values of the Mach number. The growth rate obtained for the mode is shown in Fig. 12. For the lowest Mach number, M 0,avg =0.16, the growth rate is negative for all frequencies, and there is no instability. For all Mach number above this threshold of about M 0,avg ∼0.16, there is spatial instability over some frequency range. The frequency range sustaining instability is wider for larger Mach numbers. One observes also that the maximum of the growth rate is shifting towards higher values of ω r when the Mach number is increased. Hence, the model allows to confirm the experimental observations mentioned above. The conditions taken for the analysis are similar to that of a previsouly reported experiment: same geometry, same impedance law, same real frequency for spatial analyses, same mean flow Mach number. There are necessary simplifications: in particular, the analytical velocity profile corresponds well to the flow in a channel with rigid walls. An expression for velocity profiles (that is, equivalently, an expression for the eddy viscosity) over lined wall would be desirable. Despite these simplifications, the viscous model identifies a spatial instability whose characteristics does match correctly the experimental ones, in terms of wavenumber (wavelength, and growth rate) and of eigenfunction shapes.

Figure 1 :

 1 Figure 1: Studied configuration. The notation˜is for physical quantities .

Figure 2 :

 2 Figure 2: Experimentally measured mean velocity profile; and analytical approximation obtained by integration of Eq. (3) using Cess' expression for the total viscosity, Eq. (4).

Figure 3 :

 3 Figure 3: Total viscosity, given by Eq. (4) for Re τ =2760.

Figure 4 :

 4 Figure 4: Maps of the phase-averaged axial (ũ ) and transverse (ṽ ) velocity over the height of the channel. The liner covers the bottom wall, and the origin x=0 is counted from the position where the liner starts.

Figure 5 :

 5 Figure 5: Spectrum for: M 0 =0.5, ω=2, Y =2.5-2.5i. The plug flow and Myers boundary condition are used together. Symbols: + present numerical method with N =100 Chebyshev modes; o asymptotic approximation for surface waves as given by roots of Eq. (36).

Figure 6 :

 6 Figure 6: (a) Spectrum calculated by solving (30) for ω = ω a =0.196. (b) Trajectories of the modes when the imaginary part of ω is increased from 0 to 0.41i, the real part been kept constant at ω r =ω a . The symbols + indicate the starting points for ω i =0 (the same as in (a)); the symbols o indicates the points for ω i =0.41. The gray dots are for the trajectory in between. The symbol indicates the potentially unstable mode. The line obtained with the shooting method is a further verification of the trajectory of the potentially unstable mode.

Figure 7 :

 7 Figure 7: Trajectories of one selected branch α(ω) for ω belonging to L ω contours with: ω L,i =0.10 ;ω L,i =0.11; ω L,i =0.12 ; ω L,i =0.18; ω L,i =0.41. The curve

Figure 8 :

 8 Figure 8: (a) Spectra obtained at ω = ω a =0.196: + compressible case Eq. (30) ; • incompressible case Eq. (31). (b) Trajectories of the unstable mode along several L ω contours: ω L,i =0.10 compressible;ω L,i =0.11 compressible; ω L,i =0.12 compressible;

Figure 9 :

 9 Figure 9: (a) Spectrum obtained at ω=ω a =0.196 by solving Eq. (32). The symbol indicates the convectively unstable mode. The symbol indicates the experimental unstable mode. (b) Trajectories of the unstable mode along several L ω contours: ω L,i =0 (ω real); ω L,i =0.0009; ω L,i =0.0018. The line represents the trajectory of the unstable mode when the imaginary part of ω is increased from 0 ( symbol) to 0.41i (o symbol), the real part being fixed at ω r = ω a . The symbol • is for the mode having maximal growth rate for ω real.

Fig. 9 (

 9 Fig. 9(b). Contrary to the insviscid case, no pinch has been observed when tracking mappings of L ω contours in the α-plane. Hence, taking into account viscosity does indeed change the stability from absolute instability to convec-

Figure 10 :

 10 Figure 10: Comparison of numerical ( ) and experimental ( ) eigenfunctions for the unstable mode at ω r =ω a . (a) Amplitude of the axial component; (b) Phase of the axial component; (c) Amplitude of the transverse component; (d) Phase of the transverse component.

Figure 11 :

 11 Figure 11: (a) Trajectory of the unstable mode and (b) Growth rate as a function of the section-averaged Mach number, obtained at ω=ω a =0.196 when varying the sectionaveraged Mach number from M 0,avg ∼0.09 (+ symbol) to ∼0.45 (o symbol). The symbol indicates the convectively unstable mode at the baseline value M 0,avg =0.27.

Figure 12 :

 12 Figure 12: Spatial growth rate for the spatially unstable branch as a function of ω r for several values of the section-averaged Mach number: M 0,avg =0.16 ; M 0,avg = 0.165 ; M 0,avg =0.18 ; M 0,avg =0.21; M 0,avg =0.27; M 0,avg =0.36 ; M 0,avg = 0.45. The symbol indicates the convectively unstable mode at the baseline value M 0,avg =0.27 and ω=ω a =0.196.


C De Souza

S Tarbouriech
email: sophie.tarbouriech@laas.fr

E B Castelan
email: eugenio.castelan@ufsc.br

V J S Leite

Event-triggered Dynamic Output Feedback Controller for Discrete-time LPV Systems with Constraints ⋆

Keywords: Linear parameter-varying systems, Event-triggered control, Saturation, Dynamic controller, Discrete-time systems

This paper examines an event-triggered control design approach for discrete-time linear parameter-varying (LPV) systems under control constraints. A parameter-dependent dynamic output-feedback controller with an event-triggering condition is designed to ensure the regional asymptotic stability of the closed-loop system while minimizing a quadratic cost function. Sufficient conditions are derived in terms of linear matrix inequalities (LMIs) thanks to the use of a parameter-dependent quadratic Lyapunov function. Also, convex optimization schemes are proposed using these conditions to minimize an upper bound of a given cost function or to maximize the size of the region of closed-loop stability. Examples illustrate the proposal.

INTRODUCTION

Event-triggered control (ETC) has attracted a lot of attention in recent years due to its significant benefits in terms of reduced usage of the communication and computational resources compared to traditional control systems (Sandee et al., 2005;[START_REF] Aström | Event based control[END_REF][START_REF] Heemels | An introduction to event-triggered and self-triggered control[END_REF]Postoyan et al., 2015). The idea behind the event-triggered control consists of performing control tasks after the occurrence of an event, generated by some well-designed event-triggering mechanism, rather than the elapse of a certain fixed period of time, as in conventional periodic sampled-data control. ETC is then capable of reducing effectively the control task executions, while ensuring a satisfactory closed-loop performance. Various ETC strategies can been found in the literature: see, for example, (Wu et al., 2014a;Zhang and Han, 2014;[START_REF] Khashooei | Output-based event-triggered control with performance guarantees[END_REF] for continuous-time framework and [START_REF] Eqtami | Event-triggered control for discrete-time systems[END_REF]Tallapragada and Chopra, 2012;Wu et al., 2016) for discrete-time counterpart. However, most of the prior results have been obtained for linear and nonlinear systems, without taking into account the presence of measurable parameters or control constraints.

Besides, Linear Parameter-Varying (LPV) systems represent a very important class of dynamical systems whose dynamics depends on a priory unknown, but on-line measurable time-varying parameters. Due to their effectiveness in modeling and control of nonlinear systems, LPV systems have been extensively studied in the literature (Mohammadpour and Scherer, 2012). However, a few studies have been done on event-triggered control of LPV systems, in particular in the discrete-time context. In [START_REF] Hu | Eventtriggered guaranteed cost control for uncertain discretetime networked control systems with time-varying transmission delays[END_REF] an event-triggered guaranteed cost control is proposed for systems with norm-bounded uncertainties and timevarying transmission delays. [START_REF] Li | Co-design of eventtriggered H ∞ control for discrete-time linear parametervarying systems with network-induced delays[END_REF] propose an event-triggered H ∞ control by jointly designing an eventtriggering mechanism with two given threshold (proportional, additional). In [START_REF] Li | Event-triggered control for discrete-time uncertain linear parameter-varying systems[END_REF][START_REF] Golabi | Event-triggered control for discretetime linear parameter-varying systems[END_REF] the problem of the co-design of an event-triggering condition and a state-feedback controller is addressed.

Additionally, saturation is a common phenomenon in practical control systems. However, its presence in a control loop can cause performance degradation or even instability. Recently, the ETC problem under actuator saturation was addressed, for instance, by [START_REF] Kiener | Actuator saturation and anti-windup compensation in event-triggered control[END_REF]Seuret et al. (2016); [START_REF] Liu | Event-triggered control for linear systems with actuator saturation and dis[END_REF]; [START_REF] Li | Dynamic output feedback control for systems subject to actuator saturation via event-triggered scheme[END_REF] for continuous-time framework and by Wu et al. (2014b);Zuo et al. (2016); [START_REF] Groff | Event-triggered control co-design for discretetime systems subject to actuator saturation[END_REF]Ma et al. (2019) for discrete-time counterpart. More precisely, a procedure to design a state-feedback controller under a given eventtriggering condition that maximizes the region of attraction of linear systems is proposed in (Wu et al., 2014b) and of linear piece-wise affine systems in (Ma et al., 2019). A cone complementary linearization algorithm is used in (Zuo et al., 2016) for solving the non-convex optimization procedure to treat the co-design of an event-triggering strategy and a controller under saturation. In [START_REF] Groff | Event-triggered control co-design for discretetime systems subject to actuator saturation[END_REF], the co-design of a static state-feedback gain and an event-triggering function to stabilize discrete-time LTI systems subject to actuator saturation is proposed aiming at reducing the sampling activity.

Based on the above discussion, it is clear that the eventtriggered control design remains an open issue for discretetime LPV systems under saturating actuators. Thus, the main contributions of this paper aims at filling such gap providing the following features: i) A convex procedure to design a parameter-dependent dynamic output-feedback controller with anti-windup action subject to a given eventtriggering condition; ii) The proposed methodology includes the optimization of the control performance, while transmission reductions can be configured by adapting a given parameter in the event-generator. The convex conditions are formulated in terms of LMIs applying the Lyapunov theory together with the S-procedure and the generalized sector condition, and guarantee the regional asymptotic stability of the closed-loop system. The usefulness of our design approach is illustrated by a numerical example and comparison with similar approach in the literature.

Notation: The set of real numbers is denoted by R. The set of integer numbers belonging to the interval from a ∈ Z up to b ∈ Z, b ≥ a, is denoted by I[a, b]. The set of matrices with real entries and dimensions m×n is noted by R m×n . A block-diagonal matrix A with blocks A 1 and A 2 is denoted as A = diag{A 1 , A 2 }. The ℓ th line of a vector or matrix A is indicated by A (ℓ) . The matrix 0 stands for the null matrix of appropriate dimensions and I n corresponds to the identity matrix with dimensions n × n. The symbol ⋆ stands for symmetric blocks in the matrices. • represents an element that has no influence on development.

PROBLEM FORMULATION

Consider the discrete-time system represented by:

x

(k + 1) = A(θ k )x(k) + B(θ k )sat(u(k)), y(k) = Cx(k), (1)
where

x(k) ∈ R n is the state vector, y(k) ∈ R p is the measurable output, u(k) ∈ R m is the control signal and sat(u(k)) is a symmetric decentralized saturation function given by sat(u (ℓ) (k)) = sign(u (ℓ) (k)) min(|u (ℓ) (k)|, ū(ℓ) ) (2) with ū(ℓ) > 0, ℓ ∈ I[1, m],
being the ℓ th component of the symmetric saturation level ū.

The time-varying matrices A(θ k ) ∈ R n×n and B(θ k ) ∈ R n×m belong to a polytopic set given by the convex combination of known N vertices, i.e.

[A(θ k ) B(θ k )] = N i=1 θ k(i) [A i B i ] , (3)
where θ k ∈ R n is a vector of time-varying parameters that belongs to a unitary simplex

Θ N i=1 θ k(i) = 1, θ k(i) ≥ 0, i ∈ I[1, N ] . (4)
To stabilize system (1), we adopt the following parameterdependent dynamic output-feedback controller:

x c (k + 1) = A c (θ k )x c (k) + B c (θ k )ŷ(k) -E c (θ k )Ψ(u(k)), u(k) = C c (θ k )x c (k) + D c (θ k )ŷ(k), (5)
where

x c (k) ∈ R n is the state of the controller, Ψ(u(k)) : R m → R m is the dead-zone non-linearity defined by Ψ(u(k)) = u(k) -sat(u(k))
, and ŷ(k) is a signal defined as

ŷ(k) = y(k), if y(k) is updated, ŷ(k -1), if y(k) is not updated. ( 6
)
Note that we are employing an event-triggering mechanism to determine whether or not the current measurement y(k) should be transmitted to the controller. The idea is to reduce the number of packets exchanged between these nodes, so that we can save energy and communication bandwidth while ensuring the stability of the closed-loop system with a certain performance index. The updating of the output signal is performed ever that the following event-triggering condition is verified ŷ

(k -1) -ŷ(k) 2 > σ 2 (µ y(k) 2 + (1 -µ) u(k) 2 ) (7) for a given σ > 0 and µ ∈ I[0, 1].
Additionally, we assume the following linear structures for the matrices of the controller (5):

[A c (θ k ) B c (θ k )] = N i=1 N j=i (1 + ς ij )θ k(i) θ k(j) A cij 2 B cij 2 , [C c (θ k ) D c (θ k )] = N i=1 θ k(i) [C ci D ci ] , E c (θ k ) = N i=1 θ k(i) E ci ,
with θ k ∈ Θ and ς ij = 1 if i = j and ς ij = 0 otherwise.

In this technical note, we consider not only the stability of the discrete-time saturated LPV system by using a dynamical event-triggered control, but also the satisfaction of a certain level of performance measured by a quadratic cost function. The considered cost function is defined as follows

J ∞ = ∞ k=0 J(k) = ∞ k=0 x(k) ′ Qx(k) + u(k) ′ Ru(k), (8)
where Q ∈ R n×n and R ∈ R m×m are symmetric and positive definite matrices.

Since there is a non-linearity in the loop, the global stability of the origin may not be guaranteed. In this case, the region of attraction R A , defined in terms of the augmented state vector

ξ(k) = [x(k) ′ x c (k) ′ ] ′ ∈ R 2n
, must be considered. As the exact characterization of the R A is, generally, a hard task, it is important to determine subsets with a well-fitted analytical representation (see, for example, in (Tarbouriech et al., 2011)). By denoting R E the estimated attraction region, then we are interested in computing R E ⊆ R A and in optimizing its size with respect to some criteria (Tarbouriech et al., 2011).

Hence, the problem we intend to solve can be summarized as follows. Problem 1. Given σ and µ, design a parameter-dependent dynamic output-feedback controller (5) under the eventtriggering condition ( 6)-( 7), that ensures the regional asymptotic stability of the closed-loop system, while minimizing the quadratic cost function (8).

MAIN RESULTS

Preliminary Results

Define the error variable as

e(k) = ŷ(k) -y(k). (9)
Based on the definition (7), the following inequality 10) is always satisfied. Thus, with the control law (5) and ŷ(k) = e(k) + y(k), the closed-loop system takes the form

e(k) 2 ≤ σ 2 (µ y(k) 2 + (1 -µ) u(k) 2 ) (
ξ(k + 1) = A(θ k )ξ(k) -B(θ k )Ψ(u(k)) + E(θ k )e(k) u(k) = K(θ k )ξ(k) + D c (θ k )e(k) (11) y(k) = Cξ(k) where ξ(k) = [x(k) ′ x c (k) ′ ] ′ ∈ R 2n and A(θ k ) = A(θ k ) + B(θ k )D c (θ k )C B(θ k )C c (θ k ) B c (θ k )C A c (θ k ) , B(θ k ) = B(θ k ) E c (θ k ) , E(θ k ) = B(θ k )D c (θ k ) B c (θ k ) , K(θ k ) = [D c (θ k )C C c (θ k )] and C = [C 0] .
The time-varying matrices of the closed-loop system (11) also verify

[A(θ k ) E(θ k )] = N i=1 N j=i (1 + ς ij )θ k(i) θ k(j) A ij 2 E ij 2 , [B(θ k ) K(θ k ) ′ ] = N i=1 θ k(i) [B i K ′ i ] ,
with θ k ∈ Θ and ς ij = 1 if i = j and ς ij = 0 otherwise.

To investigate the regional asymptotic stability of the closed-loop system (11), we propose the following candidate Lyapunov function

V (k) = ξ(k) ′ P -1 (θ k )ξ(k), (12)
where P (θ k ) = N i=1 θ k(i) P i , with 0 < P i = P ′ i ∈ R 2n×2n and θ k ∈ Θ. Thus, we are searching for a dynamic controller such that the difference ∆V (k) = V (k+1)-V (k) along the trajectories of the closed-loop system satisfies

∆V (k) = ξ(k + 1) ′ P -1 (θ k+1 )ξ(k + 1) -ξ(k) ′ P -1 (θ k )ξ(k) < -J(k). ( 13
)
By summing (13) up from k = 0 and k = ∞, it yields to J ∞ < ξ(0) ′ P -1 (θ 0 )ξ(0) < tr(P -1 (θ 0 )) ξ(0) 2 , (14) with J ∞ defined in (8). Thus, we can conclude that the upper bounded of the cost function J ∞ is related to the Lyapunov matrix P -1 (θ 0 ) and the initial state ξ(0). The estimation of the region of attraction of the origin for the closed-loop system is computed, in this case, by

R E = θ k ∈Θ E(P (θ k ) -1 , 1) = i∈I[1,N ] E(P -1 i , 1) (15)
with

E(P -1 i , 1) = ξ(k) ∈ R 2n : ξ(k) ⊤ P -1 i ξ(k) ≤ 1 for all i ∈ I[1, N ].
To deal with the saturating actuator, we use the following property directly derived from Tarbouriech et al. (2011).

Lemma 1. Consider a matrix G(θ k ) = N i=1 θ k(i) G i with G i ∈ R m×2n , i ∈ I[1, N ] and θ k ∈ Θ. If ξ(k) belongs to S(ū) given by S(ū) {ξ(k) ∈ R 2n : |G(θ k )ξ(k)| ≤ ū},
(16) then, for any diagonal positive definite matrix T ∈ R m×m , the non-linearity Ψ(u(k)) satisfies the following inequality.

Ψ(u(k)) ′ T (Ψ(u(k))-(K(θ k )-G(θ k ))ξ(k)-D c (θ k )e(k)) ≤ 0. ( 17
)

Theoretical conditions

Inspired by Scherer et al. (1997), we define the real matrices X, Y, W, and Z ∈ R n×n such that

U = X • Z • , U -1 = Y • W • , and Ω = Y I n W 0 .
Therefore, we have

U Ω = I n X 0 Z and Û = Ω ′ U Ω = Y ′ F ′ I n X , (18)
where, by construction

F ′ = Y ′ X + W ′ Z. ( 19
)
Furthermore, by using the partitioning P = P 11 P 12 ⋆ P 22 , we

obtain P = Ω ′ P Ω = P11 P12 ⋆ P22 (20) with P11 = Y ′ P 11 Y + W ′ P ′ 12 Y + Y ′ P 12 W + W ′ P 22 W , P12 = Y ′ P 11 + W ′ P ′
12 and P22 = P 11 . From that, the following result can be stated to design the dynamic controller (5) when the parameters of the eventgenerator (7) µ and σ are assumed to be given. Theorem 1. (Control Synthesis). Given σ > 0, µ > 0 and κ > 0, consider that there exist symmetric positive definite matrices Pi , positive definite diagonal matrix S and matrices X, Y , T , Âcij , Bcij , Ĉci , Dci and Êci of appropriate dimensions such that (21) (provided at the top of the next page) and the following LMI are feasible.

Û + Û ′ -Pi ⋆ H i(ℓ) ū2 (ℓ) > 0, i ∈ I[1, N ] ℓ ∈ I[1, m]; (22)
with

Π 1ij = ( Dci + Dcj )C Ĉci + Ĉcj , Π 2ij = Y ′ (A i + A j ) + Bcij C A i + A j + (B i Dcj + B j Dcj )C Âcij (A i + A j )X + (B i Ĉcj + B j Ĉci ) , Π 3ij = Bcij B i Dcj + B j Dci , Π 4ij = Êci + Êcj (B j + B i )S , Π 5ij = Q1 /2 2I n 2X ( Dci + Dcj )C Ĉci + Ĉcj , Q = Q 0 0 R + β 2 I m Π 6ij = Q1 /2 0 Dci + Dcj , β 1 = κσ 2 µ, and
β 2 = κσ 2 μ.
Therefore, by choosing full rank matrices W and Z such that (19) holds, we have that the LPV system (1) under the dynamic output-feedback controller (5) with matrices

                  Û + Û ′ -1 ⋆ 1 2 Π 2ij 1 2 Π 3ij -1 2 Π 4ij Pr ⋆ ⋆ C CX 0 0 0 1 β1 I p ⋆ 1 2 Π 5ij 1 2 Π 6ij 0 0 0 I n+m                   > 0, r, i ∈ I[1, N ] j ∈ I[i, N ]; (21)
D ci = Dci , C ci = ( Ĉci -D ci CX)Z -1 , B cij = (W -1 ) ′ ( Bcij -Y ′ (B i D cj + B j D ci ), A cij = (W -1 ) ′ ( Âcij -Y ′ (A i + A j + (B i D cj + B j D ci )C)X -W ′ B cij CX -Y ′ (B i C cj + B j C ci )Z)Z -1 , E ci = (W -1 ) ′ ( Êci S -1 -Y ′ B i ), (23)
is regionally asymptotically stable, and the region R E , defined in (15), is an estimation of the region of attraction of the origin for the closed-loop system.

The proof of Theorem 1 is presented in the Appendix A.

Optimization procedures

We are more particularly interested in the minimization of the upper bound of the cost function J ∞ . In this case, the objective is to design the dynamic controller (5) that minimizes the cost function J ∞ . According (14), J ∞ is upper-bounded by ξ(0) ′ P -1 (θ 0 )ξ(0). Then, considering (20), we have

J ∞ < ξ(0) ′ Ω P -1 (θ 0 )Ω ′ ξ(0) (24)
i.e.

J ∞ < x(0) x c (0) ′ Y I W 0 P -1 (θ 0 ) Y ′ W ′ I 0 x(0) x c (0) (25)
By setting x c (0) equal to zero, the above inequality results in (

J ∞ < x(0) ′ [Y I] P -1 (θ 0 ) [Y I] ′ x(0). So,
) 27
Note that if θ 0 is not known a priori, we have to check the condition (26) for all Pi with i ∈ I[1, N ]. In this case, the same kind of optimization procedure can be used to maximize the region of stability R E when no cost function optimization is taken into account.

SIMULATION RESULTS

Consider the saturated LPV system (1) with the following data:

A(ρ) = (1 + ρ) 1.0018 0.01 0.36 1.0018 , B(ρ) = (1 + ρ) -0.001 -0.184 , C = [0.5 1] ,
where the measured parameter is ρ ≤ 0.05 and symmetric saturation limit is ū = 1. This system can be cast in the polytopic form ( 3)-( 4) with A 1 = A(-0.05), B 1 = B(-0.05), A 2 = A(0.05) and B 2 = B(0.05).

We have run the optimization procedure ( 27) with κ = 15, Q = 0.2I, R = 1, and different values of σ and µ such that 0 < σ < 0.2 and 0 < µ < 1. Then, for each pair (σ, µ), we have simulated the respective closed-loop system with the initial condition ξ(0) = [0.02085 0.1203 0 0] ′ , and an θ k sequence that leads the open-loop system to have unstable modes; thus requiring a greater control effort. In each simulation we computed the cost function J ∞ and percentage of the update rate. The J ∞ varies between 0.657 and 0.723 with no evident pattern. On the other hand, the percentage of update rate, as shown in the lefthand side of Figure 1, decreases with bigger values of σ and µ. In this case, we can conclude that as µ goes to one, it relieves the weight of the control signal, reducing the number of updates. For instance, for σ = 0.18, there was a decrease of 74.51% in the update rate when we increased µ from 0.1 to 0.9.

Another similar set of tests was performed with µ = 0.7, 12.5 ≤ κ ≤ 27.5, and 0 < σ < 0.15. For each achieved controller design and the same previous initial condition, we have simulated the closed-loop and computed both the cost J ∞ and the percentage of rate update. The J ∞ is shown in the right hand-side of Figure 1, where we note the cost increasing with κ. We have observed no pattern in rate update in this case.

We illustrate the closed-loop behavior under the eventtriggering condition with a design dismissing the optimal control law, i.e., with Q = 0, R = 0. We chose µ = σ = 0.1, κ = 1 and simulated the closed-loop response with ξ(0) = [0.1508 0.5741 0 0] ′ . The states, the control signal, the events and the used parameter-varying θ k are showed in Figure 2, where we can see the closed-loop asymptotic stability despite the saturation of the control signal for 0 ≤ k ≤ 13. This test leads to an mean update rate of 65.55%. With a constant ρ = 0, replacing the output matrix C by the identity, Q = 0, R = 0, and µ = 1 we compare our approach with the one in (Wu et al., 2014b). In this case, we set κ with the same value found in (Wu et al., 2014b). So, we show in Figure 3 two projections of R E on the x k plan, one for σ = 0.1 (line) and the other for σ = 0.02 (dashed-line), illustrating how such a region can be sensitive to σ value. However, our approach (blue), with a dynamic output feedback controller, yields to almost the same R E as with the method in (Wu et al., 2014b) (green) which is based on state-feedback. Moreover, if we consider µ = 1, which is not possible in (Wu et al., 2014b), our projections become even larger (magenta).

CONCLUSION

In this paper, we investigated an event-triggered control design approach. Our proposal allows synthesizing dynamic output-feedback controllers under event-triggering conditions for discrete-time LPV systems subject to actuator saturation. The convex conditions, in the form of LMIs, ensure the regional asymptotic stability of the closed-loop system while minimizing the upper bound of a quadratic cost function. The usefulness of the proposed methods was showed through numerical tests.

  we are interested in minimizing tr [Y I] P -1 (θ 0 ) [Y I] ′ , or, similarly, in minimizing the scalar η such that

Fig. 1 .

 1 Fig. 1. Update rate (%) with κ = 15 (left-side) and cost with µ = 0.7 (right).

Fig. 2 .

 2 Fig. 2. Time evolution of event-triggered control system.

( Pi + Pj ) ⋆ ⋆ ⋆ ⋆ ⋆ 0 κI p ⋆ ⋆ ⋆ ⋆ 1 2 (H i + H j -Π 1ij ) -1 2 ( Dci + Dcj ) 2S ⋆ ⋆

⋆ This work has been supported in part by the Brazilian Agencies CAPES under the project Print CAPES-UFSC "Automation 4.0" and CNPq; and by ANR project HANDY 18-CE40-0010.

turbances. IET Control Theory & Applications, 11(9), 1351-1359. Ma, Y., Wu, W., Görges, D., and Cui, B. (2019 By supposing the feasibility of ( 21), it follows from block (1,1) that Û + Û ′ > 0, and consequently, Û has full rank. Thus, in view of ( 19), the matrix (T ′ -Y ′ X) has full rank. Then, it is always possible to choose full rank matrices W and Z such that ( 19) is verified. This shows that the gains ( 23) are well-defined.

Moreover, by considering ( 18), ( 20), and the change of variables Âc , Bc , Ĉc , Dc , Êc , Q, β 1 and β 2 in ( 23), pre-and post-multiply ( 21) by diag{(Ω ′ ) -1 , I p , (Ω ′ ) -1 , I p , I n+m } and its transpose, respectively. Next, multiply the lefthand side of the resulting inequality by θ k(r+1) , θ k(i) and θ k(j) , sum it up for r, i ∈ I[1, N ] and j ∈ I[i, N ], and replace H(θ k ) by G(θ k )U . Then, use the fact that [P (θ k )-

, pre-and post-multiply the resulting inequality by diag{(U ′ ) -1 , I p , (S ′ ) -1 , I 2n , I p , I n+m } and its transpose, respectively, and replace S -1 = T . After, apply Schur complement and rearrange the terms to obtain (28) (reported at the top of the this page).

Finally, pre-and post-multiply both sides of (28) by the augmented vector

′ and X(k) ′ , respectively, and replace A(θ k )ξ(k)-B(θ k )Ψ(u(k))+ E(θ k )e(k) by ξ(k+1), see (11), and ξ(k+1

Therefore, we conclude that the feasibility of ( 21) ensures the positivity of ( 12) and the negativity of ∆V (k). Also, by Lemma 1, we have that ξ(k) ∈ S(ū), and by inequality (10), we have that the event-triggering condition ( 7) is always satisfied. Now, by supposing the feasibility ( 22), pre-and post multiply it by diag{(Ω ′ ) -1 , I} and its transpose, respectively. Next, multiply the left-hand side of the resulting inequality by θ k(i) , sum it up for i ∈ I[1, N ] and replace H(θ k ) (ℓ) by G(θ k ) (ℓ) U . Then, use again the fact P (θ k ) -U -U ′ ≥ U ′ P -1 (θ k )U , and pre-and post -multiply the resulting inequality by diag{(U ′ ) -1 , I} and its transpose, respectively. Finally, apply Schur complement and preand post multiply the resulting inequality by ξ(k) ′ and ξ(k), to get

ū-2 (ℓ) G(θ k ) (ℓ) )ξ(k) ≤ 0, (A.2) which ensures that |G(θ k )ξ(k)| ≤ ū2 and ξ(k) ′ P -1 (θ k )ξ(k) ≤ 1, and consequently, R E ⊆ S(ū).