J M Génevaux

N Dauchez

O Doutres

Non linear damping of a plate using Faraday instability of a fluid film

Keywords: Non linear damping, Damping of panels, Faraday instability, Fluid film, Vibroacoustics PACS: 43.55.Wk, 46.40.Ff, 68.35.J, 43.25

Damping using an instability of a fluid film in contact with a vibrating structure is investigated. Waves induced in the fluid film are the source of the added damping. A model based on the theory of Faraday instability is applied to a clamped circular plate covered by a fluid film. It is shown that this original technique can provide a significant damping, as with viscoelastic or porous material treatments. It is related to the amplitude of the waves which is a non linear function of the plate acceleration. Theoretical and experimental results are compared. The model overestimates the added damping: it is four times greater than the measured one.

Introduction

This paper examines a method to reduce the vibration and therefore the emitted noise of a structure by means of a fluid film in the low frequency range. The usual techniques for noise and vibration reduction from a structure use a viscoelastic layer bonded onto the structure [START_REF] Nashif | Vibrations damping[END_REF]. In this case, the dissipation is proportional to the loss factor of the material and to the flexural strain energy of the viscoelastic layer. To be efficient, this technique requires the use of a thick layer. Its thickness, and therefore the mass adding, can be reduced by using a light and stiff constraining sheet that increases the strain energy in the dissipating layer. Optimal partial covering may also be used to reduce the added mass [START_REF] Alvelid | Optimal position and shape of applied damping material Journal of Sound and Vibration[END_REF]. These techniques are limited by viscoelastic properties that depend on frequency and temperature [START_REF] Pritz | Loses factor peak of viscoelastic materials : magnitude to width relations[END_REF]. Designed primarily for sound absorbing, porous materials such as polymer foam may also add significant damping when mounted onto a structure [4][5][START_REF] Dauchez | Investigation and modelling of damping in a plate with a bonded porous layer[END_REF]. To improve the efficiency of Fig. 1. Oscillation of a fluid film on a vibrating plate, with w(x, y, t) the transverse displacement of a point P of the plate (coordinates in the plane (x, y)), t the time, h the water level of the fluid film at rest and ξ(x, y, t) the amplitude of the waves above the point P.

passive treatments, active control techniques have also been developed [START_REF] Baz | Robust control of active constrained layer damping[END_REF][START_REF] Guyomar | Nonlinear semi-passive multimodal vibration damping: An efficient probabilistic approach[END_REF] but require more sophisticated set up. Moreover, their robustness has to be carefully demonstrated.

In this paper, damping added by a fluid film in contact with a structure is investigated (Fig. 1). When the normal acceleration of the structure is strong enough, stationary waves appear in the fluid film (Fig. 2). This phenomenon is called Faraday instability [START_REF] Faraday | On the forms and states of fluids on vibrating elastic surfaces[END_REF][START_REF] Kityk | Spatio-temporal Fourier analysis, of Faraday surface wave patterns on a two-liquid interface[END_REF]. In case of a finite area of the fluid-air interface, the boundary conditions select countable wave lengths and several stationary mode shapes are solution of the problem [START_REF] Douady | Experimental study of the Faraday instability[END_REF]. The mode shape which has the greater amplification coefficient appears [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF][START_REF] Merlen | Duality of the supercritical solutions in magnetoacoustic wave phase conjugation[END_REF]. In case of an infinite area, this amplification coefficient will select the shape of the free surface among elementary cell patterns (roll, hexagon or square) (Fig. 2). These elementary cells can be considered as oscillators distributed over the surface of the structure [START_REF] Thompson | A continuous damped vibration absorber to reduce broadband wave propagation in beams[END_REF], which damping depends on fluid flow in a cell and on the viscosity of the fluid. Moreover, the relation between the wave amplitude and the driving acceleration of the plate is nonlinear: it is necessary to determine the acceleration threshold for waves to appear and their amplitude at saturation [START_REF] Douady | Experimental study of the Faraday instability[END_REF][START_REF] Milner | Square patterns and secondary instabilities in driven capillary waves[END_REF]. Note that for high acceleration level, ejection of droplets can be observed [START_REF] Alzuaga | Motion of droplets on solid surface using acoustic radiation pressure[END_REF]. This paper focuses on the added dissipation to the structure by the Faraday instability, using the thinnest fluid layer without droplet ejections. To the author's knowledge, this technique aimed at reducing the vibration using the Faraday instability has not been previously presented.

In a first part (Sec. 2), the modelling the method of calculating the added damping is detailed. To this end, a particular geometrical configuration is chosen: a circular plate clamped at its edge.

In section 3, the corresponding experiment is designed in order to measure the added damping of the first mode of the plate. This highlights the nonlinear The results of the model are compared with those obtained by the experiment in section 4.

Model

From local to global dissipation

This section details how to calculate the global dissipation added on the structure as function of the driving acceleration and for a given mode shape of the plate. Indeed, the modal dissipation is function of the area where the instabil-ity appears and of the nonlinear relation between the amplitude of the waves and the local acceleration.

In the present paper, the first mode of a circular clamped plate of diameter d is considered. Its mode shape is given by [START_REF] Geradin | Théorie des vibrations : application à la dynamique des structures Masson[END_REF],

φ(r, θ) = I 0 (β 01 d/2) J 0 (β 01 r) -J 0 (β 01 d/2) I 0 (β 01 r) I 0 (β 01 d/2) J 0 (0) -J 0 (β 01 d/2) I 0 (0) , ( 1
)
with r the distance of the point to the centre of the plate, I 0 as the modified Bessel function, J 0 the Bessel function and β 01 = 1.015 2π d . This mode shape is normalized so that φ(0, 0) = 1.

The modal damping ratio of the plate is given as function of modal parameters,

ζ ap = c ap 2 k p m pf , ( 2
)
where c ap is the damping coefficient added on the plate, k p the modal stiffness of the plate and m pf the modal mass of the plate loaded by the fluid film. These three terms c ap , k p and m pf are calculated in the two following sections.

Modal stiffness and modal mass

The strain energy of the plate is calculated by integrating the local strain energy over the plate. This expression gives the modal stiffness k p [START_REF] Geradin | Théorie des vibrations : application à la dynamique des structures Masson[END_REF]:

1 2 k p (φ(0, 0)) 2 = 1 2 D 2π θ=0 R r=0 φ ,rr + 1 r φ ,r 2 -2(1 -ν p ) φ ,rr φ ,r r rdrdθ,
(3) with D the stiffness of the plate defined by

D = Ee 3 /(2(1 -ν 2 p )), (4)
with E the Young modulus and ν p the Poisson's ratio of the material, e the thickness of the plate. This modal stiffness does not depend on the fluid film properties.

By the same approach, the modal mass of the plate m p is given by [START_REF] Geradin | Théorie des vibrations : application à la dynamique des structures Masson[END_REF]:

1 2 m p (φ(0, 0)) 2 = 1 2 2π θ=0 R r=0 ρ p eφ(r, θ) 2 rdrdθ, (5)
with, ρ p e = ρ s e + ρ f h = ρ s e(1 + ρ), (6) the equivalent mass per unit area of the system, based on ρ s the density of the aluminum, ρ f the density of the fluid and an added mass due to the fluid layer. In this equation, the kinetic energy of the fluid due to the flow relative to the plate is neglected. Assuming that ρ p is constant along the plate and that φ(0, 0) = 1, equation (5) can be rewritten,

m p = ρ p e 2π θ=0 R r=0 φ(r, θ) 2 rdrdθ. ( 7
)
The value of ρ s is given by fitting the first resonance frequency of the bare plate,

f 1 = 1.015 2 2π π 2 (d/2) 2 D ρ s e . ( 8
)

Added damping coefficient

The added damping coefficient c ap is determined considering the dissipated energy during one period of the plate vibration [START_REF] Faraday | On the forms and states of fluids on vibrating elastic surfaces[END_REF][START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF][START_REF] Matthiessen | Ueber die Transversalschwingungen tnender tropfbarer und elastischer Flssigkeiten[END_REF][START_REF] Rayleigh | On the crispation of fluid resting upon a vibrating support[END_REF][START_REF] Maksimov | Transient processes near the threshold of acoustically driven bubble shape oscillations[END_REF],

