points per element (left) versus 512 (right) showing good matching between the vorticity fields. . . . . . . . . . . . .

3.7 Evolution of the vorticity of the symmetric optimal perturbations for horizon time T = 2 (top) and T = 5 (bottom) at times T /4, T /2, 3T /4 and T . As the baseflow boundary layer detaches, the optimal perturbation takes the form of two displacement modes: one acting on the primary vortex, the other on the secondary. Solid lines correspond to the iso-contours of vorticity (ω z /ω 0,max = ±0.1, ±0.2) of the primary vortices. . . . . . . . . . . . . . . . . . . . . . . . . . 62

3.8 Evolution of the gain in kinetic energy with the wavenumber for various horizon times T = 2, 3, 4, 5 for the two symmetries. Similar overall behaviour is observed for both symmetries where two wavenumbers seem to arise, k = 0.8 and 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.9 Axial vorticity for the symmetric optimal pertubations, k = 0.8 (left) and k = 6 (right), SYM perturbations, horizon time T = 5. Stronger ground effects are observed for k = 6: the iso contour of vorticity |ω z | 5 are observable in the vicinity of the ground plane whereas they are not for longer wavelengths. . . . 64 3.10 Time sequence of axial vorticity of symmetric optimal perturbation at k = 0.8 (top) and k = 6 (bottom), SYM perturbations, horizon time T = 5. The optimal perturbation initially acts as a displacement mode inside the primary vortex cores and when boundary layer separates, the perturbation drastically grows in the boundary layer. When the boundary layer has detached, displacement modes grow into the secondary vortices giving rise to structures similar to those that were observed in two dimensions. Vorticity magnitude ranges between 10 and 900 (top). and between 10 and 3 × 10 4 (bottom). . . . . . . . . 65

3.11 Time sequence of vorticity magnitude of the unperturbed baseflow. The flow remains two-dimensional before perturbations set in after t = 20. . . . . . . . . 66 3.12 Time sequence of vorticity magnitude of the nonlinear response to the symmetric long-wave optimal perturbation computed for T = 5 and k = 0.8. The early breakdown of structures into short waves can not be captured by the mesh resolution in the axial direction which is made for the long-wave dynamics. . . 67 3.13 Time sequence of vorticity magnitude of the nonlinear response to the symmetric short wave optimal perturbation computed for T = 5 and k = 6 with an initial amplitude of ε = 10 -3 . The three levels of iso-contours are ω/ω 0,max = 0.05, 0.1 and 0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.14 Time sequence of vorticity magnitude of the nonlinear response to the symmetric short wave optimal perturbation computed for T = 5 and k = 6 with an initial amplitude of ε = 10 -2 . The three levels of iso-contours are ω/ω 0,max = 0.05, 0.1 and 0. for T = 1 and p = 1.6. On the left, the optimal perturbation remains m = 2, on the right, the optimal perturbation is a sum of various modes. . . . . . . . . 85 4.9 p-norm objective function as a function of p. The dotted line represent "spiral" solution and the solid line represent "whip" solutions. These two branches appear after amplification of noise during the initialilization of the adjoint. . . 86 4.10 Flow past a finite span wing (purple) at Re = 1000 with AR = 6 and α = 5 optimal perturbations for a flat plate wing. Computed at T = 10, AR = 6 and α = 5. As expected the p-norm procedure localizes the perturbation. The contour represent the streamwise vorticity levels at ±0.2 of the perturbations.

The shape of the optimal perturbation do not vary when p ≥ 2, this saturation occurs when the diffusion and p-norm lengthscales are equivalent. . . . . . . . 90 4.15 Perturbation streamwise vorticity at horizon time T = 10 (red and blue). The green chevrons represent the streamwise vorticity of the initial optimal perturbation. The chevrons invert while travelling downstream. . . . . . . . . . . . . 91

Introduction

More than 500 years ago, the polymath Leonardo Da Vinci was named a "Master of Water" by the Florentine authorities. In this role he had to explore deviating the river Arno away from Pisa. In his notebook, he depicted how fluids flow in a diagram, problably without knowing that he was describing the three-dimensional nature of flowing water and the fact that turbulent flows consist of a range of co-existing eddies, varying in scale from large to small. These scales were mathematically formalized by the mathematician A. N. Kolmogorov, and compose the well known Kolmogorov cascade of energy. Energy of bigger eddies gets con- verted to that of smaller eddies. This process goes on until the eddy is small enough so that the energy dissipates directly to heat in the fluid. Or, as this was all phrased in a poem by the English mathematician Lewis Fry Richardson:

INTRODUCTION "Big whirls have little whirls, Which feed on their velocity, And little whirls have lesser whirls And so on to viscosity."

Vortices, "the sinews and muscles of fluid" as described by [START_REF] Küchemann | Report on the iutam symposium on concentrated vortex motions in fluids[END_REF], constitute the elementary structures of all fluid flows more particularly in transitional and turbulent flows.

Studying their intrinsic dynamics and interactions helps us understand the behavior of the complicated flows. Vortices are inevitable and are constantly encountered in industrial flows. During the combustion of a fuel, they are beneficial as they increase efficiency by improving the mixing the fuel with the oxydiser. Contrariwise, the vortices forming at the tip of rotating blades in axial compressors (also known as tip-leakage vortices TLV), can cause significant power loss and are linked to the stall of a compressor. In this case, vortex removal would be beneficial. In aviation, the wingtip vortices from an airplane dissipate slowly and linger in the atmosphere long after the airplane has passed, becoming a potential flight hazard.

Wake vortices : the stakes

This Ph.D. program is part of the PHYWAKE convention, between ONERA and the French Civil Aviation Authority. The aim of the convention is to improve theoretical knowledge on wake vortices, evaluate the environmental impact of condensation trails, optimize safety distances and improve detection and measurements of wakes vortices. Aircrafts remain airborn due to the presence of lift. The latter results by the difference of pressure between the lower side of the wing (intrados) and its upper side (extrados). This pressure gradient results in the rotation of the fluid about the wing tip as the fluid is pushed from below. Once the fluid rollup is complete, two counter rotating vortices appear downstream of the aircraft, as sketched in figure 3. Measuring and characterizing wake vortices features would allow a better traffic management in airports.

Flight Safety : Airport Proximity

For the past fifty years, air traffic has been in almost progress with a yearly average of 5%. This will continue as part of the democratization of air transport in many parts of the world (see, Airbus SAS (2015)). Among other issues the increase in the number of flights has led to the saturation of major airports, especially at peak hours. In that context separation distances imposed to mitigate the risk of wake vortices encountered at take-off and landing appear as an important factor to play with in order to increase frequency of arrivals/departures. The goal is to improve airport capacity without building new infrastructures, which is often impossible for technical, environmental or societal reasons. Wake vortices are robust structures and contain a lot of energy in the circular motion of the fluid. The vortex strength is characterized by its circulation, Γ. From the Kutta-Joukowski theorem, one shows that the circulation of the vortices is proportinal to the lift (hence the weight of the aircraft). This causes safety problem in airports. Indeed, an aircraft encountering wake vortices from a preceding aircraft may lead to lethal consequences. The Federal Aviation Administration and the International Civil Aviation Organisation established separation distances (see table 1) between the take-off and landing of aircraft of each category (small, medium, heavy,jumbo). This would ensure that vortices would dissipate and no longer pose safety issues. This is particularly true for wake vortices impacting the ground as the vortices persist above the runways, hence reducing the frequency of arrivals departures in airports. The need for increasing airport capacity and flight safety becomes therefore of obvious interest.

Objectives of the thesis

The thesis focuses on the linear and nonlinear dynamics of counter rotating vortices impinging a ground plane. This setup is relevant for aircrafts in the vicinity of runways. Considering the previously described safety issues and measurement goals, the following question arises: "How can one enhance vortex decay in ground proximity?".

To answer the first question, conceptual solutions based on flow control are envisaged.

One efficient way would be to promote the inviscid dynamics of the vortex pair. Indeed the presence of viscosity, due to the induced boundary layer and entrainment effect associated with secondary vorticity, results in the persistence of the vortices below the flight path. This strongly constrains the possibility for reduced separation distances between aircrafts. This naturally suggests modifying the boundary conditions applied at the ground. This is achieved using the optimal control method in order to reduce the effect of viscous friction at the wall. These results and others will be detailed.

The second approach consists in determining the optimal perturbation, that is the perturbation that maximizes the kinetic energy of the vortex system. The objective of this thesis is to determine the linear optimal perturbation of the aircraft vortex wake. Although this goal is yet to be achieved in its entirety, several successful steps are taken towards accomplishing this aim. Due to the complex spatial distribution of the optimal perturbation, applying the results obtained with linear optimal perturbaions remains challenging. To tackle that very problem, we will explore the use of so-called p-norms. Conventional optimal perturbations are obtained using the classical L 2 norm of the velocity field. Higher order norms, on the other hand use the L p norm in order to localize the perturbation in a confined region of space. This method was first used by [START_REF] Foures | Localization of flow structures using ∞-norm optimization[END_REF] and is worthwhile for industrial applications.

Organisation of the thesis

As the work described in this thesis is exclusively numerical, chapter [1] starts by describing the mathematical and numerical tools that are used throughout the manuscript. We start with a brief overview of the Spectral Element Method (SEM). We follow with a presentation of the various variational formulations of the constrained optimization problems in order to introduce the concepts of optimal control and optimal perturbations.

In chapter [2], we present the non-linear dynamics of a pair of counter-rotating vortices in ground effect. Following a review of the literature, we present an approach aiming at increasing the lateral position of the vortices using an optimal control stategy acting on the ground. The contents form an article published in the Journal of Fluid Mechanics, see [START_REF] Wakim | A vortex pair in ground effect, dynamics and optimal control[END_REF]. Chapter [3] addresses the optimal perturbation problem in the case of a two and threedimensional counter-rotating vortices in the vicinity of the ground. The goal is to reduce the lifespan of the vortices by adding optimally located perturbation in the fluid domain. Lastly, in chapter [4], we focus on the localization of the optimal perturbation by changing the classical L 2 -norm approach with higher order norms. We first focus on the case of an isolated Lamb-Oseen vortex and then investigate the three-dimensional p-norm optimal perturbation of the flow past a flat plate wing following the study of [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF].

INTRODUCTION

6

CHAPTER 1

Numerical Tools

The work described in this thesis is exclusively numerical. In this chapter, we describe the tools and methods that are used.

Spectral Element Method tool

Most of the Direct Numerical Simulations (DNS) were carried out using the incompressible Navier-Stokes (N-S) open source solver Nek5000 (Fischer et al. (2008)) well renowned for its parallel computation performance and the availability of its adjoint solver, among other features.

In the following sections, the flows are governed by the incompressible N-S equations :

∇ • u = 0 ∂ t u + (u • ∇)u = -∇p + 1 Re ∇ 2 u (1.1)
which are solved using the Spectral Element Method (SEM) implemented in Nek5000.

Space Discretization

The SEM is an approximation scheme based on the Galerkin method and was introduced by [START_REF] Patera | A spectral element method for fluid dynamics: laminar flow in a channel expansion[END_REF]. It combines the geometrical flexibility of the Finite Element Method with CHAPTER 1. NUMERICAL TOOLS the convergence properties of Spectral Methods. More details are provided in [START_REF] Deville | High-order methods for incompressible fluid flow[END_REF]; [START_REF] Abgrall | High-order methods for cfd[END_REF]. In this subsection, we briefly summarize the theoretical content on spatial discretisation.

The Galerkin approximation allows the solving of Partial Differential Equations (PDE) using their variational (weak) form. Let's consider the following PDE :

Du = f (1.2)
with D being the differential operator. The weak formulation of the previous equation is expressed as follows :

A(u, v) = F (v) ∀v ∈ V (1.3)
where A is a bilinear operator arising from the integration by part of equation 1.2 and F a linear function.

To approximate the exact solution û, we consider a finite-dimensional subspace V h ∈ V , a set of basis functions φ n and numerical coefficient u n such that :

u N (x) = N n=0 u n φ n (x) A(u N , v) = F (v) ∀v ∈ V h (1.4)
Nek5000 uses the Legendre Polynomials as basis functions. In addition to orthogonality, they provide the best approximation using the H 1 -norm. For a one-dimensional equation, the finite-length domain [a, b]

is divided into E elements. such that a = x 0 < x 1 < . . . < x E = b. Each element x i-1 < x < x i is mapped onto the interval [-1, 1]
. Pieced together over all E elements, the Legendre Polynomials of degree p form the basis φ n . The points {ξ 0 , ξ 1 , . . . , ξ p } which form the interpolation grid (within each element) are the Gauss-Lobatto-Legendre quadrature points. They are the solution of :

1 -ξ 2 i dL N dξ (ξ i ) = 0 (1.5)
where L N is the Legendre polynomial of degree N. The spectral approximation u N,e of degree N is given by :

1.1. SPECTRAL ELEMENT METHOD TOOL u N,e (ξ) = N i=0 u i,e π i (ξ), π i (ξ) = -1 N (N + 1) 1 -ξ 2 L N (ξ) (ξ -ξ i ) L N (ξ i ) 0 ≤ i ≤ N (1.6)
Spatial discretisation is therefore achieved, first, by splitting the integration domain is divided into E elements, then by increasing the order of the polynomials.

Refining a mesh can therefore be performed in two ways :

• increasing the number of elements E.

• increasing the polynomial order P .

Time Discretization

Nek5000 uses a semi-implicit BDFk/EXTk time scheme to discretize the N-S equations. Diffusion is treated implicitly. The time derivative is approximated using the k-order Backward Differentiation Formula. The whole scheme is globally second order accurate in time. For the remainder of the manuscript, we use the BDF2/EXT2 scheme. The resulting discretised N-S equations are :

∇ • u n = 0 3u n -4u n-1 + u n-2 2∆t = -∇p n + 1 Re ∇ 2 u n -( 2(u • ∇)u| t n-1 -(u • ∇)u| t n-2 ) (1.7) CHAPTER 1. NUMERICAL TOOLS
Finally the constant timestep ∆t is chosen such that the Courant-Friedrichs-Lewy number remains below 0.6 (see, Fischer et al. (2008)).

Boundary Conditions

Depending on the nature of the flow that is simulated, the boundary conditions of the calculation vary. In the case of the bi-dimensional vortex pair in ground effect, no-slip Boundary Condition (BC) is applied on the ground and symmetric BC are imposed on the other boundaries. The top and right boundaries are placed sufficiently far from the main flow to ensure that their influence is negligible. A sketch is provided in figure (2.1). In the case of a threedimensional vortex pair in ground effect, the previously described bi-dimensional mesh is extruded in the axial direction and periodic conditions are applied on the axial boundaries.

Tool Validation

The validation of this numerical setup is provided in figure (2.2) where we compare vortex centroid trajectories obtained by Türk et al. (1999) and Zheng & Ash (1996) for three distinct values of the Reynolds number and two different initial conditions.

Variational Formulation and Constrained Optimization

In this thesis, the vortex problem is assessed in two ways. The first method, presented in Chapter 2 relies on the Optimal Control Method, which consists in finding the optimal timedependant boundary condition to maximize a flow related cost function. The second method relies on the optimal perturbation concept, detailed in Chapter 3. Its goal is to determine a perturbation field which maximizes a perturbation related cost function.

Both of these methods can be described within the same mathematical framework, namely, the constrained optimization framework.

Mathematical Framework

Let us define the following variables :

• The system state is described by a velocity field and a pressure field (etc.) for every position and time that we define collectively as X.

• The control variables that we define are denoted collectively as U .

• The cost function that needs to be optimized is denoted as J (X, U ).

Usually it includes a penalization of the control as follows

J (X, U ) = E(X) + l||U || 2 ,
where || • || is a norm and l is a positive parameter. This usually allows to be a well posed problem in which the control amplitude is bounded.

• The constraints which represent the governing equations are denoted by F(X, U ) = 0.

These, in the present case, are the N-S equations.

The general idea of control is to design a cost functional J expressing the control objectives such as energy, enstrophy, drag reduction, shape optimization (etc.) and the seek for the extremum of this cost function. Depending on the chosen method, the goal is either to find the optimal boundary condition at the wall which maximizes the chosen cost function, or to find the initial perturbation field which allows the maximum energy growth over a given time interval [0, T ].

From Constrained to Unconstrained Optimization

The optimization of J is a constrained optimization problem :

Determine the optimal control and state such that J is maximized (or minimized) under the constraint F(X, U ) = 0.

A natural approach is to compute the "gradient" of the cost function with respect to the control in which we link X and U under the constraint F(X, U ) = 0. In other terms, we want to compute dJ dU (X, U ) where X depends on U . The optimal control would then satisfy dJ dU (X(U ), U ) = 0. An other (very elegant) way is to introduce the so-called Lagrange Multipliers. Let us introduce the following inner product :

f , g (1.8)
whose definition depends on the context and the following Lagrangian function :

L(X, U, X) := J (X, U ) -F(X, U ), X (1.9)
where X which will later be assossiated to the adjoint (or co-state) of X. As opposed to the initial problem, the variables are assumed to be independant from each other. The initial problem transforms into the following :

"Determine the optimal control, state and adjoint variable such that L is stationary."

CHAPTER 1. NUMERICAL TOOLS

The following question arises : How does one find the stationary points of the Lagrangian function? Adjoint-based optimization techniques are commonly designed with the intention of finding these points using efficient iterative algorithms. As optimization theory is mainly based on the concept of directional derivatives, let us have a reminder:

A function f is Gateaux differentiable if f has a directional derivative along all directions at x. This means that there exists a function g such that

g(h) = lim t→0 f (x + th) -f (x) t (1.10)
g is also known as the differential of f at x. We then define the weak form of the gradient as the unique vector ∇f such that g(h) = ∇f (x), h .

The stationary points are found by cancelling the first variation of L with respect to X, X and U . ( The prefix δ denotes arbitrary test functions mostly used in inner products. All adjoint functions are superscripted by a tilde " • ". ). The weak form of the gradient of L is defined, for any test function δ X, δX, U , as follows :

∂L ∂ X , δ X := lim ε→0 L(X, U, X + εδ X) -L(X, U, X) ε ∂L ∂X , δX := lim ε→0 L(X + εδX, U, X) -L(X, U, X) ε ∂L ∂U , δU := lim ε→0 L(X, U + εδU, X) -L(X, U, X) ε (1.11)
• First variation with respect to X :

∂L ∂ X , δ X = lim ε→0 J (X, U ) -J (X, U ) ε - F(X, U ), X + εδ X -F(X, U ), X ε ∂L ∂ X , δ X = F(X, U ), X = 0, ∀ δ X (1.12)
Therefore, it yields :

∂L ∂ X = F(X, U ) = 0 (1.13)
and we retrieve the governing equation F(X, U ) = 0.

• First variation with respect to X :

1.2. VARIATIONAL FORMULATION AND CONSTRAINED OPTIMIZATION ∂L ∂X , δX = lim ε→0 J (X + εδX, U ) -J (X, U ) ε -lim ε→0 F(X + εδX, U ), X -F(X, U ), X ε ∂L ∂X , δX = ∂J ∂X δX - ∂F ∂X δX, X (1.14)
Let us define the adjoint operator of L satisfying :

L(Y ), X = Y, L( X) + boundary terms (1.15)
Note that the boundary terms appear from successive integration by part. If X is the adjoint solutions then the boundary terms are nill, therefore

L(Y ), X = Y, L X (1.16)
If we define the adjoint operator of ∂F ∂X , we obtain:

∂F ∂X δX, X = δX, ∂F ∂X X (1.17) Therefore: ∂L ∂X , δX = ∂J ∂X , δX - ∂F ∂X X, δX (1.18)
Yielding the adjoint equation:

∂F ∂X X = ∂J ∂X (1.19)
• First variation with respect to U :

The same reasoning as above is applied, yielding the optimality condition

∂F ∂U X = ∂J ∂U (1.20)
Equations 1.19 and 1.20 form a PDE system whose solution yields the optimal control U , the optimal state variable X and the optimal adjoint field X. Note that, the adjoint equation 1.19 is linear with respect to the adjoint variable X.

Numerical Resolution

If we assume that one can solve the system of partial derivative equations formed by equations 1.13,1.19 and 1.20, then no iteration is needed.

However in various cases it is not possible to solve this system, therefore requiring an iterative procedure (see figure 1.2) described below:

Step 1 : Start with an initial guess value U

Step 2 : Solve the governing equations including the control

Step 3 : Compute the optimality condition

Step 4 : Solve the adjoint equation

Step 5 : Compute the gradient of the Lagrangian with respect to U

Step 6 : Update the control using the gradient of the Lagrangian

U

Guess value

Governing Equation This procedure based on the determination of the adjoint variable X is a particular case of optimization method where the gradient cost function is evaluated using the adjoint state (see equation 1.20) .

Adjoint

Adjoint Navier-Stokes Equations & Nek5000

As an example, we derive the adjoint N-S equations when we try to find the initial condition that maximizes the gain in kinetic energy :

E(T ) E(0) = 1 2 Ω u 2 (x, T )dΩ 1 2 Ω u 0 2 (x)dΩ (1.21)
of the flow over the time interval [0, T ] in a fluid domain denoted Ω. We recall the Navier-Stokes equations :

∇ • u = 0 ∂ t u + (u • ∇)u+∇p - 1 Re ∇ 2 u = 0 1.2.

VARIATIONAL FORMULATION AND CONSTRAINED OPTIMIZATION

or in a more compact way F(q) = 0, where q = (u, v, w, p). The initial conditions writes :

u(x, t = 0) -u 0 (x) = 0
or simply G(u, u 0 ) = 0. Then, we define the following inner products :

a, b = Ω T 0 a(x, t) • b(x, t) dT dΩ (c, d) = Ω c(x, t) • d(x, t) dΩ [f , g] = T 0 f (t) • g(t) dT
The Lagrangian function is expressed as follows :

L (q, u 0 , q, ũ0 ) = E(T ) E(0) -F (q) , q -( G(u, u 0 ) , ũ0 ) (1.22)
As previously shown, the first variation of the Lagrangian function yields the constraints, we will directly go through the derivation of the first variation of the Lagrangian with respect to the direct variables.

