J Zhang

A Hector

M Rabaud

F Moisy

Wind-wave growth over a viscous liquid

Experimental and theoretical studies on wind-wave generation have focused primarily on the airwater interface, where viscous effects are small. Here we characterize the influence of the liquid viscosity on the growth of mechanically generated waves. In our experiment, wind is blowing over a layer of silicon oil, of viscosity 20 and 50 times that of water, and waves of small amplitude are excited by an immersed wave-maker. We measure the spatial evolution of the wave slope envelope using Free-Surface Synthetic Schlieren, a refraction-based optical method. Through spatiotemporal band-pass filtering of the surface slope, we selectively determine the spatial growth rate for each forcing frequency, even when the forced wave is damped and coexists with naturally amplified waves at other frequencies. Systematic measurements of the growth rate for various wind velocities and wave frequencies are obtained, enabling precise determination of the marginal stability curve and the onset of wave growth. We show that Miles' model, which is commonly applied to water waves, offers a reasonable description of the growth rate for more viscous liquids. We finally discuss the scaling of the growth rate of the most amplified wave and the critical friction velocity with the liquid viscosity.

I. INTRODUCTION

The generation of waves by the wind, whether on a puddle, a pond, or at sea, remains a fascinating phenomenon that has challenged researchers for decades [1][START_REF] Sullivan | Dynamics of winds and currents coupled to surface waves[END_REF][START_REF] Rabaud | The Kelvin-Helmholtz instability, a useful model for wind-wave generation?[END_REF][START_REF] Pizzo | How does the wind generate waves?[END_REF][START_REF] Ayet | The dynamical coupling of wind-waves and atmospheric turbulence: A review of theoretical and phenomenological models[END_REF]. Gaining a comprehensive understanding of this fundamental problem is of first importance for wave forecasting, assessing heat and mass transfers in the ocean, and in engineering applications involving liquid and gas transport in pipes or heat exchangers [START_REF] Boomkamp | Classification of instabilities in parallel two-phase flow[END_REF].

Since the pioneering works of Phillips [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF] and Miles [START_REF] Miles | On the generation of surface waves by shear[END_REF], several decades of research have led to the following two-stage scenario. In the first stage, incoherent wrinkles of very low amplitude are excited by the turbulent pressure fluctuations of the wind (Phillips mechanism), resulting in a linear increase in wave energy over time [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF][START_REF] Zavadsky | Water waves excited by near-impulsive wind forcing[END_REF][START_REF] Li | The principal stage in wind-wave generation[END_REF], with a saturation to a finite amplitude governed by the liquid viscosity for moderate wind [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF][START_REF] Nové-Josserand | Effect of a weak current on wind-generated waves in the wrinkle regime[END_REF]. In the second stage, for sufficient wind, the airwater shear flow becomes unstable, leading to a subsequent exponential increase of the wave energy (Miles mechanism) [START_REF] Miles | On the generation of surface waves by shear[END_REF][START_REF] Plant | Growth and equilibrium of short gravity waves in a wind-wave tank[END_REF][START_REF] Shemer | Spatially evolving regular water wave under the action of steady wind forcing[END_REF]. However, this highly idealized scenario still requires experimental and numerical confirmations. In particular, there is no consensus on the nature of the transition between the Phillips' and Miles' regimes and its dependence on the various physical parameters of the problem [START_REF] Li | The principal stage in wind-wave generation[END_REF][START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF]. In this context, accurate laboratory measurements of growth rates under controlled conditions, for physical parameters distinct from the usual air-water configuration, are highly valuable.

The difficulty in modeling the wind-wave generation problem has several origins. First, the flow on the air side is turbulent, and both the characteristics of the mean velocity profile and the statistical properties of the turbulent stress fluctuations play significant roles. Second, because of the transport of wave energy, the temporal growth of waves translates into spatial growth [START_REF] Plant | Growth and equilibrium of short gravity waves in a wind-wave tank[END_REF], but a quantitative relationship between the two requires steady forcing conditions over long distances and durations, which are rarely satisfied in open sea conditions. Third, the wind forces not only waves but also currents, which may reach saturation on a timescale distinct from the growth time of the waves [START_REF] Zavadsky | Water waves excited by near-impulsive wind forcing[END_REF]. Fourth, surface contamination, which is inevitable in open sea environments and large water tanks, can significantly impact the onset of wave generation [START_REF] Ryan | Recurrent frontal slicks of a coastal ocean upwelling shadow[END_REF]. Lastly, even in controlled laboratory experiments, the presence of walls introduces additional damping effects and undesirable wave reflections [START_REF] Plant | Growth and equilibrium of short gravity waves in a wind-wave tank[END_REF][START_REF] Hidy | Wind action on water standing in a laboratory channel[END_REF][START_REF] Caulliez | The generation of the first visible wind waves[END_REF].

In this paper, we consider the problem of wind-wave generation on viscous liquids (silicon oils of viscosity 20 and 50 times greater than that of water), with the aim of circumventing some of the previous difficulties: the flow in the liquid is purely laminar, enabling accurate analytical modeling; the surface drift remains moderate, simplifying the air-liquid energy transfer problem; the impact of surface contamination is reduced by the use of a liquid of low surface tension; and interference with reflected waves is reduced by the strong viscous damping.

As long as the liquid viscosity is not too high (typically less than 100 mm 2 s -1 ), the two-stage scenario of the air-water configuration, with linearly increasing wrinkles at small wind (Phillips mechanism) and exponentially growing waves at large wind (Miles mechanism), is still relevant; larger liquid viscosities are not considered here, as they exhibit a very distinct behavior governed by the Kelvin-Helmholtz instability [START_REF] Rabaud | The Kelvin-Helmholtz instability, a useful model for wind-wave generation?[END_REF][START_REF] Francis | Wave motions and the aerodynamic drag on a free oil surface[END_REF][START_REF] Miles | On the generation of surface waves by shear flows. Part 3. Kelvin-Helmholtz instability[END_REF][START_REF] Aulnette | Wind-sustained viscous solitons[END_REF][START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF]. While previous studies using viscous liquids focused on waves naturally generated by the wind [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF], here we are interested in waves triggered by a small mechanical disturbance. Provided that their initial amplitude is sufficient, this forcing bypasses the first stage of Phillips [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], and the system directly enters the second stage of Miles [START_REF] Miles | On the generation of surface waves by shear[END_REF]. This shortcut in the two-stage wave generation problem has the advantage of directly assessing the growth rate of each wave number, at the cost of disregarding the realistic conditions required for the "true" wave generation process starting from an initially undisturbed surface. As a result, different critical wind velocities may be expected for natural and mechanically generated waves.

In the exponential growth regime, Miles demonstrated that the normalized temporal growth rate of a given wave number k is proportional to the normalized wind shear stress (u * /c) 2 , where u * is the friction velocity and c(k) the phase velocity, with a coefficient β that depends on the curvature of the mean velocity profile at the (k-dependent) critical height where the mean flow velocity matches the phase velocity. This growth is mitigated by the dissipation in the liquid, either laminar or turbulent, which defines in principle a critical friction velocity for the growth of this wave number k. However, unlike standard hydrodynamic instabilities arising from the infinitesimal disturbance of a base state, the two-stage transition here also depends on the incoherent base state from which the waves grow, which itself depends on the liquid viscosity and the wind forcing conditions (intensity, duration) [START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF]. This complexity likely contributes to the persistent challenge in defining a clear critical wind velocity for the onset of wave generation under realistic conditions.

The growth rate of wind-generated waves has been the subject of a large number of experimental works, either from their temporal growth at a fixed point [START_REF] Larson | Wind-generated gravity-capillary waves: Laboratory measurements of temporal growth rates using microwave backscatter[END_REF][START_REF] Kawai | Generation of initial wavelets by instability of a coupled shear flow and their evolution to wind waves[END_REF][START_REF] Geva | Excitation of initial waves by wind: A theoretical model and its experimental verification[END_REF], or from their spatial growth in steady conditions [START_REF] Hidy | Wind action on water standing in a laboratory channel[END_REF][START_REF] Wilson | Wind-induced growth of mechanically generated water waves[END_REF][START_REF] Mitsuyasu | Wind-induced growth of water waves[END_REF][START_REF] Grass | Measurement of the initiation and growth of surface water waves under the action of a laminar air flow[END_REF][START_REF] Tsai | On the spatial linear growth of gravity-capillary water waves sheared by a laminar air flow[END_REF][START_REF] Liberzon | Experimental study of the initial stages of wind waves' spatial evolution[END_REF] (the equivalence between these two approaches is discussed in Ref. [START_REF] Plant | Growth and equilibrium of short gravity waves in a wind-wave tank[END_REF]), but only in the air-water configuration. These measurements, initially limited to one-point probes, have made significant progress with the advent of optical methods resolved in space and time [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF][START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF][START_REF] Shemer | On evolution of young windwaves in time and space[END_REF][START_REF] Yousefi | Momentum flux measurements in the airflow over wind-generated surface waves[END_REF][START_REF] Liu | Experimental investigation on three-dimensional structures of wind wave surfaces[END_REF] and direct numerical simulations [START_REF] Li | The principal stage in wind-wave generation[END_REF][START_REF] Lin | Direct numerical simulation of wind-wave generation processes[END_REF][START_REF] Zonta | Growth and spectra of gravity-capillary waves in countercurrent air/water turbulent flow[END_REF][START_REF] Wu | Wind wave growth in the viscous regime[END_REF][START_REF] Wu | Revisiting wind wave growth with fully coupled direct numerical simulations[END_REF][START_REF] Burdairon | A combined vof-rans approach for studying the evolution of incipient two-dimensional wind-driven waves over a viscous liquid[END_REF], revisiting this old problem with new, high quality data.

In our experiments, the spatial evolution of mechanically generated waves are measured using Free-Surface Synthetic Schlieren [START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF], a refraction-based optical method offering a wave slope resolution of 3 × 10 -3 . To isolate the spatial evolution of the wave at each forcing frequency, we use a dedicated spatiotemporal band-pass filtering. This method enables us to measure the growth or damping rate for each forcing frequency, even for cases above wind onset where natural waves are also present, from which we determine the marginal stability curve of the instability with unprecedented resolution. We show that Miles' model, which is commonly applied to water waves, offers a reasonable description of the growth rate for more viscous liquids, and we propose a scaling for the growth rate of the most amplified wave with the liquid viscosity and friction velocity. We finally discuss the difference between the critical friction velocity that triggers wave growth in the case of mechanically generated waves and natural waves.