2π/ω e t=0 1 2 c ap ẇ2 (0, 0, t)dt = 2π/ω e t=0   1 2 2π θ=0 R r=0 ĉ ∂hα(r, θ) cos(ω f t) ∂t 2 rdrdθ   dt, (9)
with w(0, 0, t) = W sin(ω e t) the instantaneous transverse displacement of the centre of the plate, W the amplitude of the displacement at this point, ĉ the local damping coefficient per unit area, α(r, θ) = A(r, θ)/h (Eq. 12) the dimensionless amplitude of the waves , A(r, θ) = max t (ξ(r, θ, t) the amplitude of the waves at saturation (Fig. 2), ξ(r, θ, t) the instantaneous position of the free surface, ω f the circular frequency of the free surface waves due to the Faraday instability, ω e the circular frequency of excitation and t the time.

Taking into account the axisymmetry of the first mode shape, that ω f = ω e /2 and ẇ(0, 0, t) = W ω e sin(ω e t), the integration in time and in angle θ gives,

c ap W 2 = ĉh 2 π R r=0 α(r) 2 rdr. ( 10
)
The integration is made on the area where the instability appears. It is a disk of radius r d , so that for r > r d , α(r) = 0. Assuming that the amplitude of the waves does not depend on the acceleration of the other points of the plate (local hypothesis), r d may be defined by

ẅ(r d , θ) = ˜ c g, ( 11
)
with ˜ c the dimensionless acceleration threshold for the Faraday instability to appear and g the gravity.

For r < r d , a linear relation between α 2 and ˜ is accounted for as suggested in ref. [START_REF] Douady | Experimental study of the Faraday instability[END_REF]: the amplitude of the waves at saturation, when ∂A ∂t = 0, is an affine function of ˜ ,

α 2 = ˜ A ω f -B ω f , ( 12
)
where A ω f and B ω f are two parameters which will be determined experimentally in this paper (Sec. 3 and Fig. 7). Note that ˜ c is ˜ for α = 0 in equation ( 12):

˜ c = B ω f A ω f . ( 13
)
Equation ( 10) then writes

c ap = πĉh 2 W 2 A ω f g r d r=0 ˜ (r)rdr -B ω f r d r=0 rdr , ( 14
)
where the dimensionless acceleration ˜ (r) can be related to φ(r, θ) the normalized mode shape of the plate (Eq. 1) by

˜ (r) = W ω 2 e g φ(r, θ). ( 15
)

Equivalent damping coefficient of a fluid cell

The aim of this section is to calculate the modal damping coefficient ĉ of the fluid per unit area, which is required to calculate c ap (Eq. 14). It is given by,

ĉ(r, θ) = 2ζ 0 m f ω f , ( 16
)
with m f the modal mass per unit area and ζ 0 the damping ratio of the fluid. ζ 0 can be related to the logarithmic decrement α s of free oscillations [START_REF] Geradin | Théorie des vibrations : application à la dynamique des structures Masson[END_REF] by

ζ 0 = α s 4π 2 + α 2 s , ( 17
)
with α s = 2πδk [START_REF] Milner | Square patterns and secondary instabilities in driven capillary waves[END_REF]. Here δ = 2ν/ω f is the thickness of the viscous layer for a fluid with a kinematic viscosity ν solicited at the circular frequency ω f , and k the wave number. This wave number is solution of [START_REF] Génevaux | Gravity effects on coupled frequencies of a 2D fluid-structure problem with free surface[END_REF]]

ω 2 f = gk tanh(kh) 1 + σk 2 ρ f g , ( 18
)
with σ the surface tension of the fluid-air interface.

The modal mass per unit area m f is

m f = m f cell (λ/2) 2 , ( 19
)
with m f cell the modal mass of one cell whose area is (λ/2) 2 , with λ = 2π/k the wave length. Here, the first mode of a fluid cell for a given pattern of the free surface is considered. Several patterns appear successively on the free surface with the increase of the level of the acceleration [START_REF] Kityk | Spatio-temporal Fourier analysis, of Faraday surface wave patterns on a two-liquid interface[END_REF]: roll, hexagon, then square. In this paper, the square pattern which is present for high accelerations is chosen (Fig. 2). The modal shape can be defined according to a potential function for the flow: the Reynolds is large enough (Re = 1.7 10 6 ) for the thickness of viscous boundary layer (δ = 0.063 mm) to be smaller than the fluid depth. Let us consider a fluid cell defined by a volume of incompressible fluid whose dimensions are λ/2 in x and y directions, and h in z direction. Symmetry arguments allow the use of a flow v relative to the plate associated to Φ f (x, y, z, t) = φ f (x, y, z)g(t), with g(t) a harmonic function and

φ f (x, y, z) = β(r, θ) cos(q(2x/λ -1)) cos(q(2y/λ -1)) cosh( √ 2q, z/λ), ( 20
)
with β the amplitude of the potential function which depends on the position of the fluid cell (cylindrical coordinates r,θ), and (x,y,z) the coordinates of a point into the cell. The velocity field v is given by,

h ∂φ f ∂x g(t) = v.x, h ∂φ f ∂y g(t) = v.y, h ∂φ f ∂z g(t) = v.z. (21)
The boundary conditions at x = λ/2 and y = λ/2 give for the first mode of the fluid cell, q = π.

The kinetic energy of the fluid in the cell is used to calculate m f cell :

1 2 m f cell ∂α(r,θ) cos(ω f t) ∂t 2 = 1 2 ρ f g(t) 2 λ x=0 λ y=0 h z=0 φ 2 f,x + φ 2 f,y + φ 2 f,z dxdydz, ( 22
)
with φ f,x = ∂f ∂x . β(r, θ) is derived from α(r, θ) the amplitude of the waves, by equating the vertical velocity of a point of the free surface (z = h):

α(r, θ) ∂ cos(ω f t) ∂t = g(t) ∂φ f ∂z ∀t, ( 23
)
and using equation ( 20) it gives,

α(r, θ) = β(r, θ) 2 √ 2π λ sinh( √ 2πh/λ), ( 24
)
g(t) = -ω f sin(ω f t) (25) so that, β(r, θ) = α(r, θ)λ 2 √ 2π sinh( √ 2πh/λ) . ( 26
)
Using equation ( 26), and integrating the second member of equation ( 22), the modal mass of the 1 degree of freedom system associated to the fluid cell writes,

m f cell = ρ f λ 3 √ 2 256π 1 -exp -8 √ 2πh λ exp 4 √ 2πh λ cosh √ 2πh λ 2 -1 . ( 27
)
The modal mass per unit area m f and the modal damping coefficient by unit area ĉ may be obtained with equations ( 16) and ( 19).

Experimental results

The aim of this section is to determine experimentally the plate damping induced by the Faraday instability of the fluid film. Because the system is nonlinear, all experimental quantities are referenced to the acceleration at the centre of the plate, called driving acceleration, and denoted W a (0, 0). Four water levels and four values of the driving acceleration are considered.

Experimental set-up

The geometrical configuration is presented in figure 3. A circular aluminum plate of diameter d = 0.290 m is clamped. Its characteristics are summarized in tables 1 and 2. It is excited by a shaker connected to its centre via a force transducer giving the excitation force F . An accelerometer of mass m a = 0.0042 kg is bonded at 1 cm of the centre of the plate to get the reference acceleration W a (0, 0) = ω 2 e W . A laser vibrometer is focused on the plate through the fluid film. The laser spot area is less than gives a measured area of the order of 1 mm 2 . The frequency range is set between 40 Hz and 120 Hz to ensure Faraday instabilities for experimentally achievable accelerations. The thickness of the fluid layer is chosen to be of the order of the wavelength: this ensures a total covering of the plate and a small added mass.

The frequency response functions are obtained using a step by step harmonic excitation for a constant acceleration amplitude W a (0, 0).

Determination of the free surface velocity of the fluid

The signal given by the laser vibrometer is proportional to the apparent velocity of the plate d w dt . It is a combination of the plate and of the free surface velocities:

d w(r, θ, t) dt = dw(r, θ, t) dt - c 0 c 1 -1 d(h + ξ(r, θ, t)) dt , (28)
with w(r, θ, t) the displacement of the plate at the point P (r, θ), c 0 the light velocity in the air, c 1 the light velocity in the fluid, h(r, θ, t) the water level crossed by the beam. The time signal of the apparent speed of the plate (Fig. 4) exhibits a sub-harmonic at ω e /2. Their contributions can be easily identified because they are at different frequencies (see Fig. 4). To do so, equation ( 28) is rewritten as:

d w(r, θ, t) dt = ã1 cos(ω e t + ã2 ) + ã3 cos ω e 2 t + ã4 . ( 29
)
The parameters ã1 and ã2 are associated to the plate velocity, and the parameters ã3 and ã4 to the fluid interface movement. They are determined using a nonlinear optimization. The position of the nodes of the free surface mode shape are not stationary: ã3 = 0 if the laser beam focuses through a node, and is maximum if the laser beam focuses through an antinode. Indeed, for level of accelerations far higher the threshold, second order effects [START_REF] Martin | Drift instability of standing Faraday waves[END_REF] induce a slow drift of the positions of the antinode of the fluid air interface. Thus, the time of acquisition must be long enough to contain measurement on an antinode. The amplitude of the waves is then deduced from the greatest value of ã3 .