L (u, v, w, p, u 0 , v 0 , w 0 , ũ, ṽ, w, p, ũ0 , ṽ0 , w0 ) = Ω u 2 (T ) + v 2 (T ) + w 2 (T ) dΩ Ω u 2 0 + v 2 0 + w 2 0 dΩ - Ω T 0 p (∂ x u + ∂ y v + ∂ z w) dtdΩ - Ω T 0 ũ ∂ t u + u∂ x u + v∂ y u + w∂ z u + ∂ x p - 1 Re ∆u dtdΩ - Ω T 0 ṽ ∂ t v + u∂ x v + v∂ y v + w∂ z v + ∂ y p - 1 Re ∆v dtdΩ - Ω T 0 w ∂ t w + u∂ x w + v∂ y w + w∂ z w + ∂ z p - 1 Re ∆w dtdΩ - Ω ũ0 (u(0) -u 0 ) dΩ - Ω ṽ0 (v(0) -v 0 ) dΩ - Ω w0 (w(0) -w 0 ) dΩ (1.23)
L is derived with respect to each of the independant variables. After applying successive integration by parts, one can determine the following set of adjoint equations:

-(∂ x ũ + ∂ y ṽ + ∂ z w) = 0 -∂ t ũ -u∂ x ũ -v∂ y ũ -w∂ z ũ + u∂ x u + ṽ∂ x v + w∂ x w + ∂ x p - 1 Re ∆ũ = 0
-∂ t ṽ -u∂ x ṽ -v∂ y ṽ -w∂ z ṽ + u∂ y u + ṽ∂ y v + w∂ y w + ∂ y p -1 Re ∆ṽ = 0

-∂ t w -u∂ x w -v∂ y w -w∂ z w + u∂ z u + ṽ∂ z v + w∂ z w + ∂ z p - 1 Re ∆ w = 0 (1.24)
CHAPTER 1. NUMERICAL TOOLS or a more compact way:

-∇ • ũ = 0 -∂ t ũ -(u • ∇)ũ + (∇u) T ũ + ∇p - 1 Re ∆ũ = 0 (1.25)
and the compatibilty equation arises by matching the time integration by part:

ũ(x, T ) = 2 • u(x, T ) E(0) (1.26)
The gradient (update direction) with respect to the initial condition is obtained in a similar way :

∂L ∂u 0 = -2 • E(T ) E(0) 2 • u(x, 0) + ũ(x, 0) (1.27)
The adjoint equation (1.25) is linear with respect to the adjoint variable and as opposed to the direct N-S equation, the time derivative is preceded by a minus sign. One interpretation is to consider that the equation travels backwards in time by considering the change of variable t = Tt and the compatibility equation is the initial condition of the adjoint field.

Finite Element linear optimization tool

In the interests of speeding up the determination of linear optimal perturbations of the isolated Lamb-Oseen vortex, a one-dimensional tool is used implementing the Finite Element Solver Freefem++ [START_REF] Hecht | New development in freefem++[END_REF]. The results obtained from this tool are presented in Chapter 4. The linearized Navier-Stokes equations are solved to obtain normal mode solutions. The base flow composed of an isolated Lamb-Oseen vortex, is considered to be frozen. The velocity profile of the Lamb-Oseen vortex is U = (0, V (r), 0) where

V (r) = Γ 2πr 1 -exp(-r 2 /a 2 )
where a is the vortex dispersion radius and Γ is the circulation.

Spatial Discretization

As is the case with spectral element method, the Finite Element Method is based on a variational formulation of the equations. We use Taylor Hood finite elements of second order P 2 for the velocity components and of the first order P 1 for the pressure. These elements and their degrees of freedom are presented in figure 1.3.

FINITE ELEMENT LINEAR OPTIMIZATION TOOL

Time Discretization

We solve the linearized Navier-Stokes equations in cylindrical coordinates:

1 r ∂ r (ru) + 1 r ∂ θ v + ∂ z w = 0 ∂ t u + V r ∂ θ u - 2V r v = -∂ r p + 1 Re 1 r ∂ r (r∂ r u) + 1 r 2 ∂ 2 θ u + ∂ 2 z u - u r 2 - 2 r 2 ∂ θ v ∂ t v + ∂ r V u + V r ∂ θ v + V r u = - 1 r ∂ θ p + 1 Re 1 r ∂ r (r∂ r v) + 1 r 2 ∂ 2 θ v + ∂ 2 z v + 2 r 2 ∂ θ u - v r 2 ∂ t w + V r ∂ θ w = -∂ z p + 1 Re 1 r ∂ r (r∂ r w) + 1 r 2 ∂ 2 θ w + ∂ 2 z w (1.28)
After introducing q = qe imθ+ikz the equations read:

∂ r (rû) + imv + ikr ŵ = 0 r 2 ∂ t û + imrV û -2rV v = -r 2 ∂ r p + 1 Re r∂ r (r∂ r û) -m 2 + k 2 r 2 + 1 û -2imv r 2 ∂ t v + imrV v + r 2 ∂ r V + rV û = -imr p + 1 Re r∂ r (r∂ r v) -m 2 + k 2 r 2 + 1 v + 2imû r 2 ∂ t ŵ + imrV ŵ = -ikr 2 p + 1 Re r∂ r (r∂ r ŵ) -m 2 + k 2 r 2 ŵ (1.29)
The solution q n+1 at time t n+1 is computed as a function of q n and q n-1 . The lineat terms of the Navier-Stokes equations are handled semi-implicitly with the BDF2 scheme. The time derivative is expressed, at order 1 as:

∂ t u = u n+1 -u n ∆t (1.30)
and at order 2 as:

∂ t u = 3u n+1 -4u n + u n-1 2∆t (1.31)
CHAPTER 1. NUMERICAL TOOLS with ∆t the time step. The integration is started at order 1 and switches to order 2 after two iterations. This yields:

3 2 u n+1 -∆tL q n+1 = 2u n - 1 2 u n-1 (1.32)
with L q n+1 i the linear terms of the Navier-Stokes equations.

Normal Modes

Normal modes are solution of the linearized Navier-Stokes equations of the form q = qe σt , with q the normal mode and σ a complex frequency σ = σ R + iσ I with σ R the growth rate and σ I the frequency. After introduction into the linearized Navier-Stokes equations, one obtains:

∂ r (rũ) + imṽ + ikr w = 0

r 2 σ + imrV ũ -2rV ṽ = -r 2 ∂ r p + 1 Re r∂ r (r∂ r ũ) -m 2 + k 2 r 2 + 1 ũ -2imṽ r 2 σ + imrV ṽ + r 2 ∂ r V + rV ũ = -imr p + 1 Re r∂ r (r∂ r ṽ) -m 2 + k 2 r 2 + 1 ṽ + 2imũ r 2 σ + imrV w = -ikr 2 p + 1 Re r∂ r (r∂ r w) -m 2 + k 2 r 2 w
(1.33) These equations can be cast in matrix format Aũ = σB ũ with:

A =       -imrV + D -ν 2rV -2iνm 0 -r 2 ∂ r -r 2 ∂ r V -rV + 2iνm -imrV + D -ν 0 -imr 0 0 -imrV + D -ikr 2 r∂ r + 1 im ikr 0       (1.34) D is the viscous diffusion: D = ν r(∂ r + r∂ 2 r ) -(m 2 + k 2 r 2 ) B =       r 2 0 0 0 0 r 2 0 0 0 0 r 2 0 0 0 0 0       (1.35)
B is the mass matrix. Adjoint modes can be similarly defined using the adjoin Navier-Stokes equations:

A + ũ+ = σ + B + ũ+ with A + =       imrV + D -ν -r 2 ∂ r V -rV -2iνm 0 -r 2 ∂ r 2rV + 2iνm imrV + D -ν 0 -imr 0 0 imrV + D -ikr 2 r∂ r + 1 im ikr 0       (1.36) 1.3. FINITE ELEMENT LINEAR OPTIMIZATION TOOL
σ+ is the adjoint complex frequency. The adjoint is defined through the Hermitian dot product (u, v) = w H N u where the superscript H denotes the Hermitian conjugate and

N =       r 0 0 0 0 r 0 0 0 0 r 0 0 0 0 0       (1.37)
The direct and adjoint modes are normalized as such:

(w i , u j ) = δ ij (u i , u i ) = 1 (1.38)
The result of the adjoint simulation based on the B mass matrix needs to be mapped to the dot product based on N rather than B. This can be achieved by writing

w H i N B -1 Bu i = 1 = B -1 N w i H Bu i . After identifying the B-based adjoint field w i = B -1 N w i , the N -based
adjoint is obtained via applying w i = N -1 Bw i .

Numerical Resolution

The problem to solve is: B∂ t ũ = Aũ. Using the BDF2 scheme, the solution at time t n+1 is therefore: ũn+1 = (3B -2∆tA) -1 B (4ũ nũn-1 ) (1.39)

The same reasonning can be used for the adjoint equations. Finally, as for the optimization procedure described in 1.2 the iterative direct-adjoint process has been implemented using Freefem++.

Conclusion:

The aforementioned numerical tools are thorougly used in this thesis to describe and control the dynamics of vortices in ground effects. The following chapters focus on the use of these tools to hasten vortex decay. More details on the different algorithms and implementations are provided within each chapter.

CHAPTER

2

A Vortex Pair in Ground Effect, Dynamics and

Optimal Control

The contents of this chapter form a self-contained article, see [START_REF] Wakim | A vortex pair in ground effect, dynamics and optimal control[END_REF], available at the following hyperlink: Wakim, Arnold, Brion, Vincent, Dolfi-Bouteyre, Agnès & Jacquin, Laurent, 2020. A vortex pair in ground effect, dynamics and optimal control. Journal of Fluid Mechanics, 885.

In the following we describe and study the optimal control of a pair of counterrotating vortices in ground effect of in a two-dimensional, incompressible and laminar configuration. We present a conceptual approach to deal with aircraft separation distances in airport airspace by increasing the sideward position of the vortices, away from the runway. The dynamics and control of a vortex pair in ground effect are investigated in a planar, incompressible and laminar setting. The evolution of the vortices obtained numerically shows vortex rebound as a consequence of the separation of the boundary layer induced at the wall by the vortices. An optimal control approach is developed and employed for vortex Reynolds numbers of 200 and 1000 in order to identify the optimal Dirichlet boundary condition at the wall to counteract this rebound and allow for an increased lateral displacement of the vortex, similarly to the inviscid evolution of the flow, which features hyperbolic trajectories. The work is primarily a conceptual approach to deal with aircraft separation distances in airport airspace by moving the vortices laterally, away from the runway but may also apply to the control of coherent structures in wall bounded turbulence. The most efficient control is able to double the lateral position and yields mostly vertical in and outflow at the wall. An optimal horizon time is found, equal to 5 characteristic time units of the vortex system, beyond which control is not able to further displace the vortices. The control is shown to delay the separation of the boundary layer at the origin of vortex rebound by applying suction ahead of the vortex, and to generate a vorticity flux at the wall, leading to a pusher vortex of sign opposite to that of the primary vortex, that attenuates the effect of the no-slip boundary condition at the wall by pushing the vortex outward.

Key words: control theory, variational methods, vortex dynamics

Introduction

Air traffic has been in almost constant progress, with a yearly average increase of 5 %, during the past 50 years and will continue to grow as a consequence of the democratization of air transport in many parts of the world (Airbus SAS 2015). Among other issues, the increase in the number of flights has led to the saturation of major airports, especially at peak hours. In that context, separation distances imposed to mitigate the risk of wake vortices encountered at take-off and landing appear to be an important factor to play with in order to increase the frequency of arrivals/departures. The goal is to improve airport capacity without building new infrastructures, which is often impossible for technical, environmental or societal reasons. In the present work, we propose a new vortex mitigation strategy based on flow control applied at the ground. The action of control gives novel insight into the mechanisms affecting vortex evolution in ground effect and how these can be manipulated.

For an inviscid flow [START_REF] Lamb | Hydrodynamics[END_REF]) a pair of rectilinear counter-rotating vortices set parallel to the ground in a quiescent environment moves along hyperbolic trajectories. The initial descent is related to the downward momentum provided to the flow by the lifting surface. The subsequent lateral displacement is the consequence of the momentum induced by the image vortices about the ground. With viscosity, the phenomenon of vortex rebound occurs, as reported by [START_REF] Dee | Flight Measurements of Wing-Tip Vortex Motion Near the Ground[END_REF], i.e. the vortex leaves the hyperbolic trajectory, gains relative altitude and features little sideward motion. From an air traffic management perspective, this corresponds to a persistence of the vortices about the runway. As a consequence of the pragmatic importance of this matter and the required fundamental knowledge of vortex dynamics and turbulence near walls, vortex rebound has received much attention in the past. Following [START_REF] Harvey | Flowfield produced by trailing vortices in the vicinity of the ground[END_REF], vortex rebound results from the formation of a detached boundary layer at the ground, the former being formed due to the velocity field induced by the vortices, and the generation of opposite sign vorticity. The production of this secondary vorticity at the wall during the dipole-wall collision was in particular studied by [START_REF] Coutsias | Fundamental interactions of vortical structures with boundary layers in two-dimensional flows[END_REF] and [START_REF] Clercx | Dissipation of kinetic energy in two-dimensional bounded flows[END_REF]. [START_REF] Barker | The motion of two-dimensional vortex pairs in a ground effect[END_REF] explained vortex rebound as the elliptic deformation of the vortex as it nears the ground, with the subsequent oscillation of the ellipticity creating the vertical motion of the vortex centre. In [START_REF] Peace | A viscous vortex pair in ground effect[END_REF], a matched asymptotic expansion of the vortex induced boundary layer shows that the increase of vortex height is also to be attributed to the thickness effect associated with the developing transverse boundary layer. The rebound has been simulated by [START_REF] Orlandi | Vortex dipole rebound from a wall[END_REF], [START_REF] Coutsias | Fundamental interactions of vortical structures with boundary layers in two-dimensional flows[END_REF], Zheng & Ash (1996), Türk, Coors & Jacob (1999) and [START_REF] Clercx | Dissipation of kinetic energy in two-dimensional bounded flows[END_REF]. [START_REF] Zheng | Prediction of turbulent wake vortex motion near the ground[END_REF] used a specific vorticity-streamfunction computational method, a technique that allows for efficient calculation of the flow evolution, up to large Reynolds numbers (Zheng & Ash 1996;Türk et al. 1999).

The inviscid dynamics of the vortex pair would be beneficial to the airport problem. Indeed, the presence of viscosity, due to the induced boundary layer and entrainment effect associated with secondary vorticity, results in a dramatic persistence of the vortices below the flight path that strongly constrains the possibility for reduced separation distances between aircraft. This naturally suggests that a control aimed at removing the effect of viscous friction at the wall would be beneficial. Following this simple idea, the present work looks at how to apply blowing/suction of air at the wall to favour an improved lateral displacement of trailing vortices approaching ground. A previous work to control wake vortices at the ground is described in [START_REF] Holzäpfel | Enhanced wake vortex decay in ground proximity triggered by plate lines[END_REF] and uses a system of plate lines developed in a heuristic way.

It is noteworthy that the dynamics of vortices in ground effect generally results from multiple other effects, notably cross-wind, stratification, turbulence and three-dimensionality [START_REF] Robins | Potential hazard of aircraft wake vortices in ground effect with crosswind[END_REF]Zheng & Ash 1996;[START_REF] Gerz | Commercial aircraft wake vortices[END_REF][START_REF] Pailhas | Behaviour of trailing vortices in the vicinity of the ground[END_REF][START_REF] Harris | Instability of secondary vortices generated by a vortex pair in ground effect[END_REF][START_REF] Bricteux | Les investigation of the transport and decay of various-strengths wake vortices in ground effect and subjected to a turbulent crosswind[END_REF]Leweke, Le Dizes & Williamson 2016). Gentle side wind is in particular known to be particularly harmful for the persistence of wake vortices [START_REF] Gerz | Commercial aircraft wake vortices[END_REF]. In the current study, which targets the numerical implementation of the control and its application to the most simple case, none of these effects are accounted for. The resulting scenario, two-dimensional laminar vortices in a quiescent environment, is yet one of the worst cases for air traffic and consequently quite conservative since turbulence, stratification and three-dimensional effects all lead to an accelerated dissipation of the vortex system. The limitation to a two-dimensional physical set-up amounts to considering infinite trailing vortices, parallel to the ground.

Irrespective of aircraft wakes, vortices are elementary components of wall turbulence. The understanding of vortex/ground interaction is therefore of primary interest for all turbulent flows and flow control in this domain could play an interesting role in accomplishing objectives such as drag reduction, relaminarization, delay of flow separation and mixing enhancement. A significant part of the drag generated by wall turbulence results from the shear stress powered by coherent streamwise vortices located near the wall [START_REF] Kravchenko | On the relation of near-wall streamwise vortices to wall skin friction in turbulent boundary layers[END_REF]. The dynamics of these vortices includes sweep and lift (rebound) motions that resemble those of the large scale trailing vortices shed by airplanes. The typical Reynolds number involved is in the range O(10 2 -10 3 ). These vortices play an active part in the generation of turbulence at its early stages, by promoting longitudinal velocity streaks and the formation of secondary vorticity, whose destabilization leads the turbulence cycle. Conversely, their attenuation is able to trigger flow relaminarization [START_REF] Jiménez | The autonomous cycle of near-wall turbulence[END_REF]. In developed wall turbulence [START_REF] Lee | Application of neural networks to turbulence control for drag reduction[END_REF] and [START_REF] Choi | Active turbulence control for drag reduction in wall-bounded flows[END_REF] show that the application of a blowing/suction control strategy produces significant drag and turbulence reduction by modifying the vortex lateral and vertical motions. Along this line [START_REF] Koumoutsakos | Active control of vortex-wall interactions[END_REF] developed a sensor-based control strategy with mass exchange at the wall to manipulate wall vorticity. Vortex annihilation could even be attained this way. [START_REF] Akhavan | Turbulence control in wall-bounded flows by spanwise oscillations[END_REF] realized control by lateral blowing and oscillatory wall deformation and also observes that this shifts the vortices and in turn reduces turbulence production. In all these studies it seems that reducing vortex rebound is a way to temper turbulent production. Such a feature could also be of interest to maintaining longitudinal vortices for longer times near the wall and controlling flow separation using vortex generators [START_REF] Lin | Review of research on low-profile vortex generators to control boundary-layer separation[END_REF]. One could also find interest in such controls for turbomachinery tip leakage vortex flow. [START_REF] Storer | Tip leakage flow in axial compressors[END_REF] mentions that the tip leakage region leads to increased flow blockage and the tip vortex is a source of instabilities in the near stall regime [START_REF] Furukawa | The role of tip leakage vortex breakdown in compressor rotor aerodynamics[END_REF].

The present work adopts an optimal control strategy, as detailed by [START_REF] Bewley | Flow control: new challenges for a new renaissance[END_REF][START_REF] Choi | Instantaneous control of backward-facing step flows[END_REF][START_REF] Choi | Instantaneous control of backward-facing step flows[END_REF], to do so. In the case when a large set of degrees of variables needs to be optimized, the associated direct-adjoint technique is particularly well suited. Here, the variable to be optimized is the velocity at the wall and falls into that category. The technique has been successfully applied for boundary layers and wakes [START_REF] Walther | Optimal control of Tollmien-Schlichting waves in a developing boundary layer[END_REF][START_REF] Homescu | Suppression of vortex shedding for flow around a circular cylinder using optimal control[END_REF][START_REF] Airiau | A methodology for optimal laminar flow control: application to the damping of Tollmien-Schlichting waves in a boundary layer[END_REF][START_REF] Guégan | Optimal energy growth and optimal control in swept hiemenz flow[END_REF][START_REF] Flinois | Optimal control of circular cylinder wakes using long control horizons[END_REF]. In [START_REF] Önder | Optimal control of a transitional jet using a continuous adjoint method[END_REF], the optimal control method is used to improve the mixing of an incompressible axisymmetric jet and in [START_REF] Bewley | Dns-based predictive control of turbulence: an optimal benchmark for feedback algorithms[END_REF] flow relaminarization and drag reduction in a plane channel are demonstrated. The control algorithm for an unsteady flow and nonlinear dynamics typically requires a large amount of storage as the full flow field information is needed at each time step. As a result, sub-optimal boundary control strategies have been designed to reduce the associated cost of computation. Such a strategy was for instance applied by [START_REF] Choi | Instantaneous control of backward-facing step flows[END_REF] for the case of the recirculation bubble behind a backward-facing step.

In the present study, linear interpolation is used to reduce this storage issue.

A26-4

A. Wakim, V. Brion, A. Dolfi-Bouteyre and L. Jacquin

The paper is organized as follows. In § 2 we define the governing equations and present the numerical set-up to compute the dynamics of the vortices in ground effect. Section 3 presents the principal features of the uncontrolled flow. In § 4, we outline the optimal control approach that allows the maximization of the horizontal position of the vortices. In § 5 the strategy is applied, the various results are presented and we discuss the effect of optimal control on wake vortices. Conclusions and outlooks for application are provided in § 6.

Governing equations

In the following, the fluid domain is (x, y 0) where y = 0 is ground. The flow is started by a pair of counter-rotating Lamb-Oseen vortices, located at (x c 0 , y c 0 ) = (±0.5b, y c 0 ) with initial circulation Γ 0 . Initially the vortex dispersion radius that defines the vortex velocity profile is set to a = 0.15b. Length and velocity are normalized on vortex separation b and on the drift velocity of the vortex pair v * = Γ 0 /2πb, yielding a reference time T * = 2πb 2 /Γ 0 and a circulation-based Reynolds number Re = Γ 0 /2πν where ν is the kinematic viscosity of the fluid. The Reynolds number is considered in the range 200 Re 5000. For the remainder of the paper and unless otherwise stated, all the quantities are normalized upon the aforementioned variables. The flow is governed by the incompressible Navier-Stokes equations,

∇ • u = 0, ∂ t u + (u • ∇)u = -∇p + Re -1 ∇ 2 u, u(x, y = 0, t) = u w ,    (2.1)
where u w is the velocity at the wall, which is zero in the uncontrolled case, and optimized in the controlled case. These equations are solved using the incompressible computational tool Nek5000 (Fischer, Lottes & Kerkemeier 2008). Numerical simulations have been carried out using parallel computations with 128 cores and solving the set of (2.1) limited to the right half of the domain, denoted Ω, shown in figure 1 with a mesh containing 143 × 129 = 18 447 elements. Each element contains 8 × 8 Gauss-Lobatto points, leading to 1.2 M discretization points. Boundary conditions are symmetry conditions at all frontiers except the wall where Dirichlet boundary condition is applied. Numerical validation is carried out using comparisons of vortex trajectories against data from different simulations published by Türk et al. (1999) and Zheng & Ash (1996). Usually the vortex centroid M c = (x c , y c ) is obtained by the first moment of the normal to the plane vorticity field ω