II. WAVE GROWTH AND ATTENUATION MECHANISMS

The spatial evolution of wind-generated waves results from a combination of energy transferred from the air flow and dissipated in the liquid. We first recall the modeling of these two contributions. We start from a simple harmonic surface deformation, ζ (x, t) = ζ 0 e i(kx-ωt) , of wave number k and angular frequency ω, propagating over a liquid layer of depth h, density ρ ℓ , and surface tension σ. In the absence of dissipation and forcing, ω and k are real and are related by the dispersion relation

ω 2 = gk + σk 3 ρ ℓ tanh (kh) = gk 1 + k 2 k 2 cap tanh (kh) , (1)
where k cap = 2π/λ cap = ρ ℓ g/σ is the capillary wave number. The energy density per unit surface of the wave writes

E = ρ ℓ gζ 2 0 [1 + (k/k cap ) 2 ].
Neglecting nonlinear interactions and surface currents, the transport and growth of the wave energy are governed by [START_REF] Kahma | A laboratory study of the minimum wind speed for wind wave generation[END_REF][START_REF] Grare | Growth and dissipation of wind-forced, deep-water waves[END_REF][START_REF] Peirson | On the wind-induced growth of slow water waves of finite steepness[END_REF]

∂ ∂t + c g ∂ ∂x E = P -D, (2)
with c g = ∂ω/∂k the group velocity, P the local power injected by the wind and D the local power dissipated in the liquid by viscosity. The injected power is the work per unit time of the wind stress (pressure and viscous shear stress) applied on the liquid surface, and can be modeled in two ways, depending on the wave amplitude. For very low wave amplitude, typically smaller than a fraction of the viscous sublayer thickness, the wavy shape of the liquid surface does not modify the flow in the air. In this case, only the turbulent pressure fluctuations in the air deform the interface, resulting in an injected power P independent of the wave energy E, and hence a linear growth of E: this is the starting hypothesis of the Phillips model [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF]. On the other hand, for larger wave amplitude, a feedback of the wave on the air flow takes place, and the wind stress becomes modulated by the wave profile. The resulting injected power P is now proportional to E, yielding an exponential growth of the wave energy: this is the Miles model [START_REF] Miles | On the generation of surface waves by shear[END_REF].

In the present paper, we are interested in mechanically generated waves of initial amplitude larger than the viscous sublayer thickness, which places us directly in the second stage of wave growth. In this regime, both P and D are proportional to E, so the right-hand-side of Eq. ( 2) can be written in the form P -D = γE, with γ the temporal growth rate. In this framework, the onset of wave growth is defined by γ = 0, when the injected power overcomes the viscous dissipation.

We first focus on the injected power P . Since the flow in the air is necessarily turbulent in our problem, the wave-induced wind stress variations (pressure and viscous shear stress) are conveniently scaled with the applied mean shear stress τ = ρ a u * 2 , with u * the friction velocity and ρ a the air density. In a developing boundary layer, this friction velocity u * depends on the free-stream wind U a and the fetch x (distance along which the stress is applied) [START_REF] Schlichting | Boundary Layer Theory[END_REF].