Influence of the driving acceleration and water levels

The experiments are performed with water to test the influence of the fluid level on the frequency response functions. The water levels remain below 0.0133 m and above 0.004 m for all the experiments. Below this lower limit, the wetting mechanism prevents the water from covering the whole plate. The ratio of added mass per unit area is:

ρ = ρ f h ρ s e . ( 30
)
Near the resonance frequency (75 Hz), stationary waves are observed on the free surface (Fig. 2). Their amplitude is higher at the centre of the plate. This coincides with the maximum amplitude of the first mode shape. They appear as square patterns within a circle at the centre of the plate. Their wavelength, which are consistent with capillarity waves on the surface, increases as the frequency decreases.

To determine the acceleration threshold for Faraday instability, the frequency response function is measured for several amplitudes of the driving acceleration for water levels from 0 to 0.0133 m. Nonlinear effects are observed (Fig. 5) when increasing the amplitude of acceleration:

• the existence of a threshold for waves to appear,

• an increase in the area of the circular surface on which waves are present,

• a slight increase in the frequency resonance,

• a reduction of the amplitude of the peak (up to 13 dB).

In figure 5, the threshold is between 6 m s -2 and 13 m s -2 . It corresponds to a sharp increase in damping. Figure 6 shows damping ratio ζ pf which is estimated using the half-power method. The measured damping for this non optimized configuration can be half the value of the damping measured in the case of unconstrained-layer treatment, and is greater than the damping added by foam or fiber layer. This increase is not observed for the bare plate (continuous line). The thinner the fluid layer, the stronger the damping is.

Comparison between theoretical and experimental results

Comparison between theoretical and experimental values of the plate damping induced by the fluid film is now evaluated. Water with the thinner thickness is used to induce the stronger damping. The numerical values of the parameters of the model are given in tables 1 and 2 for this configuration.

The added damping ratio depends on the local dissipation of the fluid layer and on the area of the instability area.

This local dissipation is a function of the relation between the local acceleration of the plate and of the wave amplitude (Fig. 6). The relation between the wave amplitude and the motion of the plate is quantified by the coefficients A ω f and B ω f (Eq. ( 12)) which are identified on the experimental data by a least-squares regression method (Fig. 7). The uncertainties on the position of the regression line are calculated with a confidence of 95%. Note that three points are not taken into account: the noise on the signal and the low spatial stability of the waves do not allow a correct evaluation of the amplitude of the waves for these three acceleration levels. In the following, above the threshold ˜ c , the dimensionless amplitudes of the waves lie in the following boundaries: 0.0135(˜ -0.82) < α 2 < 0.0163(˜ -0.82).

The wave length λ = 0.004 m given by the model using σ = 0.072 N m -1 the surface tension of the fluid-air interface, is near the measured experimental wave length (0.00626 m< λ < 0.00634 m, Sec. 3.3). The experimental average value λ = 0.0063 m is used in the model to evaluate the damping.

Moreover, the equivalent damping coefficient c ap added to the structure depends on the area on which the damping occurs. The first mode shape of the plate induces a circular limit of instability of radius r d where, W a (r d , θ) = 8.04 m s -2 . The relationship between this radius and the acceleration in the centre of the plate is given by the mode shape of the plate φ(r, θ) (Eq. 1) and the instability threshold. The evolution of c ap , as function of W a at the centre of the plate, can then be calculated taking into account the uncertainty of the amplitude of the waves (Eq. 31). From the evolution of c ap , the damping ratio added by the fluid ζ ap can be plotted (Fig. 8). For the experimental point of view, the damping ratio added by the fluid can be deduced by the comparison of the damping ratio observed with the fluid ζ pf and the damping ratio of the bare plate ζ p (Fig. 6). The damping ratio of the bare plate is,

ζ p = c p 2 k p m p , ( 32
)
with c p the modal damping coefficient of the plate, m p the modal mass of the plate. With the fluid layer the damping ratio is given by,

ζ pf = c p + c ap 2 k p m pf = c p 2 k p m pf + ζ ap . ( 33
)
By neglecting the kinetic energy of the waves compared to the kinetic energy of the loaded plate, the relation m pf = m p (1 + ρ) gives,

ζ pf = c p 2 k p m p (1 + ρ) + ζ ap . ( 34
)
Thus, the experimental value of the added damping ratio is The comparison (Fig. 8) between experimental and theoretical values of ζ ap shows that above the threshold the amplitude of the theoretical damping is four times greater than the experimental damping.

ζ ap = ζ pf -ζ p 1 √ 1 + ρ . ( 35)
The source of this discrepancy can be sought in the high sensitivity of the model to the amplitudes of the waves. The relation between the wave amplitude and the acceleration has been determined assuming a locally reacting behaviour of each fluid cell. This hypothesis is valid where the acceleration is uniform: this is the case in the vicinity of an antinode (centre of the plate for the first mode). This assumption may not be valid where an acceleration gradient is present. Moreover, in this transition region, the cell shape is no longer square (Fig. 2).

Nevertheless, note that the damping added by the fluid film is comparable to thoses obtained with viscoelastic or porous material treatments [START_REF] Nashif | Vibrations damping[END_REF][START_REF] Dauchez | Investigation and modelling of damping in a plate with a bonded porous layer[END_REF].

Conclusion

Damping of a vibrating structure by means of a heavy fluid film subjected to Faraday instability has been studied. It is shown that added damping may be

Fig. 2 .

 2 Fig. 2. Stationary waves at the fluid-air interface (top view of the plate). Driving acceleration W a (0, 0) : 13.8 m s -2 behaviour of the instability of the fluid film and the influence of the parameters governing the phenomenon.

Fig. 3 .

 3 Fig. 3. Experimental set-up.

Fig. 4 .

 4 Fig. 4. Signals (V) of the acceleration of the plate (dotted line 31.6 m s -2 / V) and its apparent velocity measured by the laser vibrometer (bold line 0.025 m s -1 / V) maximum velocity of the plate : 0.039 m s -1 ; water level 0.0046 m ; ρ = 0.88 ; frequency of excitation 65 Hz.

Fig. 5 .Fig. 6 .

 56 Fig.5. Influence of the driving acceleration W a (0, 0) on the frequency response function (accelerometer/force) for a water level of 0.0046 m, ρ = 0.88: continuous line, W a (0, 0) = 3.16 m s -2 ; dashed line, W a (0, 0) = 6.32 m s -2 ; dotted line, W a (0, 0) = 12.7 m s -2 ; dashed dotted line, W a (0, 0) = 25.3 m s -2 .

Fig. 7 .

 7 Fig. 7. Experimental evolution of the dimensionless oscillation amplitude α = A/h with the acceleration above the threshold for a given frequency (75 Hz) with water.The two lines corresponds to the uncertainty with a confidence of 95% of the mean square line which fits the data.

Fig. 8 .

 8 Fig. 8. Evolution of the damping ratio added by the fluid versus the driving acceleration W a (0, 0): theoretical values (between the two lines: the uncertainties on the measurement of the water waves induce two curves of dimensionless damping with 95% of confidence), experimental values (crosses).

Table 1

 1 Numerical inputs of the model.

	Plate
	Diameter of the plate (m)	d	0.290
	Thickness of the plate (m)	e	0.001
	Material of the plate		aluminum
	Young's Modulus (Pa)	E	7 10 10
	Poisson's ratio	ν p	0.3
	Frequency of the first resonance
	of the bare plate (Hz)	f 1	100
	Frequency of excitation (Hz)	f e	75
	Fluid
	Gravity (ms -2 )	g	9.81
	Fluid layer thickness (m)	h	0.0046
	Kinematic viscosity of water (m 2 s -1 )	ν	10 -6
	Density of the fluid (kg m -3 )	ρ f	1000
	Wave length of the free surface (m)	λ	0.0063
	First amplitude of the wave coefficient	A ω f [0.0135, 0.0163]
	Second amplitude of the wave coefficient B ω f [0.0111, 0.0134]