x c = Ω xω dΩ Ω ω dΩ ; y c = Ω yω dΩ Ω ω dΩ . (2.2a,b)
However, in the present configuration, M c does not match the primary vortex locus M pv = (x pv , y pv ) because of the boundary layer and secondary vortices. Instead, the vortex locus M pv is taken as the barycentre of the vorticity field ω pv given by Zheng & Ash (1996) are respectively represented with the dashed, dash-dotted and solid lines and the trajectories obtained in the simulations are marked with squares, triangles and circles.

ω pv (x, y) = ω(x, y) • [ω(x, y) > 0] • [Q(x, y) > 0], ( 2
∇u [START_REF] Hunt | Eddies, streams, and convergence zones in turbulent flows[END_REF]Kolář 2007). The barycentre leading to M pv is then obtained by the same volume integration as in relation (2.2), with ω replaced by ω pv .

The comparison against Türk et al. (1999) is conducted for Re 0 = Γ 0 /ν = 330 and 3300 and that against Zheng & Ash (1996) for Re 0 = Γ 0 /ν = 7650 (note the use of Re 0 = 2πRe). The flow is initialized by a pair of Lamb-Oseen vortices localized at ±(x c 0 , 2x c 0 ) and (±x c 0 , 3.8x c 0 ) to match the settings of Türk et al. (1999) and Zheng & Ash (1996), respectively. The ensemble of trajectories, shown in figure 2, provides comparisons that validate the present numerical approach.

Uncontrolled flow

In potential flows, vortices initially located in (±x c 0 , y c 0 ) follow a hyperbolic trajectory (x c (t), y c (t)) described by

1 x 2 c (t) + 1 y 2 c (t) = 1 x 2 c 0 + 1 y 2 c 0 . ( 3

A26-6

A. Wakim, V. Brion, A. Dolfi-Bouteyre and L. Jacquin ). Note the existence of an average cross-wind (u crosswind ∼ -1 m s -1 ) in the lidar case, not taken into account in the simulated data.

With viscosity, as seen previously in figure 2 and now in figure 3 which plots the vortex trajectories for Re = 500 and 1000, the dynamics is modified when the vortices near the ground. This so-called rebound results in an oscillation of the altitude of the vortex and a reduction of its sideward motion. It is noteworthy that such an observation is coherent with practical, high Reynolds number cases of aircraft wake vortices impacting ground, as is shown by the lidar (light detection and ranging) measurements also plotted in figure 3. The data come from a set of Doppler lidar measurements performed at Toulouse-Blagnac airport [START_REF] Hallermeyer | Traitement du signal dun lidar doppler scannant dédié á la surveillance aéroportuaire[END_REF]. From the simulated data it can be remarked that with a larger Reynolds number the vortex follows the inviscid path up to a larger distance. Rebound yet remains practically unchanged, given that both the current simulations and the real life lidar measurements show identical overall kinematics. Rebound is illustrated in figure 4 with plots of the vorticity field at various instants of the vortex evolution. The figure also features a set of velocity profiles in the near wall region and the wall-wise pressure gradient. The velocity profiles show the flow at the wall that is induced by the momentum brought by the vortices. As initially described by [START_REF] Harvey | Flowfield produced by trailing vortices in the vicinity of the ground[END_REF] and due to the large Reynolds number, this subsequently forms a boundary layer. The vortex also imposes a positive pressure gradient at its sides because of the low pressure present in its core. As a consequence of this (adverse) pressure gradient, the boundary layer separates. This occurs when the vortex gets closer to the ground (figure 4b). The separated flow then rolls into opposite sign vorticity, soon embarked around the primary vortex (c). Additional detachments of the separated boundary layer later form other vortices, which themselves induce an opposite boundary layer flow (d).

The question of how the vortex circulation evolves in this context can be analysed using the circulation of the right part of the flow Γ tot = Ω ω dΩ, whose evolution is to be inferred from the equation of the vorticity ω

dΓ tot dt = 1 Re ∂Ω ∇ω • n∂Ω = - 1 Re ∂Ω 1 ∂ω ∂y dx + ∂Ω 4 ∂ω ∂x dy . (3.2)
Since ∂ω/∂x > 0 at ∂Ω 4 , the left boundary, and ∂ω/∂y > 0 at ∂Ω 1 , the bottom boundary, total circulation is expected to decrease. Physically the second and first terms on the right-hand side relate to friction with the left vortex and ground, respectively. It is important to mention that total circulation Γ tot is not equal to Γ 0 initially because of the presence of the boundary layer, which features negative vorticity.

To characterize the vortex circulation, Γ pv is introduced, which is defined as The configuration imposes Γ pv = Γ 0 initially. The zone Ω pv defines the spatial extent of the primary vortex. Here, to discriminate it from the several vortical zones featured by the flow, Ω pv is defined as the zone of positive vorticity about the primary vortex locus M pv = (x pv , y pv ). Note that the simple integration of ω pv over the half-domain Ω, as used to define M pv in (2.3), does not yield the vortex circulation as desired, as, due to the Q-criterion constraint, it misses some part of the vortical area associated with the vortex. The zone Ω pv cannot be easily described and therefore a specific procedure is employed to automatically track its shape in the snapshots of the flow. The first step is to locate the primary vortex centre using procedure (2.3). Once the vortex centre is known, the method consists in using a cylindrical coordinate system about M pv and to locate at all azimuthal positions θ the radial extent R where ω pv (R, θ ) = 0. In practice, an ensemble of N θ = 128 azimuthal positions θ i is considered, and the radial direction is discretized into R i with a spacing of δ r = 0.01. The ensemble of (θ i , R i ) forms a polygonal area that is taken as Ω pv .