To first order in wave slope, since the velocity of the surface is essentially vertical, the work per unit time of the horizontal component of the stress (shear stress) can be neglected, and only the vertical component (pressure) contributes to the injected power. This contribution writes P = ℜ⟨p † ζ⟩, where p † is the complex conjugate of the oscillating component of the pressure along the surface and ζ the vertical velocity of the surface [1,[START_REF] Miles | On the generation of surface waves by shear[END_REF]. It therefore depends on the phase relation between the surface elevation and the wave-induced pressure oscillation. Following Miles [START_REF] Miles | On the generation of surface waves by shear[END_REF], we can write the oscillating pressure in the form p = (α + iβ)ρ a u * 2 kζ to first order in wave slope, with α + iβ a nondimensional parameter coding for this phase shift [START_REF] Wu | Revisiting wind wave growth with fully coupled direct numerical simulations[END_REF]. Only the pressure contribution in phase with the wave slope (β ̸ = 0) provides work, yielding a growth rate γ = P/E in the form [START_REF] Miles | On the generation of surface waves by shear[END_REF][START_REF] Plant | A relationship between wind stress and wave slope[END_REF]

γ = βω ρ a ρ ℓ u * c 2 . (3)
The normalized phase velocity c/u * is often referred to as the "wave age", because young waves are short and have small phase velocity c. The key result of Miles theory [START_REF] Miles | On the generation of surface waves by shear[END_REF][START_REF] Lighthill | Physical interpretation of the mathematical theory of wave generation by wind[END_REF] is that β is proportional to the curvature -U ′′ (z c ), with z c the elevation of the critical layer where the wind velocity U (z c ) matches the phase velocity of the considered wave. This indicates that too small or too large wavelengths, for which z c is such that U ′′ (z c ) vanishes (in the viscous sublayer or outside the boundary layer), cannot be amplified by this mechanism. The main difficulty in the determination of β therefore relies in the detailed knowledge of the velocity profile in the air.

Since the original work of Miles, the experimental verification of the scaling (3), and the question of the extent to which the parameter β can be considered a constant, have been the subject of numerous studies. For naturally growing waves, β corresponds to the most amplified wave number, but it is possible to examine its dependence with k using mechanically generated waves. Compilation of field and laboratory experiments suggests that the scaling (u * /c) 2 approximately holds for the air-water interface over a significant range of wave ages (1 < c/u * < 20), with an average value β ≃ 32 ± 16 [1,[START_REF] Plant | A relationship between wind stress and wave slope[END_REF]. This large variability can be ascribed to a number of parameters, such as the dependence with wave number, the presence of currents, the incorrect modeling of wave dissipation for nonlaminar flow in the liquid, nonlinear effects (dependence with wave slope [START_REF] Peirson | On the wind-induced growth of slow water waves of finite steepness[END_REF]), etc.

We now describe the dissipation in the liquid, focusing on the laminar dissipation mechanisms, which are relevant for the viscous liquids considered here. Viscosity introduces a negative imaginary part in the dispersion relation, ω = ω r + iω i , yielding a temporal damping rate γ = -D/E = 2ω i < 0 of the wave energy. The full dispersion relation for arbitrary viscosity (without side effects) is given by Lamb [START_REF] Lamb | Hydrodynamics[END_REF] in infinite depth, later improved by LeBlond and Mainardi [START_REF] Leblond | The viscous damping of capillary-gravity waves[END_REF] for arbitrary depth. For moderate viscosity (including that considered in the present work), the real part ω r remains close to the inviscid dispersion relation (1), i.e. the phase velocity c = ω r /k is marginally reduced by viscosity. The damping rate γ = 2ω i can be approximately split in two contributions: bulk dissipation due to the internal shear stress,

γ bulk = -4ν ℓ k 2 , (4)
and dissipation at the bottom wall due to the shear stress in the Stokes boundary layer of thickness δ = ν ℓ /ω r (provided that δ is much smaller than h) [START_REF] Lighthill | Waves in fluids[END_REF],

γ BW = - √ 2ν ℓ ω r k sinh 2kh . ( 5
)
The bulk dissipation is dominant at small wavelength, while the bottom-wall dissipation is dominant for wavelength much larger than liquid depth. Finally, for waves in a channel of finite width W , a third contribution arises from the friction with the side walls [START_REF] Tsai | On the spatial linear growth of gravity-capillary water waves sheared by a laminar air flow[END_REF][START_REF] Hunt | Viscous damping of waves over an inclined bed in a channel of finite width[END_REF][START_REF] Peirson | Rain-induced attenuation of deep-water waves[END_REF],

γ SW = - √ 2ν ℓ ω r W . (6)
In the following, we use the numerically resolved LeBlond-Mainardi solution, noted γ LM , which has the bulk and bottom-wall approximations (4) and ( 5) as asymptotic solutions for large and small k, respectively, to which we add the side-wall contribution γ SW [START_REF] Boomkamp | Classification of instabilities in parallel two-phase flow[END_REF]. The relative magnitude of these contributions, for the specific geometry of our experimental setup, is discussed in the Appendix. To summarize, the net temporal growth rate γ accessible to experiments is the sum of the positive contribution (3) and a combination of negative contributions (4) and [START_REF] Boomkamp | Classification of instabilities in parallel two-phase flow[END_REF]. In most practical configurations (including the present paper), waves are forced under permanent conditions and develop spatially, so the energy transport equation ( 2) reduces to

d dx E = γ s E,
with γ s = γ/c g the spatial growth rate [START_REF] Gaster | A note on the instability between temporally-increasing and spatially-increasing disturbances in hydrodynamic stability[END_REF]. Assuming homogeneous forcing and dissipation (independent of x), this yields an exponential variation of the energy profile, E(x) = E 0 e γsx , with γ s > 0 for growing waves and γ s < 0 for damped waves. In the following, we perform measurements of γ s by fitting E(x) in the range of exponential growth or decay for various k and friction velocities u * , and discuss the relevance of the various growth and damping rate contributions in the case of viscous liquids.

III. EXPERIMENTS AND MEASUREMENTS

A. Experimental setup

The experimental setup is composed of a rectangular tank of length L = 1.5 m and width W = 296 mm, located at the bottom of a wind tunnel; see Fig. 1 and Refs. [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Aulnette | Wind-sustained viscous solitons[END_REF][START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF] for details. The wind velocity U a , measured in the center of the cross-section, ranges from 0 to 7 m s -1 . The tank, of depth h = 35 mm, is filled with silicon oil (Rhodorsil ® oil 47 V20 and V50). Two kinematic viscosities are used: ν ℓ = 20.5 mm 2 s -1 and ν ℓ = 50.1 ± 0.1 mm 2 s -1 . The experiments are carried out at temperature of 24.0 ± 1.0 • C, giving an oil density ρ ℓ = 956 kg m -3 and surface tension σ = 20.7 mN m -1 . The corresponding capillary wavelength and capillary-gravity frequency are λ cap = 9.33 mm and f cap = 18.3 Hz, and the minimum phase velocity is c min = (4gγ/ρ) 1/4 = 0.171 m/s.

Waves are generated using an immersed wave maker, made of a vertically oscillating stainless steel cylinder of diameter 5 mm and length matching the width of the tank. The wave maker is located 125 mm after the beginning of the tank and defines the origin x = 0 of our measurements. It is hung by a rectangular frame and oscillated using an electromagnetic shaker (Sinocera JZK-20) via a crank rod system [Fig. 1(b)]. The upper part of the cylinder is at z(t) = z 0 + A cos(ωt), with z 0 = -14 mm the mean depth and A the stroke amplitude. In order to minimize wave reflection, two inclined perforated planes of length 85 mm and slope 30 • are added at both ends of the tank to absorb the wave energy.

The forcing frequency f = ω/2π ranges from 3 to 12 Hz. This range is below the capillary-gravity crossover frequency f cap , indicating that excited waves are in the gravity regime. The corresponding wavelengths range from 150 mm to 15 mm. The frequency range is limited at low frequency by the fact that the small diameter of the cylinder cannot excite efficiently large wavelengths, and at large frequency by the fact that the wavelength and the attenuation length become too small to be measurable. The Reynolds numbers Re = λ cap c min /ν ℓ , which characterize gravito-capillary waves of wavelength λ cap propagating at the minimum phase velocity c min , are Re = 78 and 32 for ν ℓ = 20.5 mm 2 s -1 and ν ℓ = 50.1 mm 2 s -1 , respectively. Since the wavelengths considered in the following are significantly larger than λ cap , these values represent lower bounds for the actual Reynolds numbers, confirming that the waves are in the weakly damped regime [START_REF] Leblond | The viscous damping of capillary-gravity waves[END_REF].

The stroke amplitude of the wave maker is kept fixed to A = 1 mm in all this study. This forcing induces waves of initial amplitude ζ 0 between 0.1 and 0.4 mm depending on the forcing frequency. This initial amplitude is small enough to satisfy the linear wave approximation (the maximum wave slope near the wave maker remains below kζ 0 ≃ 0.02). On the other hand, this initial amplitude is large enough to be directly in the Miles' exponential growth regime: the characteristic surface deformation corresponding to the transition from the Phillips' incoherent wrinkles to the Miles' regular waves is estimated to ζ ≃ (0.11 -0.14)δ ν in Refs. [START_REF] Li | The principal stage in wind-wave generation[END_REF][START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF], with δ ν = ν a /u * the thickness of the viscous sublayer, which lies in the range 40 -100 µm.

B. Friction velocity

In order to check the scaling of the growth rate (3), it is necessary to determine accurately the friction velocity u * for a given wind velocity U a . We deduce u * from the stationary current U s that develops at the liquid surface because of stress continuity ρ a u * 2 = ρ ℓ ν ℓ ∂u x /∂z| z=0 . The flow in the liquid being laminar (the maximum Reynolds number U s h/ν ℓ is less than 70), a stationary Couette-Poiseuille flow of zero flow rate develops, with a surface current U s given by [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF]

U s = ρ a u * 2 h 4ρ ℓ ν ℓ . ( 7
)
Close to the onset of wave growth, this surface current is of the order of 0.02 m s -1 , which is less than 10% of the typical phase velocity. In Aulnette et al. [START_REF] Aulnette | Wind-sustained viscous solitons[END_REF], we show that the friction velocity u * deduced from U s at various x locations using Eq. ( 7) is in good agreement with the empirical relation for a developing turbulent boundary layer -see, e.g., Eq. (21.12) in Schlichting [START_REF] Schlichting | Boundary Layer Theory[END_REF] u

* 2 (x) U a 2 ≃ C (x + x 0 )U a ν a -0.2 , (8)
with x 0 = 475 mm the distance between the beginning of the flat plate at the end of the windtunnel convergent and the wave maker, and C = 0.029. This decrease of u * along x is related to the development of the boundary layer along the tank, with a thickness that increases as δ(x) ≃ x 0.8 [START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF].

In the following, we use Eq. ( 8) to compute the friction velocity u * from the prescribed wind velocity U a close to the wave-maker, in the range x ≃ 100 -300 mm where the growth rate is measured. The uncertainty on u * is of the order of 5% [START_REF] Aulnette | Wind-sustained viscous solitons[END_REF].

C. Wave measurement

The surface deformations are determined using Free-Surface Synthetic Schlieren (FS-SS) [START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF]. This method is based on the analysis of the refracted image of a random-dot pattern through the interface. The dots are 1 mm in diameter (corresponding to 2.3 pixels) with a density of 60 dots per cm 2 . The pattern is printed on a transparent film located at the bottom of the liquid tank and illuminated from below by a LED panel. Images of the pattern are taken through the transparent upper-wall of the wind tunnel. To reduce parallax errors, images are acquired through a 45 o mirror located above the wave tank by a camera located 3 m from the liquid surface [Fig. 1(b)]. Images are acquired using a high-speed camera (Photron FASTCAM Mini WX50), fitted with a Canon 85 mm f /1.8 macro lens. Images of dimension of 890 × 320 mm 2 are recorded at a frame rate at least 10 times larger than the wave maker frequency f , for a recording duration of at least 20 wave periods.

The surface slope field is obtained by computing the apparent displacement field δr of the dot pattern induced by the surface deformation using an image correlation algorithm [START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF]. In the limit of weak slopes and within the paraxial approximation, the displacement field is proportional to the surface height gradient, δr = -h * ∇ζ, with h * = 14 mm an effective distance that includes the surface-pattern distance (h p = 48.2 mm) and the refraction indices of the liquid and intermediate layers between the surface and the pattern [START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF]. Although the surface height field ζ(r, t) can be reconstructed by inverting the ∇ operator, as done in Refs. [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF], here we work directly with the refraction-induced displacement field δr: for a quasimonochromatic wave propagating in the x direction, ζ(x, t) = ζ 0 (x) cos(kx -ωt), where ζ 0 (x) is a slowly varying amplitude, the displacement is simply proportional to the wave height,

δx = -h * ∂ζ ∂x ≃ h * kζ 0 (x) sin(kx -ωt),
so the spatial growth rate of the wave amplitude ζ 0 (x) can be obtained directly from that of the displacement δx.

The image correlation algorithm is based on interrogation windows of size 12 pixels (5.2 mm) with 50% overlap, yielding a spatial resolution of 2.6 mm. This is approximately 6 times smaller than the smallest wavelength considered here, λ ≃ 15 mm for f = 12 Hz. The resolution in the apparent displacement δx is 0.1 pixel (0.04 mm), yielding a resolution in surface slope kζ ≃ 0.003.

Figures 2(a) and (b) show two snapshots of the magnitude of the wave gradient field ∇ζ = -δr/h * measured at a wind velocity U a = 5.0 m s -1 in the case of an amplified wave (f = 5 Hz, left) and damped wave (f = 9.5 Hz, right), for the oil of viscosity ν ℓ = 20.5 mm 2 s -1 . The slight tilt observed in the wave crests originates from the slightly nonhomogeneous wind in the y direction, of the order of 3% [START_REF] Aulnette | Ondes non linéaires générées par le vent à la surface d'un liquide visqueux[END_REF]. From the wave slope in the x direction, we compute the mean-square envelope, which (c,d) Snapshots of the filtered surface gradient |∇ ζ| using the spatiotemporal filter [START_REF] Li | The principal stage in wind-wave generation[END_REF]. (e,f) Wave slope profiles ∂ζ/∂x (red dashed curves: raw data; black continuous curves: filtered data). (g,h) Mean square envelopes ⟨(∂ζ/∂x) 2 ⟩ (raw and filtered). The two vertical lines show the range used for the exponential fit [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF].

writes for a quasi-monochromatic wave as

⟨(∂ζ/∂x) 2 ⟩ ≃ 1 2 k 2 ζ 2 0 (x), (9)
with ζ 0 (x) the slowly varying wave amplitude and ⟨•⟩ the average over time and y. Examples of wave slope profiles ∂ζ/∂x at fixed y and their corresponding mean-square envelopes are shown as red dashed lines in Figs. 2(e)-2(h). In the amplified case, the envelope shows a well defined exponential growth at moderate fetch, followed by a saturation at large fetch [Fig. 2(g)]. This saturation, which results from a combination of nonlinear interactions with other waves and the small decrease of the wind shear stress along the developing boundary layer, is not considered in the following. The damped case is more complex [Fig. 2(h)]: after a short range of exponential decay, the envelope saturates or even slightly increases at large fetch, because of the presence of other waves of frequencies different from the forcing frequency. These unforced waves are the amplified waves naturally present for a wind velocity above the critical threshold, and may hide the underlying decay of the forced wave, which makes it difficult to accurately measure the damping rate.

To extract the forced waves from the sea of ambient waves present in the field, we perform a spatiotemporal band-pass filtering of the displacement field δx = -h * ∂ζ/∂x. Figure 3(a) shows the power spectrum |δ x| 2 (k, ω) (averaged in the y direction) in the case of a damped wave. The displacement field being real, its Fourier transform is symmetric under (k, ω) → (-k, -ω). The quadrants (+, +) and (-, -) correspond to waves propagating in the x > 0 direction, while the quadrants (+, -) and (-, +) correspond to waves propagating in the x < 0 direction. In addition to the injected energy around the forcing ±(k 0 , ω 0 ) (red ellipses), energy is also visible in the reflected waves ±(-k 0 , ω 0 ), in harmonics, and along the viscous dispersion relation (white dashed lines). Energy along the dispersion relation is systematically found for a wind velocity above the onset, and is due to the natural waves growing from initial disturbances and amplified by the wind.

To filter out the undesired wave components, we compute for each y the space-time Fourier transform δ x(k, ω), convolute it with a Butterworth filter kernel

G (k, ω) = H(kω) 1 + k-k0 ∆k 2 + ω-ω0 ∆ω 2 N/2 , (10)
and compute the inverse Fourier transform to reconstruct the filtered displacement field δ x(x, t), and hence the filtered wave slope ∂ ζ/∂x. We choose a filter kernel centered around the forcing ±(k 0 , ω 0 ), with a spectral bandwidth (∆k, ∆ω) and a decay characterized by the exponent N , chosen here equal to 6. A Heaviside operator H(kω) is included to remove the reflected waves (kω < 0). A narrow temporal bandwidth ∆ω is chosen to select precisely the forced waves, but a wider spatial bandwidth ∆k is chosen to correctly measure the spatial decay of the wave. The temporal bandwidth is chosen as ±3 spectral points π/T acq (with T acq the acquisition duration) surrounding the forcing frequency ω 0 . The spatial bandwidth ∆k is chosen equal to k 0 , so the smallest measurable damping length is of the order of the wavelength. The resulting filtered spectrum |G(k, ω)δ x| 2 (k, ω) in Fig. 3(b) shows that the energy of the forced wave is correctly selected by these parameters.

The snapshots for the filtered instantaneous wave gradient ∇ ζ = -δr/h * , shown in Figs. 2(c)-2(d), are very close to the raw data in the amplified case, but are significantly modified in the damped case, indicating that the naturally amplified waves are efficiently filtered out by our procedure. The corresponding wave slope profiles and envelopes are also shown as continuous black lines in Figs. 2(e

)-2(h).

To validate our procedure over the full range of wind velocity, we show in Fig. 4 a series of meansquare filtered wave slopes at a fixed frequency (f = 4.1 Hz) for increasing U a , illustrating the gradual transition from damped to amplified waves. Close to the wave maker, the growth or decay is well described by an exponential function, from which we fit the spatial growth rate γ s ,

⟨(∂ ζ/∂x) 2 ⟩ ∝ e γsx . ( 11
)
We choose the range x ∈ [1.8, 4]λ for the fits [shown by the two vertical lines in Figs. 2(g) and 2(h) and 4], for which a robust exponential fit can be performed at all f and U a . The overall uncertainty on γ s is less than 8%. ). The vertical lines show the range used to fit the spatial growth rate γs using Eq. [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF].

IV. RESULTS

A. Phase velocity

To compare the spatial growth rates γ s determined experimentally to the temporal growth rates γ predicted theoretically, it is necessary to convert the measured spatial variations of the wave energy to the corresponding temporal variations. In the limit of small viscosity, the ratio between the two growth rates is given by the group velocity [START_REF] Gaster | A note on the instability between temporally-increasing and spatially-increasing disturbances in hydrodynamic stability[END_REF], which cannot be directly measured in our system. We must therefore use the group velocity obtained from the dispersion relation of the free waves in a viscous liquid, i.e. including viscous effects but without wind, as discussed in Sec. II. To check to what extent we can rely on this viscous free-wave dispersion relation, we first measure here the phase velocity of the waves and compare it with this prediction.

Figures 5(a) and 5(b) show the phase velocity c exp = ω/k for the two liquid viscosities as a function of the wave number for various wind velocities U a (here k is simply measured from the distance between wave crests, ensuring high accuracy). The data are compared to the inviscid prediction (1) in solid line, and to the numerically computed viscous prediction of LeBlond and Mainardi [START_REF] Leblond | The viscous damping of capillary-gravity waves[END_REF] including depth effects in dashed line. These predictions are very close in our range of frequencies (3 -12 Hz), within 2%, so we will simply consider the inviscid prediction in the following. At low frequency, the inviscid phase velocity tends to the nondispersive shallow water limit √ gh ≃ 0.58 m/s, while the viscous prediction falls to 0 when the thickness of the Stokes boundary layer becomes of the order of the liquid depth (this cutoff is well below the frequencies relevant to our study).

For our range of frequencies, the measured phase velocity without wind, c 0-exp , is almost indistinguishable from the prediction (either inviscid or viscous). In the presence of wind, it still remains close to the prediction, to within 3% at ν ℓ = 20.5 mm 2 s -1 and 8% at ν ℓ = 50.1 mm 2 s -1 . We can therefore use with reasonable accuracy the predicted group velocity computed from the inviscid dispersion relation (1) to convert spatial growth rates to temporal ones in the following.

Before further proceeding, it is interesting to analyze in more detail the small differences between the measured phase velocity with and without wind (c exp and c 0-exp , respectively). The difference c exp -c 0-exp is plotted in Fig. 5(c) and 5(d), with uncertainties shown as shaded areas. It can be explained by a combination of two antagonistic effects [START_REF] Valenzuela | The growth of gravity-capillary waves in a coupled shear flow[END_REF], the importance of which depending on the liquid viscosity. The first one corresponds to an increase of the phase velocity by the surface current induced by the wind shear stress, which dominates at small viscosity, and the second one corresponds to a decrease the phase velocity due to the wave-induced aerodynamic pressure, which dominates at large viscosity.

The increase in phase velocity by surface current is simply given by U s for short wavelengths [see Eq. ( 7), with U s ≃ 0.02 m s -1 close to the wave onset], but is less pronounced for long wavelengths because the perturbation they induce penetrates deeper into the liquid, where the current is lower. This transport effect therefore depends on the mean velocity profile in the liquid. For the Couette-Poiseuille profile considered here [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF], the modification in the phase velocity has been computed by Lilly, and can be found in the appendix of the paper by Hidy and Plate [START_REF] Hidy | Wind action on water standing in a laboratory channel[END_REF],

∆c w (k) = U s 1 - 1 + 2 cosh 2kh kh sinh 2kh + 3 2k 2 h 2 . ( 12
)
The second effect is the slowdown of the phase velocity by the variations of the aerodynamic pressure over the wave profile, and can be explained qualitatively as follows. The phase velocity of gravity waves is governed by gravity, i.e. by the weight of the liquid deformation. In the presence of wind, a fraction of this weight is carried by the wave-induced aerodynamic pressure (Bernoulli suction over the wave crests and over-pressure on the wave troughs), of the order of ρ a U 2 a kζ. The surface deformation therefore experiences an effective reduced gravity, yielding a reduced phase velocity, (c,d) Difference between the phase velocity with and without wind. Solid lines are the predictions including the (positive) surface drift correction [START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF] and the (negative) aerodynamic pressure correction [START_REF] Nové-Josserand | Effect of a weak current on wind-generated waves in the wrinkle regime[END_REF], with uncertainties shown as color shaded areas.

which can be written as [START_REF] Plant | Phase speeds of upwind and downwind traveling short gravity waves[END_REF]

c(k) = c 0 (k) 2 - ρ a ρ ℓ U (z k ) 2 , ( 13
)
with U (z k ) the wind velocity evaluated at a characteristic elevation that depends on the wave number k. We evaluate this effect by considering a classical logarithmic profile U (z) = u * [κ -1 log(zu * /ν a )+B] (with κ = 0.41 and B = 3.5) computed at the elevation z k = ϵ/k with ϵ ≃ 0.23 [START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF]. The correction to the phase velocity, combining the surface current effect [START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF] and the aerodynamic pressure effect [START_REF] Nové-Josserand | Effect of a weak current on wind-generated waves in the wrinkle regime[END_REF], is plotted in Figs. 5(c) and 5 (d). It provides a reasonable description of the data: At low viscosity and large wave numbers, the waves are essentially transported by the surface current, yielding a phase velocity larger than the inviscid prediction. On the other hand, waves at low wave numbers are systematically slower than the free-wave celerity, as a result of the dominant aerodynamic pressure effect.

To summarize, the phase velocity of the waves in presence of wind can be correctly described by the finite-depth inviscid phase velocity, by including the corrections due to the surface drift and the aerodynamic pressure. However, in the range of liquid viscosities and wind velocities considered here, these corrections remain small: for a wind velocity close to the wave onset, the measured phase velocity matches the free-wave inviscid prediction to within 10%. By extension, we can consider that the group velocity is also only marginally affected by the surface current and aerodynamic pressure, which allows us to use the free-wave inviscid prediction c g derived from Eq. (1) to infer the temporal growth rate γ from the measured spatial growth rate γ s .

B. Spatial growth rates

Figures 6(a) and 6(b) shows the measured spatial growth rate γ s as a function of the measured wave number, for different wind velocities. As expected, without wind γ s is negative for all k. The corresponding damping rate is compared to the bulk prediction γ s = -4ν ℓ k 2 /c g (solid line) and to the numerically solved LeBlond-Mainardi prediction [START_REF] Leblond | The viscous damping of capillary-gravity waves[END_REF] (including bulk and bottom-wall dissipation) completed with side-wall dissipation (6) (dashed line). The large wave numbers (k > 70 m -1 ) are equally well described by the two models, but the smaller k are much better described by the refined model, especially at large viscosity, highlighting the significant effect of the bottom dissipation at small k (see Appendix V).

As the wind velocity is increased, the waves become gradually less damped and, above a critical wind u * c , the growth rate becomes positive in a finite range of wave numbers. We recall that the friction velocity u * is deduced from the free-stream velocity U a by using Eq. ( 8) for x close to the wave maker. From the set of growth rate curves γ s (k) at various velocities, we can compute the marginal stability curve u * (k), defined as the minimum friction velocity required for a given wave number k to become unstable. It is plotted in Figs. 6(c) and 6(d) for the two liquid viscosities. The lowest point (k c , u * c ) of the marginal stability curve defines the critical wave number and critical friction velocity, which are summarized in Table I for the two liquid viscosities. The critical wave numbers correspond to wavelengths λ c ≃ 60 -80 mm, much larger than the capillary wavelength λ cap = 9.33 mm.

As the liquid viscosity is increased, the critical friction velocity strongly increases, but the critical wave number is marginally modified. Increasing u * shifts the most amplified wave number (shown as blue triangles) to larger values, and rapidly widens the domain of unstable wave numbers. The low-k branch of the stability curve, for which the bottom-wall and the bulk dissipation are of comparable importance (see the Appendix), shows a weak dependence with u * . On the other hand, the large-k branch, which is essentially governed by the bulk dissipation, strongly varies with u * . We note that this higher branch is defined with high accuracy, because γ s sharply crosses 0, but that the lower branch is more difficult to measure, because γ s shows little variations with u * for small k.

C. Comparison with Miles' scaling

We now investigate to what extent Miles' model, usually applied to the air-water interface, also applies to the case of the more viscous liquids considered here. Following Eq. (3), we plot in Fig. 7(a) and 7(b) the temporal growth rate γ = γ s c g normalized by the wave frequency ω as a function of the inverse wave age u * /c (only amplified waves are considered here). In this traditional representation, "young" waves (high k, low c) are on the right and "old" waves (small k, large c) are on the left. However we note that, strictly speaking, there is wave aging in our experiments: each data point corresponds to a single wave excited at a given frequency.

Although individual curves for each u * fail to follow the (u * /c) 2 scaling, the upper envelope of the set of curves is compatible with this scaling. The strong departure at small and large wave age results from a combination of dissipation effects (dominated by the bulk dissipation for large u * /c and bottom-wall dissipation for small u * /c) and the selective amplification of the Miles' mechanism: wave numbers k far from the most amplified wave numbers correspond to a height where the curvature of the mean velocity U ′′ (z), and consequently the growth rate, becomes zero.

To further test the (u /c) 2 scaling, is therefore necessary to subtract the contribution from viscous dissipation. In order not to on specific dissipation model, we use the damping rate γ 0 without wind as the best estimate for the dissipation. We therefore plot the normalized FIG. 8. Dimensionless net growth rate (γ -γ0)/f (with γ0 the damping rate without wind) vs the inverse wave age, showing our data (red circles: ν ℓ = 20.5 mm 2 s -1 , blue triangles: ν ℓ = 50.1 mm 2 s -1 ) superimposed to data obtained in water and compiled by Plant [START_REF] Plant | A relationship between wind stress and wave slope[END_REF] (see also [1,[START_REF] Komen | Dynamics and modelling of ocean waves[END_REF][START_REF] Melville | The equilibrium dynamics and statistics of gravity-capillary waves[END_REF]). Open symbols: field data; filled symbols and crosses: laboratory data; solid line: Miles's prediction [START_REF] Miles | On the generation of surface waves by shear flows. Part 2[END_REF].

net growth rate (γ -γ 0 )/ω in Figs. 7(c) and 7(d), which is now consistent with the scaling (u * /c) 2 , at least for large wind velocity and large wave numbers (large u * /c). The best fit with Eq. (3) gives the average value β ≃ 45 with an uncertainty range β ≃ 23 -75 (shown in the shaded area).

The significant scatter in the dimensionless growth rates is not specific to our experiments performed with viscous oils, and is also present in experiments in the air-water case. In Fig. 8 we superimpose our normalized net growth rates to the classical Plant's compilation of laboratory and field measurements [START_REF] Plant | A relationship between wind stress and wave slope[END_REF], reproduced from Janssen [1]. The correct agreement between our data in viscous oils and the data in water supports the robustness of the scaling (u * /c) 2 at large u * /c, with a comparable scatter. The solid line shows the prediction of Miles [START_REF] Miles | On the generation of surface waves by shear flows. Part 2[END_REF], which gives the asymptotic scaling (u * /c) 2 at large u * /c, relevant to our "young" (slow) waves, but a lower growth rate for "older" (faster) waves when the phase velocity becomes of the order of the wind velocity.

D. Maximum growth rate

We now turn to the maximum growth rate γ max (growth rate of the most amplified wave number) and its dependence on the liquid viscosity. Figure 9(a) shows γ max as a function of the friction velocity u * for the two viscosities. By definition γ max increases with u * and crosses 0 at the critical friction velocity u * c (given in Table I). We note that negative values of γ max can be defined for u * < u * c because slightly below u * c the growth rate shows a local (negative) maximum, corresponding to the least damped wave number.

Given the various sources of dissipation, with distinct scaling laws with respect to ν ℓ (bulk dissipation in ν ℓ and wall dissipation in ν 1/2 ℓ ), a unique scaling for γ max with ν ℓ cannot be derived in general. I). (b) Viscosity-weighted maximum temporal growth rate γmaxν 3 ℓ as a function of the friction velocity u * . The line corresponds to Eq. ( 15) with β = 45, and the shaded area shows the range β = 23 -75.

However, an approximate trend can be proposed, at the cost of a number of approximations: (1) we ignore the k dependence of the parameter β in Eq. ( 3), which is acceptable close to the most amplified k, and consider the average value β = 45; (2) we consider that waves are in the deep-water gravity regime, ω = √ gk, which again is acceptable because k c /k cap < 0.15 and tanh(k c h) ≃ 0.99; and (3) we consider the bulk dissipation γ bulk = -4ν ℓ k 2 as the only source of dissipation. This third approximation is the most questionable in our experiments, with a bulk-bottom transition at k ≃ 70 m -1 (see the Appendix), which is close to the critical wave number k c . Under these approximations, the temporal growth rate can be modeled as

γ ≃ βs u * 2 k 3/2 g 1/2 -4ν ℓ k 2 , (14)
with s = ρ a /ρ ℓ the density ratio. Solving for ∂γ/∂k = 0 yields the maximum growth rate

γ max = 3 3 2 14 β4 s 4 u * 8 g 2 ν 3 ℓ . (15)
To verify this scaling, we plot in Fig. 9(b) the combination γ max ν 3 ℓ as a function of u * . The data for the two viscosities approximately gather along the power law u * 8 , with variations of the same order as the uncertainty on the parameter β found previously (the black line gives the average value β = 45, and the shaded area shows the same range of β as in Fig. 7). Although the data for the two viscosities do not overlap, because of the very limited range of accessible u * (at most twice the critical value u * c ), the high-ν ℓ curve reasonably extends the low-ν ℓ curve, which is consistent with the scaling law [START_REF] Shemer | Spatially evolving regular water wave under the action of steady wind forcing[END_REF].

The variation of γ max with u * has been the subject of many studies, but only in the air-water case [START_REF] Larson | Wind-generated gravity-capillary waves: Laboratory measurements of temporal growth rates using microwave backscatter[END_REF][START_REF] Kawai | Generation of initial wavelets by instability of a coupled shear flow and their evolution to wind waves[END_REF][START_REF] Creamer | Surface films and wind wave growth[END_REF][START_REF] Tsai | Stability analysis on the initial surface-wave generation within an air-sea coupled shear flow[END_REF][START_REF] Zeisel | Viscous effects on wave generation by strong winds[END_REF]. Approximate power laws γ max ∼ u * n are reported, but with significantly smaller exponents n, in the range 1.5-3.5. The asymptotic scaling n = 8 may be difficult to access experimentally in water experiments because nonlinear effects (wave saturation, frequency downshift, wave breaking) rapidly appear as u * is increased, even close to the onset u * c .

E. Critical friction velocity

We finally consider the influence of the liquid viscosity on the critical friction velocity u * c for the growth of mechanically generated waves. This critical friction velocity cannot be directly inferred from Eq. ( 14), based on the sole bulk dissipation, because it predicts a range of unstable k starting from k = 0 as soon as u * > 0. This directly follows from the approximation β(k) ≃ β and the omission of other dissipation sources. In our experiments, the dissipation being given by a combination of bulk and bottom-wall contributions (see the Appendix), we can modify Eq. ( 14) by still assuming a constant coefficient β, while now including the Leblond-Mainardi damping rate

γ ≃ βs u * 2 k 3/2 g 1/2 + γ LM (k). (16)
The critical friction velocity u * c , obtained by solving numerically γ = 0, is plotted in Fig. 10 as a function of ν ℓ . The numerical solution, shown again for the average value β = 45 (black line) and its uncertainty range (shaded area), is in good agreement with the measured u * c for the two values of ν ℓ . Interestingly, the measured u * c and the numerical solution are well described by the approximate power law u * c ∼ ν 1/2 ℓ . Although this power law cannot be directly derived from Eq. ( 16), it is compatible with the fact that the most unstable wave number k c shows little variation with ν ℓ . Considering k c constant in Eq. ( 14) and solving for γ = 0 indeed yields

u * c ≃ 4ν ℓ ω c βs , (17)
with ω c = ck c the frequency of the critical wavenumber. This is equivalent to the relation proposed by Wu and Deike [START_REF] Wu | Wind wave growth in the viscous regime[END_REF], u * /c ≃ a √ ν l k 2 T , with T = 2π/ω the wave period and a a numerical factor. The best fit to our data using this relation yields a ≃ 3.0 ± 0.2, which is a factor 2 below the value a ≃ 6.76 found in the numerical simulations of Ref. [START_REF] Wu | Wind wave growth in the viscous regime[END_REF]. In spite of this difference, probably related to the difference in the velocity profile (which is linear in the numerics), this suggests that the scaling [START_REF] Hidy | Wind action on water standing in a laboratory channel[END_REF] for the critical friction velocity in the case of mechanically generated waves is robust.

We close this discussion by comparing this critical friction velocity and the one obtained for natural waves, without mechanical forcing. In Fig. 10 we add u * c for natural waves as measured in the same setup in Refs. [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF], using silicon oils and water-glycerol mixtures. Here u * c is defined as the transition from the small-amplitude wrinkles (Phillips regime) to the exponentially growing waves (Miles regime). This natural-wave threshold is larger than the forced-wave threshold, but shows a shallower dependence with the liquid viscosity, u * c ∼ ν 0.2 ℓ , with a very similar numerical prefactor (within 5%) for the two types of fluids.

The two-stage scenario provides a natural framework to interpret this difference: Mechanicallygenerated waves grow from a disturbance of initial amplitude sufficient to disturb the air flow and trigger the Miles instability mechanism. On the other hand, natural waves grow from a noisy base state defined by the low-amplitude incoherent wrinkles excited by the turbulent pressure fluctuations in the air, so they require a larger wind velocity for the wrinkle amplitude to disturb the flow and trigger the instability. More specifically, for u * < u * c , the wrinkle amplitude results from the balance between the power injected by the turbulent pressure fluctuations and the viscous dissipation in the liquid, yielding ζ w ≃ sδ(u * 3 /gν ℓ ) 1/2 , with δ the boundary layer thickness. Assuming a Phillips-Miles transition when the wrinkle amplitude ζ w becomes of the order of the viscous sublayer thickness δ ν = ν a /u * [START_REF] Li | The principal stage in wind-wave generation[END_REF][START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF], we obtain u * c ∼ ν

1/5
ℓ , in good agreement with the data in Fig. 10. The strong difference between these two scalings, u * c ∼ ν 1/2 ℓ for forced waves and u * c ∼ ν 1/5 ℓ for natural waves, emphasizes the difficulty in defining a unique onset for the wave generation problem, and the interest of varying the liquid viscosity to gain more insight into this problem. •, Natural waves in water-glycerol mixtures [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. △, Natural waves in silicon oils [START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF]. □, Forced waves in silicon oils (present data). The black line shows the numerical solution of γ = 0 using Eq. ( 16) for β = 45, and the shaded area shows the range β = 23 -75.

V. CONCLUSION

In this paper we investigated the wind-induced growth of mechanically generated waves over liquids of viscosity 20 and 50 times larger than that of water, with the aim to gain insight into the scaling of the growth rate in the exponential regime. This growth rate has been extensively studied in the literature but only on the air-water configuration. Varying the liquid viscosity offers an interesting opportunity to test existing theories and shed light on this problem.

Using Free-Surface Synthetic Schlieren with a spatio-temporal filtering provides wave slope measurements with a very good resolution, and the possibility to access the growth and damping rates of the forced waves even when masked by other naturally growing waves. From these measurements we reconstruct the marginal stability curves with unprecedented accuracy. Once the viscous dissipation is subtracted, our measurements show correct overall agreement with Miles' scaling, with a significant scatter comparable to that of experiments performed in water. Our analysis suggests a maximum growth rate (growth rate of the most amplified wave number) that scales as u * 8 ν -3 ℓ , in reasonable agreement with our data.

Our work raises the question of the origin of the critical friction velocity u * c for wave growth, and its dependence with liquid viscosity ν ℓ . For mechanically generated waves, provided that the initial wave amplitude is not too small, the system is directly in the second stage of wave growth (exponential amplification of the wave energy). In this case, the threshold is governed by the balance between the power injected by the wind and the dissipation in the liquid. On the other hand, for natural waves, a larger threshold u * c is found, because the wave instability is triggered from the low-amplitude incoherent wrinkles excited by the turbulent pressure fluctuations in the air [START_REF] Li | The principal stage in wind-wave generation[END_REF][START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF]. Previous works [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF][START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF] for mechanically generated waves, in qualitative agreement with recent numerical results [START_REF] Wu | Wind wave growth in the viscous regime[END_REF]. Extending these results to arbitrary liquid depth and to a larger range of viscosities, including the case of water, would be valuable to clarify this difference.

APPENDIX: RELATIVE IMPORTANCE OF THE VARIOUS CONTRIBUTIONS TO THE WAVE DISSIPATION

The damping rate of waves at the surface of a viscous liquid in a container, provided that the thickness of the boundary layers is much smaller than the container dimension, can be split into three contributions: bulk (γ bulk ), bottom-wall (γ BW ) and side-wall (γ SW ), given by Eqs. ( 4), [START_REF] Ayet | The dynamical coupling of wind-waves and atmospheric turbulence: A review of theoretical and phenomenological models[END_REF], and (6), respectively. The first two contributions are contained in the finite-depth damping rate γ LM for arbitrary viscosity derived by LeBlond and Mainardi [START_REF] Leblond | The viscous damping of capillary-gravity waves[END_REF], which has γ bulk and γ BW as asymptotic solutions at low viscosity for large and small k respectively.

To estimate which contributions dominate in our experiments (tank depth h = 35 mm, width W = 296 mm), we plot in Fig. 11 the spatial damping rates -γ s = -γ/c g from the numerically resolved Leblond-Mainardi dispersion relation (γ LM ) and from the three separate contributions (γ bulk , γ BW , γ SW ) for the two liquid viscosities ν ℓ = 20.5 and 50.1 mm 2 s -1 . In the range of accessible wave numbers, shown as gray bars, the bulk dissipation γ bulk (4) is the dominant contribution for k > 70 m -1 , while the bottom-wall contribution γ BW becomes significant for lower k. The side-wall contribution γ SW , shown by the dashed blue line, is never dominant, but it contributes approximately 20% to the wall dissipation at low k.

[1] P. Janssen, The interaction of ocean waves and wind (Cambridge University Press, 2004).

FIG. 1 .

 1 FIG. 1. Sketch of the experimental setup. (a) Side view; (b) cross sectional view. L = 1500 mm, h = 35 mm, H = 105 mm, W = 296 mm. The wave maker is a cylinder of 5 mm diameter, located at Lwm = 125 mm from the beginning of the tank and -14 mm below the liquid surface.

FIG. 2 .

 2 FIG. 2. Spatial evolution of the instantaneous and mean-square wave slope at wind velocity Ua = 5 m s -1 and oil viscosity ν ℓ = 20.5 mm 2 s -1 . Left column: f = 5 Hz, for which the wave is amplified; right column: f = 9.5 Hz, for which the wave is damped. (a,b) Snapshots of the raw surface gradient |∇ζ|.(c,d) Snapshots of the filtered surface gradient |∇ ζ| using the spatiotemporal filter[START_REF] Li | The principal stage in wind-wave generation[END_REF]. (e,f) Wave slope profiles ∂ζ/∂x (red dashed curves: raw data; black continuous curves: filtered data). (g,h) Mean square envelopes ⟨(∂ζ/∂x) 2 ⟩ (raw and filtered). The two vertical lines show the range used for the exponential fit[START_REF] Paquier | Viscosity effects in wind wave generation[END_REF].

  FIG. 3. (a) Power spectrum |δ x| 2 (k, ω) of the displacement field in logarithmic scale for Ua = 4.2 m s -1 , f = 7.3 Hz and oil viscosity ν ℓ = 20.5 mm 2 s -1 . The white dashed line shows the viscous dispersion relation [48]. The energy of the forced wave ±(k0, ω0) is selected by the spatiotemporal filter (10) (the red ellipses show the isovalue G = 0.05). (b) Resulting filtered spectrum |G(k, ω)δ x| 2 (k, ω).

FIG. 4 .

 4 FIG.4. Mean square envelopes ⟨(∂ ζ/∂x) 2 ⟩ (filtered data) at various wind velocities Ua, for ν ℓ = 50 mm 2 s -1 and f = 4.1 Hz (k = 69 mm -1 ). The vertical lines show the range used to fit the spatial growth rate γs using Eq.[START_REF] Paquier | Viscosity effects in wind wave generation[END_REF].

FIG. 5 .

 5 FIG.5. Phase velocity as a function of the wave number for various wind velocities and for the two liquid viscosities, (a) ν ℓ = 20.5 mm 2 s -1 and (b) ν ℓ = 50.1 mm 2 s -1 . Error bars correspond to the standard deviation of the data. The black lines show the phase velocity for free waves, either inviscid (solid line) or with viscosity (dashed line). (c,d) Difference between the phase velocity with and without wind. Solid lines are the predictions including the (positive) surface drift correction[START_REF] Perrard | Turbulent windprint on a liquid surface[END_REF] and the (negative) aerodynamic pressure correction[START_REF] Nové-Josserand | Effect of a weak current on wind-generated waves in the wrinkle regime[END_REF], with uncertainties shown as color shaded areas.

  FIG. 6. (a,b) Spatial growth rate γs as a function of the wave number for various wind velocities, for (a) ν ℓ = 20.5 mm 2 s -1 and (b) ν ℓ = 50.1 mm 2 s -1 . Black solid line: bulk dissipation (-4ν ℓ k 2 /cg); black dashed line: numerically solved LeBlond-Mainardi dissipation (bulk and bottom-wall dissipation) combined with side-wall dissipation. Color lines: local polynomial regression ("loess" smoothing). (c,d) Marginal stability curve (green circles) and most amplified wave number (blue triangles). Thin dashed lines are guides for the eyes.

FIG. 7 .

 7 FIG. 7. Dimensionless growth rate γ/ω as a function of the inverse wave age u * /c, for (a,c) ν ℓ = 20.5, and (b,d) 50.1 mm 2 s -1 . The solid line shows Eq. (3) for β = 45. (a, b) Raw growth rates, without viscous correction. (c, d) Growth rates corrected by the viscous contribution γ0 measured without wind. The shaded area depicts the range β = 23 -75.

u* 8 FIG. 9 .

 89 FIG. 9. (a)Maximum temporal growth rate γmax as a function of the friction velocity u * for ν ℓ = 20.5 mm 2 s -1 (red) and ν ℓ = 50.1 mm 2 s -1 (blue). Vertical dash-dotted lines show the critical friction velocities u * c (TableI). (b) Viscosity-weighted maximum temporal growth rate γmaxν 3 ℓ as a function of the friction velocity u * . The line corresponds to Eq. (15) with β = 45, and the shaded area shows the range β = 23 -75.

5 FIG. 10 .

 510 FIG. 10. Comparison between the critical friction velocity u *c for natural and forced waves. •, Natural waves in water-glycerol mixtures[START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. △, Natural waves in silicon oils[START_REF] Aulnette | Kelvin-Helmholtz instability and formation of viscous solitons on highly viscous liquids[END_REF]. □, Forced waves in silicon oils (present data). The black line shows the numerical solution of γ = 0 using Eq. (16) for β = 45, and the shaded area shows the range β = 23 -75.

suggested a dependence u * c ∼ ν 1/ 5 ℓ

 5 for natural waves, while the present work, available for two liquid viscosities only, suggests a dependence u * c ∼ ν 1/2 ℓ

2 FIG. 11 .

 211 FIG. 11. Various contributions to the spatial damping rates -γs = -γ/cg > 0 for (a) ν ℓ = 20.5 and (b) 50.1 mm 2 s -1 . LM indicates the numerically resolved Leblond-Mainardi solution, which includes the bulk dissipation and the bottom-wall dissipation. The gray horizontal bar shows the experimentally accessible range of wave numbers k.

TABLE I .

 I Critical parameters at the onset of wave generation for the two liquid viscosities.

	ν ℓ (mm 2 s -1 )	Uac (m s -1 )	u * c (m s -1 )	kc (m -1 )	fc (Hz)
	20.5	3.18 ± 0.04	0.173 ± 0.005	80 ± 5	4.35 ± 0.1
	50.1	5.92 ± 0.04	0.303 ± 0.005	100 ± 6	4.60 ± 0.1

ACKNOWLEDGMENTS

We are grateful to F. Charru, J. Magnaudet, R. Mathis, F. Burdairon, M. Aulnette and S. Perrard for fruitful discussions, and to P. Balondrade and J. Sant-Anna for assistance with measurement and data processing. We thank A. Aubertin, L. Auffray, J. Amarni and R. Pidoux for experimental help. This work was supported by the project "ViscousWindWaves" (ANR-18-CE30-0003) of the French National Research Agency.


Sebastien Kolb
email: sebastien.kolb@ecole-air.fr

Flight dynamics analysis and trajectory optimization of some gliding phases

As far as the gliding phases are concerned, the specificities of the aircraft flight dynamics and the surrounding aerology must be taken into account seriously, especially since no propulsion is at disposal to cope with possible problems. In this study, several issues are examined. On the one side, pilot induced oscillations of the space shuttle command channel are expertized and on the other side, a performance (range, mean velocity, altitude loss) is maximized for a sailplane in an environment with a given repartition of thermals (or wind). Even if the study is more applied mathematics oriented, the results are interesting for engineers and pilots since it helps to take in-flight decisions and to make pre-flight preparations.

Finally for such studies, one hope is surely to help pilots to make wise in-flight decisions and to prepare as well as possible the gliding phases taking into consideration the amount of data and their level of precision at disposal but also to compare different scenarios and to avoid unscheduled off-field landing and more generally to diminish soaring risks. The mathematical analysis, optimization and underlying physics are explained in theoretical and practical terms in order to be valuable for mathematicians and aerospace specialists.

The analysis of some issues coming from nonlinear flight dynamics is accomplished here by means of the bifurcation theory and continuation algorithms. Some nonlinear features are actually due to couplings with pilot reaction here and may imply unexpected events and it is then important to evaluate whether there are dangerous or not. The main topic of this part concerns pilot-aircraft couplings for a configuration with forced oscillations as inputs. [START_REF] Duc | Effect of actuator saturation on pilot-induced oscillation: A nonlinear bifurcation analysis[END_REF] The application concerns the space shuttle which meets some issues in its gliding phase. [START_REF] Smith | Analysis of a longitudinal pilot-induced oscillation experienced on the approach and landing test of the space shuttle[END_REF] The analysis of the command channel with a rate-limited actuator reveals the presence of jumps (flying qualities cliffs), indeed a little variation of the input pulsation may imply large variation (of amplitude and phase). According to the describing function method (DFM), a continuation algorithm allows to predict the amplitude (and the phase shift) at the entry of the actuator thanks to the resolution of an harmonic balance equation (HBE) which involves the transfer functions of the different parts and the describing function of the nonlinearity (saturation or rate limiting). Besides another topic is the flight trajectory optimization of a sailplane. When (the mean velocity is) calculated in a stationary environment with only one vertical thermal, it is associated to the so-called MacCready speed, whereas configurations with several thermals raise some questions such as the optimal speeds and the best strategy to adopt (which speed to take or which thermals to target, in which order, etc). Another strategy might be to realize a dolphin-style gliding and keep on flying (in thermals). In this case, optimal control can furnish the best speeds (and altitudes) to follow during a longitudinal flight (in the vertical plane) depending the magnitude distribution of the vertical wind.

Introduction

Gliding phases are specific in the sense that it is not possible to use thrust. It implies a necessity of carefull handling of the other controls and an obligation to take seriously into account the aerology (and aerodynamics). In this study, two themes are dealt with, the types of aircraft are not the same and the mathematical tools employed neither but all of them present interesting features.

Concerning the pilot-induced oscillations of the space shuttle during its landing, the trigger is a rate-limited actuator. As a consequence, the classical linear criteria are not sufficient and efficient to assess the susceptibility to these aircraft-pilot couplings successfully. Several methods are at disposal to analyze these phenomena and to attenuate or suppress it. Indeed major accidents of fighters (Gripen, YF22, YF16, etc) occurred for example during the nineties and have forced the aerospace industry and research community to develop methods in order to predict and control the flight dynamics during dangerous phases which demand an important pilot effort like landing, take-off, refueling, etc, thus exacerbating rate limiting effects. Unfortunately the linear theory which was widely employed in the engineering process was not sufficient, therefore some tools and concepts coming from the research world needed to be taken into account and transferred to the industry. Sometimes these last ones might be too conservative in order to be employed successfully for a design process, nevertheless they always present a clear interest in understanding the phenomenon and trying to master it. Describing function method, harmonic balance and bifurcation theory are powerful methodologies, since they allow to diagnose existence, stability of periodical orbits and their changes of characteristics or multiplicity. They are the mathematical core of this study and are extensively used. Nevertheless other methodologies are also relevant and some new ideas are trying to take better into account the transitory characteristics of nonlinear flight dynamics.

As far as the optimization of sailplane flight is concerned, it is always interesting to find a compromise between mathematics and concrete piloting. Indeed on the one side, optimization can furnish precise results but on the other side the pilot must be able to apply them successfully in concrete situations. Especially the pilot does not know perfectly the aerological environment quite often and the in-flight instruments have an inherent delay before displaying the correct measure. Consequently, in order to optimize its flight, the pilot must have an idea of the diverse possible strategies and some knowledge about their mathematical and physical foundations, but also a methodology in order to take quick decisions which may be partly intuitive. For example, when estimating the ascending speed in a thermal, the pilot experience may play a role but also when choosing a strategy amongst several ones. A focus is put in the use of optimal control, since it permits to consider a large variety of flight types, performances to optimize and aerologic situations. Some care must nevertheless be taken in the mathematical formulation in order to facilitate the calculation process.

Even if both themes are different, they each present huge interest and some common mathematical methodologies and flight experiences can contribute to improve the pilot ability and the engineering design process.

Pilot induced oscillations of the space shuttle

The space shuttle met pilot induced oscillations during its approach and landing tests. Nasa technical reports furnishes well documented details on the command channel as well as possible improvements in its technical memoranda. [START_REF] Smith | Analysis of a longitudinal pilot-induced oscillation experienced on the approach and landing test of the space shuttle[END_REF][START_REF] Smith | Design of a nonlinear adaptive filter for suppression of shuttle pilot-induced oscillation tendencies[END_REF] In order to evaluate the susceptibility to such aircraft -pilot couplings, several methodologies exist. A direct attempt consists in applying forced oscillations with various amplitudes and pulsations.

Forced oscillations

In order to evaluate the susceptibility of the space shuttle command channel, forced oscillations may be applied as pilot stick inputs like described in figure 1. Depending on the reaction of the aircraft to these direct solicitations, it is already possible to have an idea of the potential issues and to their hazardousness. Since the source of pilot induced oscillations is a rate limited actuator, a nonlinear analysis is required (linear considerations cannot catch such behaviours). The describing function method is here employed, thus we try to find solutions to equation 1 thanks to a pseudo-arclength continuation. [START_REF] Dhooge | Matcont: A matlab package for numerical bifurcation analysis of odes[END_REF][START_REF] Ravnas | Continuation and bifurcation software in matlab[END_REF] Ae jϕ = Filter( jω)

1 -N(A, ω)Actuator( jω)Air f rame( jω)Feedback( jω)Filter( jω) Feed f orward( jω)Input (1)
The results (for the worst case) with a time delay T = 0.4 s and an input amplitude of 15 • are shown in figure 2. Some bifurcations of periodical orbits are diagnosed and are responsible for sudden jumps. From the flight dynamics viewpoint, these last ones may take the pilot by surprise and lead to an accident. Another issue is that once the branch with higher amplitudes is reached by a little change of input frequency, a huge decrease in the frequency of the inputs is necessary so as to jump back to the branch with lower amplitudes. This hysteresis phenomenon may render the situation quite difficult to tackle. In figure 3, the time simulations realized for different input pulsations (ω = 1.2 rad/s and ω = 1.3 rad/s) are converging towards the limit cycle of low amplitude or high amplitude. If this jump occurs, the pilot may be astonished. But the time simulations show that once a second limit cycle exists, the switch to the high amplitudes can happen anytime. This may create a dangerous situation and a linear analysis is far from predicting such an unpleasant occurrence.

Other methodologies may also be employed in order to analyze the behaviour. Some fresh ideas are brought by [START_REF] Duc | Effect of actuator saturation on pilot-induced oscillation: A nonlinear bifurcation analysis[END_REF] with (equivalent) nonlinear bode diagrams. Besides continuation of limit cycles is a classical attempt to use bifurcation theory. [START_REF] Dhooge | Matcont: A matlab package for numerical bifurcation analysis of odes[END_REF] Nevertheless since fold bifurcations of limit cycles are the only critical elements, employing a sinusoid with a sweeping frequency (increasing and decreasing) as input seems to be enough in order to catch the two branches of periodical orbits and to have an idea on the underlying dynamics and especially the fold bifurcation occurrences. On the time simulations of the figure 4, the amplitudes at the entry of the actuator saturation are suddenly jumping at frequencies whose values correspond approximately to the ones coming from the analysis with the describing function method. The inputs are sinusoids with a sweeping frequency from 0.2 rad/s to 2 rad/s (increasing and decreasing on 1000 s).

After accomplishing an analysis and diagnosis of pilot induced oscillations, we will try next to attenuate or suppressing this issues thanks to anti-windup techniques.

Anti-windup

When trying to cope with PIO problems due to rate limiting, the most direct way consists in adapting the sensitivity of the stick and thus to adapt the shaping function. A first proposal is made at the end of the technical memorandum [START_REF] Smith | Analysis of a longitudinal pilot-induced oscillation experienced on the approach and landing test of the space shuttle[END_REF] and another one is exposed in the report. [START_REF] Smith | Design of a nonlinear adaptive filter for suppression of shuttle pilot-induced oscillation tendencies[END_REF] Moreover, it is also possible to add a feedback block so as to re-use the PID tuning know how.

For example, when adding a PID block before the lead-lag filter, we may hope to change sufficiently the dynamics by a fine tuning so as to modify the reaction to strong solicitations in the stick. In the next configuration a PID under the form PID(s) = k p + k i s + k d s T f s+1 is selected. A comparison is made between the nominal configuration and the one with this modification for forced oscillations of amplitude 5 deg as input (and a time delay of 0.5s) and is illustrated in the figures 5. A direct anti-windup strategy consists in measuring the difference of the input and the output of the rate-limited element, than to filter the signal with a low-pass filter and to apply the little correction at the entry like described in the modified command channel 6. The effect can be direct and are illustrated in the time simulations 7. By adapting equation 1 so as to employ the describing function method and a continuation algorithm, it is possible to determine the (first harmonic) characteristics of the periodical orbits at the entry of the actuator saturation in function of the input pulsation for this specific case with a sole feedback loop. This last simple anti-windup technique has bad reputation for the high frequencies and is thus often extended with a bypass. [START_REF] Rundqwist | Rate limiters with phase compensation[END_REF] But here it seems to offer sufficient effects to tackle this issue.

The interest of these nonlinear analysis tools are clearly shown, since they permit to know in advance whether a configuration might be hazardous or not toward PIO or PIO-like. Once the methodology and the numerical tools are adapted and implemented, it is possible to compute systematically the different control systems, the limitations (numerical values or types of command channels) being known in advance by experience.

The first topic of this study dealt with the pilot induced oscillations of the space shuttle. An analysis was accomplished and strategies to cope with it were evoked. The next part of this study is now devoted to soaring in an environment with a specific aerology.

Sailplane flight in an environment with vertical wind

Several strategies are explored in order to obtain the best performances for a sailplane facing an environment with some thermals during a longitudinal flight. The classical method is first described for a unique thermal before trying to optimize the flight in more general wind configurations.

MacCready speed-to-fly

When competing in the presence of an (ascending) thermal, it is possible to choose an optimal airspeed, the so-called MacCready speed as exposed in textbooks [START_REF] Gudmundsson | General Aviation Aircraft Design[END_REF][START_REF] Reichmann | Cross Country Soaring[END_REF] on sailplane. It aims at achieving cross-country flight as quickly as possible. During the first phase of the flight in calm air, the glider flies with a faster airspeed and may thus loose some excessive altitude but in the thermal, it is possible to return to the initial altitude in a reasonable amount of time. This permits to develop an optimal strategy in the context of a competition.

Graphically the speed-to-fly can be read thanks to a tangent line to the speed polar curve crossing the point 0, V Z 1 where V Z 1 is the estimated vertical speed of the glider in the thermal.

As mathematical foundation to this theory with simplified hypotheses (no time lost in stabilizing the glider in the thermal for example), the demonstration leads to the equality of the ratio

V mean V a = V Z 1 V Z 1 -V Z (V a ) ( 2
)
where V mean = D T total is the mean speed required to perform the overall flight that is to say the distance (to the thermal) divided by the total amount of time T total which englobes a gliding phase in still air and an ascending phase in the thermal. Browsing all the points of the speed polar V Z (V a ) allows the optimisation of V mean and leads to the best airspeed V a = V speed-to-f ly . The figure 9 plots a speed polar calculated (with data taken in flight dynamics literature 1 ) and its tangent line passing through 0, V Z 1 . The ascending speed in the thermal is here equal to V Z 1 = 2m/s and the speed-to-fly is calculated to be V speed-to-f ly = 33.7m/s.

After presenting and employing a standard way to tackle soaring with thermals, a more mathematical approach is employed and uses optimal control theory. This will be the next evoked subject.

Optimal Control

Soaring in a environment with a certain wind distribution can be treated thanks to specific mathematical tools. The model and some results exploiting optimal control theory can be found in [START_REF] Bion | Minimum altitude-loss soaring in a sinusoidal vertical wind distribution[END_REF][START_REF] Bion | Minimum-time soaring through a specified vertical wind distribution[END_REF] and are valuable since writing the equations and formulating the optimization problem correctly is important in order to be able to employ successfully such tools. For optimal control calculations, the code ICLOCS2 5 is here employed. Aircraft data are here the ones of the Slingsby Dart T51 (taken in 1 ). But identification of the speed polar in the flight manual may be possible for other gliders. A polynomial model of reduced order (2 for the drag polar and 3 for the speed polar) may be sufficient even if at the highest speed, the error may sometimes begin to be significant.

The equations 6 describe the aircraft motion in the aerodynamic frame R b and include the vertical wind distribution W h .

dh dt = V a sin γ + W h (3)
dV a dt = - ρV 2 a C X S 2m -V a cos γ dW h dX dX dt + g sin γ (4)
dγ dt = ρV a C Z S 2m -cos γ dW h dX dX dt + g V a cos γ (5) (6)
Concretely, re-writing all the time derivatives with distance derivatives seems to simplify greatly the computation process. That's why we will use the following equations for the numerical resolution, exploiting the fact that dx dt = V cos γ.

dh dx = tan γ + W h (x) V a cos γ ( 7
)
dV a dx = - ρV 2 a C X S 2m -g sin γ / (V a cos γ) - dW h (x) dx sin γ (8)
dγ dx = ρV a C Z S 2m - g V a cos γ / (V a cos γ) -cos γ dW h (x) dx /V a (9)
When the objective consists in flying as far as possible, the performance to optimize (altitude loss in fact) can be re-written

J = t f 0 (W h + V sin γ) dt = x f 0 W h + V sin γ V cos γ dx ( 10
)
A calculation is performed in the case of two vertical thermals of different ascending airspeeds. Here the ascending speeds of the two thermals are equal to 2m/s and 3m/s (respectively between 20 km and 25 km and between 40 km and 45 km). The speed to adopt by the glider is plotted in figure 10 and also the corresponding altitude during the flight. As constraint, the lift aerodynamic coefficient cannot exceed its maximum value C Z max .

Figure 10: Best speeds to adopt in order to minimize the altitude loss of a glider in a calm environment with two stationary ascending thermals. Results obtained thanks to optimal control of dolphin soaring.

Another performance to optimize may be the time needed so as to cover a certain range. The mathematical objective is here the following one.

J = t f 0 dt = x f 0 1 V cos γ dx (11)
As constraints, the lift aerodynamic coefficient cannot exceed its maximum value C Z max and the certified maximum airspeed implies the existence of a minimal value C Z min . With the same aerology as before, calculations based on optimal control are accomplished for two ranges of 50 km (with a final altitude of 5000 m) and 100 km (with a final altitude of 1000 m). According to figures 11, the velocities to adopt proves to be different.

Figure 11: Best speeds to adopt in order to minimize the time needed to cover a certain range for a glider in a calm environment with two stationary ascending thermals. Results obtained thanks to optimal control of dolphin soaring for two ranges of 50 km (bottom) and 100 km (top).

After addressing the classical issue of minimizing the time and the altitude loss when flying in an aerology containing some fixed thermals, it is also possible to examine other situations which may also happen.

A more general vertical wind profile

The distribution of the wind profile for which the gliding is optimized can be more general. The practical knowledge of the aerology comes most of times from an estimation of the ascending speed in the forecoming thermals by observing the other gliders or by collecting information on them. It can also be the experience of the pilot who has an idea of the properties of a thermal under a cloud for example. Here another situation is examined. The wind profile may be estimated by other means (like computation). The mathematical frame remains described by equations 6 and 10 but the wind profile W h is more general (it depends just on the position).

For example, in the case of a sinusoidal wind profile of periodic length 5000 m and of amplitude 2 m/s, the optimal trajectory is computed. The figures 12 present the airspeed to follow and the altitude in function of the distance covered.

Figure 12: Best speeds to adopt in order to minimize the altitude loss of a glider in an environment with sinusoidal ascending wind. Results obtained thanks to optimal control of dolphin soaring.

The indications and concrete optimal values are interesting, but such an optimal trajectory is maybe more difficult to follow as far as sailplane piloting is concerned.

Conclusion

In this study of gliding phases, two main topics were dealt with. An analysis of the pilot induced oscillations of the space shuttle due to rate limiting was first accomplished thanks to the describing function method and the associated curves were computed by means of a continuation algorithm. It revealed to be a pertinent method to diagnose the hazardousness of a configuration before and after potential modifications of the command channel. The practical experience consisted in applying direct oscillatory inputs and observing (or not) some jumps corresponding to handling qualities cliffs. Several anti-windup strategies were tested in order to cope with this issue. The other topic focused on the trajectory optimization of a sailplane with a distribution of vertical wind. After describing the mathematical model employed, some practical configurations, situations were considered with isolated thermals or with a periodic wind distribution. The loss of altitude was optimized but other performances such as the time needed could be sought as well.

In both parts, mathematical treatments allowed to determine valuable results and pieces of information which can help pilots to know how to react adequately or in an optimal manner. In the design process, engineers can also apply these numerical tools to fix in advance potential issues or to estimate optimal solutions.

Figure 1 :

 1 Figure 1: Command channel of the space shuttle.

Figure 2 :

 2 Figure 2: Amplitude A and phase ϕ at actuator entry in function of the pulsation ω of the pilot inputs according to DFM. Diagnosis of jumps due to a rate limited actuator in the space shuttle command channel.

Figure 3 :

 3 Figure 3: Time simulations at actuator entry for pilot inputs of amplitude 15 • and of pulsation ω = 1.2 rad/s (left) and ω = 1.3 rad/s (right). Stabilization to the branch of low amplitude or to the one of high amplitude.

Figure 4 :

 4 Figure 4: Time simulations at actuator entry with an increasing (left) and decreasing (right) pulsation ω of the pilot inputs. Diagnosis of jumps in the space shuttle command channel response.

Figure 5 :

 5 Figure 5: Influence of a PID in a command channel without (left) or with it (right). Presence or absence of jumps for the space shuttle command channel with forced inputs as shown by the diagrams presenting amplitude A (of the periodical orbit) in function of pulsation ω.

Figure 6 :

 6 Figure 6: Low-pass filter correction applied to a command channel and taking into account the diffence between the input and the output of the rate-limited actuator.

Figure 7 :

 7 Figure 7: Influence of a low-pass filter correction. Time simulations for an input of amplitude 15 • and pulsation ω = 2 rad/s. Disappearance of jumps with the feedback correction (left) in comparison with the original command channel (right).

Figure 8 :

 8 Figure 8: Diagrams presenting amplitude A and phase ϕ (towards forced oscillations inputs) of the periodical orbit at actuator entry in function of pulsation ω for the pitch attitude control system with a low-pass filter correction according to DFM.

Figure 9 :

 9 Figure 9: Speed polar of the Slingsby Dart T51 (calculated) and tangent line for the determination of the MacCready speed-to-fly.

  YY pitching moment of inertia C X aerodynamic coefficient of drag C Z aerodynamic coefficient of lift Cm aerodynamic coefficient of pitching moment R X drag force R Z lift force R a aerodynamic frame R b body-fixed frame R 0 Earth frame S(main wing) reference area V a airspeed W h vertical wind distribution

		Nomenclature
	α	angle-of-attack (AoA)
	δ e	elevator angle (pitch control)
	γ	flight-path angle
	ρ	air density
	θ	pitch angle
	c	chord
	g	Earth's gravitational acceleration
	h	altitude
	m	mass
	q	pitch rate
	I

DOI: 10.13009/EUCASS2023-864