Table 2

 2 Numerical outputs of the model.

	Loaded plate
	Stiffness of the plate (N m)	Eq. (4)	D	8.82
	Dimensionless added mass	Eq. (30)	ρ	0.88
	Density of the plate with air loading (kg m -3 )	Eq. (8)	ρ s	5230
	Modal stiffness (Nm -1 )	Eq. (3)	k p	1.08 105
	Equivalent mass of the plate with fluid (kgm -3 )	Eq. (6)	ρ p	7215
	Modal mass (kg)	Eq. (7)	m p	0.0875
	Fluid
	Threshold of instablility (ms -2 )	Eq. (13)	˜ c g	8.04
	Local damping ratio of the fluid	Eq. (17)	ζ 0	14.5%
	Modal mass of a fluid cell (kg)	Eq. (27) m f cell 2.23 10 -6
	Modal mass per unit area (kg m -2 )	Eq. (19)	m f	2.25 10 -1
	Equivalent damping coefficient per unit area (kg m -2 s -1 ) Eq. (16)	ĉ	9.69
	Damping ratio added by the fluid	Eq. (2)	ζ ap	Fig. 6

ACKNOWLEDGMENTS

Thanks to Lazhar Benyahia from the Laboratoire Polymères, Colloides, Interfaces of Université du Maine (Le Mans, France) for our productive discussions. Thanks to Denis Ritter for the English corrections.

List of Figures

Oscillation of a fluid film on a vibrating plate, with w(x, y, t) the transverse displacement of a point P of the plate (coordinates in the plane (x, y)), t the time, h the water level of the fluid film at rest and ξ(x, y, t) the amplitude of the waves above the point P. Influence of the driving acceleration W a (0, 0) on the frequency response function (accelerometer/force) for a water level of 0.0046 m, ρ = 0.88: continuous line, W a (0, 0) = 3.16 m s -2 ; dashed line, W a (0, 0) = 6.32 m s -2 ; dotted line, W a (0, 0) = 12.7 m s -2 ; dashed dotted line, W a (0, 0) = 25.3 m s -2 . 11 6 Influence of the driving amplitude W a (0, 0) on the plate damping for several water levels: continuous line, 0.0 m ; dashed line, 4.6 mm ; dotted line, 8.0 mm ; dashed dotted line, 13.3 mm. 12 7 Experimental evolution of the dimensionless oscillation amplitude α = A/h with the acceleration above the threshold for a given frequency (75 Hz) with water. The two lines corresponds to the uncertainty with a confidence of 95% of the mean square line which fits the data. [START_REF] Milner | Square patterns and secondary instabilities in driven capillary waves[END_REF] 8 Evolution of the damping ratio added by the fluid versus the driving acceleration W a (0, 0): theoretical values (between the two lines: the uncertainties on the measurement of the water waves induce two curves of dimensionless damping with 95% of confidence), experimental values (crosses) 16

List of Tables 1 Numerical inputs of the model. [START_REF] Merlen | Duality of the supercritical solutions in magnetoacoustic wave phase conjugation[END_REF] 2 Numerical outputs of the model. 14


Juliette Pottier

Maria Mendez Real

Sébastien Pillement

Sébastien Pillement First

First analysis and protection of the micro-architecture of a RISC-V core

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

FirstWP 4 -•

 4 INSTITUT D'ÉLECTRONIQUE ET DE TÉLÉCOMMUNICATIONS DE RENNES SEC-V Project: Methodology and organization WP 1 -Dynamic code transformation unit On-the-fly decoding modification/alteration Dynamic instrumentation Instruction set tailoring/customization/adaptation WP 2 -Micro-architectural modifications Alternative approaches to traditional caches (scratchpads, TCM) Dynamic cache management Inserting execution noise (access instructions for example) stage pipeline partially out-of-order (Execute Stage) o Single issue Micro-architecture analysis • Vulnerabilities targeted by side-channel attacks and covert channels: • Kocher, Paul, et al. "Spectre attacks: Exploiting speculative execution. "Communications of the ACM"(2020) • Ge, Qian, et al. "A survey of microarchitectural timing attacks and countermeasures on contemporary hardware." Journal of Cryptographic Engineering (2018) Prototype and evaluation Inclusion in the CVA6 core of the OpenHW Group Assessment (indicators and metrics) of security levels Security strategies: ➢ Monitoring : detection of contexts favorable to side-channels attacks and/or covert channels ➢ Dynamic management : micro-architectural defenses and micro-decoder ➢ Deployment of a complete solution on target, while preserving performance • Data cache :


Chunyao Ma

Sameer Alam

Qing Cai
email: qcai@ntu.edu.sg

Daniel Delahaye
email: daniel.delahaye@enac.fr

Dynamic Air Traffic Flow Coordination for Flow-centric Airspace Management

Keywords: traffic flow coordination, flow-centric operation, air traffic flow prediction, transformer neural networks

The air traffic control paradigm is shifting from sectorbased operations to cross-border flow-centric approaches to overcome sectors' geographical limits. Under this paradigm, effective air traffic flow coordination at flow intersections is crucial for efficiently utilizing available airspace resources and avoiding inefficiencies caused by high demand. This paper proposes a dynamic air traffic flow coordination framework to identify, predict, assess, and coordinate the evolving air traffic flows to enable more efficient flow configuration. Firstly, nominal flow intersections (NFI) are identified through hierarchical clustering of flight trajectory intersections and graph analytics of daily traffic flow patterns. Secondly, spatialtemporal flow features are represented as sequences of flights transiting through the NFIs over time. These features are used to predict the traffic demand at the NFIs during a given future period through a transformer-based neural network. Thirdly, for each NFI, the acceptable flow limit is determined by identifying the phase transition of the normalized flight transition duration from its neighboring NFIs versus the traffic demand. Finally, when the predicted demand exceeds the flow limit, by evaluating the available capacity at different NFIs in the airspace, the flow excess is alternated onto other NFIs to optimize and re-configure the air traffic demand to avoid traffic overload. An experimental study was carried out in French airspace using the proposed framework base on the ADS-B data in December 2019. Results showed that the proposed prediction model approximated the actual flow values with the coefficient of determination (R 2 ) above 0.9 and mean absolute percentage error (MAPE) below 20%. Acceptable flow limit determination showed that for above 68% NFIs, the flight transition duration increases sharply when the demand exceeds a certain level. The flow excess at an NFI whose demand was predicted to exceed its limit was coordinated, and the potential increase in the flight transition duration caused by the flow excess was avoided.

I. INTRODUCTION

The scalability limit of traditional sector-based Air Traffic Control (ATC) services, i.e., difficulty in subdividing heavily loaded sectors, is becoming a barrier to the sustainable growth of air traffic [START_REF] Gianazza | Learning air traffic controller workload from past sector operations[END_REF]. Researchers have started examining and testing the concept of sectorless ATC, which views the airspace as a whole instead of the current practice of dividing the airspace into small sectors. One primary practice of sectorless airspace is flow-centric operation [START_REF] Zeki | Business models for flight-centric air traffic control[END_REF], which relies on controlling and monitoring flow-based formation and evolution of air traffic, i.e., the management of dynamic flow corridors [START_REF] Reitmann | Advanced quantification of weather impact on air traffic management[END_REF]. It opens the opportunity to distribute air traffic more efficiently in the airspace without being constrained by sector boundaries [START_REF] Undertaking | European atm master plan: the roadmap for delivering high performing aviation for europe: executive summary[END_REF].

Despite the benefits of the flow-centric concept, its implementation has been limited. One primary challenge is the efficient coordination of air traffic flow at the intersections to avoid inefficiencies that may jeopardize flight safety [START_REF] Gerdes | From free-route air traffic to an adapted dynamic mainflow system[END_REF]. Research focusing on the traditional sector-based air traffic coordination, such as sector traffic prediction and flow optimization for workload balancing between sectors [START_REF] Chen | A network based dynamic air traffic flow model for en route airspace system traffic flow optimization[END_REF], no longer adapts the flow-centric operations where coordination is primarily used to avoid potential inefficiencies or conflicts between the intersecting air traffic flows [START_REF] Undertaking | European atm master plan: the roadmap for delivering high performing aviation for europe: executive summary[END_REF]. Therefore, for safe and efficient airspace management under flow-centric operations, it is crucial to develop a flow-centric-based air traffic coordination framework that can dynamically coordinate air traffic flow in advance based on the collaborative identification, prediction, and inefficiency assessment of the evolving air traffic flows. For instance, traffic flow can be re-routed in advance when the predicted flow demand exceeds the acceptable flow limit at the intersections [START_REF] Undertaking | European atm master plan: the roadmap for delivering high performing aviation for europe: executive summary[END_REF].