Γ pv = Ω pv ω dΩ pv . ( 3
Figure 5 shows the evolution of Γ tot and Γ pv as a function of time for the two Reynolds numbers Re = 200 and 1000. Total circulation Γ tot is seen to decrease significantly up to t = 4 and to level off after that. This can be analysed in the light of figure 4, which corresponds to the case Re = 1000. Indeed the development of the boundary layer, its separation and the formation of the secondary vortex (see figure 4a-c) explain the initial strong decrease of Γ tot in figure 5 while the mild decrease of total circulation after t = 4 stems from the end of these processes and corresponds to the slow viscous interaction between the different vorticity zones.

It is observed that the evolution of Γ tot hardly changes between Re = 200 and 1000, suggesting a weak Reynolds number dependency of this quantity. The formulation of circulation written below

Γ tot = ∂Ω u • dl, (3.4)
can help understand why this is so. In (3.4) dl is the tangent vector along the closed contour ∂Ω. In the no-slip case, the contour integration reduces to the symmetry boundary ∂Ω 4 . Indeed, the velocities at the top and right boundaries are very small and their contribution can thus be neglected. Equation (3.4) can hence be rewritten as where ∂Ω 4 is oriented downward. The vertical velocity at the symmetry plane v(0, y) essentially evolves as a function of the lateral distance of the vortex and one expects little influence of the no-slip boundary condition applied at the ground. As a check, we introduce Γ inv = -∂Ω 4 v inv (0, y) dy, the contribution of the left symmetry frontier to the total circulation in the inviscid case. The velocity v inv results from the Biot-Savart law applied to the point vortex along the hyperbolic trajectory and its images. This circulation can be calculated following

Γ tot = - ∂Ω 4 v(0,
Γ inv = ∂Ω 4 Γ 0 2π 2x c x 2 c + (y -y c ) 2 - 2x c x 2 c + (y + y c ) 2 dy. (3.6)
Crossed markers in figure 5 show the time evolution of total circulation using (3.6).

The evolution of Γ inv follows closely that of Γ tot provided by (3.2) up to t = 3. Departure from the inviscid evolution appears afterwards, indicating the impact of no slip at the ground on the lateral position of the vortex. However the difference between Γ inv and Γ tot remains small, which concludes on the weak influence of Reynolds number on Γ tot . The circulation of the primary vortex Γ pv shows a different evolution than Γ tot , with three distinct phases and a net effect of Reynolds number. The first phase, that corresponds to the initial descent of the vortices, up to t = 2.5, preserves circulation, i.e. Γ pv Γ 0 . The second phase, up to t = 5, features a sudden drop of circulation and coincides with the rebound of the vortex. In this step the vortex experiences strong friction at the wall because ∇ω • n is increased by the presence of opposite sign vorticity at the wall (boundary layer). The decay rate then slightly levels off, corresponding to the third phase of evolution. Vortex decay in the two last phases shows a clear dependence upon the Reynolds number. The start of the third phase (t 5) corresponds to the slow decay phase of the total circulation Γ tot , and to the quasi-disappearance of the pressure gradient at the wall, which goes along with the quasi-suppression of the boundary layer, see figure 4. The vortex system altitude is now large due to the previous rebound and the friction at the wall is not as significant.

Optimal control approach

We consider the effect of either blowing or suction applied at the ground to counteract the stagnation of the vortex system and promote increased lateral displacement. Practically, a modification of the boundary condition at ∂Ω 1 in (2.1) is applied with u ∂Ω 1 = u w (x, t) where u w (x, t) is free to vary in intensity and orientation. The control is restricted to a zone of finite extent ∂Ω 1,w = {0 x 10}. Figure 6 provides a schematic of the control set-up.

The optimal control strategy requires an appropriate objective function J to optimize u w and maximize the lateral position of the vortex system at a given horizon time T. A choice could have been the first moment of vorticity Ω xω dΩ as it evaluates, once divided by Γ tot (see (2.2)), the horizontal component of the barycentre of vorticity.

However, because of the different vorticity areas, which yield opposite signs, this optimization would lead to maximizing the lateral position of the positive vorticity (primary vortex) and minimizing that of the negative vorticity (separated boundary layer). Therefore, the maximization of the lateral position of the vortex seems not to be guaranteed. A way to shape the objective function to our goal is to consider instead the square of vorticity ω 2 . The benefit is that it mechanically embarks the vorticity 885 A26-10 A. Wakim, V. Brion, A. Dolfi-Bouteyre and L. Jacquin

(u w , √ w )(x, t) ™Ø 1,w FIGURE 6
. Schematic of the control scheme applied along the portion ∂Ω 1,w of the wall. u w and v w denote the flow velocity imposed at the wall.

of the primary vortex and the opposite vorticity coming from the separated boundary layer alike.

The implementation hence proceeds from an evaluation of the horizontal position of the vortex based on the horizontal component of the barycentre of vorticity squared, that is

J (u; u w ; T, l) = J 1 (u; u w ; T, l) -lJ 2 (u w ; T, l), (4.1) with J 1 (u; u w ; T, l) = Ω xω 2 (T) dΩ Ω ω 2 (T) dΩ , (4.2)
and

J 2 (u w ; T, l) = 1 2 T 0 ∂Ω 1 u 2 w dx dt. (4.3)
The term J 2 (spatio-temporal norm of the control) penalizes J 1 by acting as a cost function for control. The mathematical parameter l > 0 (which holds no physical interpretation) prevents u w from reaching an unphysically large amplitude and increases the convergence of the optimization algorithm. As such, it favours the existence and unicity of a solution as well as eases convergence of the optimization loops. The case l → ∞ amounts to the uncontrolled case i.e. J m (•, ∞) = J 1 (•, ∞) since lim l→∞ J 2 = 0 (a validation is provided further down, see § 5.1). The maximum J m (T, l) of J in the space of u w for a given T and l is found using the Lagrange multiplier method in which the constrained optimization problem is transformed into an unconstrained one. To this end, the cost functional and the constraints are assembled into the Lagrangian function L, in the form

L(u, u + , p, p + , u w , u + w ) = J (u w ) -∂ t u + (u • ∇)u + ∇p - 1 Re ∇ 2 u, u + -[∇ • u, p + ] -(u -u w , u + w ), (4.4)
where [•, •] and (•, •), denote the following inner products The Lagrangian function L depends on the direct (u, p, u w ) variables and Lagrange multipliers (u + , p + , u + w ), also known as adjoint variables. The optimal control is found by cancelling the first variation of the Lagrangian with respect to its variables. For a complete derivations of these equations the reader is referred to [START_REF] Choi | Instantaneous control of backward-facing step flows[END_REF].

[a(x, t), b(x, t)] = T 0 Ω a(x, t) • b(x, t) dΩdt, (c(x), d(x)) = ∂Ω 1 c(x) • d(x) dx,          (4.
Here, only the most important relations are provided. The first variation of L with respect to the adjoint variables yields the equations of the constraints and the first variation with respect to the direct variables provides a set of equations for the adjoint variables,

∇ • u + = 0, -∂ t u + -(u • ∇)u + + (∇u) u + = -∇p + + 1 Re ∇ 2 u + , (4.6)
where denotes the transpose operator. The minus sign in front of the time partial derivative can be interpreted as a change of variable t -→ Tt. As a consequence, the adjoint equation is solved backward in time from t = T to t = 0. By matching terms of the integration by part of the partial time derivative, one obtains the following compatibility equation (which is derived in appendix A) for t = T

u + (x, T) = 2 I 0 (T) (x∇ × ω(T)e z + ∇x × ω(T)e z ) - 2I 1 (T) I 0 (T) 2 (∇ × ω(T)e z ), (4.7)
where

I 0 (t) = Ω ω 2 (t) dΩ, (4.8)
and

I 1 (t) = Ω xω 2 (t) dΩ. (4.9)
Relation (4.7) initializes the adjoint variable u + (x, T). The gradient of the Lagrangian with respect to u w

∇ u w L = -lu w + u + w = -lu w + p + (x, y = 0, t) - 1 Re ∂u + ∂y (x, y = 0, t) e y , (4.10)
eventually provides the direction to follow in order to update the control variable. This optimal system is based on a set of nonlinear partial differential equations. While this nonlinearity could lead to difficulties in converging an optimal solution, in the present case the dynamics of the vortex to be controlled follows a deterministic evolution, so that similarly to [START_REF] Önder | Optimal control of a transitional jet using a continuous adjoint method[END_REF], a steepest descent algorithm with a generic backtracking line search is successfully used in order to perform the optimization. The whole procedure is sketched in figure 7. Starting from an initial u w,0 = 0, the algorithm follows two loops. The main loop goes towards the optimal velocity at the wall by computing the gradient and updating u w with a factor α k of the gradient. The second loop concerns the line search to adapt the value of α k when the chosen value does not optimize. The starting value is 0.5. Whenever the chosen value does not lead to optimization, α k is halved and the direct simulation is repeated until optimization is achieved. The iterative algorithm is eventually stopped when (J k+1 -J k )/J k ε = 10 -4 . In the current algorithm the cost of storing the control variable when integrating the adjoint is decreased by using linear interpolation over uniform time intervals

[t i , t i+1 ] of the horizon time T u w (t k ) = β k u w (t i ) + (1 -β k )u w (t i+1 ), (4.11)
where β k = (t i+1 -t k )/ t and δt i = t i+1 -t i is chosen such that δt i 5 %. The convergence of the results is provided in appendix B.1 and convergence in δt i is shown in appendix B.2.

Controlled flow

Optimal parameters

Optimizations have been carried out for T ∈ [2, 10] and l ∈ [1, 100], and Reynolds numbers of 200 and 1000.

In figure 8(a) the cost of control J 2 is plotted as a function of l. It shows, as expected, that growing l returns to the uncontrolled case J 2 → 0. Low values of l have been tested and have shown that when l < 1 the optimization breaks down: the intensity of u w becomes too high for the capability of the mesh and time step that are used (which correspond to the initial simulation of the uncontrolled case). Figure 8(b) shows the optimized vortex barycentre J 1 normalized upon its value for infinite l (no control) as a function of l for Re = 200 and 1000. It shows that the gain increases as l is lowered, reaching a maximum of 1.5 and 2 for l = 1 at the two Reynolds numbers, respectively. Increasing Reynolds number yields a favourable effect for augmenting control capacity. Plots of the vorticity field at final time t = T = 5 for selected values of l ∈ [1, ∞] and the two values of Reynolds number are displayed in figure 9. They show the effect of control on the separation pattern and timing of the boundary layer and the intensity, location and number of the secondary vortices that rotate around the primary one.

x 0 1 3 2 x 0 1 3 2 x 0 1 3 2 x 30 0 -30 l = 1 l = 10 l = 100 l = ∞ l = 1 l = 10 l = 100 l = ∞ (a) (b) (c) (d) (e) (f) (g) (h)
The effect of horizon time T is shown in figure 10 for Re = 200 and Re = 1000 at l = 1, the best value within the range under consideration. Figure 10(a) plots the absolute position of the vortex J 1 (T, l) for the uncontrolled and controlled cases. The increase in lateral distance provided by control is apparent from the offset between the two curves. Above T = 5, the lateral position of the vortex shows a limit at approximately x = 3 and no decisive increase beyond that time. Since l = 1 provides the strongest control effect, this shows that optimization at a larger time is not required: T = 5 appears to be quasi-optimal for the present parameter space.

For Re = 200, it is found that optimization cannot be carried out beyond T = 5. Indeed, the routine produces irregular vorticity at the end of the domain ∂Ω 1,w , which is attributed to the finite length of the control area at the ground. Coincidentally, the vortex is also sucked by the control at the wall and disappears. Calculations beyond T = 5 are not conducted for Re = 200.

Figure 10(b) shows the amplification of the full functional J m between the controlled and uncontrolled cases. It shows that the optimal for T = 5 is also the time for which the dynamics can be optimally balanced between the possibility of increased lateral displacement and cost of control. This shows that an optimal time exists over the parameter range investigated here, at approximately 5.

Physical mechanisms

The physics behind the action of control is now analysed. The time sequence displayed in figure 11, which is obtained for T = 5 and corresponds to the same Reynolds number (Re = 1000) as the uncontrolled one shown in figure 4, is used in

A26-14

A. Wakim, V. Brion, A. Dolfi-Bouteyre and L. Jacquin 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 3.0 2.5 2.0 1.5 1.0 0.5 0 No data are consequently available for larger time.

3 2 1 T T j 1 (T, 1) j m (T, l = 1)/ j m (T, ∞) (a) (b)
0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 2 1 2 1 0 30 0 -30 x y y x x x t = 0.25 t = 1 t = 1.25 t = 2.5 t = 3 t = 3.75 t = 4 t = 5 (a) (b) (c) (d) (e) (f) (g) (h)
FIGURE 11. Vorticity field in the controlled case, computed for (Re, T, l) = (1000, 5, 10), showing the sequence up to the optimization time t = 5. Thicker borders refer to the times used in figure 4 showing the uncontrolled case at the same Reynolds number, with the same iso-contours. Note however the different spatial extensions of the plots.

preamble. Some of the time indices used in figure 4 are also repeated in figure 11, to ease the comparison. The main effect of control which can be seen in this time sequence is the delay in boundary layer separation, which is only apparent near the final time t > 4, while in the uncontrolled case (figure 4) it is already effective at t = 2.5. The doubling of the lateral distance, presented in the previous section, is also directly visible here, changing from x 1.2 (see figure 4d) to x 3 when control is active. Interestingly, the vortex kinematics under the effect of control bears a strong similarity to the vortex evolution obtained with the v-control scheme developed by [START_REF] Choi | Active turbulence control for drag reduction in wall-bounded flows[END_REF] (figure 22b). The evolution of the velocity at the wall imposed by control is displayed in figure 12 with some of the same time stamps as in figure 11. It is observed that control initially applies blowing at the centre of the two vortex system x ∈ [0, 1]. After t = 2, suction starts to act, essentially below the vortex, with increasing strength up to

0 1 2 3 x 0 1 2 3 x 0 1 2 3 x 0 1 2 3 x 2 0 -2 2 0 -2 u w , √ w u w , √ w t = 0.25 t = 1 t = 1.25 t = 2.5 t = 3 t = 3.75 t = 4 t = 5 (a) (b) (c) (d) (e) (f) (g) (h)
FIGURE 12. Evolution of the optimal control u w along the control area for different time steps. The thicker borders refer to the same times as those used in figure 4. Solid (respectively dashed) lines represent v w (respectively u w ).

0 1 2 3 4 5 2 0 -2 u w √ w t FIGURE 13
. Evolution of the optimal control velocity at the wall as a function of time at x = 1. Solid (respectively dashed) lines represent v w (respectively u w ).

t = 3. At this time, the primary vortex has reached its lowest altitude (see figure 11e). Subsequent control essentially reduces to slight blowing at the left of the primary vortex. After t = 2, suction starts to act, essentially below the vortex, with increasing strength up to t = 3. At this time, the primary vortex has reached its lowest altitude (see figure 11e). Subsequent control essentially reduces to slight blowing at the left of the primary vortex. Figure 13 shows the time trace of the control at x = 1, that is, under the descent path of the vortex. The suction then blowing strategy of the control described above is grasped immediately using this plot, along with the dominance of the vertical velocity component upon the lateral component. Two main factors have been identified to explain how control physically achieves increased lateral displacement. One is related to (i) delayed separation, which reduces the vortex updraft by secondary vorticity, hence maintaining lower vortex altitude and increased induction effect by image vorticity about the ground, and the other one to (ii) reduced effective friction at the ground to decrease the opposition to the vortex lateral displacement.

The factor (i) relates to the suction that is applied by control in the region x ∈ [1, 1.5], slightly to the right of the primary vortex. Suction is a generic principle to delay boundary layer separation [START_REF] Lachmann | Boundary Layer and Flow Control: Its Principles and Application[END_REF]. Its interest in the present The second factor (ii) is vortex induction. It can be identified using figure 11(d,e), which indicates that a specific vorticity pattern is formed at the wall by control between t = 2 and 3, this pattern being absent in the flow at the same time when control is absent (figure 4b). The field of Q-criterion between the uncontrolled and controlled cases at t = 2 and focused in the region between the primary vortex and the wall is shown in figures 14(a) and 14(b), respectively, to highlight the associated structure. The plot shows that a rotational area is brought about by control which, referring to figure 11(e), amounts to negative vorticity. Albeit the zone does not present a pressure minimum (at least not apparently), for simplicity's sake we coin it as a vortex. Considering Biot-Savart induction, it can be stated that the effect of this vortex upon the primary vortex favours the motion to the right of the latter, and effectively pushes it along the wall. In the following, this vortex is conveniently referred to as the pusher vortex.

The formation of this pusher vortex can be directly linked to the velocity imposed by the control at the wall. Indeed figure 12 shows that the vertical velocity component generates a strong gradient ∂v/∂x below the primary vortex for time t = 2 to 3 (figure 12d,e) and this is to be analysed as creating the ω associated with the pusher vortex. In effect ω ∂v/∂x since |u w | |v w | and characteristic length scales for u w and v w are similar.

An interesting consequence of control concerns circulation, which is found to be significantly diminished with the application of control. The behaviours of Γ tot and Γ pv are shown in figure 15(a) for Re = 200 and 1000, T = 5 and l = 1. Comparison can be made with the uncontrolled case in figure 5. For Re = 200, the circulation drops by more than 60 % at t = 5 compared to 50 % after 15 time units in the uncontrolled case. For Re = 1000, the decrease is approximately 15 % compared to less than 10 % in the uncontrolled case, also at t = 5.

Figure 15(b) sums up the previous analysis by showing the trajectory of the controlled vortex against the uncontrolled and inviscid trajectories, with equal time stamps to enable the comparison of the vortex progress in each case. The controlled vortex is initially moved above the uncontrolled trajectory as a consequence of the initial blowing imposed by control. Then, as the uncontrolled vortex rebounds, it moves laterally and close to the ground, reaching an altitude lower than the inviscid ˝p√ (200,5,1) ˝p√ (1000, 5, 1) ˝tot (200,5,1) ˝tot ( 1000 hyperbolic path. At the final time, the controlled vortex achieves a more outward position than the inviscid vortex, concluding the efficiency of the blowing/suction strategy.

Control robustness

In this subsection, we discuss the robustness of the optimal control approach with respect to the initial height, radius and activation time of the control law. It is noteworthy that assessing the robustness of the method with respect to the initial vortex circulation would imply only varying the circulation-based Reynolds number. Previous results comparing Re = 200 and 1000 show no major modification of the present strategy with Reynolds number. Three distinct simulations H (for height), R (for radius) and DAT (for delayed activation time) have been carried out where we modify the initial height, dispersion radius and activation time separately. In simulation H, all parameters are kept the same, except for the initial height y c 0 = 5 (instead of 2.5) and the horizon time (T = 10 instead of 5). Note that doubling the descending distance justifies doubling of the horizon time in order to compare with the case at (1000, 5, 1) and y c 0 = 2.5. Simulation R keeps all the same parameters as for the case at (1000, 5, 1) except for the initial dispersion radius a = 0.1 (instead of 0.15). In simulation DAT, we use the optimal control law of the optimal case (1000, 5, 1) for the situation when the vortex initially departs from y c 0 = 5. Control is then only active between t a = 2.5 and t = 7.5. Table 1 summarizes the objective function results and compares them to the reference case (1000, 5, 1), y c 0 = 2.5 and a = 0.15. The performance of the control, when looking at the value of J 2 , is moderately affected for cases H and R, with a relative 18 % decrease, and is unchanged for case DAT. This suggests that the control does not suffer much from changes in the parameters of the present model. The final vorticity fields and time trace of the optimal control at x = 1 are shown in figure 16(a,b) for cases H and R. The case DAT leads to negligible changes compared to the reference in figures 11(h) and 13 and is thus not shown in figure 16. Similar behaviours and amplitudes of the cost function are observed between the optimal reference case (1000, 5, 1) and cases H, R and DAT. In particular, the time traces shown in figure 16 resemble those in figure 13(c,d). To be compared to the optimal case (1000, 5, 1), see figures 11 and 12.

5.4. Vorticity flux at the wall The physical analysis can be backed by looking theoretically at the evolution of I 1 (t), the numerator of J 1 (t), introduced in § 4 (see (4.8)). Following the calculations given in appendix C, the time derivative of I 1 yields

İ1 = - 2 Re Ω x|∇ω| 2 dΩ + 1 Re ∂Ω 1 x(∇ω 2 ) • n dx + Ω uω 2 dΩ = İ11 + İ12 + İ13 . (5.1)
This equation shows that I 1 is reduced by the dissipative term ( İ11 ), and increased by the vorticity flux at the bottom boundary ( İ12 ) which also can be written as

İ12 = - 2 Re ∂Ω 1 xω ∂ω ∂y dx, (5.2)
and by the advection of the vorticity field ( İ13 ).

Figure 17 plots the evolution of İ11,12,13 for the controlled and uncontrolled cases at Re = 1000. The term İ13 is almost unaffected by control, unlike the terms İ11 and İ12 , which are increased by nearly 3 orders of magnitude when reaching their maximum values (shortly after t = 2). The term İ12 becomes larger than İ11 , which implies the net increase of I 1 and also means that the vorticity flux provoked by control at the wall, which results in the negative pusher vortex being present there, is the dominant mechanism at Re = 1000. As explained previously, this pusher vortex creates a horizontal induction to the primary vortex which minimizes the impact of zero velocity at the wall. This way, it allows the primary vortex to slide along the wall and move to larger lateral distances.

Conclusion

This work shows the effective control of vortex rebound at the ground in the situation of a laminar set of two opposite symmetric vortices, in a quiescent environment. The physics of vortex rebound results from the interplay between the vortex and the boundary layer it induces. An optimal control approach is employed with mass flow exchange at the wall to counteract this rebound and allow for a larger lateral displacement of the vortex, a concept initially devised for the interest of wake vortices in airport airfields. At the same time, the optimization offers novel insights into the dynamics of vortices impacting the ground and the effective physical leverages for its modification. The optimal blowing/suction strategy is shown to target a delayed separation of the boundary layer formed at the wall by a predominant suction of the flow ahead of the vortex, and the formation of a pusher vortex, of opposite sign vorticity compared to the primary vortex, to reduce the no-slip constraint at the wall. This helps the lateral movement of the vortex by mutual induction by the inserted vortex. Importantly, it is found that blowing laterally to the wall is essentially not required. In the most efficient control case, the lateral movement of the vortex is almost doubled, reaching a larger lateral distance than without viscosity. At the lower Reynolds number Re = 200 envisaged in this study, the control is also found to suppress the vortex. From a broader point of view, the findings of the present analysis show perspectives for vortex control of the coherent structures embedded in wall turbulence, in line with the conclusions of [START_REF] Akhavan | Turbulence control in wall-bounded flows by spanwise oscillations[END_REF], [START_REF] Choi | Active turbulence control for drag reduction in wall-bounded flows[END_REF] and [START_REF] Lee | Application of neural networks to turbulence control for drag reduction[END_REF], all the more so as the Reynolds number range targeted here is coherent with this type of flow.

40

For the problem of airport and airport terminal management (ATM), the effectiveness of the control potentially makes it an appealing strategy to apply ahead of runways to reduce vortex encounter. From a technical point of view, blowing and sucking at the ground could be achieved by compressed air supply and sucking devices (e.g. vacuum chambers or ejector suction devices) wrapped into a web of holes on the ground, vertically oriented, and a controller system for the multiple valves involved. This would still represent a great technical challenge. Moreover, in spite of the apparent generality of the mechanisms exposed in the present study, there are also physical questions on the variety of situations that could be tackled by such a system. Further analysis of the effect of turbulence, presence of wind, thermal updrafts and terrain effects to name a few would be necessary to fully qualify the robustness of the control strategy. Some of these questions could be answered by carrying out higher Reynolds number simulations of the control laws developed here, including possibly variations in the boundary conditions to incorporate, e.g. side wind. Note, however, that the effect of Reynolds number on the early stages of vortex dynamics is small, which justifies the choice of low Reynolds analysis performed here, but becomes important for the development of turbulence in the boundary layer flow at the ground and turbulent interactions with the vortex. As it is, the current algorithm stands as a proof of concept that shows both the physical interest in optimizing what happens at the ground and the working of the control in a simple configuration.

Another side question on the application of the concept is that of parallel runways. Indeed, as the vortices are pushed to the sides, they may interfere with the nearby runway. The lateral distances realized in the present work are of the order of 3 to 4 wingspans at most (we neglect here the ratio between vortex separation and wingspan, as it is close to unity). A non-exhaustive list of some 86 airports that operate parallel runways is available in libraries dedicated to flight safety (see for instance [START_REF] Skybrary | Parallel runway operation[END_REF]). For airports with parallel runways far apart, such as London Heathrow airport, the matter is not a question. For airports with nearby runways, such as Philadelphia International airport, the question becomes important as the separation distance, in this case approximately 400 m, ranges between 10 wingspans for small airliners to 5 for heavy ones. In the present analysis however, it is found that as the vortex control stops, the vortex rebounds and almost stops its sideward motion. This means that the control concept could certainly cope with an additional constraint of maximum lateral distance for the vortex system.

At this stage certainly one of the most interesting discussions is on the cost of the control. Using the normalization on Γ 0 and b used in the present study, the power P required to apply control can be evaluated by integrating the flux of kinetic energy at the wall, that is

P = T 0 ∂Ω 1 |u w | 2 v w dx dt. (6.1)
This quantity scales with ρΓ 3 0 /(16π 3 b 2 ) hence, for a typical large aircraft (b = 50 m, Γ 0 = 500 m 2 s -1 ), this would yield P 100 W m -1 (the integrated quantity is order one). The set-up being two-dimensional, control must be applied over an a priori infinite distance. In practice, the control area would be finite yet the questions of what length of application would be necessary is open. There would also be different constraints when considering departing or arriving aircraft. This would deserve a detailed analysis and especially require three-dimensional effects to be looked into.

One certainly interesting use of the proposed control would be to exploit the capacity to move the vortex normal to its axis to induce three-dimensional waves along the vortex. This could be obtained by modulating the control in the axial direction or even implementing only crenelated patches of control. Besides drastically decreasing the required power, this would help three-dimensional instabilities that naturally grow in the wake [START_REF] Williamson | Phenomena, dynamics and instabilities of vortex pairs[END_REF]. Indeed, the range of horizon times T = 0..10 contains the characteristic times of Crow and elliptic instabilities, respectively T cr 6 [START_REF] Williamson | Phenomena, dynamics and instabilities of vortex pairs[END_REF] and T el 4.9 (Widnall, Bliss & Tsai 1974;[START_REF] Leweke | Cooperative elliptic instability of a vortex pair[END_REF][START_REF] Sipp | Widnall instabilities in vortex pairs[END_REF]. Such an opportunity could be usefully investigated in future works with the question of the optimal wavelength to be applied in mind.

Other perspectives are opened up by this work. One interesting outlook concerns the control of the shape of the wall instead of that of the boundary condition for the velocity field. Using the same technique that has been applied in the present study, the optimization of the shape of the wall to achieve the increased lateral displacement of the vortex or more simply to target its reduction could be envisaged. As a first step toward such an approach the present results can be used to derive such wall deformations since (u w , v w ) = (∂x w /∂t, ∂y w /∂t) can relate the movement of the wall to the velocity imposed at its surface. However, this approach would not take into account the effect of surface deformation on the flow and would hence be suited to very small surface deformations. For large displacements, new developments would be required to incorporate the deformation of the numerical domain (mesh) in the procedure. Like the blowing/suction strategy, adaptive walls could be an issue in practice due to implementation complexity. A way to deal with this would be to consider a static shape to control the vortices, in a way similar to what has been achieved at Amsterdam airport to counter the noise of aircraft [START_REF] Sorvig | Quietly defend silence[END_REF]. The optimization procedure in this case would not require the use of the adjoint technique since the shape of the wall could be described by a small set of parameters. Another way could be the use of compliant walls. The subject has gained strong interest in recent years to deal with hydro and aerodynamic problems such as drag reduction [START_REF] Endo | Direct numerical simulation of turbulent flow over a compliant surface[END_REF], flow instabilities [START_REF] Yeo | The three-dimensional stability of boundary-layer flow over compliant walls[END_REF][START_REF] Davies | Instabilities in a plane channel flow between compliant walls[END_REF][START_REF] Hoepffner | Mechanisms of non-modal energy amplification in channel flow between compliant walls[END_REF] and turbulence [START_REF] Luhar | A framework for studying the effect of compliant surfaces on wall turbulence[END_REF]. In particular, [START_REF] Endo | Direct numerical simulation of turbulent flow over a compliant surface[END_REF] shows that passive compliant walls tend to reduce the intensity of vortex structures at the wall. where f and g are functions of space. By cancelling the first variation of the Lagrangian function (4.4) with respect to u(T), one obtains ∂L ∂u(T)

, δu(T) = ∂J ∂u(T)

, δu(T) -u + (T), δu(T) = 0, (A 2)

where δu(T) is an arbitrary vector. To compute the term (∂J /∂u(T)), δu(T) , one needs to compute (∂I m /∂u(T)), δu(T) where I m=0,1 is defined in (4.8) and (4.9). For the sake of clarity, let us rewrite I m with the following expression:

I m (u(T)) = V x m (∇ × u(T)) 2 dΩ. (A 3)
This allows us to compute the directional derivative of I m with respect to u(T) as follows:

∂I m ∂u(T) , δu(T) = lim →0 1 (I m (u + δu(T)) -I m (u)), (A 4) yielding ∂I m ∂u(T) , δu(T) = Ω 2x m ω(T) • (∇ × δu(T)) dΩ. ( A 5)
We can use the following formula ∇

• (a × b) = b • (∇ × a) -a • (∇ × b)
and integration over Ω using the Green-Ostrogradski theorem, which yields

∂I m ∂u(T) , δu(T) = Ω ∇ × (2x m ω) • δu(T) dΩ. ( A 6)
Remark that the previous equality holds true for any arbitrary vector δu(T), thus

∂I m ∂u(T) = ∇ × (2x m ω) = 2∇(x m ) × ω + 2x m ∇ × ω. (A 7)
Hence, one can evaluate the initialization of the adjoint variable in two dimensions

u + (T) = ∂J ∂u(T) = 2 I 0 (T) (x∇ × ω(T)e z + ∇x × ω(T)e z ) - 2I 1 (T) I 0 (T) 2 (∇ × ω(T)e z ), (A 8)
and find (4.7).

Appendix B. Convergence analysis

B.1. Optimization This appendix provides examples of the convergence of the algorithm described in § 4. The assessment of numerical convergence is evaluated using the following residual

r n = J n -J n-1 J n . (B 1)
Figure 18 plots r n as a function of iteration number and displays two typical convergence results: approximately 20 iterations are required to obtain a wellconverged optimal (r n 10 -3 ).

B.2. Influence of the interpolation time interval δt i

This appendix discusses the convergence as a function of δt i . One example is considered: (Re, T, l) = (1000, 5, 1). Convergence results are provided in table 2. The value δt i = 0.025 is chosen as a consequence of the satisfactory convergence it enables.

Appendix C. Temporal evolution of the cost function

This appendix details the derivation of relation (4.4). We consider the general case of I m = Ω A m dV with A m = x m ω 2 and m = 0, 1. The time derivative of I m is obtained as follows:

İm = d dt Ω A m dV = Ω ∂A m ∂t dV + ∂Ω A m (u • n) dS = Ω ∂A m ∂t dV + Ω ∇ • (A m u) dV. (C 1) The last term of the previous equation develops into A m ∇ • (u) + u • ∇(A m ) = u • ∇(A m ) = mx m-1 uω 2 + x m (u • ∇)ω 2 .
The vorticity equations in two dimensions are written as

∂ω ∂t + (u • ∇)ω = 1 Re ω. (C 2)
This is used to simplify (C 1) The term x m ω ω can be worked out using the identity ∇

İm = Ω 2ωx m ∂ω ∂t + (u • ∇)ω + mx m-1 uω 2 dV = Ω 2 Re x m ω ω + mx m-1 uω 2 dV. ( C
• ( f ∇(g)) = ∇( f ) • ∇(g) + f g which yields İm = 2 Re Ω ∇ • (x 2 ω∇ω) -∇(x m ω) • ∇ω + mx m-1 uω 2 dV. (C 4)
Then, the middle term of (C 3) can be developed and integrated. This gives

İm = 2 Re Ω ∇ • (x 2 ω∇ω) -x m ∇(ω) 2 + mx m-1 uω 2 dV. (C 5)
Eventually, using the Green-Ostrogradski theorem, we get

İm = 2 Re Ω -x m |∇ω| 2 dV + 1 Re ∂Ω x m (∇ω 2 ) • n dS + Ω mx m-1 uω 2 dV. ( C 6)
In figure 19 the time evolution of the cost functional along with its derivative for the uncontrolled case (Re, T, l) = (1000, 15, ∞) are shown. The derivative is obtained either by using (C 6) leading to J = ( İ1 I 0 -I 1 İ1 )/I 2 0 or by numerically differentiating J (t). The two approaches to the differentiation match well.

CHAPTER 2. A VORTEX PAIR IN GROUND EFFECT, DYNAMICS AND OPTIMAL CONTROL

Future studies could focus on the development three-dimensional instabilities of counter-rotating vortices impinging the ground after applying the optimal control strategy, however because of time constraints, they are out of the scope of this thesis.

A starting point would be the study of bi and tree-dimensional optimal perturbations of counter-rotating vortices in ground effect, as addressed in the next chapters.

CHAPTER 3

Optimal Perturbations of a Counter-Rotating

Vortex Pair in Ground Effect

The reader may refer to Chapter 2 in which the dynamics of a vortex pair in ground effect has been presented, along with its control by flow actuation at the ground. One of the main motivations for studying optimal perturbation of vortex pairs in ground effect is that these perturbations are expected to hasten the transition to turbulence, by triggering three-dimensional mechanisms. Thereby, reducing the danger of a vortex encounter.

Vortex Instabilities & Ground Effect

Late in the 19-th century, Kelvin (1880) analyzed wave motions in a column of uniform vorticity surrounded by irrotational motion (Rankine vortex model). It became a starting point for the study of various vortical flows. Many of them rely on the Linear Stability Analysis (LSA). LSA is a powerful tool used to study the asymptotic (t → +∞) dynamical behavior of a fluid system subjected to small perturbations. In the long term, the least stable mode is expected to dominate the evolution of the flow. However, many flow mechanisms lead to the development of strong perturbations from initially small amplitude perturbations over a finite time interval that cannot, therefore, be predicted by LSA. This process, known as transient growth of perturbation energy, relates to the non-normality of the Linearized Navier-Stokes operator [START_REF] Farrell | Optimal excitation of perturbations in viscous shear flow[END_REF]; [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF]). As a result of this non-normality, the eigenmodes of the homogeneous problem are non-orthogonal, making it possible to form an initial perturbation (from a combination of different eigenmodes) for which the time derivative of perturbation energy is positive, albeit the linear eigenmode are individually stable. During the transient phase of both stable and unstable flows liable to non-normality and short term disturbance amplification, the perturbations growing in amplitude may lead to further linear development (secondary instabilities) and non-linear interactions. This is able to drive the dynamics away from the initial flow state. The plane Couette flow is an example of LSA predicting stability for all Reynolds numbers while experiments show that turbulence is observed for relatively low Reynolds numbers (Re ∼ 300). Flow non-normality and transient growth mechanisms lie behind this surprising observation. The term optimal perturbations has been defined by [START_REF] Farrell | Optimal excitation of perturbations in viscous shear flow[END_REF] while studying the perturbation that would lead to the greatest growth over a finite period of time [0, T ]. Mathematically, the optimal perturbation is the perturbation of the base-flow that optimizes a given cost function. Kinetic energy is commonly used to quantify the growth of a perturbation to the flow. However other quantities can also be considered depending on the physics at play and mechanism to target. For instance, enstrophy is also an interesting choice as it targets the capacity of the flow to dissipate its momentum through turbulence. It is also particularly adapted to flows dominated by vortical motion.

In the present chapter we describe the transient growth potential in vortices interacting with a wall. As an introduction we first review some of the single and double vortex instabilities.

Most studies have focused on the Rankine model (uniform and circular patch of vorticity) both for its good repreentation of most vortex flows and also for its analytical simplicity, the Lamb-Oseen (L-O) vortex (gaussian vorticity profile) for its good agreement with experimental data and the Batchelor vortex (gaussian vorticity including axial velocity). Therefore we will go over the range of these different vortex models in the course of this review, first looking at the single vortex case and then the vortex pair. In a second step, we introduce the mathematical formalism for the transient growth analysis and lastly we present the results for the wall interacting vortices.

Stability and optimal perturbations of an isolated vortex

The optimal perturbation for an isolated vortex has been studied mainly within the linear framework [START_REF] Antkowiak | Transient energy growth for the lamboseen vortex[END_REF]; [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF]; [START_REF] Mao | Continuous spectra of the batchelor vortex[END_REF][START_REF] Hecht | New development in freefem++[END_REF]). The eigenmodes of the L-O vortex were determined by [START_REF] Fabre | Kelvin waves and the singular modes of the lamb-oseen vortex[END_REF] and are composed of singular modes and Kelvin waves. However, due to the non-normality of the Navier-Stokes operator applied to the L-O vortex, strong transient growth are known to be possible. The linear optimal perturbations of the L-O vortex were described by [START_REF] Antkowiak | Transient energy growth for the lamboseen vortex[END_REF][START_REF] Kolá Ř | Vortex identification: new requirements and limitations[END_REF] and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF]. Pertubation on an isolated axisym-

VORTEX INSTABILITIES & GROUND EFFECT

metric vortex can be decomposed into azimuthal modes of m periodicity. The axisymmetric m = 0 optimal perturbation takes the form of superimposed filaments of alternating positive and negative azimuthal velocity around the vortex core (see [START_REF] Antkowiak | On vortex rings around vortices: an optimal mechanism[END_REF]). These structures evolve into counter-rotating rings of azimuthal vorticity through a mechanism of anti-lift-up. This term is coined in reference to the process of lift-up, because it proceeds reversely to the process of the lift-up mechanism commonly observed in shear flows (see [START_REF] Ellingsen | Stability of linear flow[END_REF]; [START_REF] Landahl | Wave breakdown and turbulence[END_REF]).

The nonaxisymmetric optimal perturbations of the L-O vortex were described by [START_REF] Antkowiak | On vortex rings around vortices: an optimal mechanism[END_REF] for m = 1 and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] up to m = 4. In these cases, perturbation growth is achieved through the Orr mechanism and resonance between structures located outside the core and eignemodes of the vortex core. The shear mechanism is fully described by the [START_REF] Orr | The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. part ii: A viscous liquid[END_REF] energy equation [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF]):

dE dt = d t E = - V u • ∇U • u dV - 1 Re V ∇u : ∇u dV = P -D (3.1)
where D is the viscous dissipation rate and The second mechanism consists of a resonance phenomenon between the perturbation outside the vortex and core modes. Both [START_REF] Antkowiak | Transient energy growth for the lamboseen vortex[END_REF] and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] showed that the optimal perturbation tends to select a radial location for the perturbation that rotates at the frequency of the least stable mode. Radial velocity u