Effective air traffic flow identification is the cornerstone for flow-centric-based practices regarding traffic flow analysis, prediction, and coordination [START_REF] Gerdes | From free-route air traffic to an adapted dynamic mainflow system[END_REF]. In the literature, air traffic flow has been identified and described in accordance with the airspace configuration, such as the behaviors of groups of flights transiting through area control centers [START_REF] Sridhar | Modeling and optimization in traffic flow management[END_REF], waypoints [START_REF] Murca | Identification, characterization, and prediction of traffic flow patterns in multi-airport systems[END_REF], sectors [START_REF] Ma | Sector entry flow prediction based on graph convolutional networks[END_REF], and airways [START_REF] Zhang | A hierarchical heuristic approach for solving air traffic scheduling and routing problem with a novel air traffic model[END_REF]. Such a characterization of air traffic flow serves the traditional air traffic operations where air traffic control units are geographical sectors and flights have to fly along airways consisting of a set of fixed waypoints. However, flow-centric air traffic management focuses on the flow of air traffic in the airspace from a holistic view disregarding the fixed airways and sectors. It requires air traffic flow identification methods to explore the underlying flow patterns, such as the spatial-temporal evolution of flow locations and structures [START_REF] Olive | Identifying anomalies in past en-route trajectories with clustering and anomaly detection methods[END_REF].

In addition to air traffic flow identification, effective air traffic flow coordination requires the constant viewing of air traffic demand according to the available capacity at the flow intersections. Air traffic flow prediction models in the literature primarily take the time series information of the number of flights at single or multiple geographical locations in the airspace as the input and predict the future number of flights through Linear Dynamic System Models (LDSM) or neural networks such as Long Short-Term Memory (LSTM) [START_REF] Gui | Machine learning aided air traffic flow analysis based on aviation big data[END_REF] and Convolutional Neural Networks (CNNs) [START_REF] Lin | Deep learning based short-term air traffic flow prediction considering temporal-spatial correlation[END_REF]. Such time-series-based methods mainly predict future traffic demand by analyzing the past time series [START_REF] Dalmau-Codina | A machine learning approach to predict the evolution of air traffic flow management delay[END_REF]. Important spatial-temporal flow dynamics were not incorporated, limiting the modeling and prediction accuracy, which is critical for flow-centric operations where real-time decision-making based on accurate data is necessary. Moreover, the traditional sector-based capacity estimation method lacks consideration of the air traffic flow features, such as the Monitor Alert Parameter(MAP) model regarding the hand-off service workload [START_REF] Marr | Controller workload-based calculation of monitor alert parameters for en route sectors[END_REF] and the sector merge/split-based model [START_REF] Gianazza | Learning air traffic controller workload from past sector operations[END_REF]. Therefore, developing flow prediction and acceptance limit identification methods at the flow intersections concerning the evolving flow features is the building block for flow-centric operations.

Information on future air traffic demand and the corresponding acceptable flow limit at the flow intersections enables reconfiguring the air traffic demand transiting through the flow intersections in advance so that air traffic flow can be restricted within a level that does not overload the system excessively [START_REF] Lin | From aircraft tracking data to network delay model: A datadriven approach considering en-route congestion[END_REF]. In traditional sectors, air traffic overload is managed by sector operations such as merging and splitting [START_REF] Liu | A framework for strategic online en-route operations: Integrating traffic flow and strategic conflict managements[END_REF]. Under the flowcentric paradigm, traffic flow density and complexity change over time, rendering static flow control operations underloaded or overloaded during the day. Analogically, dynamic flow coordination according to the evolving airspace conditions, such as traffic flow merge/split/re-routing, gives flow-centric airspace an option to address the anticipated flow excess without compromising the flow demand and overloading the acceptable limit at the flow intersections. Therefore, developing dynamic flow coordination algorithms depending on the time-varying traffic flows is an important enabler of efficient flow management decisions and optimal utilization of airspace resources.

In view of the above analysis, this paper proposes a dynamic air traffic flow coordination framework to identify, predict, assess, and coordinate the evolving air traffic flows to enable more efficient flow-centric airspace management. Firstly, nominal flow intersections (NFI) are identified through hierarchical clustering of flight trajectory intersections and graph analytics of daily traffic flow patterns. Secondly, based on the identified NFIs, spatial-temporal flow features are represented as sequences of flights transiting through the NFIs over time. These features are used to predict the future traffic demand at the NFIs through a transformer-based neural network. Thirdly, for each NFI, the acceptable flow limit is determined by identifying the phase transition of the normalized flight transition duration from its neighboring NFIs versus the traffic demand during different periods. Finally, when the predicted future flow exceeds the flow limit, by evaluating the available capacity at different NFIs, the flow excess is alternated onto other NFIs to re-configure the air traffic demand to avoid traffic overload at the NFIs.

II. PROBLEM DESCRIPTION

This paper focuses on the problem of dynamic air traffic flow coordination at the Nominal Flow Intersections (NFIs) for efficient Flow-centric Airspace Management. This problem can be further decomposed into four sub-problems as shown in Fig. 1: 1) identification of NFIs through graph analysis of air traffic flow patterns; 2) flow representation and prediction based on the spatial and temporal dynamics feature of air traffic flows at the NFIs; 3) flow acceptance limit identification at the NFIs through flow transition efficiency analysis; 4) Flow coordination, i.e., flow-excess re-routing, based on the capacity availability and predicted flow demand at the NFIs.

1) NFI Identification: Given the freedom of airspace users taking direct or user-preferred flight routes, one may argue that flow-centric implies insignificant spatial structure and temporal patterns of the air traffic flow. However, taking into account that all flights have to depart and land at airports, the positions of the airports and the scheduled flights between airports restrict the air traffic to an appropriate pattern of main flows. A traffic flow analysis of French airspace, where the free route airspace (a potential coupled working method to flow-centric operations) has been implemented in nearly 50% of the airspace above flight level 195, shows that the majority of flight trajectories are aggregated as major flows connecting major traffic hubs [START_REF] Ma | Air traffic flow representation and prediction using transformer in flow-centric airspace[END_REF]. Thus, the first sub-problem is identifying the locations and the flow inter-connections of the NFIs through constructing and analysing spatial and temporal patterns of air traffic flows.

2) Flow Representation and Prediction: Most existing methods in the literature represent the spatial-temporal air traffic flow features as a time series of the number of flights at different locations without considering the spatial and temporal dynamics feature within the air traffic flow, such as the spatial and temporal distributions of flights in the airspace. The spatial distribution of flights is vital for determining the air traffic complexity and density, while the temporal distribution is essential for describing the dynamic evolution of air traffic [START_REF] Delahaye | Air traffic complexity map based on linear dynamical systems[END_REF]. They are the primary influencing factors of the future air traffic at the NFIs, without which the prediction model may fail to learn the causal relations between the input air traffic feature and the prediction target, limiting the model's accuracy in predicting future air traffic flows. Therefore, the second sub-problem is the representation of the spatial-temporal flight distribution and accurate prediction of the future number of flights transiting through the NFIs.

3) NFI Flow Acceptance Limit Identification: Due to the increasing demand for air travel, an NFI may be overloaded by excessive air traffic flows to be coordinated. In a traditional sector-based ATC system, the sector capacity is usually quantified as the maximum number of flights that may enter a sector per hour averaged over a sustainable period and is used to manage a safe, orderly, and efficient traffic flow [START_REF] Welch | Sector workload model for benefits analysis and convective weather capacity prediction[END_REF]. Similarly, in flow-centric ATC systems, if the demand at an NFI exceeds an acceptable limit, the efficiencies of air traffic flow transiting through the NFI can be degraded. When an NFI is overloaded, air traffic congestion can be induced, and it will take a longer time for air traffic transiting to the overloaded NFI on the flight paths due to regulatory measures such as vectoring and speed control [START_REF] Baumgarten | The impact of hubbing concentration on flight delays within airline networks: An empirical analysis of the us domestic market[END_REF]. Therefore, the third sub-problem is to identify a flow acceptance limit at each NFI above which extra traffic demand will degrade the flow transition efficiency, i.e., the flight transition duration from neighboring NFIs, to the NFI significantly.

4) NFI Flow Coordination: Through solving the above three sub-problems, the locations, the inter-connectivity, the future traffic demand, and the flow acceptance limit of the NFIs can be determined. With this information in hand, the fourth subproblem is to come up with a flow-centric air traffic excess re-routing algorithm in advance to re-configure the air traffic demand transiting through the NFIs, i.e., re-routing the excessive air traffic flow onto NFIs with spare flow acceptance capability, to avoid exceeding the flow acceptance limit at the NFIs which can cause inefficiencies.

III. METHODOLOGY

A. NFI Identification