P = - V u • ∇U • u dV = -

CHAPTER 3. OPTIMAL PERTURBATIONS OF A COUNTER-ROTATING VORTEX PAIR IN GROUND EFFECT

is generated within the core through the advection of perturbation axial vorticity ω by the baseflow.

Linear Optimal Perturbation of a Vortex Pair

The stability of a vortex pair was studied in many earlier works [START_REF] Crow | Stability theory for a pair of trailing vortices[END_REF]; Widnall et al. (1974); [START_REF] Tsai | The stability of short waves on a straight vortex filament in a weak externally imposed strain field[END_REF]; Pierrehumbert (1980)). Vortices in a pair can become unstable to short-wave and long wave instabilities. These instabilities are described extensively in reviews by [START_REF] Widnall | The structure and dynamics of vortex filaments[END_REF] and Leweke et al. (2016). Asymptotic optimal perturbations of vortex pairs, in the absence of ground, was determined numerically by [START_REF] Brion | Optimal amplification of the crow instability[END_REF], for the wavelength of the Crow instability. Jugier ( 2016) investigated the optimal perturbation to the vortex pair for a pure planar flow, further accounting for stratification effects, along with [START_REF] Ortiz | Three-dimensional instabilities and optimal perturbations of a counter-rotating vortex pair in stratified flows[END_REF]. In the particular case of the Crow wavelength, the perturbation takes the form of a pair of opposite-signed vorticity sheets situated close to the plane separating the vortices. The evolution of the perturbation goes first through an advection toward the leading hyperbolic point of the flow, where the perturbation vorticity is then amplified and transfered to the vortex core through an induction mechanism with the velocity field of the base flow (Biot-Savart). An acceleration of the Crow instability is obtained through this process, which leads to a reduction of approximatey 2.5 normalized time units τ = 2πb 2 /Γ of the time required to reach the Crow instability, for Re = Γ/ν = 3600. The investigation performed by [START_REF] Jugier | Stabilité bidimensionnelle de modèles de sillage d'aéronefs[END_REF] of the linear optimal pertubations for a counterrotating pair of vortices in two dimensions provided results showing that the transient growth can lead to higher energy gain than the most unstable mode. [START_REF] Ortiz | Three-dimensional instabilities and optimal perturbations of a counter-rotating vortex pair in stratified flows[END_REF] investigated a similar bi-dimensional configuration with stratification. The latter is observed to significantly increase the instability as the two vortices are pushed together, resulting in an enhanced straining field, and hence a more rapid growth of the elliptic instability.

Instabilities in Ground Effect

In terms of linear stability, ground proximity has a notable influence on the development of long and short-wave instabilities in vortex pairs. The modification to the long-wavelength instability was investigated by [START_REF] Crow | Stability theory for a pair of trailing vortices[END_REF]; [START_REF] Asselin | Influence of a wall on the three-dimensional dynamics of a vortex pair[END_REF]. Several regimes are observed depending on the height of the vortex system above the ground. If the initial height of the vortex pair above the ground plane is sufficiently large, then the Crow instability and the transformation into vortex rings will occur before surface interaction. Howerver, at relatively smaller initial heights, the growth of the Crow instability is found to be inhibited by the presence of the ground. What occurs is that while the long-wave instability develops, regions of the vortex pair closest to the wall interact with the boundary layer first leading to an increase in local pressure, and then pushing the flow axially away from these regions. This dynamics tends to modify the long wave motion of the vortex cores which supports the long wave Crow instability of the pair, and results in a weaker growth.

Following this influence of vortex height, experiments by [START_REF] Asselin | Influence of a wall on the three-dimensional dynamics of a vortex pair[END_REF] further detail that there exists three regimes of vortex-surface interaction. If the ratio h 0 /b 0 (h 0 is the initial height and b 0 is the initial separation) is above 9 then Crow instability prevails (note that according to [START_REF] Spalart | Airplane trailing vortices[END_REF] the typical time for the development of Crow instability is 5 (vortex pair time units), which is coherent with the value 9 corresponding to this first regime). If the vortex pair is generated below a critical height (3 < h 0 /b 0 < 9), the wall tends to inhibit the Crow instability and two other modes, of shorter wavelength, are observed. For smaller intermediate heights (3 < h 0 /b 0 < 6), the primary vortices interact with the secondary vortices to form vertically oriented rings (see figure 3.2). These 'vertical rings' are the typical structures observed when three-dimensional deformations of the flow are accounted for in the wall interaction, on top of the rebound phenomenon that we could observe in the section dedicated to the two-dimensional situation. For completeness, it must be noted that the Crow instability also develops before wall interaction for intermediate heights (6 < h 0 /b 0 < 9). However, as the perturbation grows larger, the reduction of circulation at the trough (were the vortex tubes are closer together) leads to a higher pressure compared with the peaks of the vortices, and triggers the collapse of the primary vortices as a result of the strong axial flows from the troughs towards the peaks. This yields "horizontal rings modes". The structures found by [START_REF] Asselin | Influence of a wall on the three-dimensional dynamics of a vortex pair[END_REF] are comparable to those observed in vortex-ring impingement at a wall, studied experimentally by [START_REF] Lim | An experimental study of a vortex ring interacting with an inclined wall[END_REF] and [START_REF] Cheng | Numerical study of a vortex ring impacting a flat wall[END_REF].

Similar results are found by [START_REF] Dehtyriov | Optimal growth of counter-rotating vortex pairs interacting with walls[END_REF] while studying the transient growth of a counter-rotating vortex pair impinging a wall with the particular setting h 0 /b 0 = 6. These authors show that the linear growth of the Crow instability is inhibited by the wall. In addition, the evolution results in the suppression of the secondary vortices and a strong reduction of the rebound phenomenon. Furthermore, the linear short-wave instability outgrows the longwaves modes, thus clearly underlining the importance of the elliptic mode in the presence of a plane surface. At last, experiments of a bi-dimensional vortex pair impinging on a wavy wall, carried by [START_REF] Morris | Impingement of a counter-rotating vortex pair on a wavy wall[END_REF], show the generation of rebounding vortex rings and the acceleration of the decay of the primary vortex pair.

For the rest of this thesis, the single vortex and vortex pairs are modelled, as for the initial condition, by a single L-O vortex or a superposition of L-O vortices. Recall that the radial velocity profile of the L-O vortex is given by

V θ (r, t) = Γ 2πr 1 -exp - r 2 r 2 c (3.3)
where Γ and r c (t) = 4νt + r 2 c (0) respectively denote the vortex circulation and dispersion radius.

METHOD

The following sections are dedicated to the study of a spatially periodic (in the axial direction) linear optimal perturbations of a counter-vortex pair impinging on a surface plane and initially located at h 0 /b 0 = 2.0 with an initial vortex radius a = 0.2b 0 . We first discuss the pure two-dimensional case and then the three-dimensional case. We detail the mathematical formulation of the problem in section 3.2.1. The presentation and discussion of the results are given in section 3.3.2. Finally section 3.4 shows the results of the nonlinear evolution of the baseflow disturbed by a selection of linear optimal perturbations with finite initial amplitude.

.

Method

Optimization approach

The flow is decomposed into a base flow U and a perturbation u. The amplitude of the latter is considered to be small in comparison with the base flow such that one can write u tot = U + εu and p tot = P + εp where ε 1. One can rewrite the Navier-Stokes (N-S) equations as follows :

∇ • (U + εu) = 0 ∂ t (U + εu) + ((U + εu) • ∇) (U + εu) = -∇ (P + εp) + 1 Re ∇ 2 (U + εu) (3.4)
By introducing this decomposition, one retrieves the Navier-Stokes equations applied to the base flow at order ε 0 :

∇ • U = 0 ∂ t U + (U • ∇) U = -∇P + 1 Re ∇ 2 U U (y = 0) = 0 (3.5)
and the linearised evolution of the perturbation at order ε 1 when the second order term ε 2 (u • ∇) u sufficiently small (that is when u has a small amplitude). This yields the first order linearized N-S equations:

∇ • u = 0 ∂ t u + (U • ∇) u + (u • ∇) U = -∇p + 1 Re ∇ 2 u u(y = 0) = 0 (3.6)

METHOD

For the methodology and the detailed derivation of the equations, the reader may refer to Chapter 1 and papers by [START_REF] Farrell | Optimal excitation of perturbations in viscous shear flow[END_REF] and [START_REF] Corbett | Optimal perturbations for boundary layers subject to stream-wise pressure gradient[END_REF]. Reaching the optimal perturbation is achieved through the iterative procedure described in figure 1.2. The unit energy constraint is enforced by the geometric update technique proposed by [START_REF] Douglas | On gradient adaptation with unit-norm constraints[END_REF] to iterate over q 0 . First, the component of the Lagrangian function that is normal to q k 0 is scaled so as to satisfy the energy constraint, as follows :

N k = E(0) 1/2 ∇ q k 0 L ⊥ ∇ q k 0 L ⊥ (3.13)
where k is the current iteration and ⊥ denotes the normal component. The update is thus

q k+1 0 = q k 0 cos(α k ) + N k sin(α k ) (3.14)
where α k is the step size with initial value α 0 = 1. The line search is achieved by dividing α by 2 when the gain doesn't increase. The procedure stops when the residual

r n = G k -G k-1 G k-1 (3.15)
is inferior to 10 -6 (value allowing gain convergence inferior to 0.01%). Finally, the first initial guess is set as random white noise. The whole procedure is sketched in figure 3.3 q 0,k q 0 α 0 = 1 Direct equations

Adjoint equations

Compatibility equation

q k+1 0 = q k 0 cos(α k ) + N k sin(α k ) q k+1 0 = q k 0 cos( α k 2 ) + N k sin( α k 2 ) G k+1 ? ≥ G k r k ≤ε NO NO yes : exit

Definition of the computational configuration

The computations are carried out at Re = 500, in the laminar regime. The Reynolds number effect is small in this range of rather large Reynolds number and laminar flow conditions, and therefore only this value is considered.

CHAPTER 3. OPTIMAL PERTURBATIONS OF A COUNTER-ROTATING VORTEX PAIR IN GROUND EFFECT

Due to the simplicity of the flow configuration that is considered, the horizon time and the axial wavelength are retained as the two main parameters of the study. The range of horizon time is T ∈ [2, 5] and the axial wavenumber varies from 0 to 10.

The initial height of the vortex pair is constant throughout the study. No parametric study over this initial height is considered given that, as reviewed above, its effect is mainly to delay the occurence of the influence of the ground. Here the vortex pair initial height, at h 0 /b 0 ≥ 2, is lower than the one considered in [START_REF] Asselin | Influence of a wall on the three-dimensional dynamics of a vortex pair[END_REF]; [START_REF] Dehtyriov | Optimal growth of counter-rotating vortex pairs interacting with walls[END_REF] where h 0 /b 0 ≥ 4, which leads to a shorter time for the vortex pair to experience wall influence. This choice allows for a consistent physics while reducing the computational cost of the simulation and optimization.

The baseflow being symmetric with respect to the central plane x = 0, one may show that two orthogonal subsets of perturbations can be considered : symmetric and anti-symmetric perturbations, respectively denoted SYM and ASYM throughout the rest of the manuscript. Symmetric perturbations are the one verifying mirror symmetry with respect to x = 0, that is:

u(-x, y, z, t) = -u(x, y, z, t) v(-x, y, z, t) = v(x, y, z, t)
w(-x, y, z, t) = w(x, y, z, t) (3.16) conversely anti-symmetric perturbations are the one that verify :

u(-x, y, z, t) = u(x, y, z, t) v(-x, y, z, t) = -v(x, y, z, t) w(-x, y, z, t) = -w(x, y, z, t) (3.17)
For the upcoming simulations, the simulation domain is either two-dimensional or threedimensional.

2D domain:

The simulation domain consists in a rectangle of size 20b × 10b composed of N x ×N y = 222×111 spectral elements with polynomial order P = 5 (which yields 6 gridpoints per element in each direction).

3D domain:

The simulation domains is the two-dimensional plane extruded in the zdimension. Hence it becomes a box of size 20b × 10b × λ. The associated mesh is made of N x × N y × N z = 222 × 111 × 3 = 73926 spectral elements. The polynomial order is P = 5. Each spectral element is composed of L x × L y × L z = (5 + 1) 3 = 216 grid points. All in all, the mesh comprises approximately 6.3 × 10 7 degrees of freedom.

RESULTS

A mesh convergence analysis has been performed, whose results are provided in Appendix A. The chosen grid size allows for well converged results and therefore, for the remainder, the results relie on a polynomial order P = 5. Due to the homogeneity of the base flow in the axial direction perturbations can be decomposed into axially periodic components. A Fourier modal decomposition is therefore implemented

[u, p, ũ, p] = +∞ -∞ [u, p, ũ, p]e ikz dk (3.18)
with the simulation domain axial dimension being equal to the spatial wavelength λ. Due to the linearity of the problem the Fourier modes with different wavenumbers k = 2π/λ are treated separately. In practice, to enforce a unique wavelength λ = 2π/k and prevent the development of unwanted disturbance walength compatible with the computational domain (harmonics of the principal wavelength specified by the domain axial length), a projection onto the space (cos(kz), sin(kz)) is executed regularly during the time integration of the direct and adjoint equations. The perturbation equation being linear, this projection operation removes undesirable wavelengths arising from numerical errors, keeping the selected wavelength untouched.

Wrapping up the previous details, in the following we consider the set of parameters given by:

• wavenumber k ∈ [0, 10]

• horizon time T ∈ 2, 5

• symmetry of the perturbation.

to conduct the analysis of the ground interaction dynamics.

Results

Two-dimensional case

This subsection describes and analyses the two-dimensional optimal perturbations.

Figure 3.4 shows the optimal gain as a function of time horizon T , for the symmetry and anti-symmetric sets. The figure suggests that the optimal gain is an increasing function of the horizon time regardless of the symmetry. Furthermore it is found that symmetric perturbations dominate for lower values (T < 3.7) of horizon time whereas anti-symmetric perturbations dominate for larger horizon times, although from a global point of view the range of amplification remains quite similar for the two symmetries. This relatively low incidence of flow symmetry on transient amplification growth puts forward the importance of the ground effect which is not determined by symmetry. In order to understand the mechanisms that lead to transient growth, one may refer to equation (3.1) which writes the instantaneous growth of kinetic energy. In this equation the second term represents viscous dissipation and the first term corresponds to the production of kinetic energy. The latter occurs by the alignment of the perturbations with the proper directions of the strain rate tensor of the flow (baseflow). As is evident from equation (3.2), the region of the flow where the shear production per unit volume is negative contributes to the growth of the perturbations.

To understand where, in the vortex pair flow interacting with the wall, this production region lies, we show in figure (3.5) the norm tr( ∇U • ∇U ) 1/2 of the strain rate tensor at times t = 0, 1.25, 2.5, 3.75.

At initial time, it is found that the regions contributing to the growth are predominantly located inside the Kelvin oval of the vortex pair (fig. 3.5(a)). However, when time increases and the vortices approach the plane of the ground, the strength of the strain rate tensor shifts to the boundary layer that concurrently forms at the wall (fig. 3.5(b-c)). When the primary vortices rebound (fig. 3.5(d)), the intensity of the strain rate tensor partially shifts back to the vortex cores. This can be interpreted as a consequence of the relaxation of the effect of the ground due to the increased distance of the vortices above it (the shearing effect due to the vortices reduces). The consequence of this is shown in figure 3.6 where the initial vorticity of ASYM and SYM optimal perturbations is plotted for various values of the horizon time T . When the horizon time is low (T 2), the perturbation is principally located in the vicinity of the Kelvin oval of the primary vortices. As horizon time grows (T = 3 → 5), the spatial distribution of the 3.3. RESULTS Figure 3.5: Distribution of the norm of the strain rate tensor of the baseflow for time t = 0, 1.25, 2.5 and 3.75 respectively corresponding to subfigures a), b) c) and d). Green solid lines correspond to the iso-contours of vorticity (ω z /ω 0,max = ±0.1, ±0.2) of the primary vortices. initial perturbation yields stronger structures in the vicinity of the wall and in the vortex cores. For the lower values of the horizon time T , the optimal perturbation is principally Figure 3.6: Initial vorticity distribution of ASYM (top row) and SYM (bottom row) optimal perturbation computed for different horizon times 2, 3, 4 and 5 (from left to right). Increasing horizon time implies a shift of the optimal perturbation shape towards the ground, pointing out the influence of horizon time on the boundary layer effect. PAIR IN GROUND EFFECT located inside the Kelvin oval, at the center plane x = 0 and near the leading and trailing hyperbolic points. These region are identified by the norm of the strain rate in figure 3.5(a) although they do not correspond to the zones where the maximum production rate occurs. In fact lower optimal time yields a better match to the instantaneous optimal production zone (not shown). The interesting feature is that as horizon time increases, the optimal perturbation becomes stronger near the wall, indicating increased sensitivity to perturbation there.

The optimal perturbation at horizon time is shown in 3.7, in terms of axial vorticity. The final shape of the optimal perturbation for T = 2 is a displacement mode located inside the vortex cores, see figure 3.7(d). For T = 5, the baseflow boundary layer separates to form secondary vortices. From figures 3.7f) and g), it appears that the perturbation grows in the secondary vortices where it forms two displacement modes. One displacement mode acts on the primary vortices (dotted lines in fig. 3.7-f), the second on the secondary vortices.

Figure 3.7: Evolution of the vorticity of the symmetric optimal perturbations for horizon time T = 2 (top) and T = 5 (bottom) at times T /4, T /2, 3T /4 and T . As the baseflow boundary layer detaches, the optimal perturbation takes the form of two displacement modes: one acting on the primary vortex, the other on the secondary. Solid lines correspond to the iso-contours of vorticity (ω z /ω 0,max = ±0.1, ±0.2) of the primary vortices.

Concluding on two-dimensional optimal perturbation, it must be reminded that in the absence of the ground the linear optimal perturbations of a vortex pair consists in vortex filaments in a spiral arrangement around the vortex core and that its evolved state has a quadripolar structure near the vortex center [START_REF] Navrose | Optimal perturbation for two-dimensional vortex systems: route to non-axisymmetric state[END_REF]). In addition they have shown that the mechanism for the largest gain in perturbation energy for a vortex pair is the same as that of an isolated vortex [START_REF] Navrose | Optimal perturbation for two-dimensional vortex systems: route to non-axisymmetric state[END_REF], figure 20). In the presence of ground, we find here that the linear optimal perturbation is a displacement mode acting on the primary vortices (see, figure 3.7) for lower horizon times (T = 2) and that for higher 3.3. RESULTS values of the horizon time T , the optimal perturbation grows into displacement modes in the secondary vortices as well as in the primary vortices. In this case, the flow is sensible about the initial vortex pair footprint but also in the region near the ground. for the two symmetries. Overall, it appears that long wavelength instabilities prevail for short horizon time, replaced by short wavelengths for large horizon times. These results agree with those obtained by [START_REF] Asselin | Influence of a wall on the three-dimensional dynamics of a vortex pair[END_REF]. That is, as explained in the introduction, the presence of the wall tends to inhibit the longwave instabilities and to promotes shorter wavelengths. We find two characteristic wavenumbers, k = 0.8 and k = 6.0. Wavenumber k = 0.8 corresponds to the Crow wavelength. Wavenumber k = 6.0 yields the most critical optimal perturbation at larger horizon times, regardless of the symmetry.

Three-dimensions case

For what follows, we focus on symmetric optimal perturbations as they yield higher gains, more particularly at T = 5. The two wavenumbers that will be investigated are k = 0.8 and k = 6. Figure 3.9 shows iso-contours of vorticity magnitude for these specific wavenumbers (k = 0.8 in the left panel and k = 6 in the right panel). One can see that the effect of the shape of the optimal perturbation ressembles the two dimensional optimals previously shown in figure 3.6, modulated in the z-direction. In addition low values of the contours ||ω|| 5 are more present near the ground for k = 6 compared to k = 0.8 certainly due to the fact that shorter axial wavelength are more compatible with the small dimension of the boundary layer. This would indicate that a possible explanation for the domination of the shorter wavelengths is PAIR IN GROUND EFFECT that they may take better advantage of the dynamics of the boundary layer and separated flow that forms during its evolution. The sheared flow region at the frontier of the separated flow area offers strong potential for disturbance growth, both for 2D (Orr mechanism) and 3D situations. During the first stage of the primary vortex descent (t = 0 → t = 2), the effect of the wall is weak and the perturbations grow into longwave deformations. As the boundary layer of the baseflow separates at the ground the perturbations appear to grow strongly in this region of the boundary layer (t = 2 → t = 3). At later times the perturbation then acts principally on the secondary vortices (t = 4 → t = 5). and k = 6 (bottom), SYM perturbations, horizon time T = 5. The optimal perturbation initially acts as a displacement mode inside the primary vortex cores and when boundary layer separates, the perturbation drastically grows in the boundary layer. When the boundary layer has detached, displacement modes grow into the secondary vortices giving rise to structures similar to those that were observed in two dimensions. Vorticity magnitude ranges between 10 and 900 (top). and between 10 and 3 × 10 4 (bottom).

Nonlinear response to the linear optimal perturbation

The potential for anticipated destruction of a counter-rotating vortex pair using the linear optimal perturbation of the vortex pair in ground effect is assessed. Direct numerical simulation is used to study the development of instabilities and the subsequent evolution of the flow up to 25 characteristic times (t = 25). To achieve this goal, the previously computed optimal perturbations are superimposed to the unperturbed flow at initial time such that the initial velocity field for each DNS described in this section is expressed as u ε (0) = U (0)+ε||U (0)||u 0

NONLINEAR RESPONSE TO THE LINEAR OPTIMAL PERTURBATION

that the mesh contains only three elements in the z-direction, each composed of 6 gridpoints. For computational cost purposes, as the simulations would need mesh refinement in the axial direction (thus drastic increase in gridpoints would ensue), long-wave dynamics is not investigated here, and interest is devoted to shorter wavelengths, as they are the linearly dominant ones.

Figure 3.12: Time sequence of vorticity magnitude of the nonlinear response to the symmetric long-wave optimal perturbation computed for T = 5 and k = 0.8. The early breakdown of structures into short waves can not be captured by the mesh resolution in the axial direction which is made for the long-wave dynamics.

Short wave optimal perturbations: k = 6, ε = 10 -3 Figure 3.13 shows the evolution of the symmetric optimal perturbation computed for T = 5 and k = 6 superimposed to the baseflow with an initial amplitude of ε = 10 -3 . The first effects of the perturbations are felt early during the descent of the vortex, around t = 2 where we can observe the deformation of the boundary layer as the secondary vortices form. The perturbations grows inside the secondary vortices, causing the secondary vortices to undergo strong stretching at midplane z = λ/2 (cf. time t = 3 → 5). This stretching causes the vortices to link in the middle plane x = 0 until the formation of a vertical ring (cf. time t = 5 → 6). This ring stretches outwards until the identified iso-levels break (cf. t = 7 → 10) and provoke a slow but steady decay of the primary vortices. This decay can be inferred from the decreasing size of the considered vorticity iso-level. Short wave optimal perturbations: k = 6, ε = 10 -2 Figure 3.14 shows the evolution of the symmetric optimal perturbation computed for T = 5 and k = 6 superimposed to the baseflow with an initial amplitude of ε = 10 -2 . In the early stages, the flow ressembles the U-shaped vortical structure of the linear optimal perturbation above the boundary layer at time t = 1, see figure 3.10b. These U-shaped structures rotate around the baseflow (t = 2) and link with the secondary vortices while being stretched (t = 3 → 4). The secondary vortices slowly lose strength and get advected away from the ground and the primary vortices (cf. t = 5 → 8). The primary vortices decay slowly but faster than the case ε = 10 -3 , for t > 8.

Short wave optimal perturbations: k = 6, ε = 10 -1

The structures of the perturbation in the early stage ressemble the ones in the case ε = 10 -2 but stronger. The U-shaped structures link with the boundary layer sooner, at times t = 0.5 → 1 and we also notice that the rebound phenomenon is reduced. The primary vortices undergo strong deformations in the axial direction (cf. t = 2 → 5) and decay rapidly thereafter (t = 6 → 8). The vortex system is almost completely suppressed beyond t = 10.

Kinetic energy: k = 6

The temporal evolution of the volume integrated kinetic energy writes as where the volume Ω is the entire computational domain.

dE dt = -2µ Ω d : d dxdydz (3.19)
Figure 3.16 shows the time evolution of the kinetic energy that results from the unpertur-PAIR IN GROUND EFFECT bated and perturbated simulations. In the first case (no perturbations) the plot shows how the kinetic energy first decays at a rather steady rate, with a marked increased after about 17 time units, finally reaching 20% of its initial value after 25 time units. The initial decrease of the kinetic energy results from the work of the dissipative term. The increased reduction of the kinetic energy after 17 time units seems to correspond to the appearance of small scale structures, which probably arise due to the discretisation errors of the computational solver. These structures, precursor of the turbulent state, reinforce the dissipation process by intensifying the deformations of the flow. From the reference state, increasing the amplitude of the initial amplitude of perturbation hastens the occurence of the vortex decay, with a breakdown time reduced to t = 5 for ε = 0.001 and 0.01 and almost immediate for ε = 0.1. In the evolution of the flow with finite values of ε (see figures 3.13 3.14 3.15) it seems that the anticipated breakdown relates to the first appearance of the deformations of the secondary flow about the primary vortices, which only after transfer to them. More precisely, the deformations start in the boundary layer and continue then in the rolling-up vortex sheet that forms the secondary vortices. This stresses the important role of the boundary layer flow in the evolution of the interacting vortices. Figure 3.17 and 3.18 show the evolution of the circulation present in sections of the flow at two axial positions, z = λ/2 and λ/4. Only a half of the section and the axial vorticty is considered for the integration of the circulation. No distinction is made of the primary vortices thus these plots provide the total circulation of the half flow. The surprising result is that there is no overall effect of on the evolution of the total circulation, all evolutions showing a rapid decrease of the circulation a the two section considered. We eventually look at the evolution of enstrophy in figure 3.19 in the section of the flow. Enstrophy is a good indicator of the appearance of small vortical structures in the flow and viscous dissipation. The figure shows that the simulations made at the various values of yield a peak of enstrophy at initial time. In the reference test case = 0 the first enstrophy peak corresponds to the rebound and formation of secondary vorticity from the formation of the boundary layer and its detachment. In the simulations with finite values, this first peak is complement by additional formation of smaller scale structures due to the initial perturbation. In the = 0.001 case this causes a second peak right after the rebound peak while in the case = 0.01 the first and second peaks occur together, generating a large and higher single peak. In the most perturbed case = 0.1 the rebound peak is anticipated by the strong initial deformations due to the initial perturbation and the peak due to rebound is only secondary in the overall dynamics. The enstrophy levels in this last case reduce strongly and become negligible after about 10 time units, which corresponds to the disappearance of the flow (see figure 3.15). In this chapter we have investigated the potential for vortex decay in the vicinity of the ground by perturbating the system at initial time using linear optimal perturbations. The vortices show to be principally affected by displacement modes for both two and three-dimensional vortices. The gain in kinetic energy was much higher in the case of the three-dimensional perturbations. The ground effect promotes short wavelengths at short times and then small scale perturbations arise in the flow. The separation of the boundary layer provides great potential for destabilization. Non-linear DNS of the total flow show noticeable vortex decay as compared to the unperturbed flow.

On p-norm Optimal Perturbations

Variational formulation based on a direct-adjoint optimization mostly relies on the kinetic perturbation density defined as

E(t) = V 1 2 u(x, t) 2 dV
which corresponds to the common L 2 norm of the perturbation field. Using kinetic energy-based techniques often yields widespread structures (cf. chapter [3]). Taking a more local approach to identify the zones of optimal energy growth is of great interest, especially for vortex systems. To tackle the aforementioned application issue, which is also of fundamental interest, one may consider using the following p-norm,

E p (t) = V e p (x, t)dV 1/p
, where e(x, t) = u(x, t) 2 2 introduced by [START_REF] Foures | Localization of flow structures using ∞-norm optimization[END_REF].

Note that in the present section, the initial kinetic energy is written E(0) instead of E 0 to prevent any misunderstanding with the p index.

The vast majority of studies on stability are based on a measure of the kinetic energy density defined as e(x, t) = (u 2 + v 2 + w 2 )/2. A common choice to quantify the growth of perturbations is the L 2 norm of e defined as

E 1 (t) = V e(x, t)dV
where V is the fluid domain. Using such energy based approaches yield global structures that are widespread throughout the entire fluid domain as no spatial restrictions are explicitly enforced. With the growing interest driven by industrial flows (pipe, junction, vortex flows, etc.), locating the regions responsible for optimal energy growth is central. Recently, constrained optimization of the p-norm of the disturbance energy defined as:

E p (t) = V e p (x, t)dV 1/p (4.1)
has been introduced by [START_REF] Foures | Localization of flow structures using ∞-norm optimization[END_REF] with the aim of promoting the localization of the optimal perturbation. The methodology was applied to the case of a simple two-dimensional channel flow. Locating and identifying 'hotspots' for energy growth can be achieved by choosing E p for the optimization process. We show below how, for a finite dimension space, the p-norm tends to the maximum value as p increases, that is, considering a vector

x = [x 1 , x 2 , . . . , x n ] lim p→∞ x p = x ∞ = max 1≤i≤n {x i } (4.2)
The heuristic proof is obtained by bounding the p-norm as follows

max |x i | ≤ (max |x i | p ) 1/p ≤ n i=1 |x i | p 1/p ≤ n i=1 max |x i | p 1/p = n 1/p max |x i |
Taking the limit for large p yields the desired result. A similar proof exists for continuous functions but is not developped here. This shows that the p-norm targets the part of the perturbation with maximum possible growth. Such a localized regions is coined as the 'hottest spot'.

The reader can refer to section (3.2.1) in chapter [3] for the derivation of the variational framework when p = 1 (that is, L 2 ). When it comes to arbitrary values of p, there are two minor differences for finding p-norm optimal perturbations as opposed to p = 1:

• The cost function transforms into G(T ) = E p (T )/E 1 (0) where E 1 (0) = E(0) is the kinetic energy.

• The initialization of the adjoint field, or the so-called compatibility equation, becomes:

ũ(x, T ) = e(x, T ) p-1 E p (T ) p-1 u(, x, T ) (4.3)
where e(x, t) = u(x, t) 2 /2 denotes the perturbation kinetic energy density. Note that one retrieves the result derived from Chapter [3] with p = 1.

Equation 4.3 is the key localization step in the procedure. If we consider the limit of large p the adjoint initial condition approaches a Dirac function centered on the maximum values of the direct energy density field, see figure 4.1. When advected backwards in time, this will therefore lead to spatially localized gradient making the procedure converge toward a localized initial perturbation. Hence, depending on the initialization of the optimal perturbation, one can expect to have multiple optmimal solutions. For p = 1, [START_REF] Foures | Localization of flow structures using ∞-norm optimization[END_REF] have recovered traditional energy-based stability analysis results where the initial optimal perturbation is an array of vortices aligned against the mean flow shear (see Foures et al., 2013, figures 2b, 3 and 7a). For large values of p 1 they have observed 'hotspots' where significant energy growth can be expected. Two different regions of energy growth have arisen: the center of the channel (see Foures et al., 2013, figure 7b) and at the walls (see Foures et al., 2013, figure 7c) of the domain each obtained at T = T opt = 10.4 (their characteristic time is based on the average flow velocity in the channel and channel width). For either solution, a saturation of E ∞ is observed, that occurs much before the limit p → ∞ is reached. For instance it is observed for p = 50 when the Reynolds number is equal to 4000. The reason for this is that, increasing p reduces the localization length scale (see figure 4.1). There is a competition between this length scale and that of diffusion. Specifically, when the localization length scale becomes smaller than that of diffusion, the effect of p disappears, causing the observed saturation. Therefore the saturation effect of p depends upon the value of the Reynolds number. In a subsequent work, [START_REF] Farano | Subcritical transition scenarios via linear and nonlinear localized optimal perturbations in plane poiseuille flow[END_REF] also took a look into the potential of p-norm optimization, specifically considering:

• p-normed objective functions within the linear framework

• 1-normed objective function within the nonlinear framework in order to analyze the structure of localized optimal perturbations in a three-dimensional plane Poiseuille flow for subcritical conditions. They found that 1-normed optimal perturbations can be grouped in four families depending on the amplitude of the initial perturbation kinetic energy and horizon time (at finite horizon time, the energy of perturbation can become sufficient to trigger non-linear effects):

• Linear (LOP), for low values (E(0) ≤ 3 • 10 -7 ) of the initial energy for all horizon times.

• Weakly NonLinear (WNLOP) for horizon times T ≤ 20 and E(0) ≤ 1.5 • 10 -6 .

• Highly NonLinear (HNLOP) for horizon times T ≥ 20 and E(0) ≥ 4 • 10 -7 .

• hairpin-like optimal perturbations obtained using large values E(0).

The LOP is composed of arrays of vortices parallel to the streamwise direction. The WNLOP are composed of alternated vortices, inclined with respect to the streamwise direction. The HNLOP are highly localized, do not have a preferred direction in space and do not show any symmetry. The reader may refer to figure 2 of [START_REF] Farano | Subcritical transition scenarios via linear and nonlinear localized optimal perturbations in plane poiseuille flow[END_REF] for further description of the optimal perturbations. The progressive localization of the nonlinear optimal perturbations has been achieved by progressively increasing E(0) and T . For linear p-normed optimal perturbations, they found that p = 50 is the value that maximizes E p . As opposed to [START_REF] Foures | Localization of flow structures using ∞-norm optimization[END_REF], who have found multiple solutions (center and wall) depending on the initialization of the optimal perturbations, the linear p-norm optimal perturbations led to one type of solution only that varies slightly with p. The linear p-norm optimal perturbations are characterized by vortices inclined upstream, (Farano et al., 2016, figure 6). Though the nonlinear perturbation outperform the linear p-norm optimal perturbations in terms of hastening the transition to turbulence when used in the initialization of direct numerical simulations, one conclusion of the study is that linear p-norm solutions provide a good compromise between efficiency and computation cost. In summary, the physical mechanisms associated with energy growth often rely on large structures able to extract energy from the baseflow. In L 2 norm the initial perturbation has no specific restrictions for its spatial distribution. In L p norm with large enough p, the constraint for thinner localization targets regions of the flow with the most potential for growth, yet with overall lower perturbation growth. As such p-norm techniques can be very useful for designing experimental setups, in which perturbations can be injected only locally using control devices, either due to technical constraints, accessibility or cost of control devices, which necessarily limit the distribution of such devices. Also control devices generally operates at or near walls, with limited reach within the flow, making much more attractive theoretical optimals localized in the wall vicinity. Therefore constraining more localized structures and choosing those occuring near walls help to achieve the realization of the control strategies, 4.1. APPLICATION TO THE ISOLATED LAMB-OSEEN VORTEX although with sub-optimality compared to less restrained optimal that would not have to deal with such realization constraints.

We use the same algorithm as chapter [3] with the modification due to the definition of the cost function and compatibility equation given above. We apply the algorithm for an isolated vortex and for the case of the flow past a finite wing. In the isolated vortex case, because the problem is axisymmetric, we further use a one dimensional optimization tool which discretizes the perturbation equations in radial coordinates, with the p-norm constraints, see Chapter [1].

Application to the Isolated Lamb-Oseen Vortex

In this part we optimize E p for an isolated Lamb-Oseen vortex whose velocity profile in cylindrical coordinates is given by U = 0, Γ 2πr 1exp(-r 2 /a 2 ) , 0 . The radius of the vortex a is set equal to 1. The circulation based Reynolds number Re = Γ/2πν = 5000. The horizon times considered are 2T * , 5T * and 10T * where T * = 4π 2 a 2 /Γ is the characteristic time of rotation of the isolated vortex. Throughout the chapter the baseflow is considered to be frozen, meaning that the vortex does not undergo viscous diffusion. This approximation can be applied when the time evolution of the perturbative part of the flow given by T * is small compared to the viscous time scale T ν = 2πa 2 /ν. The ratio of the two time scales T * /T ν is equal to the inverse of the Reynolds number and is therefore sufficiently small.

As a consequence of the axysymmetric base flow, perturbations can be decomposed as: u = û(r, t) exp(mθ) with m the azimuthal wavenumber. Modes of different m are independant as the equations are linear. We further consider planar perturbations with no axial wavelength.

4.1.1 p = 1 case (L 2 norm)

In order to validate our p-norm optimization tool, we refer to figure 4.2 to ensure the correct retrieval of values when p = 1. The latter provides the linear optimal gains for varying horizon time for m = 1, 2, 3 at Re = 5000 as determined by [START_REF] Bisanti | Linear and nonlinear optimal perturbation analysis of vortices in incompressible flow[END_REF], also retrieved by [START_REF] Johnson | Nonlinear dynamics of wake vortices[END_REF]. As noted by [START_REF] Antkowiak | Transient energy growth for the lamboseen vortex[END_REF] the growth of the m = 1 displacement mode increases linearly with the horizon time at large T . This linear behaviour is preceded by a stage where modes of higher azimuthal wavenumber m = 2 and m = 3 prevail. As presented in chapter 3, the inviscid linear growth described by [START_REF] Antkowiak | Transient energy growth for the lamboseen vortex[END_REF] and [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] can be explained upon writing the perturbation kinetic energy As the strain S = r∂ r (V /r) is negative in a LO vortex, kinetic energy production occurs when uv > 0 (positive Reynolds stress). Hence, the linear optimal perturbations of an isolated vortex takes the form of spiral as displayed in figure 4.3. For short times, this inviscid mechanism is equivalent of the the Orr mechanisms in plane shear flows. [START_REF] Antkowiak | Transient energy growth for the lamboseen vortex[END_REF]; [START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] have shown that the optimal perturbations select a radial location for the initial perturbation that would progressively induce a core mode within the vortex through a resonance-driven mechanism. By deriving the linearized vorticity equation for perturbation:

∂ω ∂t + V r ∂ω ∂θ + u ∂Ω ∂r = 0 (4.5)
one shows that advection of the perturbation vorticity by the baseflow induces radial velocity u within the vortex core according to the Biot-Savart law.

Effect of varying p

Several optimizations are carried out for T = 5 when varying the value of p. Figure 4.4 provides the variation of energy gain and ∞-norm gain versus the value of p for T = 5 and m = 2. When p is increased the energy gain decreases, on the other hand, the ∞-norm increases but saturates for large values of p. A similar behavior has been found also for larger target times, the ∞-norm converging toward an asymptotic value for p ≥ 50. Therefore, the The 1, 10, 50norm optimal perturbations computed for m = 1, 2, 3 and T = 5 are shown in figure 4.5. The optimal perturbations are characterized by spirals similar to the case p = 1. For m = 1, increasing the value of p yields a more localized vorticity in the vortex core. Although being less visual, similar core localization is observed for m = 2 and 3. Figure 4.7 shows the initial vorticity amplitude as a function of the radius for θ = 0. For m = 1, the perturbation switches from a peripheral to a core perturbation. By recalling equation (4.5) we notice that the p-norm optimal perturbation directly acts on the production term (second term on the left-hand side) and deactivates the induction term for increasing values of p. This also stands for short horizon times in L 2 -norm (p = 1), where only the production term participates to the growth of the perturbation. Regarding the final optimal perturbation, the shapes at horizon time are almost unchanged despite an increase of p. Finally, the value of E p (T ) converges towards the ∞-norm gain as pointed out by equation 4.2.

Short time multiple modal states

The previously described pnorm optimization procedure has been implemented in a 2D simulation setup of the vortex flow without hypothesizing the azimuthal decomposition of the perturbation field. The methodology is implemented in the Nek5000 solver. We thus perform p-norm optimization in the cartesian coordinates system and keep the Lamb-Oseen vortex as the baseflow.

The optimization are carried out for relatively short horizon times (T ≤ 1). We focus on this short time horizon because we observe an interesting behavior of the pertubration, due to the p-norm optimization, that we describe below. For this range of short horizon times, figure 4.2 shows that in L 2 norm optimization (p = 1), modes m = 1, 2, 3 yield similar energy gains. The fact that the instability of these azimuthally periodic solutions do not differ much from each other creates the possibility for non selectivity at these short time, due to the optimization process. In p-norm this indefiniteness is amplified by equation (4.3), whose e p (T ) numerator dictates the adjoint initial perturbation and, upon the noise inherently present in the simulation, can potentially promote non periodic, more localized solutions. We could observe the occurence of this effect when initializing the optimization process at short horizon time with an isolated mode m, and playing on the value of p. The procedure led to the ap-

Flow past a wing

In order to apply the p-norm procedure to a more realistic case, we present below, the p-norm optimal perturbations for the flow past a finite span wing, following the work of [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF].

L 2 optimal perturbations

In their study, [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF] investigated the L 2 linear optimal perturbations in the steady and fully developped flow past a finite aspect ratio wing. The parameters that were explored are the horizon time T , the Reynolds number, the aspect ratio AR (span/chord), and the angle of attack α. The computational set-up consisted of a rectangular wing placed in a cubiform domain with a uniform inlet velocity. The streamwise length of the domain is 64 times the chord length (c): 4c upstream and 60c downstream from the leading edge of the wing. In the vertical direction, the position of the wing is such that the leading edge is at an equal distance of 10c from the bottom and top boundaries. Due to the symmetry of the flow about the midplane, the computations were carried out in the half-domain with the symmetry boundary condition enforced on the mid-plane y = 0. The reader may refer to figure 4.10 to visualize the computational setup. For Re = 1000, with an aspect ratio AR = 6 and α = 5 • of angle of attack, [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF] have found that the gain in kintetic energy increases optimal perturbation is mainly located at the wing and its vicinity. When increasing the aspect ratio, the optimal gain increases. Moreover, the computations have shown that the 2-D optimal gain is higher than that of the finite span wing. When α increases, the optimal gain increases.

Changing the shape of the airfoil to a NACA0012 increases the gain by an order of magnitude as compared to the flat plate. Regardless of the gain being affected by the changes of the parameters, [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF] have shown that the shape of the optimal perturbation is akin over the parameter space. When instantiating linear and nonlinear direct numerical simulations using the linear optimal perturbation, it is observed that the perturbation evolves as a wavepacket travelling at nearly the freestream velocity. In addition its energy increases in the near wake region following the Orr mechanism. It is also noted that the finite amplitude, nonlinear perturbation provokes an elliptical motion of the vortex core whereas the linear evolution does not affect the locus of the vortex.

The fact that the linear optimal perturbation is located near the wing (which is desirable from a practical control point of view) encourages the investigations of p-norm optimal perturbations.

p-norm optimal perturbations