This paper identifies the NFIs based on the flight trajectories analysis using ADS-B data, including intersection points clustering, daily flow pattern representation, and graph analysis.

1) Intersection Points clustering & Flow Pattern Representation:

The NFIs in this paper are defined as the positions in the airspace where air traffic flows intersect. This paper proposes to determine the NFIs by clustering the intersection points of flight trajectories to extract the natural groupings of the intersection points. Hierarchical clustering relies on the hierarchical decomposition of the data based on group similarities to find a multilevel hierarchy of clusters, where clusters at one level are joined as clusters at the next level [START_REF] Cai | Hierarchical clustering of bipartite networks based on multiobjective optimization[END_REF]. Air transportation networks are commonly a nodal hierarchy that follows the spoke-hub structure in which air traffic flows range from regional feeders to international hubs [START_REF] Seabra | Determinants of brazilian international flights: The role of hub-and-spoke and infrastructure variables[END_REF]. Therefore, this paper adopts the hierarchical clustering algorithm to discover the NFIs for flow representation.

To determine the optimal number of clusters that should be identified by the clustering algorithm, this paper proposes to model the daily air traffic flow patterns as graphs based on the clustering outcome. The constructed graphs across different days are further investigated and compared to determine the best-fitted number of NFIs for flow pattern representation. The daily air traffic flow pattern is represented as a weighted graph G = (V, E), where V is the set of nodes denoting the air traffic flow components, i.e., the NFIs. The flow connectivity between the nodes can be described by the weighted edges E. The weight of each edge is quantified as the air traffic volume transiting through it, i.e., the number of flights whose trajectories consecutively cross the two nodes connected by the edge.

2) Graph Analysis: The consistent and dependable performance of airspace users is an essential requirement for improving ATM system predictability [START_REF] Murc ¸a | Flight trajectory data analytics for characterization of air traffic flows: A comparative analysis of terminal area operations between new york, hong kong and sao paulo[END_REF]. Therefore, the representation of traffic flow patterns should describe the behavioral consistency of air traffic flows in the airspace even though there are daily alternations in the geographical positions of the NFIs and the traffic flow connectivity between NFIs. Under this consideration, this section proposes to determine the optimal number of clusters by modeling and analyzing the consistency of daily flow patterns versus the changes in the number of clusters. The flow pattern consistency is evaluated from two perspectives: a) geographical consistency in node (NFI) locations; b) structural consistency in flow connectivity.

a) Geographical Consistency in NFI locations: To measure the geographical consistency in the daily node locations, a nearest neighbor-based analysis is conducted to match the nodes on the temporal horizon. Given a number of n nodes in the graph Based on the distances, the nearest neighbouring node v ai k+1 of node v i k for i = 1, ..., n can be determined. Similarly, the nearest neighbouring node v bj k of node v j k+1 for j = 1, ..., n can be determined with the same process. b j is the index of the identified nearest node of v j k+1 in V k . Therefore, two sets of matched node pairs between the two graphs can be obtained:

G k constructed for day D k , let V k = {v 1 k , v 2 k , ..., v i k , ..., v n k } represent the set of nodes. Similarly, let V k+1 = {v 1 k+1 , v
S 1 : (v 1 k , v a1 k+1 ), (v 2 k , v a2 k+1 ), ..., (v n k , v an k+1 ) and S 2 : (v b1 k , v 1 k+1 ), (v b2 k , v 2 k+1 ), ..., (v bn k , v n k+1 ).
Then, the geographical consistency of node locations is quantified as the number of mutually matched nodes divided by the total number of nodes, which is formulated as

gc 1 = |S 1 ∪ S 2 |/n,
gc 2 = l-1 i=1 l j=i+1 min(w i j c k , w i j c k+1 ) n-1 i=1 n j=i+1 (w i j k + w i j k+1 ) (1)
Fig. 2 shows a diagram of the graph analysis to determine the node (NFI) location and flow structure consistency. The first column of the figure shows the graphs G 1 and G 2 constructed for the air traffic flow on day 1 and day 2, respectively. The second column of the figure shows the enlarged portions of G 1 and G 2 . The nodes, i.e., V 1 and V 2 , are marked by solid gray circles. The third column of the figure depicts the mutually paired nodes C 1 and C 2 in the node location consistency analysis of graphs G 1 and G 2 . C 1 and C 2 are marked by the red solids, and the paired nodes are labeled by the same number. The fourth column of Fig. 2 displays the two sub-graphs of G 1 and G 2 formed by C 1 and C 2 . The numbers on the edges are the edge weight W c 1 and W c 2 . The fifth column of Fig. 2 shows the determination of mutual connections in the two sub-graphs, i.e., min{W c 1 , W c 2 }. min{W c 1 , W c 2 } represents taking the smaller value of the weights on the corresponding edges in the two sub-graphs. For instance, on day 1, the weight on the edge connecting nodes 44 and 166 is 380, meaning there are 380 flights transiting through this edge. While on day 2, the edge weight is 416. Thus, the mutual flow connection on this edge during the two days is considered 380. Eventually, the flow structure consistency is calculated by summing up the mutual connections in the two sub-graphs and dividing this value by the total flow connections in G 1 and G 2 .

Calculating the geographical consistency gc 1 and the structural consistency gc 2 versus the varying number of clusters, the "saddle point" on the curves can be adopted as the optimal number of clusters for traffic flow representation. Consequently, the NFIs can be determined by hierarchical clustering based on the determined number of clusters.

B. Flow Representation and Prediction

The above section illustrates the methodology for identifying the NFIs which characterize the nominal air traffic flow intersections across different days of traffic. This section will describe the representation of air traffic flow features and the air traffic flow prediction model. This paper proposes to describe air traffic flow dynamics using the spatial-temporal flights distribution in the airspace. More specifically, the traffic flow feature in the airspace is represented as a text paragraph describing the sequences of flights transiting through the NFIs over time. The proposed method for flow feature representation and prediction consists of three steps: trajectory registration, flow dynamics description using spatial-temporal (S-T) flow distribution, and flow prediction at NFIs.

1) Trajectory Registration: With the identified NFIs, a flight trajectory can be represented as a sequence of the NFIs. The objective of trajectory registration is to search for a sequence of NFIs that can optimally approximate the original trajectory. This objective is formulated as finding the minimum dissimilarity between the original trajectory and a representative trajectory constituted by a subset of the NFIs. More details regarding the trajectory registration can be found in [START_REF] Ma | Air traffic flow representation and prediction using transformer in flow-centric airspace[END_REF].

2) Flow Dynamics Description Using S-T Flight Distribution: After trajectory registration, we can understand at what time a flight has transited through which NFI. The traffic flow dynamics in the airspace is represented as the spatial-temporal flight distribution, described by "paragraph" whose "sentences" are the sequences of flights transiting through the NFIs. In this paper, a flight is referred to by its callsign. Even though the flight callsign may be changed when airlines introduce updated schedules, in most cases, a flight callsign remains fixed for a particular flight operating on the same route regularly during a relatively long period (such as a season or a few months).

Use t 0 and t 1 to denote the start and end times of a period. Let F P k denote the traffic flow at NFI v k during t 0 and t 1 . F P k can be described as a sequence of flights (callsigns) according the time flights passing p k , i.e., F P k :

f 1 k , f 2 k , ..., f i k , ..., f m k k with f i
k denoting the i th flight in the sequence and m k denoting the total number of flights passing v k . If there is no flights transited through a NFI, the callsign sequence for this NFI will be replaced by the phrase: "No flights". Analogizing F P k as a "sentence" depicting the flow context at v k , the combinations of "sentences" during t 0 and t 1 constitute a "paragraph" description of the flow context in the entire airspace. The paragraphs for various periods will be used as inputs to the transformer-encoder-based model, introduced in the next section, to learn the contextual relations between air traffic flow and make predictions about the future air traffic flow.

3) Flow Prediction at NFIs: The prediction model is based on the neural network structure developed in by the authors in [START_REF] Ma | Air traffic flow representation and prediction using transformer in flow-centric airspace[END_REF], whose components are tokenization, embedding, transformer encoder blocks, and a fully-connected layer to produce the prediction outcome.