Despite being suboptimal in terms of pure kinetic energy, using p-norm optimization can be applied for more realistic flows. In this subsection, we discuss the early investigations of p-norm optimal perturbations of the flow past a flat plate wing. The computational setup is identical to the one presented in the previous subsection. For more detail on mesh convergence, the reader may refer to table 1 from the original paper by [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF]. Unless otherwise stated, the selected parameters are: Re = 1000, T = 10, AR = 6 and α = 5 • . perturbation yields a gain in kinetic energy E 1 (T = 10)/E(0) = 1720 whereas the L 2 optimal perturbation yields E 1 (10)/E(0) = 2206. The 22% decrease in energy is expected as the localization step reduces the perturbation support i.e. the subset of the domain containing the elements which are not mapped to zero. For higher values of p, the kinetic energy gain decreases and suddenly drops for p ≥ 50. Figure 4.15 plots the shape of the optimal perturbation at horizon time T = 10 for p = 5. The chevron pattern is inverted compared to the one observed in figure 4.14 in which the perturbation is oriented downstream. This is the consequence of the action of the Orr mechanism. The importance of the Orr mechanism points out the two-dimensional nature of the perturbation growth in the near wake. As a reminder, [START_REF] Navrose | Transient growth in the near wake region of the flow past a finite span wing[END_REF] have shown that the gain following optimal perturbation was the highest for an infinite aspect ratio wing. The question that would need to be answered is whether, as it does when nonlinearity is triggered within the evolution of the L 2 -norm optimal, is whether such horizon time perturbation can transfer to the vortex core and provokes its displacement. Numerical simulation up to larger time would be worthwhile to carry out in future researches. Throughout this chapter, we have explored the use of p-norm optimal perturbation in the case of the single vortex and the flow past a plate. We have shown how the higher order exponent reduces the spatial extent of the perturbation field. The achieved growth rates in terms of kinetic energy remains at significant levels in spite of the p-norm constraint. It is found that initial optimal perturbation retains a complex spatial distribution, however support is shrunk.

Conclusion

Summary of the main results and outlook on future works

Real trailing wakes vortices persist for long durations, typically between one to three minutes. This duration is likely longer in the worst case scenario when crosswind and atmospherical turbulence, which are among the mechanisms for accelerated turbulent decay, are absent.

In the vicinity of the ground this persistence is complicated by the rebound effect, which provokes the stagnation of the vortices above the runway area. The rebound is a result of the viscosity acting on the flow provoked by the vortices at the ground. In the absence of viscosity the vortices follow an hyperbolic motion due to the image effects of the ground.

The viscous action at the ground generates a boundary layer that separates in the counter pressure gradient provoked by the vortices, which leads them to sustain entrainment by secondary, detached vorticiy. Within their hovering above the ground, the vortices are found to loop multiple times. In real air trafic, the lengthy presence of the vortices constraints the take-off and landing frequencies of aircrafts. The risky potential for hazardous wake vortex is greatest in this phase of an aircraft mission, compared to cruise, and has deserved much of the focus of the regulations on vortex separation distances. The way to circumvent these rules and gain take-off and departure frequencies, without compromising safety, has been tackled in recent works mostly by gaining increased knownledge on vortex dynamics in ground effect, and by developping control strategies, for instance by selecting ideal configurations of obstacles pattern at the ground or by activating long-wave instabilities ideal wavelength by modulated ground altitude. Vortex control has been at the heart of this thesis, with the objective in mind to find new and especially optimal ways to reduce vortex lifetime in and out of ground proximity. We have explored several theoretical control strategies, based on optimal control and optimal perturbation techniques.

We started in chapter [2] by describing the vortex dynamics in ground effect. Although largely discussed and presented in the literature before, we could verify and analyse the cause for vortex rebound and looping with precise two-dimensional direct numerical simulations.

The flow was then employed to calculate an optimal control strategy based on blowing and suction at the ground. The motivation behind this theoretical setting was based on the a priori simplicity of such a system, although in the end an attempt to evaluate the cost of such a control showed the applicative difficulty. We discuss in the outlook how, still, the sum of the thesis work could be combined efficiently for a reduced cost of such a control. The optimal control strategy was very rich from the point of view of the physics learned. First we showed that the vortex rebound effect could be countered, with consequence that the controlled vortex is able to follow an inviscid like trajectory, that is, a hyperbolic line first down the ground and then parallel to it, outward. The mechanism of the control was found to be the stabilization of the boundary layer formed at the ground, with dominant suction to prevent flow separation. In effect the control strategy was able to double the lateral distance of the vortices. An optimal horizon time of the control was also found, meaning that acting beyond it would not provide more benefit to the lateral displacement. The analysis ended by evaluating the energy required to practically apply the optimal control in at a real airport runway.

In the next chapter, chapter [3], we have explored the potential for vortex mitigation in ground effect through acting on the intrinsic stability of the vortices. Due to the peculiar kinematics of the vortices in ground proximity (described above : rebound and looping) we have employed a linear optimal perturbation strategy to describe the linear stability of the vortices, first in a two-dimensaionl setting and then in three-dimensions. The methodology was fully based on the flow initiated from a pair of Lamb-Oseen vortex. After the numerical implementation of the method and its validation, the results of the two-dimensional analysis first showed the importance of the boundary layer at the ground in the overall dynamics of the interaction. The vortices were shown to be affected by displacement modes, as could be anticipated from previous study of the two-dimensional dynamics and optimal perturbation results obtained for instance by [START_REF] Jugier | Stabilité bidimensionnelle de modèles de sillage d'aéronefs[END_REF]. The interest was greater when looking at three-dimensional perturbations. First the amplification potential of three-dimensional perturbations is much larger than in two-dimensions. Second, after a normalized time of about 2 to 3, the short wavelengths become dominant among the spectrum, showing that ground effect promotes the rise of small scale perturbations into the flow. This was attributed to the consequence of the boundary layer, which appears as a shorter scale compared to the vortices, and the separation of which provides strong mechanisms for flow destabilization. Yet at short time, the optimal wavelength is found to be that of the long-wavelength typical of the vortex dynamics in free space. We concluded this linear analysis by conducting non-linear DNS simulations of the flow initialized by optimal perturbation of finite amplitudes. Initial amplitudes of 0.1 to 10% were considered. When applied with an amplitude of 1% of the total energy, the system is found to be completely mitigated, developping turbulence in a short amount of time. Comparison was taken against the uncontrolled system which, in the same time, remained coherent. This validated the strong response of the vortices to short wavelength initial perturbations.

Eventually chapter [4] questions the localization of the optimal perturbations, such as those found in the previous chapter, however with a limited application to the case of a single vortex case, without ground effect, and for the flow past a wing, following some previous works of the literature. The motivation is driven by the application of optimal perturbation results in realistic configurations. Usually the spatially wide distribution of the optimal perturbation structure is a difficulty when one searches for control devices to generate them. One way to relax the issue is to optimize more localized perturbations and this can be achieved using suboptimal optimization, like p-norm optimization. p-norm replaces the usual kinetic energy of the perturbation with a p-norm evaluation of the velocity field. With high value of the p exponent, the norm tends to select only the area of the flow where the maximum growth is achieved. Although the initial optimal perturbation retains a complex spatial distribution, the support is thinned and possibly realized using control devices in fewer places. The application of the p-norm optimal shows the localization effect for the single vortex, as well as for the flow past a flat plate. The achieved growth rates in terms of kinetic energy remains at significant levels in spite of the p-norm constraint.

Outlook and Perspectives

Various aspects of the dynamics of vortices have been investigated in the course of this thesis. Behind the focus of ground effect, and the motivation of flow control for accelerated vortec decay, fundamental results showing how vortices behaves and respond to perturbations and the presence of hard surfaces placed around them have been obtained.

The results that were obtained form a partially coherent set that provides possible outlook for future works and investigations. The most interesting is to take a global look on optimal control described in chapter [2] with the optimal perturbation in three-dimensional space detailed in chapter [3]. Indeed, the optimal control can be seen as a way to trigger three-dimensional perturbation at the optimal wavelength of the vortex dynamics. If the optimal control were distributed axially at the wavenumber k 6 which is found for large horizon time, given that part of the optimal perturbation is located at the ground for such time horizons, then one could expect to trigger such flow departure from coherence that are observed in the DNS results at the end of chapter [3]. This objective of flow decay would then replace the objective that was chosen in the two-dimensional setting, of optimizing the lateral position of the vortices. A future study could replicate the optimal control algorithm with an adapted objective for perturbation growth and a three-dimensional domain. Although the computational cost would be high, it should be affordable with current ressources.

To reduce the burden of implementing the blowing/suction strategy in a real configuration, investigations towards the optimization of simpler control means would be helpful. Following the work of [START_REF] Stephan | Aircraft wake-vortex decay in ground proximity-physical mechanisms and artificial enhancement[END_REF] and the efficiency of plate lines devices, it could be envisaged to optimize fix patterns of geometric variations of the ground surface. Such geometries could be defined by a limited set of parameters thereby allowing for simpler optimization algorithms than the adjoint technique used in the present work.

Another research thread is to investigate the non-linear fate of perturbations and taking non-linearity into account from the optimization stage. With the increased compute capability and storage, adding the effect of non-linearities to the optimal perturbations in ground effect would be achievable, and provide worhtwhile results. [START_REF] Navrose | Optimal perturbation for two-dimensional vortex systems: route to non-axisymmetric state[END_REF] have shown that a nonlinear optimal perturbation can outperform the linear optimal perturbation for a 2-D isolated vortex and counter-rotating vortex pair. The complex kinematics of the vortices in ground effect would legitimate the use of a non-linear model because it is expected that perturbations to the secondary structures of the flow may not be insignificant to them at their early stage of formation (separation of the boundary layer for instance). Retroactions of the perturbation upon the flow are to be expected.

Eventually it would be worthwhile to attempt experimentally some of the theoretical results obtained in this thesis. The question whether the highlighted mechanisms are efficient in a realistic environnement with imperfections of symmetry, wall roughness, boundaries, and account of the way to generate the initial vortices would be a powerful incentive to dig into a broader domain of the vortex dynamics. Successful results would open new perspectives for regulating air trafic with regards to wake hazard.

Long French Résumé

Les tourbillons de sillage générés derrière les avions de transport peuvent persister sur des longues périodes de temps, généralement comprises entre une à trois minutes. Cette durée est encore plus longue dans le cas le plus défavorable, lorsque le vent de travers et la turbulence atmosphérique, qui font partie des mécanismes de décroissance accélérée des tourbillons, sont absents. Au voisinage du sol, cette persistance est amplifiée par l'effet de rebond, qui provoque la persistance des tourbillons au-dessus de la zone de la piste. Le rebond résulte de la viscosité agissant sur l'écoulement provoqué par les tourbillons au niveau du sol. En l'absence de viscosité, les tourbillons suivent une trajectoire hyperbolique qui s'explique par les effets d'image dû au sol. L'action visqueuse au niveau du sol génère une couche limite qui se sépare dans le gradient de pression adverse provoqué par les tourbillons, ce qui les conduit à subir l'entraînement par de la vorticité secondaire détachée. Au cours de leur mouvement au-dessus du sol, les tourbillons effectuent de multiples boucles. Dans le trafic aérien réel, cette présence prolongée des tourbillons limite les fréquences de décollage et d'atterrissage des avions. C'est dans cette phase de la trajectoire d'un avion, par rapport à la phase de croisière, que le risque de tourbillon de sillage dangereux est le plus élevé et c'est pourquoi les réglementations sur les distances de séparation des tourbillons ont fait l'objet d'une grande attention. Le contrôle des tourbillons a été au coeur de cette thèse, avec pour objectif de trouver de nouveaux moyens, particulièrement optimaux, d'assurer la sécurité des vols. L'objectif de cette thèse a été de trouver des moyens optimaux pour réduire la durée de vie des tourbillons dans et hors de la proximité du sol. Pour cela, nous avons exploré plusieurs stratégies théoriques de contrôle, basées sur des techniques de contrôle optimal et de perturbation optimale.

Nous avons commencé dans le chapitre 2 par décrire la dynamique des tourbillons en effet de sol. Bien que largement discutée et présentée dans la littérature auparavant, nous avons pu vérifier et analyser la cause du rebond et le bouclage des tourbillons avec des simulations numériques directes en deux dimensions. L'écoulement a ensuite été utilisé pour calculer une stratégie de contrôle optimale basée sur le soufflage et la succion. La motivation derrière ce A.2. DIRECT NUMERICAL SIMULATIONS des simulations DNS non linéaires de l'écoulement initialisées par une perturbation optimale d'amplitudes finies. Des amplitudes initiales de 0,1 à 10% ont été considérées. Lorsqu'il est appliqué avec une amplitude de 1% de l'énergie totale, le système est trouvé est complètement atténué, développant la turbulence en peu de temps. La comparaison a été faite par rapport au système non contrôlé qui, dans le même temps, est resté cohérent. Ceci valide la forte réponse des tourbillons aux perturbations initiales de courte longueur d'onde.

Enfin, le chapitre 4 questionne la localisation des perturbations optimales, telles que celles trouvées dans le chapitre précédent, avec toutefois une application limitée au cas d'un seul tourbillon, sans effet de sol, et pour l'écoulement devant une aile, suivant certains travaux antérieurs. La motivation est expliquée par l'application des résultats de perturbation optimale dans des configurations réalistes. Habituellement, la large distribution spatiale de la structure de la perturbation optimale représente une difficulté lorsque l'on recherche des résultats surtout lorsqu'on cherche des dispositifs de contrôle pour les générer. Une façon de résoudre ce problème est de chercher des perturbations plus localisées, ce qui peut être réalisé en utilisant une procédure sous-optimale, comme l'optimisation en norme p. La norme p remplace l'énergie cinétique habituelle de la perturbation par une évaluation en norme p du champ de vitesse. Avec une valeur élevée de l'exposant p, la norme tend à sélectionner uniquement la zone de l'écoulement où la croissance maximale est atteinte. Bien que la perturbation optimale initiale conserve une distribution spatiale complexe, le support est amoindri et éventuellement réalisé en utilisant des dispositifs de contrôle en moins d'endroits. L'application de la norme p montre l'effet de localisation pour le tourbillon unique, ainsi que pour l'écoulement autour une plaque plane. Les taux de croissance atteints en termes d'énergie cinétique restent à des niveaux significatifs malgré la contrainte de la norme p. Abstract : Real trailing wakes vortices persist for long durations, typically between one to three minutes. This duration is likely longer in the worst case scenario when crosswind and atmospherical turbulence, which are among the mechanisms for accelerated turbulent decay, are absent. In the vicinity of the ground this persistence is complicated by the rebound effect, which provokes the stagnation of the vortices above the runway area. The rebound is a result of the viscosity acting on the flow provoked by the vortices at the ground. In the absence of viscosity the vortices follow an hyperbolic motion due to the image effects of the ground. The viscous action at the ground generates a boundary layer that separates in the counter pressure gradient provoked by the vortices, which leads them to sustain entrainment by secondary, detached vorticiy. Within their hovering above the ground, the vortices are found to loop multiple times. In real air trafic, the lengthy presence of the vortices constraints the takeoff and landing frequencies of aircrafts. The risky potential for hazardous wake vortex is greatest in this phase of an aircraft mission, compared to cruise, and has deserved much of the focus of the regulations on vortex separation distances. The way to circumvent these rules and gain take-off and departure frequencies, without compromising safety, has been tackled in recent works mostly by gaining increased knownledge on vortex dynamics in ground effect, and by developping control strategies. Vortex control has been at the heart of this thesis, with the objective in mind to find new and especially optimal ways to reduce vortex lifetime in and out of ground proximity. We have explored several theoretical control strategies, based on optimal control and optimal perturbation techniques.

Institut Polytechnique de Paris 91120 Palaiseau, France

1

  Leonardo Da Vinci, Studies of water passing obstacles and falling, c. 1508-9. Wikimedia Commons (cropped image) . . . . . . . . . . . . . . . . . . . . . . 2 Cross-section of the flow in a typical room, showing different eddy sizes. This was made visible using a laser, a glass rod (to refract the light into a plane) and a fog machine. Credit : YouTube channel 3Blue1Brown / Visualizing Turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Aircraft wakes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Sketch of a bi-dimensional grid forming an element with N = 4. . . . . . . . . . 1.2 Sketch of the iterative procedure aproximating the optimal solution. . . . . . . 1.3 Finite element P 1 and P 2 discretization. Degrees of freedom are denoted q i . . . 3.1 Positive-tilt (left) and negative-tilt (right) spirals, from (Pradeep & Hussain, 2006) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Visualization of the secondary vorticity associated with the vertical rings mode (h 0 /b 0 = 5); (a-d) show a top view and (e-h) show a side view of the vortex pair impinging on a flat surface. Figure taken from Asselin & Williamson (2017). 3.3 Sketch of the optimization algorithm, Douglas et al. (2000). The line search parameter is α k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Optimal gain a function of time for ASYM and SYM optimal perturbations. . . 3.5 Distribution of the norm of the strain rate tensor of the baseflow for time t = 0, 1.25, 2.5 and 3.75 respectively corresponding to subfigures a), b) c) and d). Green solid lines correspond to the iso-contours of vorticity (ω z /ω 0,max = ±0.1, ±0.2) of the primary vortices. . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Initial vorticity distribution of ASYM (top row) and SYM (bottom row) optimal perturbation computed for different horizon times 2, 3, 4 and 5 (from left to right). Increasing horizon time implies a shift of the optimal perturbation shape towards the ground, pointing out the influence of horizon time on the boundary layer effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES

  2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 iv LIST OF FIGURES 3.15 Time sequence of vorticity magnitude of the nonlinear response to the symmetric short wave optimal perturbation computed for T = 5 and k = 6 with an initial amplitude of ε = 10 -1 . The three levels of iso-contours are ω/ω 0,max = 0.05, 0.1 and 0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Time evolution of the kinetic energy for various values of the initial perturbation amplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Time evolution of the circulation γ z at z = λ/2 . . . . . . . . . . . . . . . . . . 3.18 Time evolution of the circulation γ z at z = λ/4 . . . . . . . . . . . . . . . . . . 3.19 Time evolution of the enstrophy E. . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Example of localization step (see eq. 4.3) with the mock function: u(x, y) = 1 1 + x 2 + y 2 , using p = 1 (left), p = 100 (middle) and p = 1000 (right). The adjoint field approaches the dirac delta function centered around the maximum value of u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Linear optimal gain as a function of the horizon time T for m = 1, 2, 3 at Re = 5000, see Bisanti (2013). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Top: Perturbation energy growth over time for the linear optimal perturbation of a Lamb-Oseen vortex for m = 2, p = 1, T = 5 and Re = 5000. Bottom left: z-vorticity at t = 0. Bottom middle: z-vorticity at t = 5. Bottom right: z-vorticity at t = 10 showing inverted spirals. The dotted line indicates the vortex dispersion radius. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Result of the p-norm optimization for m = 1 and T = 5 showing the effect of the optimization on the differnt norms as a function of p. (left) Energy gain E 1 /E(0), (middle) p-norm evaluation E p /E(0) and (right)and ∞-norm gain. 4.5 Result of the optimization in p-norm showing the shape of the initial p-norm perturbation for m = 1, 2, 3 (top to bottom) and p = 1, 10, 50 (left to right). The contour levels are linearly spaced between -10 and 10 with 512 distinct levels. The dotted lines represent the circle of radius the dispersion radius of the base Gaussian vortex. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Result of the optimization in p-norm showing shape of the final perturbation for m = 1, 2, 3 (top to bottom) and p = 1, 10, 50 (left to right). The contour levels are linearly spaced between -10 and 10 with 512 distinct levels. The dotted lines represent the circle of radius the dispersion radius of the base Gaussian vortex. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Vorticity amplitude ω z (r, θ = 0) for m = 1, 2, 3 (left, middle, right) showing the increased amplitudes for r → 0 when p increases. This phenomenon is in agreement with the observations flowing from equation (4.5). . . . . . . . . . . v LIST OF FIGURES 4.8 Spiral (left) vs. Whip (right) : Contours of vorticity magnitude. Computed

  • . The grey contour represent |ω x | = 0.1. For this very illustration, the contour surface has been shown for the full wing by reflecting the flow at the mid-plane y = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.11 Optimal perturbation for a flat plate: Isocontours ω y = ±0.2 for T = 10 of the linear optimal perturbation. The wing root lies on the plane of symmetry of the geometrical set-up. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.12 Optimal perturbation for a NACA0012 rectangular wing. Computed at T = 10, AR = 6 and α = 5. The contour represent the spanwise vorticity at different chord sections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.13 Evolution of E(T )/E(0) (left), E p (T )/E(0) (middle), E ∞ (T )/E(0) (right) as a function of p in the case of the flow past a wing at Re = 1000, α = 5 and AR = 6. The evolutions are similar to that observed in figure 4.4. The values of p correspond to 1, 2, 5, 10, 50 and 100 from left to right. . . . . . . . . . . . 89 4.14 Comparison between the 1, 2, 5 and 100-norm (left to right, top to bottom)

Figure 1 :

 1 Figure 1: Leonardo Da Vinci, Studies of water passing obstacles and falling, c. 1508-9. Wikimedia Commons (cropped image)

Figure 2 :

 2 Figure 2: Cross-section of the flow in a typical room, showing different eddy sizes. This was made visible using a laser, a glass rod (to refract the light into a plane) and a fog machine. Credit : YouTube channel 3Blue1Brown / Visualizing Turbulence.

Figure 3 :

 3 Figure 3: Aircraft wakes.

Figure 1 . 1 :

 11 Figure 1.1: Sketch of a bi-dimensional grid forming an element with N = 4.

Figure 1 . 2 :

 12 Figure 1.2: Sketch of the iterative procedure aproximating the optimal solution.

Figure 1 . 3 :

 13 Figure 1.3: Finite element P 1 and P 2 discretization. Degrees of freedom are denoted q i .

  www.cambridge.org/core/terms.

FIGURE 1 .FIGURE 2 .

 12 FIGURE 1. Schematic (not to scale) of the physical set-up. Dimensions are indicated.The right part of the domain is the computational domain and is shaded in grey. The disks represent the two Lamb-Oseen vortices at time t = 0. Symmetric conditions are applied at the right, upper and left boundaries (respectively denoted ∂Ω 2,3,4 ) and Dirichlet is prescribed at the ground denoted ∂Ω 1 .

FIGURE 3 .

 3 FIGURE 3. Trajectories of the right vortex for Re = 1000 (dashed line), for Re = 5000 (solid line) and inviscid flow (dotted line). Simulation time is t = 15. Circle markers correspond to Doppler lidar measurements achieved at Toulouse-Blagnac airport (European Commission 2015; Hallermeyer 2017) during 50 s (in this case b 40 m and circulation is Γ 0 500 m 2 s -1). Note the existence of an average cross-wind (u crosswind ∼ -1 m s -1 ) in the lidar case, not taken into account in the simulated data.