Given an input sequence, i.e., including the time and the flight callsign sequences at the NFIs, word-tokenization [START_REF] Sun | Spectral-spatial feature tokenization transformer for hyperspectral image classification[END_REF], which splits the data based on natural breaks and meaning, such as time (number), callsigns (words), and sequence separations (delimiters), is applied to convert the input into a list of integers that can be embedded into a vector space. After tokenization, the token embedding layer converts the list of integers into a list of vectors. For the model to use the order of the elements in the sequence, positional embeddings, which contain information on the relative or absolute position of the elements, are added to the token embeddings as the input to the transformer encoders. The embeddings are trained jointly with the rest of neural network. Back-propagation is carried through all the network layers up to the embeddings that are updated as other parameters. After embedding the elements in the input sequence, each of them flows through the transformer encoder blocks to encode the feature into a meaningful context tensor representation. The transformer encoder is composed of a stack of N e encoder blocks. The output from the stack of N e encoders is then forwarded to a fully connected layer to obtain the flow prediction results for different NFIs.

C. NFI Flow Acceptance Limit Identification

To identify the flow acceptance limit at an NFI, this paper proposes to use the flight transition duration from neighboring NFIs as the overload indicator. When the demand at an NFI is above the acceptable limit during a period, air traffic congestion can happen, and it will take a significantly larger time cost for air traffic transiting to the overloaded NFI on their flight paths due to regulatory measures such as vectoring and speed control.

Assume there is a number of N v NFIs identified, denoted as: v 1 , v 2 , ..., v i , ..., v Nv . With the time information of flights reaching the NFIs during the s th period, the number of flights transiting to NFI v i from NFI v j can be determined, denoted as n s ij . The set of transition duration for the n s ij flights can be denoted as H s ij = {h s k ij |k = 1, 2, ..., n s ij }, where h s k ij represents the transition duration of the k th flight among the n s ij flights. Considering air traffic flows can evolve daily in terms of locations of the NFIs, the traffic flow structure, or unexpected events such as mechanical issues and air traffic control disruptions, the duration H s ij is normalized by the daily minimum duration to reduce the effects of daily fluctuations in air traffic flow:

N H s ij = H s ij H Ds ij ( 2
)
where H Ds ij represents the minimum flight transition duration from v j to v i on the day of the period s. By doing so for all NFIs that are connected to v i , the normalized transition duration to v i for the period s can be obtained:

N H s i = {N H s ij |j ∈ N i },
where N i is the neighbor set of v i . The total number of flight transiting to v i can consequently be determined as n s i = n s ij , j = 1, 2, ..., N i . For all periods s ∈ {1, 2, ..., N T } in the traffic data, the set of flight transition duration N H i = {N H s i |s ∈ {1, 2, ..., N T }} as well as the set of corresponding traffic demand n i = {n s i |s ∈ {1, 2, ..., N T }} can be determined through the above procedures. By fitting the demand values to the transition duration values, a set of demand values X = {1, 2, ..., N } and the corresponding transition duration value Y = {y 1 , y 2 , ..., y N } can be obtained. The flow acceptance limit is determined by identifying the abrupt changes on the curve of demand versus transition duration, which is formulated as [START_REF] Killick | Optimal detection of changepoints with a linear computational cost[END_REF]:

obj = arg min k (k)var([y 1 , ..., y k ]) +(N -k)var([y k+1 , ..., y N ]) (3)
By minimizing Eq. 4, the flow acceptance limit l i at v i can be determined as k, flow demand above which will lead to abrupt increase in the transition duration of flights to the NFI v i .

D. NFI Flow Coordination

To re-route the flows at an NFI, a pre-requisite is the knowledge of the main flow paths crossing the NFI, i.e., the main flows transiting through the NFI and their corresponding paths. This paper identifies the main flow paths from the flight data based on the typical sequences of NFIs flown by the flights. Upon identifying the main flows transiting through an NFI, the next step is to predict the number of flights in each main flow during a future period. The prediction model in this step utilizes the same input features and neural network architecture as described in Section III-B, except that the model output is the main flow values at the target NFI instead of the flow demand values of different NFIs.

With the knowledge of the path and predicted demand of each main flow transiting through the NFI, the third step of flow coordination is to re-configure the flow excess by re-routing the main flows so that overload at the NFIs can be avoided. Let f i represent the predicted flow demand at the NFI v i . Assume that the predicted flow demand f * at the NFI v * exceeds the acceptance limit l * . Let A be the adjacency matrix of G, with its entry a ij = 1 if NFIs v i and v j are connected by traffic flows, otherwise, 0. Given that v * is predicted to be overload and excessive flow need to be re-routed to the other NFIs, a i * and a * j are set as 0. Assume there is a number of T main flows identified for v * . Let D = {d 1 , ..., d t , ..., d T } be the vector for the number of re-routed flights in the T flows, i.e., d t is the number of re-routed flights in the t th flow. Assume that for the t th flow, there are maximum R t accessible paths respectively denoted by P 1 t , P 2 t , ..., P r t , ..., P Rt t . The flow on path P r t is represented by f tr . Define a function δ (tr,ij) , with δ (tr,ij) = 1 if the link between v i and v j is on the path P r t , otherwise, δ (tr,ij) = 0. Let d ij and x ij represent the great circle distance and the amount of flow between v i and v j respectively. With all the above notations and definitions, the flow coordination model is formulated as:

obj = min X={xij } i,j=1,...,Nv Nv i=1 Nv j=1 a ij x ij d ij s.t. f tr ≥ 0 Rt r=1 f tr = d t T t=1 d t = l * -f * Nv j=1 a ij x ij ≤ l i -f i T t=1 Rt r=1 δ(tr, ij)f tr = x ij (4)
Optimizing the above function allows the traffic flow excess at v * to be re-routed so that minimum flight distance is occurred without exceeding the acceptance limit of all the NFIs.

IV. EXPERIMENTAL STUDY

To verify the efficacy of the proposed framework, an experimental study has been carried out on the French airspace using one-month ADS-B data from December 1 to December 31, 2019, including a number of 158856 flights. This study focuses on the en-route air traffic above 10,000 ft. The target in this experimental study is set as follows: identify the flow acceptance limits on the NFIs during a 30-minute interval, predict the number of flights that will transit through the NFIs in the future 30 minutes, and re-route the flow excess to avoid overloading the NFIs.

A. NFI Identification

By finding intersections of flight trajectories in the French airspace and clustering the intersection on a daily basis, a graph representation of the daily air traffic flow pattern can be obtained.

By calculating the geographical consistency gc 1 and the structural consistency gc 2 against a set of different cluster numbers ranging from 100 to 1500. A "saddle point" is observed for gc 1 and gc 2 around cluster number 605. Therefore, this paper takes the value 605 as the number of clusters to be identified by the hierarchical clustering algorithm. The centers of identified clusters are determined as the NFIs. Fig. 3 shows the graph representation for the one-month air traffic structure with 605 NFIs. We can see that this graph is able to depict the nodal hierarchy of air traffic flows ranging from regional feeders to international hubs, such as Paris and Geneva. The en-route air traffic flows are organized as a series of "spokes" connecting the traffic hubs or connecting outlying areas to a hub area. Fig. 6 presents the traffic flow prediction result on eight example NFIs in the airspace from 00:00 to 23:59 on Dec 24, 2019. The solid blue lines show the true number of flights passing the NFIs, while the red lines show the predicted value using the proposed prediction method. We can observe from Fig. 6 that the proposed method sustainably gives forecasts in close proximity to the actual flow value for different NFIs in the future 30-minute horizon. Furthermore, the prediction method can capture sharp changes in the air traffic demand, as can be observed from the figure that when the actual demand increases or decreases abruptly, the prediction model is still able to give a close prediction of the true value. Overall, the proposed prediction model approximates the actual flow values on the identified NFIs with the coefficient of determination (R 2 ) above 0.9 and mean absolute percentage error (MAPE) below 20%.

D. NFI Flow Excess Re-routing

Based on the identified flow acceptance limit and the predicted flow during the 30-minute interval, we can observe when an NFI is overloaded and what the flow excess is. It can be observed from Fig. 6 that, for the NFI whose prediction and acceptable flow limit are shown on the bottom-left, the flow is predicted to exceed the acceptable limit from 12:10 to 12:40 on the day. The value of the flow excess is five, meaning five flights should be rerouted to avoid overloading the NFI. This NFI is shown in Fig. 7 and marked by the red node. Three main flows are transiting through the NFI. The green, black, and blue lines visualize flight trajectories in the three main flows. The predicted flights in the three main flows t = 1, t = 2, and t = 3, are 5, 4, and 8, respectively.

With the knowledge of the number of predicted demand, f * = 23, the acceptable flow limit, l * = 18, and the number of flights in each main flow, this step will re-route the traffic in the main flows onto alternative routes consisting of other NFIs with spare capacity. Although exceptions exist due to traffic congestion or weather conditions, flights in the en-route phase usually prefer the routes with shorter flight distance due to fuel consumption and en-route charges [START_REF] Soler | En-route optimal flight planning constrained to pass through waypoints using minlp[END_REF]. Therefore, the alternative paths of each main flow, i.e., P r t , f or t = 1, 2, 3 in this experiment are set as the two shortest paths, i.e., R t = 2, f or t = 1, 2, 3 between its origin NFI and destination NFI without transiting through the overloaded NFI. The paths are identified based on the connectivity among the NFIs of the main flow. Therefore, in this case the Eq. 4 can be specified as follows:

a ij x ij d ij s.t. f tr ≥ 0 2 r=1 f tr = d t 3 t=1 d t = 5 605 j=1 a ij x ij ≤ l i -f i 3 t=1 2 r=1 δ(tr, ij)f tr = x ij (5)
Fig. 8 and Table I show the flow excess re-routing result. Five flights in the third main flow marked by solid blue lines are re-routed to its shortest alternative path, i.e., f 1 3 = 5, shown by the blue dashes. Note that the predicted demand is 23, and the flow limit is 18. The normalized flight transition duration, scaled between 0 and 1, to the NFI under a demand of 18 flights is 0.27, while under a demand of 23 is 1. Thus, the five flight excess in the traffic flow can potentially lead to a 270% increase in the flight transition duration. Through re-configuring the flow demand at the NFIs, the excess is reduced to 0 from 5, and the anticipated flow overload at the NFI is avoided in advance without causing abrupt increases in the flight durations with the utilization of the spare capacity of other underloaded NFIs. the optimal number of clusters. A graph analysis was proposed to identify the clustering outcome, which optimally sustained the consistency in the flow patterns. Secondly, based on the identified NFIs, the air traffic flow dynamics features were represented by the spatial-temporal flight distribution, characterized by a textual paragraph recording the time and the sequences of flights passing each NFI. Then, a transformer-based model was adopted to learn the text-enriched flow features and predict future air traffic at the NFIs. Thirdly, for each NFI, the acceptable flow limit was determined by identifying the transition point of the normalized flight transition duration from its neighboring NFIs versus the traffic demand. Finally, based on the identified acceptable flow limit and the predicted demand, the flow excess at an NFI was alternated onto other NFIs to optimize and reconfigure the air traffic demand to avoid sharp increases in flight durations. An experimental study was carried out in French airspace using the proposed framework base on one-month ADS-B data in December 2019. Results showed that the proposed prediction model approximated the actual flow values with the coefficient of determination (R 2 ) above 0.9 and mean absolute percentage error (MAPE) below 20%. Moreover, acceptable flow limit determination showed that for above 68% NFIs, the flight transition duration increased sharply when the demand exceeded a certain level. Flow overload was avoided by rerouting the detected flow excess to other NFIs with spare capacity.

Research in this paper provides a basic framework for flowcentric air traffic flow coordination. In the future, constraints from both perspectives of air traffic control services and airspace users, such as Special Use Airspace (SUA), weather conditions, airline schedules, and users' preferred routes, can be incorporated to refine this framework for more effective and efficient coordination of air traffic flow. For instance, SUA and weather conditions can be adapted to produce a more accurate and dynamic estimation of the flow acceptance limit. At the same time, airline schedules and users' preferred routes can assist the proposed flow coordination algorithm to re-configure the flow demand to the preferences of airspace users and airlines.

Fig. 1 :

 1 Fig. 1: Conceptual diagram of the proposed dynamic air traffic flow coordination framework, including: 1) NFI identification through graph analysis of air traffic flow patterns; 2) flow representation and prediction based on the spatial-temporal flow distribution; 3) NFI flow acceptance limit identification through flow transition efficiency analysis; 4) Flow coordination, i.e., flow-excess re-routing, based on the capacity availability and predicted flow demand at the NFIs.

  where |S 1 ∪ S 2 | represents the number of node pairs in the union of S 1 and S 2 . It will be denoted as l in the rest of the paper. b) Structural Consistency in Flow Connectivity: Upon determining the geographical consistency in NFI locations, the next step is to quantify the structural consistency in the daily air traffic flow connectivity between the NFIs. Let S = S 1 ∪ S 2 represent the set of paired nodes from graph G k and graph G k+1 based on NFI location consistency. While C k = {c 1 k , c 2 k , ..., c l k } denotes the nodes in S from graph G k and C k+1 = {c 1 k+1 , c 2 k+1 , ..., c l k+1 } denotes the corresponding paired nodes from G k+1 . Let e i j c k represent the edge connecting nodes c i k and c j k and w i j c k represent the weight on it. Let W c k represent the entire set of edge weights {w 1 2 c k , w 1 3 c k , ..., w 1 j c k , ..., w 1 l c k , ..., w 2 3 c k , w 2 4 c k , ..., w 2 j c k , ..., w 2 l c k , ..., w i (i+1) c k , w i (i+2) c k , ..., w i j c k , ..., w i l c k , ..., w (l-1) l c k }. The structural consistency in the air traffic flow patterns is measured by the mutual flow connectivity in the two graphs. More specifically, with the set of paired nodes C k and C k+1 from G k and G k+1 respectively, the flow structure consistency is evaluated by the ratio of the mutual flow connections between the graphs characterized by nodes C k and C k+1 compared to the union of flow connections in G k and G k+1 . It is formulated as:

Fig. 2 :

 2 Fig. 2: The diagram for graph analysis of node (NFI) location consistency and flow structure consistency. G 1 and G 2 are the graphs constructed for the air traffic flow on the day 1 and day 2, respectively. The nodes are marked by solid gray circles. The red solids, i.e., C 1 and C 2 , in the graphs denote the nodes showing location consistency across the two days. W c 1 and W c 2 are the edge weights in the two sub-graphs of G 1 and G 2 formed by C 1 and C 2 . The structural consistency is characterized by the mutual connections in the two sub-graphs, i.e., min{W c 1 , W c 2 }, compared to the total flow connections in G 1 and G 2 .

Fig. 3 :

 3 Fig. 3: Graph representation for the one-month air traffic flow using 605 NFIs.

Fig. 4 :

 4 Fig. 4: Comparison of the identified NFIs for two different days.

Fig. 3

 3 Fig. 3 presents the NFIs from a spatial structural perspective, and Fig. 4 depicts the temporal dynamics in the NFIs in describing air traffic flows during different days. Fig. 4 shows the identified NFIs for two days, yellow dots for day one and red dots for day two. The dots' sizes are proportional to the traffic volume transiting through the corresponding NFIs. It can be observed that the NFIs show consistent patterns in geographical distribution and traffic volume across different days, although

Fig. 5 :

 5 Fig. 5: Flight transition duration versus the flow demand on eight example NFIs. The blue circles show the observations from the traffic data, the solid red lines show the third-degree polynomial fitting of the observations, and the pink dashes bound the 95% confidence intervals of the fitting. The solid black line indicates the identified acceptance limit of the NFIs.

Fig. 6 :

 6 Fig. 6: Flow prediction result on Dec 24, 2019 for eight example NFIs. The solid blue lines show the true number of flights passing the NFIs, while the red lines show the predicted value.

Fig. 7 :

 7 Fig. 7: Three main flows, depicted by the green, blue, and black lines, transiting through the overloaded NFI (red node).

Fig. 8 :

 8 Fig. 8: Flow excess re-routing results. Five flights in the original flow marked by solid blue lines are re-routed to the path shown by the blue dashes.

  2 k+1 , ..., v j k+1 , ..., v n k+1 } represent the set of nodes in the graph G k+1 constructed for day D k+1 . For each node v i k in V k , the proposed algorithm searches for its nearest node v ai k+1 in V k+1 . a i is the index of the identified nearest node of v i k in V k+1 . By representing node v i

	k 's latitude and k+1 's as (φ i k ) and node v i k , λ i longitude as (φ i k+1 , λ i k+1 ), the nearest neighbour of node v i k in V k+1 is identified according
	the great circle distance.

TABLE I .

 I Flow excess re-routing result.

	flow	f t0	before excess	f t1	f t2	dt	f t0	after excess
	t = 1	5			0	0	0	5
	t = 2	4	5		0	0	0	4	0
	t = 3	8			5	0	5	3
				605	605
	obj =	min
		X={xij } i,j=1,...,605	i=1	j=1

ACKNOWLEDGMENT

This research is supported by the National Research Foundation, Singapore, and the Civil Aviation Authority of Singapore, under the Aviation Transformation Programme. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not reflect the views of National Research Foundation, Singapore and the Civil Aviation Authority of Singapore. The authors would like to thank professor Patric Senac, ENAC, for hosting the first author at ENAC, Toulouse, where part of this research was conducted.