FIGURE 4 .

 4 FIGURE 4. Uncontrolled vorticity field computed for Re = 1000 at t = 0.25, 2.5, 3.75, 5, 10, 15. Thick lines represent boundary layer profiles are given at locations x = 0.5, 1, 1.5, 2.0, 2.5. Dashed lines represent the pressure gradient along the wall (at y = 0). (a) Descending vortex t = 1.25. (b) Boundary layer and detachment of a secondary vortex t = 2.5. (c) Rebound of the vortex t = 3.75. (d) Entrainment of the secondary vortex t = 5. (e) Entrainment of the secondary vortex t = 10. ( f ) Entrainment of the secondary vortex t = 15.

FIGURE 5 .

 5 FIGURE 5. Evolution of Γ pv and Γ tot at Re = 200, 1000 and T = 15 for the uncontrolled case. Vertical dotted lines correspond to times t = 1.25, 2.5, 3.75, 5, 10 referring to figure 4.

  and a, b are functions of space and time and c, d are functions of space alone.

FIGURE 7 .

 7 FIGURE 7. Sketch of the optimization algorithm. The line search parameter is α k .

  https://doi.org/10.1017/jfm.2019.1002 Downloaded from https://www.cambridge.org/core. /www.cambridge.org/core/terms.

FIGURE 9 .

 9 FIGURE 9. Vorticity field at horizon time T = 5 for Re = 200 (a-d) and 1000 (e-h) for increasing l. Dashed lines show the locations of the vorticity-squared barycentres, J 1 (5, l).

  https://doi.org/10.1017/jfm.2019.1002 Downloaded from https://www.cambridge.org/core. /www.cambridge.org/core/terms.

FIGURE 10 .

 10 FIGURE 10. Effect of horizon time. (a) Evolution of the vortex position J 1 (T, l) for Re = 200, 1000, respectively in dashed and dotted lines, for the uncontrolled case l = ∞ (black) and for l = 1 (red). (b) Evolution of the ratio J m (T, 1)/J m (T, ∞), computed for Re = 200 (dashed line), Re = 1000 (solid lines) and l = 1. For this value of l and Re = 200, increasing the horizon time, beyond T 5, leads to the complete absorption of the primary vortex. No data are consequently available for larger time.

  www.cambridge.org/core/terms.

FIGURE 14 .

 14 FIGURE 14. Plot of the Q-criterion at time t = 2 for (a) (1000, 5, ∞) and (b) (1000, 5, 1). Creation of a rotational area in the boundary layer (positive values of the Q-criterion) that corresponds to the pusher vortex.

  www.cambridge.org/core/terms.Vortex pair in ground effect, dynamics and optimal control 885 A26-17

  FIGURE 15. (a) Evolution of primary vortex and total circulation for (200, 5, 1) (black) and (1000, 5, 1) (red). For times beyond t = 5, control has been turned off. (b) Evolution of the centroid trajectories for 0 t 5 for inviscid, controlled (1000, 5, 1) and uncontrolled flow (1000, 5, ∞). Circled markers correspond to the centroid position at t = 0, 1, 2, 3, 4, 5.

  www.cambridge.org/core/terms. 38 885 A26-18 A. Wakim, V. Brion, A. Dolfi-Bouteyre and L. Jacquin(1000, 5, 1) H : y c0 = 5 R : a = 0.1 DAT : y c0 = 5, t a = 3

  FIGURE 16. (a,b) Final vorticity field and (c,d) time trace of the control at x = 1 achieved by the control. Panels (a,c) (respectively b,d) refer to simulation H (respectively R). To be compared to the optimal case (1000, 5, 1), see figures 11 and 12.

FIGURE 17 .

 17 FIGURE 17. Time evolution -İ11 , İ12 , İ13 from left to right for the uncontrolled case (dashed lines) and the controlled case (solid lines) for (Re, T, l) = (1000, 4, 10).

  www.cambridge.org/core/terms.

  www.cambridge.org/core/terms.

  www.cambridge.org/core/terms.

  www.cambridge.org/core/terms.

  . The strain S = r∂ r (V /r) of a L-O vortex is negative, thus energy production requires a positive Reynolds stress uv. This is accomplished by flow structure in the form of spirals coiled around the vortex core with positive tilt. The transient growth is naturally stopped as the differential rotation by the base flow transforms the positive-tilt spirals into negative-tilt spirals. A schematic defining positive and negative tilt spirals is shown in figure3.1 (source:[START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF]).

Figure 3 . 1 :

 31 Figure 3.1: Positive-tilt (left) and negative-tilt (right) spirals, from[START_REF] Pradeep | Transient growth of perturbations in a vortex column[END_REF] .

Figure 3 . 2 :

 32 Figure 3.2: Visualization of the secondary vorticity associated with the vertical rings mode (h 0 /b 0 = 5); (a-d) show a top view and (e-h) show a side view of the vortex pair impinging on a flat surface. Figure taken from Asselin & Williamson (2017).

Figure 3 . 3 :

 33 Figure 3.3: Sketch of the optimization algorithm, Douglas et al. (2000). The line search parameter is α k .

Figure 3 . 4 :

 34 Figure 3.4: Optimal gain a function of time for ASYM and SYM optimal perturbations.

Figure 3 .

 3 Figure 3.8 plots the evolution of the perturbation kinetic energy gain as a function of the wave number k for various horizon times and symmetries. Similar behaviours are observed

Figure 3 . 9 :

 39 Figure 3.9: Axial vorticity for the symmetric optimal pertubations, k = 0.8 (left) and k = 6 (right), SYM perturbations, horizon time T = 5. Stronger ground effects are observed for k = 6: the iso contour of vorticity |ω z | 5 are observable in the vicinity of the ground plane whereas they are not for longer wavelengths.

Figure 3 .

 3 Figure 3.10 shows the time sequence of vorticity magnitudes contours of symmetrical optimal perturbations k = 0.8 (3.10a) and k = 6 (3.10b) for t = 0, 1, 2, 3, 4, 5 (for horizon time T = 5).During the first stage of the primary vortex descent (t = 0 → t = 2), the effect of the wall is weak and the perturbations grow into longwave deformations. As the boundary layer of the baseflow separates at the ground the perturbations appear to grow strongly in this region of the boundary layer (t = 2 → t = 3). At later times the perturbation then acts principally on the secondary vortices (t = 4 → t = 5).

Figure 3

 3 Figure3.10: Time sequence of axial vorticity of symmetric optimal perturbation at k = 0.8 (top) and k = 6 (bottom), SYM perturbations, horizon time T = 5. The optimal perturbation initially acts as a displacement mode inside the primary vortex cores and when boundary layer separates, the perturbation drastically grows in the boundary layer. When the boundary layer has detached, displacement modes grow into the secondary vortices giving rise to structures similar to those that were observed in two dimensions. Vorticity magnitude ranges between 10 and 900 (top). and between 10 and 3 × 10 4 (bottom).

Figure 3 . 13 :

 313 Figure3.13: Time sequence of vorticity magnitude of the nonlinear response to the symmetric short wave optimal perturbation computed for T = 5 and k = 6 with an initial amplitude of ε = 10 -3 . The three levels of iso-contours are ω/ω 0,max = 0.05, 0.1 and 0.2

Figure 3 . 14 :

 314 Figure3.14: Time sequence of vorticity magnitude of the nonlinear response to the symmetric short wave optimal perturbation computed for T = 5 and k = 6 with an initial amplitude of ε = 10 -2 . The three levels of iso-contours are ω/ω 0,max = 0.05, 0.1 and 0.2

Figure 3 . 16 :

 316 Figure 3.16: Time evolution of the kinetic energy for various values of the initial perturbation amplitude.

Figure 3 .

 3 Figure 3.17: Time evolution of the circulation γ z at z = λ/2

Figure 3 . 19 :

 319 Figure 3.19: Time evolution of the enstrophy E.

Figure 4 . 1 :

 41 Figure 4.1: Example of localization step (see eq. 4.3) with the mock function: u(x, y) = 1 1 + x 2 + y 2 , using p = 1 (left), p = 100 (middle) and p = 1000 (right). The adjoint field approaches the dirac delta function centered around the maximum value of u.

Figure 4 . 2 :

 42 Figure 4.2: Linear optimal gain as a function of the horizon time T for m = 1, 2, 3 at Re = 5000, see Bisanti (2013).

Figure 4 . 3 :

 43 Figure 4.3: Top: Perturbation energy growth over time for the linear optimal perturbation of a Lamb-Oseen vortex for m = 2, p = 1, T = 5 and Re = 5000. Bottom left: z-vorticity at t = 0. Bottom middle: z-vorticity at t = 5. Bottom right: z-vorticity at t = 10 showing inverted spirals. The dotted line indicates the vortex dispersion radius.

Figure 4 . 4 :

 44 Figure 4.4: Result of the p-norm optimization for m = 1 and T = 5 showing the effect of the optimization on the differnt norms as a function of p. (left) Energy gain E 1 /E(0), (middle) p-norm evaluation E p /E(0) and (right)and ∞-norm gain.

Figure 4 . 5 :

 45 Figure 4.5: Result of the optimization in p-norm showing the shape of the initial p-norm perturbation for m = 1, 2, 3 (top to bottom) and p = 1, 10, 50 (left to right). The contour levels are linearly spaced between -10 and 10 with 512 distinct levels. The dotted lines represent the circle of radius the dispersion radius of the base Gaussian vortex.

Figure 4 . 6 :

 46 Figure 4.6: Result of the optimization in p-norm showing shape of the final perturbation for m = 1, 2, 3 (top to bottom) and p = 1, 10, 50 (left to right). The contour levels are linearly spaced between -10 and 10 with 512 distinct levels. The dotted lines represent the circle of radius the dispersion radius of the base Gaussian vortex.

Figure 4 . 7 :

 47 Figure 4.7: Vorticity amplitude ω z (r, θ = 0) for m = 1, 2, 3 (left, middle, right) showing the increased amplitudes for r → 0 when p increases. This phenomenon is in agreement with the observations flowing from equation (4.5).

Figure 4

 4 Figure 4.9: p-norm objective function as a function of p. The dotted line represent "spiral" solution and the solid line represent "whip" solutions. These two branches appear after amplification of noise during the initialilization of the adjoint.

Figure 4 . 10 :

 410 Figure 4.10: Flow past a finite span wing (purple) at Re = 1000 with AR = 6 and α = 5 • . The grey contour represent |ω x | = 0.1. For this very illustration, the contour surface has been shown for the full wing by reflecting the flow at the mid-plane y = 0.

Figure 4 . 11 :

 411 Figure 4.11: Optimal perturbation for a flat plate: Isocontours ω y = ±0.2 for T = 10 of the linear optimal perturbation. The wing root lies on the plane of symmetry of the geometrical set-up.

Figure 4 .

 4 Figure 4.12: Optimal perturbation for a NACA0012 rectangular wing. Computed at T = 10, AR = 6 and α = 5. The contour represent the spanwise vorticity at different chord sections.

Figure 4 .

 4 13 shows the evolution of the perturbation kinetic gain, the p-norm cost function, and the ∞-norm gain as a function of p. The trend of the curves follow that of figure 4.4. The saturation in infinite norm gain occurs

Figure 4 . 13 :

 413 Figure 4.13: Evolution of E(T )/E(0) (left), E p (T )/E(0) (middle), E ∞ (T )/E(0) (right) as a function of p in the case of the flow past a wing at Re = 1000, α = 5 and AR = 6. The evolutions are similar to that observed in figure 4.4. The values of p correspond to 1, 2, 5, 10, 50 and 100 from left to right.

  Figure4.14 compares the L 2 and p-norm optimal perturbation shapes for p = 1, 2, 5 and 100. Similar overall shapes are observed for the various perturbations, except that the higher order norm localizes the perturbation closer to the mid plane y = 0. The 5-norm optimal

Figure 4 . 14 :

 414 Figure 4.14: Comparison between the 1, 2, 5 and 100-norm (left to right, top to bottom) optimal perturbations for a flat plate wing. Computed at T = 10, AR = 6 and α = 5. As expected the p-norm procedure localizes the perturbation. The contour represent the streamwise vorticity levels at ±0.2 of the perturbations. The shape of the optimal perturbation do not vary when p ≥ 2, this saturation occurs when the diffusion and p-norm lengthscales are equivalent.

Figure 4 .

 4 Figure 4.15: Perturbation streamwise vorticity at horizon time T = 10 (red and blue). The green chevrons represent the streamwise vorticity of the initial optimal perturbation. The chevrons invert while travelling downstream.

Titre:

  Contr ôle des Paires de Tourbillons en Effet de SolMots cl és : tourbillons, contr ôle d' écoulement, optimisation R ésum é : Les tourbillons de sillage persistent pendant de longues dur ées, g én éralement entre une et trois minutes. Cette dur ée est encore plus longue dans le cas le plus d éfavorable, lorsque le vent de travers et la turbulence atmosph érique, qui font partie des m écanismes de d écroissance acc él ér ée des tourbillons, sont absents. Au voisinage du sol, cette persistance est amplifi ée par l'effet de rebond, qui provoque la stagnation des tourbillons au-dessus de la zone de la piste. Le rebond r ésulte de la viscosit é agissant sur l' écoulement provoqu é par les tourbillons au niveau du sol. En l'absence de viscosit é, les tourbillons suivent une trajectoire hyperbolique d û aux effets d'image du sol. L'action visqueuse au niveau du sol g én ère une couche limite qui se s épare dans le gradient de pression adverse provoqu é par les tourbillons, ce qui les conduit à subir l'entraînement par de la vorticit é secondaire d étach ée. Au cours de leur mouvement au-dessus du sol, les tourbillons effectuent de multiples boucles. Dans le trafic a érien r éel, la pr ésence prolong ée des tourbillons limite les fr équences de d écollage et d'atterrissage des avions. C'est dans cette phase de la trajectoire d'un avion, par rapport à la phase de croisi ère, que le risque de tourbillon de sillage dangereux est le plus élev é et c'est pourquoi les r églementations sur les distances de s éparation des tourbillons ont fait l'objet d'une grande attention. Le contr ôle des tourbillons a ét é au coeur de cette th èse, avec pour objectif de trouver de nouveaux moyens, particuli èrement optimaux, d'assurer la s écurit é des vols. L'objectif de cette th èse a ét é de trouver des moyens optimaux pour r éduire la dur ée de vie des vortex dans et hors de la proximit é du sol. Pour cela, nous avons explor é plusieurs strat égies th éoriques de contr ôle, bas ées sur des techniques de contr ôle optimal et de perturbation optimale.Title : Control of Vortex Pairs in Ground EffectKeywords : vortex, flow control, optimization

Table 1 :

 1 Current FAA standards for aircraft separation. In nautic miles (NM).

	INTRODUCTION

TABLE 1 .

 1 Effect of the initial height, radius and activation time on the objective function results. Left column corresponds to the reference optimal case (1000, 5, 1) already shown in figure11. The second, third and fourth columns correspond to simulations H, R and DAT respectively, which yield similar results.

Department of Aerodynamics Aeroelasticity Acoustics, ONERA, 8 rue des Vertugadins, 92190 Meudon, France

Departement of Optics and Associated Techniques, ONERA, 8 Chemin de la Hunière, 91120 Palaiseau, France

Acknowledgements

This work was granted access to the HPC resources IDRIS under the allocation 2016-100746 made by GENCI. The work is primarily supported by the PHYWAKE project dedicated to vortex dynamics and funded by the French Directorate for Civil Aviation (DGAC). The careful reading and support of Navrose have been deeply appreciated and were decisive to finalize this work.

Acknowledgements

Dans le chapitre suivant, le chapitre 3, nous avons exploré le potentiel d'atténuation des tourbillons en effet de sol en agissant sur la stabilité intrinsèque des tourbillons. En raison de la cinématique particulière des tourbillons en effet de sol (décrite ci-dessus : rebond et bouclage), nous avons employé une stratégie de perturbation optimale linéaire pour décrire la stabilité linéaire des tourbillons, d'abord dans un cadre bidimensionnel, puis en trois dimensions. La méthode était entièrement basée sur l'écoulement initié à partir d'une paire de tourbillons Lamb-Oseen. Après la mise en oeuvre numérique de la méthode et sa validation, les résultats de l'analyse bidimensionnelle ont d'abord montré l'importance de la couche limite au sol dans la dynamique globale de l'interaction. Les tourbillons se sont révélés être affectés par les modes de déplacement, comme on pouvait l'anticiper à partir de l'étude précédente de la dynamique bidimensionnelle et des résultats de perturbation optimale obtenus par exemple par Jugier (2016). L'intérêt était plus grand lorsqu'on s'est intéressé aux perturbations tridimensionnelles. Premièrement, le potentiel d'amplification des perturbations tridimensionnelles est beaucoup plus important qu'en deux dimensions. Deuxièmement, après un temps normalisé d'environ 2 à 3, les courtes longueurs d'onde deviennent dominantes parmi le spectre, ce qui montre que l'effet de sol favorise la l'augmentation des perturbations à petite échelle dans l'écoulement. Ceci a été attribué à la conséquence de la couche limite, qui apparaît comme une échelle plus courte par rapport aux tourbillons, et dont la séparation fournit de forts mécanismes de destruction de l'écoulement. dont la séparation fournit des mécanismes forts pour la déstabilisation de l'écoulement. Pourtant, à court terme, la longueur d'onde optimale optimale est celle de la grande longueur d'onde typique de la dynamique des vortex en espace libre. espace libre. Nous avons conclu cette analyse linéaire en effectuant

Declaration of interests

The authors report no conflict of interest.

Appendix A. Compatibility equation

The derivation leading to the compatibility (4.7) is here derived. Without loss of generality let us consider three-dimensional fields. We introduce the following inner spatial product: The goal is to find an initial perturbation u 0 that maximizes the gain of perturbation kinetic energy over a given horizon time T . The perturbation kinetic energy is defined as follows:

and the gain in kinetic energy with respect to the initial perturbation kinetic energy at horizon time T reads as:

Hence, the optimal perturbation optimizes the value of G at the horizon time. To solve such a problem, we use the Lagrange multipliers (presented in Chapter 1). The Langrangian function L, as in equation 1.22, is defined as

where F (q) = 0 represents the governing equation (3.6), H the initial condition operator , q and λ the Lagrange multipliers (or adjoint variables). E 0 constrains the perturbation kinetic energy at t = 0 to be equal to a fixed value. (In this linearized framework, E 0 can be set to 1 since G(T ) is independant of the initial amplitude of the perturbation). This constraint is applied using a geometric update technique discussed later in this section. In order to reach the optimal state, the gradient of the Lagragian functional with respect to all the variables must be equal to zero. Doing so with respect to :

• q, yields equation (3.6)

• q, yields the adjoint equations

• q(T ), gives the compatibility equation

• q 0 gives the optimality condition

with u 0 the normalized optimal perturbation velocity field and ε the initial amplitude. For what follows, ε varies in [10 -3 , 10 -1 ] and we evaluate the flow dynamics for a selection of ε values in this range. The direct numerical simulations convergence properties are provided in Appendix (A.2).

Unperturbed flow ε = 0 :

In chapter [2] we have described the dynamics of a vortex pair in ground effect without perturbation. The evolution of the vorticity magnitude, in three-dimensional iso-contours, is shown in figure 3.11, as a reference for the finite values of ε that follow. The baseflow is robust and remains two-dimensional during multiple rebounds and finally becomes subject to perturbation (t > 22). This betrays the unavoidable arising and growth of 3D perturbations in the simulation. The fact that these appear at late time in the simulation shows the good quality of the numerical resolution of the equations offered by the Nek5000 solver. Long wave optimal perturbations: k = 0.8, ε = 10 -3 Figure 3.12 shows the evolution of the symmetric optimal perturbation computed for T = 5 and k = 0.8 superimposed to the baseflow with an initial amplitude of ε = 10 -3 . Quickly, the flow seem to develop sub-wavelengths due to the non linearity of the Navier-Stokes equations and the simulation becomes under-resolved. From figure 3.8, the short wave gain in kinetic energy is four orders of magnitude higher than the long wave ones. We remind the reader,

SHORT TIME MULTIPLE MODAL STATES

parition of what we call "spiral" and "whip" perturbations for values of p ∼ 1.6, see figure 4.8. On the left of figure 4.8, for p = 1.6, we observe a perturbation that is strictly m = 2. It is obtained by restarting the optimization with a previously computed m = 2 mode. Repeating this procedure for p < 1.6 yields m = 2 modes. On the right of figure 4.8 the computations are restarted with random white noise and the optimal perturbation resembles a whip and is a sum of multiple different modes. For p ≥ 1.7, the very same procedure exclusively yields similar whip modes. These results are summarized in figure 4.9 in which we plot the cost function as a function of p for T = 1. The whip (resp. spiral) branch is represented with a solid (resp. dotted) line. It exhibits a critical value, p = 1.6, that separates the dynamics into a well defined m azimuthal wavenumber type of optimal p-norm perturbation and a region beyond this critical value that exhibits a dissymetric perturbation not following an m azimuthal decomposition. This behavior occurs at small horizon time. It is found to disappear when increasing this horizon time. A similar observation was reported by [START_REF] Foures | Localization of flow structures using ∞-norm optimization[END_REF] where "C-branch" solutions were traced to values of p as low as p = 1.5 where they progressively merged with "periodic" solutions. The mechanism for this behavior seems to relate to the adjoint initialization. Following equation (4.3) the short time optimal is sensitive to low amplitude noise (from discretization, convergence, machine precision, etc.) that naturally arises in the numerical simulation and that is subjected to strong amplitifcation by the p-norm exponent, causing the observed breaking of azimuthal symmetry. However the fact that this whip mode is robust to the type of initialilization of the optimization procedure may be indicative of a more physical dynamics.

APPENDIX

A

Optimal Perturbations : Mesh Convergence Mesh refinement using spectral elements is convenient as it can be achieved in two ways. The first consist in increasing the total number of elements. The other consists in increasing the polynomial order, this results in exponential convergence of residuals. In Chapter 3, the polynomial order is P = 5 amounting to 6 3 = 216 grid points per element (gdp/elem).

A.1 Linear Optimal Perturbations

Increasing the polynomial order to 7 marginally modifies the results. Table A.1 summarized the results on mesh convergence for 2 extreme cases (respectively A and B) : A has a relatively low horizon time T = 1.25 and B is the most optimal perturbation T = 5, k z = 6.0.

Parameters

A.2 Direct Numerical Simulations

The case T = 5.0, k z = 6.0 has been considered to assess mesh convergence when running Direct Numerical Simulations.