V Meyer

L Maxit

Y Renou

C Audoly

Experimental investigation of the influence of internal frames on the vibroacoustic behavior of a stiffened cylindrical shell using wavenumber analysis

Keywords: Stiffened shell, Vibroacoustics, Non-axisymmetric, Laser vibrometer, wavenumber analysis

The understanding of the influence of non-axisymmetric internal frames on the vibroacoustic behavior of a stiffened cylindrical shell is of high interest for the naval or aeronautic industries. Several numerical studies have shown that the non-axisymmetric internal frame can increase the radiation efficiency significantly in the case of a mechanical point force. However, less attention has been paid to the experimental verification of this statement. That is why this paper proposes to compare the radiation efficiency estimated experimentally for a stiffened cylindrical shell with and without internal frames. The experimental process is based on scanning laser vibrometer measurements of the vibrations on the surface of the shell. A transform of the vibratory field in the wavenumber domain is then performed. It allows estimating the far-field radiated pressure with the stationary phase theorem. An increase of the radiation efficiency is observed in the low frequencies. Analysis of the velocity field in the physical and wavenumber spaces allows highlighting the coupling of the circumferential orders at the origin of the increase in the radiation efficiency.

Introduction

Many works can be found in the literature dealing with the vibroacoustic behavior of cylindrical shells reinforced by circumferential stiffeners [START_REF] Arnold | Flexural vibrations of the walls of thin cylindrical shells having freely supported ends[END_REF][START_REF] Galletly | On the in-vacuo vibrations of simply supported, ring-stiffened cylindrical shells[END_REF][START_REF] Hoppmann | Some characteristics of the flexural vibrations of orthogonally stiffened cylindrical shells[END_REF]. As a matter of fact, many structures in the aeronautical or naval industry can be modeled by stiffened cylindrical shells. Williams et al. [START_REF] Williams | Broadband nearfield acoustical holography for vibrating cylinder[END_REF] measured the nearfield pressure on a submerged cylinder excited by a point force. They deduced the vibration and radiation from the cylinder thanks to a technique called the generalized nearfield acoustical holography. In a subsequent study, the authors analyzed the results in the wavenumber domain [START_REF] Williams | Experimental investigation of the wave propagation on a point-driven, submerged capped cylinder using k-space analysis[END_REF]. This analysis shows a great interest to give a physical insight of the vibroacoustic behavior of the cylindrical shell. More particularly, wavenumber diagrams can be plotted and the waves contributing to the far-field radiation can be identified. The technique is later used by Photiadis et al. to investigate the response of a stiffened cylinder in the wavenumber domain [START_REF] Photiadis | Wave-number space response of a near periodically ribbed shell[END_REF]. The wavenumber diagrams show more complex patterns, due to the stiffeners periodicity and the propagation of Bloch-Floquet waves. Slowness surfaces, which are given by the normal to the curves of large amplitude and large radius of curvature on the wavenumber diagrams, are used to predict the dominant directions of wave propagation.

The studies mentioned previously are limited to axisymmetric systems. However, in industrial applications, internal structures such as floor partitions or engine foundations can have an influence on the vibroacoustic behavior of the stiffened cylindrical shell. Experimental work has been carried out for a stiffened cylindrical shell coupled to numerous oscillators, and the influence on the scattered or radiated pressure has been highlighted [START_REF] Photiadis | The effect of internal oscillators on the acoustic response of a submerged shell[END_REF][START_REF] Photiadis | Resonant response of complex shell structures[END_REF]. The results show that adding internal degrees of freedom reduces scattering from flexural Bloch-Floquet waves and gives rise to a noisy speckle pattern. Nevertheless, in the case of a point mechanical excitation, only a few studies based on numerical simulation have investigated the effect of non-axisymmetric internal frames. Some of them use discretization methods such as Finite Element Method and Boundary Element Method to predict the radiated pressure by an immersed shell with and without internal masses or structures and show the importance of flexible mounts [START_REF] Marcus | The effect of internal point masses on the radiation of a ribbed cylindrical shell[END_REF][START_REF] Peters | Effects of internal mass distribution and its isolation on the acoustic characteristics of a submerged hull[END_REF]. Other studies use substructuring approaches based on the admittance method to include internal structures in an axisymmetric cylindrical shell [START_REF] Ettouney | The interaction of a submerged axisymmetric shell and three-dimensional internal systems[END_REF][START_REF] Meyer | Prediction of the vibroacoustic behavior of a submerged shell with non-axisymmetric internal substructures by a condensed transfer function method[END_REF]. The main results are that more circumferential orders tend to play a role in the response of the cylindrical shell when it is coupled to a nonaxisymmetric internal frame. If the problem is decomposed on a modal basis, it is equivalent to say that the non-axisymmetry tends to couple the circumferential orders [START_REF] Rebillard | Influence of an embarked springmass system and defects on the acoustical radiation of a cylindrical shell[END_REF]. As the structural waves having low wavenumbers are the main contributors to the far-field radiated pressure, the coupling with low circumferential numbers tends to increase the radiation efficiency. Another effect of the internal frame is given by the fuzzy internal theory [START_REF] Soize | Probabilistic structural modeling in linear dynamic analysis of complex mechanical systems, I-theoretical elements[END_REF][START_REF] Pierce | Resonant-frequency-distribution of internal mass inferred from mechanical impedance matrices, with application to fuzzy structure theory[END_REF]. This theory states that if the internal structure has a high modal density and if the impedance at the junction with the cylindrical shell is adapted, the addition of internal degrees of freedom can give rise to an apparent damping effect.

The aim of the present work is to discuss experimental results of a pointdriven stiffened cylindrical shell with and without a non-axisymmetric internal frame. The effect of the internal structure on the vibroacoustic behavior of the cylindrical shell, and more especially on the radiation efficiency is investigated. The main obstacle for estimating accurately the radiation efficiency lies in the estimation of the radiated power. Precision methods using sound pressure require a large number of microphones and are difficult to set up [START_REF] Iso | Acoustics -Determination of sound power levels of noise sources using sound pressure -Precision methods for anechoic and semi-anechoic rooms[END_REF]. Measurement by scanning using sound intensity [START_REF] Iso | Acoustics -Determination of sound power levels of noise sources using sound intensity -Part II: Measurement by scanning[END_REF] is widely used to measure the radiated power from machines, but it can be difficult to define an enclosing surface around the system. To tackle these issues, the methodology of the present work is based on the use of the stationary phase theorem. It relates the far-field radiated pressure to the spectral velocity field of the shell (i.e. the velocity field in the wavenumber domain). An advantage of this approach over the measurement with microphones or intensity probes is that the directivity can also be estimated. The spectral velocity field can be obtained with a scanning laser vibrometer and a 2D Fourier transform. Attention is paid to the definition of the scanning grid in order to avoid aliasing effect. The paper is organized as follows:

• Background on the estimation of the far-field radiated pressure from the spectral displacements of the shell is given in section 2.

• The system and the experimental setup are presented in section 3.

• Operational dispersion curves plotted from the spectral velocity field are discussed in section 4.

• In section 5, the stationary phase theorem is used to determine the far-field radiated power and to investigate the influence of the non-axisymmetric frame on the radiation efficiency.

• Velocity fields are compared for the axisymmetric and non-axisymmetric cases in the real space and in the wavenumber domain in section 6.

• Finally, conclusions are drawn in section 7.

Spectral displacements and far-field pressure

Let us consider a cylindrical shell of radius R 0 immersed in an acoustic domain. The shell can either be of infinite length or of finite length extended by a cylindrical baffle. Considering the cylindrical coordinates (x, r, θ), x being the axis of the shell, W (x, θ) is the radial shell displacements at a given angular frequency ω.

A two dimensional Fourier transform can be used to transform the spatial variations of a vibratory field into the wavenumber domain. The 2D Fourier transform of W (x, θ) is given by

W (x, θ) → W (k x , n) = R 2π θ=0 W (x, θ)e -j(kxx+nθ) dxdθ (1)
where k x is the axial wavenumber and n the circumferential order. As the cylindrical shell is 2π-periodic along the circumferential coordinate θ, then n is a relative integer.

The wavenumber analysis gives an opportunity to analyze the wave propagation on the structure and the role played by the contribution of each circumferential order. The wavenumber analysis also gives an insight on the spectral components contributing to the radiated pressure, through the radiation circle defined by the points of coordinates ( kx , n) in the wavenumber domain

n R 0 2 + k2 x = k 2 0 ( 2
)
where k 0 = ω c0 is the acoustic wavenumber and c 0 = 340 m.s -1 is the sound velocity in air. For a plate, the radiation circle represents the supersonic acoustic domain, in which the components radiate to the far-field. It remains a good indicator for localizing the radiating waves on a cylindrical shell when the effect of the curvature is low [START_REF] Williams | Experimental investigation of the wave propagation on a point-driven, submerged capped cylinder using k-space analysis[END_REF].

The vibratory field in the wavenumber domain can be used to estimate the far-field radiated pressure. For R in the acoustic far-field, i.e. R 2π k0 , the asymptotic expression for the Hankel function and the stationary phase theorem can be used [START_REF] Junger | Sound, structures, and their interaction[END_REF]. If ξ is the angle from the point to the normal axis and θ is the circumferential angle, the far-field pressure p can be written as:

p(R, ξ, θ) = +∞ n=-∞ jρ 0 ω 2 2π 2 Rk 0 cos ξ W (-k 0 sin ξ, n) H (2) n (R 0 k 0 cos ξ) e -jRk0+jn(θ+ π 2 ) (3)
where ρ 0 = 1.2 kg.m -3 is the density of air and H

(2) n is the derivative of the Hankel function of second kind and order n. For a cylinder of finite length, the stationary phase theorem is based on the assumption of a baffled cylinder. However, the influence of the baffle on the radiation due to radial vibrations is weak when the shell is long compared to the acoustic wavelength [START_REF] Sandman | Fluid-loading influence coefficients for a finite cylindrical shell[END_REF]. The stationary phase theorem may then be used to predict the radiated pressure from the cylindrical shell. One of the advantages of the wavenumber analysis is that it allows the calculation of the radiated pressure from vibrations measurements only.

The far-field radiated power P a is obtained by integrating the radiated pressure over the surface of the sphere of radius R:

P a = 1 ρ 0 c 0 π ξ=0 2π θ=0 p 2 (R, ξ, θ)R 2 sin θdξdθ (4)
As the acoustic domain is conservative, any value of R in the far-field can be chosen to yield P a . The radiation efficiency σ of a cylindrical shell is defined by the ratio between the actual far-field radiated power P a and the power that would be radiated by a cylindrical shell of the same radiating surface S = 2πR 0 L, having an uniform radial velocity equal to the mean quadratic velocity

< Ẇ >: σ = P a ρ 0 c 0 S < Ẇ > 2 (5)
The mean quadratic velocity is defined as follows:

< Ẇ > 2 = 1 S L x=0 2π θ=0 | Ẇ (x, θ)| 2 dxR 0 dθ (6)
with Ẇ (x, θ) the point radial velocity at the point of coordinates (x, R 0 , θ).

Presentation of the experimental setup

Axisymmetric stiffened cylindrical shell

A stiffened cylindrical shell, which geometrical characteristics are given in table 1, is shown in Fig. 1. The cylindrical shell is made of standard nonalloy steel, for which common values may be considered: Young's modulus E = 2.1.10 11 Pa, density ρ s = 7800 kg.m -3 and Poisson's coefficient ν = 0.3. It has 72 identical stiffeners of rectangular cross-section, which have been machined from a thicker cylinder. The cylindrical shell is divided in five sections with three different stiffeners spacings, as shown in the sketch in Fig. 1c. 15 mm thick caps are screwed at the ends of the cylinder. This stiffened cylindrical shell is referred to as the axisymmetric case. The ring frequency f r , for which the longitudinal wavelength is equal to the circumference, is defined by the following equation [START_REF] Junger | Sound, structures, and their interaction[END_REF]:

f r = 1 2πR 0 E ρ s (1 -ν 2 ) = 8657 Hz. ( 7
)
The coincidence frequency, for which the acoustic wavenumber k 0 is equal to the flexural wavenumber in an equivalent plate of thickness h, can be defined as [START_REF] Junger | Sound, structures, and their interaction[END_REF]:

f c = c 2 0 2πh 12ρ s (1 -ν 2 ) E = 7812 Hz. ( 8
)
In general, the radiation efficiency defined in Eq. ( 5) presents peaks at the ring frequency and at the coincidence frequency [START_REF] Manning | Radiation properties of cylindrical shells[END_REF]. It can be noticed that for this model the ring frequency and the coincidence frequency have close values. The radiation efficiency is expected to be low below these frequencies and to tend to 1 in the supersonic domain above the coincidence frequency.

Non-axisymmetric internal frame

In the second configuration, a non-axisymmetric internal frame is added inside the cylinder described in section 3.1. As shown in Fig. 2, it is a floor which is 869 mm long, 64.3 mm high and 170.6 mm wide. It consists in a 1.5 mm thick plate, for which the ends have been cut along a circular arc of radius 98.5 mm, before being bended to form right angles. The contact between the cylindrical shell and the floor is only on these two circular arcs. A square tube of 20 mm wide and 2mm thick in steel is glued under the horizontal part and reinforces the structure lengthwise. The floor is placed inside the cylindrical shell, with line contacts along the circumference at x = 0.315 m and x = 1.185 m. Several constraints need to be taken into account to assemble the two subsystems:

• the displacements continuity should be verified (i.e. rigid connection) and no damping (i.e. energy dissipation) should be introduced;

• the junctions should be strong enough to support the floor when the cylindrical shell is hanged;

• the access inside the cylindrical shell is difficult because of the small radius;

• the process should not damage the cylinder;

• the internal structure could be dismounted once the measurements done.

Tests are done on samples made of two steel plates perpendicularly assembled (not presented here for sake of brevity). Several types of glue, as well as brazing are compared. An acrylic glue with activator (Loctite AA 330) is chosen for this application. A picture of the assembled system is shown in Fig. 3.

Experimental setup

As shown in Fig. 4, the system is hanged vertically on an arm. A strap is used in order to limit the transfer of vibrations between the system and the arm. A rotating platform allows rotating the cylindrical shell around its axis and identifying precisely the angles. A 12.5 mm diameter patch is glued on the outer surface of the cylindrical shell, and screwed to a shaker. The shaker is hanged on the rotating platform in order to excite always the same point of coordinate (x, r, θ) = (1.195, 0.1, 0) in the cylindrical system (with x = 0 at the bottom of the cylindrical shell). This excitation point is located at the bottom of a stiffener, and lies at 10 mm from one of the floor junctions in the non-axisymmetric configuration. Two elements come in between the patch and the shaker: an impedance head that measures acceleration and force at the excitation point, and a threaded rod that limits the moments and allows assuming a radial excitation in the cylindrical system. The measurements are conducted in a semi-anechoic room (8.6 m long, 6.4 m wide, 4 m high and 80 Hz cut-off frequency). The floor is made of concrete and a 175 mm space separates it from the bottom of the cylinder. These acoustic conditions are supposed to weakly influence the shell vibrations as the bulk modulus of air is low (i.e. weak coupling between the fluid and the structure). The laser vibrometer is placed 2.68 m far from the cylindrical shell's surface as seen in Fig. 5. The vibrometer scans the cylinder from down to up along a single generating line. The experiment consists in keeping the vibrometer at a fixed position while rotating the set cylindrical shell-shaker.

The aim is to measure the radial vibrations on the outer surface of the cylindrical shell for the two configurations, that means with and without the non-axisymmetric internal frame. The measurements are made up to 16 kHz. The excitation is a chirp (swept-frequency sine), and the output is the transfer function with regard to the excitation force, calculated on 20 averages.

Scanning grid definition

The grid definition in the physical space determines the definition of the wavenumber domain (see the properties of discrete Fourier transform [START_REF] Maxit | Wavenumber space and physical space response of a periodically ribbed plate to a point drive: A discrete approach[END_REF]). In order to examine the vibrations results in the wavenumber domain, the scanning grid is defined as a function of the flexural wavenumber k f of an equivalent plate (in steel, of thickness h = 1.5 mm) at the highest frequency [START_REF] Junger | Sound, structures, and their interaction[END_REF]:

k f = ω h 12 (1 -ν 2 ) ρ s E 1/2 (9)
To avoid aliasing effect, the maximum axial wavenumber is k max x = k f 207 m -1 and the maximum circumferential number is n max = k f R 0 20. For a cylindrical shell, as the system is stiffer than a plate because of the curvature, the natural wavenumbers of the propagative waves are lower than k f . The criteria are thus stricter than they could be. The maximum distance d max between two consecutive measurement points is then:

d max = 2π 2k max x (10)
The numerical application yields d max = 15 mm. A scanning line with N x = 101 points along the cylinder is thus defined. The maximum angle θ max between two consecutive measurements is given similarly by:

θ max = 2π 2n max (11)
The present case gives θ max = 9 • and results in N θ = 41 lines to scan. The system being symmetric, it is sufficient to measure between 0 and 180 • . The experiment consists then in 21 linear scans, where the system is successively rotated by 9 • .

Similarly, the resolution in the wavenumber domain δk x is linked to the length of the cylindrical shell:

δk x = 2π L 4.2 (12)
The axial wavenumber resolution δk x can be enhanced using zero padding, but results show that the resolution is fine enough to observe the variations with the axial wavenumber.

Operational dispersion curves for the axisymmetric system

The 2D Fourier transform defined in Eq. ( 1) is applied to evaluate experimentally the velocity of the shell in the wavenumber domain. In practice, as the velocities have been measured on a regularly spaced grid, the integrals of Eq. ( 1) are approximated by the sum of all the points on the grid, weighted by the exponential functions:

W (k x , n) Nx i=1 N θ j=1 Ẇ (x i , θ j )e -j(kxxi+nθj ) d max θ max (13)
where (x i , θ j ) is the coordinate of a point on the scanning grid.

Operational dispersion curves are constructed by plotting the amplitude of the spectral radial velocity on the surface of the axisymmetric cylindrical shell as a function of the frequency and the axial wavenumber, for a given circumferential order. The plots in Fig. 6 show the dispersion curves for the circumferential orders n = 0, n = 1, n = 4 and n = 6 for the axisymmetric case. Theoretical models are presented in the left column and compared to the experimental results in the right column. The theoretical model uses the Fluegge equations for an infinite cylindrical shell with stiffeners included in the dynamic model using smeared theory [START_REF] Caresta | Acoustic signature of a submarine hull under harmonic excitation[END_REF][START_REF] Rosen | Vibrations of axially loaded stiffened cylindrical shells[END_REF]. An average stiffeners spacing of 20 mm is considered.

On the theoretical curves, the solid white line represents the flexural wavenumber for an equivalent plate projected on the circumferential order n:

k xf (n) = k 2 f - n R 0 2 (14)
where k f is the flexural wavenumber for a plate defined in Eq. ( 9). The dashed white line represents the acoustical wavenumber projected on the circumferential order n:

k x0 (n) = k 2 0 - n R 0 2 (15)
The area under this dashed curve represents the supersonic domain, where the structural waves are strongly coupled to the acoustic domain and radiate efficiently. For the circumferential orders n = 0 and n = 1, the dashed-dotted line respresents the longitudinal wavenumbers projected on the circumferential order n given by: where k l = ω ρ s (1 -ν 2 )/E is the longitudinal wavenumber. For the circumferential orders n = 4 and n = 6, k xl is not real and is not represented on the plots.

k xl (n) = k 2 l - n R 0 2 ( 16
) (a) (b) (c) (d) (e) (f) (g) (h)
For the theoretical model in Fig. 6a, the radial breathing mode n = 0 is plotted. As expected, one can see that the model for a plate is not adapted at low frequencies (the solid white line does not fit the curve), meaning that the curvature of the cylindrical shell has an important effect on the vibroacoustic behavior. The trend for the plate model fits better at higher frequencies (above 12 kHz), where the curvature can be neglected [START_REF] Heckl | Vibrations of point-driven cylindrical shells[END_REF]. The breathing mode appears around 8.5 kHz, at the ring frequency [see Eq. ( 7)]. The coincidence frequency for the equivalent plate is given by the intersection of the flexural wavenumber curve with the acoustical wavenumber [see Eq. ( 8)] and reads around 7.9 kHz. This value is however not relevant for the present case because of the curvature. Indeed, as the operational dispersion curve from the cylindrical shell model lies always under the dashed curve, it shows that the breathing mode always radiates efficiently. In Fig. 6a and6c, the curves seen close to the axis k x = 0 are propagating with high velocities. They are well approached by the dasheddotted lines and can be identified as longitudinal waves.

For the circumferential orders n = 1 and n = 4 in Fig. 6c and 6e respectively, the effect of curvature can also be seen. The dispersion curve have a S-shape instead of the square-root shape for the flexural waves in a plate. The coincidence frequency for the cylindrical shell model can then be read around 2.7 kHz, instead of 8 kHz for the plate. Cut-on frequencies can also be seen on the operational dispersion curves. For instance, it can be seen that the waves cannot propagate for the circumferential order n = 4 below 1500 Hz. It can be checked that the cut-on frequency increases with the circumferential order.

For the circumferential order n = 6 in Fig. 6g, the cut-on frequency is around 3700 Hz. For the circumferential order n = 6, the effect of curvature is less important and the dispersion curve is closer to the plate model. However, it can be seen that the flexural waves are faster than the plate model shown by the white solid line, due to added stiffness because of the curvature and the stiffeners. Consequently, the whole dispersion curve is in the supersonic domain and this circumferential order always radiates efficiently.

Comparing the experimental and theoretical operational dispersion curves in Fig. 6 shows a good agreement between the curves. On the experimental plots, spots of high amplitudes can be seen instead of continuous lines. These spots correspond to axial resonances due to the fact that the experimental cylindrical shell is of finite length. However, the comparison also shows that it is difficult to isolate the contribution from only one circumferential order in the experimental case. Additional dispersion curves can indeed be seen, even if their levels are lower than the main contribution. Three explanations are considered:

• As the cylindrical shell has sections with three different stiffeners spacings (see Fig. 1c), it could explain why several dispersion curves can be seen for one circumferential order. However, numerical tests show that the influence of the stiffeners spacing is weak on the dispersion curves and can only be seen at high frequencies.

• Uncertainties in the position of the points on the scanning grid can occur. These uncertainties can break the regularity of the measured velocity field and couple the circumferential orders. However, these uncertainties have been numerically tested on an infinite shell model and a equally distributed law on the points of the scanning grid (not shown here for the sake of conciseness). It can be seen that their influence is weak on the coupling of the circumferential orders.

• Although the particular experimental care, the system and the excitation are not perfectly axisymmetric. Indeed, it is really hard in practice to be sure that the shaker only excites the cylindrical shell in a perfectly radial direction. Besides there may be imperfections on the cylindrical shell geometry because of the machining (0.1 mm tolerance on the shell thickness for instance) and the damages during transportation. It is known that non-axisymmetries tend to couple circumferential orders [START_REF] Rebillard | Influence of an embarked springmass system and defects on the acoustical radiation of a cylindrical shell[END_REF].

This suggests that the additional dispersion curves seen in the right column of Fig. 6 are a coupling of the circumferential orders due to imperfections in the axisymmetry of the system.

Influence of the non-axisymmetric internal frame

Mean quadratic velocity

The mean quadratic velocity defined in Eq. ( 6) is written in a discrete form as

< Ẇ > 2 = 1 S Nx i=1 N θ j=1 | Ẇ (x i , θ j )| 2 d max R 0 θ max (17)
The mean quadratic velocity is plotted as a function of the frequency in Fig. 7a for the axisymmetric case. Only the frequency range between 0 and 2500 Hz is shown in this plot in order to clearly see the resonance peaks. By looking at the velocity repartition on the surface of the cylindrical shell for the resonance peaks, one can notice that the circumferential modes appear by groups. This phenomenon is linked to the fact that there is a cut-on frequency below which the circumferential orders cannot appear, as explained in [START_REF] Van De Loock | étude de la source de rayonnement acoustique des rsonances gnres dans une coque cylindrique partiellement immerge[END_REF] and on the dispersion curves in section 4. Before placing the non-axisymmetric internal frame in the cylinder, it was hanged with flexible strings and an accelerometer was placed at an arbitrary spot in order to measure the frequency response to an impact hammer. The transfer function is plotted in Fig. 7b. This plot allows estimating qualitatively the number of resonances of the isolated non-axisymmetric internal frame, which appears high in comparison with the number of resonances of the axisymmetric cylindrical shell in the frequency range considered (cf. Fig. 7a).

The comparison of the mean quadratic velocities for the two configurations (with and without the non-axisymmetric internal frame) is plotted as a function of the frequency in Fig. 7c. The first observation is that the vibration level is lower in the non-axisymmetric case, in particular above 500 Hz. The cylindrical shell and the internal structure have the same thickness and are made from the same material. We can expect that the mechanical impedances of the two structures are similar and then adapted to exchange vibratory energy from one 330 structure to another. Besides, the modal density of the internal structure is higher than the one of the cylindrical shell. The combination of these two elements leads to an apparent damping effect on the shell because of energy being trapped in the internal frame. Moreover, new resonance peaks can be observed, in particular at low frequency between 200 and 500 Hz. The peaks 335 frequencies (at 213, 234 and 273 Hz for instance) are directly linked to the resonance of the non-axisymmetric internal frame seen in Fig. 7b.

Radiated power and radiation efficiency

The stationary phase theorem given in Eq. ( 3) is used to calculate the pressure radiated in the far-field using the vibrations measurements on the cylindrical shell in the wavenumber domain. In practice, the sum is calculated for n between -n max and n max . The far-field radiated power is deduced by discretization of Eq. ( 4):

P a = 1 ρ 0 c 0 N ξ i ξ =1 N θ j θ =1 p 2 (R, ξ i , θ j )R 2 sin θ j δξθ max (18)
where (ξ i ) i=1,2,...,N ξ is the vector of azimuthal angles and δξ the step between two consecutive angles. For this example ξ i is varying between -187 • and +187 • with a step of δξ = 3 • The radiated power calculated from the spectral displacements and the stationary phase theorem is plotted as a function of the frequency for the axisymmetric and non-axisymmetric case in Fig. 8a. The influence of the nonaxisymmetric internal frame is mainly seen in the low frequencies (below 1 kHz).

The radiation efficiency [see Eq. ( 5)] is plotted in Fig. 8b. The coincidence frequency for the equivalent plate (h = 1.5 mm) is plotted with a vertical dotted line. It is clear that the non-axisymmetric internal frame tends to increase the radiation efficiency below the coincidence frequency. Differences up to 10 dB can be seen at some frequencies on this configuration. Above the coincidence frequency, it can be said that the structure radiates efficiently with or without the non-axisymmetry. To explain the increase in the radiation efficiency, a focus is done on the operational dispersion curves in the frequency range up to 3000 Hz, where the radiation efficiency increases in the non-axisymmetric case. The experimental operational dispersion curves for this frequency range are plotted for the axisymmetric and non-axisymmetric case in Fig. 9 for the circumferential orders n = 0 and n = 4. It can be seen that in the non-axisymmetric case, the contribution of the circumferential order n = 0 is higher than in the axisymmetric case, while the amplitude of the waves of circumferential order n = 4 is lower.

365

It means that there is a coupling between the circumferential orders because of the internal frame, and that energy is transmitted from higher circumferential orders to lower ones. Besides, it has been shown in Fig. 6a that the circumferential order n = 0 is efficiently radiating, even in the low frequency range. This explains why the cylindrical shell radiates more efficiently when it is coupled to 370 a non-axisymmmetric internal frame. It confirms the observations from numerical simulations [START_REF] Meyer | Prediction of the vibroacoustic behavior of a submerged shell with non-axisymmetric internal substructures by a condensed transfer function method[END_REF]. In the next section, a closer look is taken at the velocity field to understand the influence of the non-axisymmetric internal frame on the repartition of the vibratory energy.

Maps of the radial vibratory field

The velocity maps on the outer surface of the cylindrical shell are plotted in Fig. 10 for the axisymmetric and the non-axisymmetric configurations, at f = 1980 Hz and f = 3580 Hz. These frequencies are chosen because they show an increase in the radiation efficiency in the non-axisymmetric case (2.8 dB at 1980 Hz and 1.9 dB at 3580 Hz) while having a drop in the radiated power (from 82.8 dB to 66.5 dB at 1980 Hz and from 82.5 dB to 80.4 dB at 3580 Hz). Besides, energy transfer through the internal structure can be observed. Looking at the frequency response in Fig. 7a, it can be said that the axisymmetric cylindrical shell has a resonance at 1980 Hz and has no resonance at 3580 Hz (out of the range in Fig. 7a, but not on a resonance peak). The plots in Fig. 10 are functions of the tangential coordinate θ in abscissa and the axial coordinate x in ordinate. The excitation is placed in θ = 0 and x = 1195 mm. A spot of low velocities can be seen around the excitation, because the vibrometer cannot scan the surface behind the shaker. At f = 1980 Hz, the mode (m, n) = (8, 4) is observed for the axisymmetric case, where m is the longitudinal modal index and n the circumferential order. Adding the non-axisymmetric internal structure has two main influences on the cartography at this frequency: the amplitude decreases and the nodes of vibrations are modified.

In Fig. 10c it can be seen that the point force mainly excites the area close to the excitation point. One of the possible explanations is that the cylindrical shell has sections with different stiffeners spacings. This means that there is a mechanical impedance discontinuity between two consecutive sections, which explains the localization effect seen at f = 3580 Hz. In the non-axisymmetric case (Fig. 10d), the energy can be transmitted through the internal structure to the other end of the cylinder, where the stifeners spacing is the same than in the section which is excited (see Fig. 1c).

To see the role played by the low circumferential orders, the radial velocities are plotted in the wavenumber space for the two configurations and the two frequencies in Fig. 11. In each plot, the amplitude has been normalized to the highest value in order to compare the repartition of the radial velocities in the wavenumber space. At f = 1980 Hz, the amplitude maxima appear in the axisymmetric case at the wavenumbers corresponding to the modal wavenumbers. The circumferential order n = 4 is playing a leading role in this case. In the non-axisymmetric case, more circumferential orders play a role, especially in the low circumferential orders. The same observations can be made at f = 3580 Hz in Fig. 11c and11d.

The radiation circle defined by Eq. ( 2) is plotted in white on the plots in Fig. 11. The radiation circles are elliptical due to the plotting aspect ratio. It can be noticed that the components inside the circle, which contribute to the far-field radiated pressure, have a higher relative level in the non-axisymmetric case. The radiation efficiency is thus expected to increase at these frequencies, as shown in section 5.2.

Conclusions

Experiments have been conducted on a stiffened cylindrical shell to investigate the effect of a non-axisymmetric internal frame. A scanning laser vibrometer has been used to measure the radial velocities on a grid on the outer surface of the point-driven cylindrical shell, with and without a non-axisymmetric internal frame. A wavenumber analysis of the vibrations of the shell is used to discuss the vibroacoustic behavior of the shell. It shows that more circumferential orders tend to play a role when the cylindrical shell is coupled to the non-axisymmetric internal frame. More particularly energy is transmitted from high to low circumferential orders. As the low wavenumbers are in the supersonic domain, the structure tends to radiate more efficiently. This phenomenon explains the increase in radiation efficiency below the coincidence frequency seen in the non-axisymmetric case.

Further work consists in conducting the same experimental process at higher frequencies. It will allow studying the influence of the non-axisymmetric internal frame on the propagation of Bloch-Floquet waves due to periodic stiffeners. Besides, if the model is submerged in water, the coincidence frequency is about 20 times higher. In this case, it is likely that the effect of the internal frame is 435 seen on a wider frequency range. Moreover, the strong coupling with the fluid can be investigated.

Figure 1 :

 1 Figure 1: Pictures of the axisymmetric stiffened cylindrical shell: (a) hanged on an arm and (b) interior view. (c) Sketch of the cylindrical shell with the sections of different stiffeners spacings.

Figure 2 :

 2 Figure 2: Pictures of the non-axisymmetric internal frame: (a) top view and (b) bottom view.

Figure 3 :

 3 Figure 3: Interior view of the stiffened cylindrical shell with the non-axisymmetric internal frame.

Figure 4 :

 4 Figure 4: Sketch of the experimental setup.

Figure 5 :

 5 Figure 5: Picture of the experimental setup

Figure 6 :

 6 Figure 6: Operational dispersion curves plotted from the spectral radial velocities (dB ref 1 m.s -1 ) in the axisymmetric case. (a) Theoretical, n = 0. (b) Experimental, n = 0. (c) Theoretical, n = 1. (d) Experimental, n = 1. (e) Theoretical, n = 4. (f) Experimental, n = 4. (g) Theoretical, n = 6. (h) Experimental, n = 6.

Figure 7 :

 7 Figure 7: (a) Mean quadratic velocity as a function of the frequency for the axisymmetric stiffened cylindrical shell. (b) Transfer function between an impact point on the top of the floor and an accelerometer placed on the side. (c) Comparison of the mean quadratic velocities for the two configurations.

Figure 8 :

 8 Figure 8: Estimation as a function of the frequency for the two configurations using the stationary phase theorem of (a) the radiated power (dB ref 1e-12 W) and (b) the radiation efficiency (dB ref 1).

Figure 9 :

 9 Figure 9: Operational dispersion curves plotted from the spectral radial velocities (dB ref 1 m.s -1 ) in the low frequency range. (a) Axisymmetric case, n = 0. (b) Non-axisymmetric case, n = 0. (c) Axisymmetric case, n = 4. (d) Non-axisymmetric case, n = 4

Figure 10 :

 10 Figure 10: Map of radial velocities (dB ref 1 m.s -1 ) on the outer surface of the cylindrical shell. (a) Asymmetric case at f = 1980 Hz. (b) Non-axisymmetric at f = 1980 Hz. (c) Asymmetric case at f = 3580 Hz. (d) Non-axisymmetric at f = 3580 Hz.

Figure 11 :

 11 Figure 11: Map of radial velocities (dB ref 1 m.s -1 ) in the wavenumber domain (kx, n). (a) Asymmetric case at f = 1980 Hz. (b) Non-axisymmetric at f = 1980 Hz. (c) Asymmetric case at f = 3580 Hz. (d) Non-axisymmetric at f = 3580 Hz.

Table 1 :

 1 Stiffened cylindrical shell dimensions (mm).

	Parameter	Notation Value
	Radius	R0	100
	Length	L	1500
	Shell thickness	h	1.5
	Stiffeners height	hw	5
	Stiffeners thickness	tw	1
	Stiffeners spacing 1	α	24
	Stiffeners spacing 2	β	20
	Stiffeners spacing 3	γ	17

Acknowledgment

This work was funded by DCNS and performed within the framework of the LabEx CeLyA of Université de Lyon, operated by the French National Research

Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007). The authors are grateful to Patrick Blachier for his great contribution to the experimental setup.


D Chatain

C Mariette

V S Nikolayev
email: vadim.nikolayev@espci.fr

D Beysens

Quench cooling under reduced gravity

Keywords: numbers: 47.55.np, 68.03.Fg, 47.15.gm

We report quench cooling experiments performed with liquid O2 under different levels of gravity as provided by a magnetic gravity compensation setup. A copper disk is quenched from 270K to 90K. It is found that when gravity is zero, the cooling time is abnormally long in comparison with any other gravity level. This phenomenon can be explained by the insulation effect of the gas surrounding the disk. The liquid subcooling is shown to drastically improve the heat exchange, thus reducing the cooling time (about 20 times). The effect of subcooling on the heat transfer is analyzed at different gravity levels. It is shown that such type of experiments cannot be used for the analysis of the critical heat flux (CHF) of the boiling crisis. The minimum heat flux (MHF) of boiling is analyzed instead.

I. INTRODUCTION

The rocket engines for spatial applications need to be able to restart under low gravity conditions to place satellites on different orbits. The solid fuel engines are currently used for this purpose. The cryogenic fuel like hydrogen and oxygen is however much more advantageous energetically and is currently under study for next generations of the rocket engines. There is however an important limitation. After some period of inactivity, the temperature of the fuel injectors of such engines may rise above the temperatures at which the fuel may be in the liquid state because of the solar radiation. The engine cannot function before the injectors are not chilled down. The injectors are usually cooled by putting them into the contact with the cryogenic fluids (i.e. with the fuel), i.e. by the quench cooling. One thus needs to master the cooling time as a function of different parameters. The understanding of the latter phenomenon under reduced gravity conditions thus becomes important. In the present article we study the quenching of a hot metal piece in cryogenic fluids as a function of several parameters like the gravity level.

The prediction of cooling time necessitates the knowledge of the heat transfer coefficient which depends on temperature. To measure the influence of the gravity level on the value of the heat transfer coefficient, we have performed experiments in a liquid oxygen-filled cell. A copper disc was quickly immersed into this bath and cooled down (quenching method). The Earth gravity was compensated by the magnetic force in the OLGA (Oxygen Low Gravity Apparatus) facility located at CEA/Grenoble in France.

The objective of this article is three-fold. We aim to show that (i) one cannot measure the critical heat flux of the boiling crisis in a quenching experiment and (ii) the subcooling improves greatly the cooling efficiency in microgravity that is otherwise very poor. Finally, we measure the dependence of the minimum flux of boiling (i. e. the threshold of the departure from film to nucleate boiling) as a function of gravity. Consider what happens when a hot sample is dipped into a liquid. Fig. 1 presents a typical temperature evolution of a copper disk initially at about 270K dipped into a liquid O 2 bath at 90K under Earth gravity.

At a first stage (a-b zone), the temperature decreases slowly because the solid surface is completely covered by a layer of vapor which insulates the solid from the liquid. This is the film boiling regime, during which the vapor bubbles are generated at the vapor-liquid boundary. When the surface temperature becomes low enough, the liquid may come into direct contact with the sample and the bubbles start nucleating directly at the solid surface. This signifies the beginning of the nucleate boiling regime which provides much more efficient heat transfer and the temperature decreases very quickly (b-c zone). At the final stage, the solid temperature is not any more sufficiently high to induce boiling; the solid temperature approaches that of the liquid O 2 bath due to convection. The evolution of heat transfer regimes can be traced on the Nukiyama diagram (Fig. 2) presenting the heat flux q from the heater as a function of the sample temperature T . The system evolves from the right (film boiling regime) to the left in Fig. 2. The moment where the direct solid-liquid contact begins to occur corresponds to the minimum heat flux (MHF). From this moment, the transition from film to nucleate boiling starts. The line along which the system evolves depends on the thermal inertia of solid. If it was vanishingly small, a transition to the nucleate boiling regime would occur instantaneously (leftwise directed arrow in Fig. 2). In reality, a finite time is necessary to cool down the solid and the system evolves along the dotted curve. The heat flux attains a maximum value at a point where no part of the solid is covered by the persisting vapor film. The nucleate boiling regime occurs from that moment on. It is very important to distinguish the quenching curve from the boiling (solid) curve. The maximum heat flux (called critical heat flux, CHF) can only be recovered in the stationary boiling experiments (see, e.g. [1]) performed at fixed heat flux during each experiment. Till 1960s, the quenching experiments were considered to be suitable for the CHF determination [2,3]. With the development of the quenching stud-ies for the metallurgy (steel quenching), it became clear that only MHF may be obtained by this method [4], the maximum flux value depending strongly on the sample thermal inertia. The actual CHF may only be obtained for the samples of the large enough thermal inertia. However the CHF may hardly be studied in that case, which is evident from the following consideration.

During the cooling, the heat Q transferred from the sample to the liquid bath during the time ∆t may be expressed as

Q = mC P (T ) (T (t) -T (t + ∆t)) , (1)
where m is the mass of the disk, T its temperature, and t is the time. The specific heat of copper C P depends on T ,

C p (T ) = -1.355 × 10 -7 T 4 + 1.303 × 10 -4 T 3 -4.798 × 10 -2 T 2 + 8.331T -217.964
The transferred heat may also be expressed with the heat transfer coefficient h,

Q = hS (T -T L ) ∆t, ( 2
)
where S is the total surface of the sample and T L is the temperature of the liquid bath. By equating these expressions, one finds

h = mC P [T (t)][T (t) -T (t + ∆t)] S∆t[T (t) -T L ]
and the heat flux q = Q/(S∆t). This method of measurements is valid only if the temperature of the sample is uniform during the cooling. For this, the Biot number Bi = hL c λ must be lower than 0.1 [START_REF] Kreith | Transmission de la chaleur et thermodynamique[END_REF]. Here L c is the characteristic length (a half thickness of the disk); λ is the copper thermal conductivity. It has been checked that the criterion Bi < 0.1 holds for all described experiments. One can see now that the requirements of small Bi < 0.1 and large thermal inertia of the sample are contradictory. This explains why the quench method is hardly suitable for the CHF study. The objective of this article is the analysis of the MHF rather than of CHF. Besides we study the cooling dynamics and the impact of subcooling on it at various gravity levels.

II. EXPERIMENTAL APPARATUS

Experiments were performed in the OLGA facility (Fig. 3). It is equipped with a superconductive coil and creates a magnetic force strong enough to compensate the buoyancy force in O 2 . Oxygen is contained in the experimental cell (Fig. 4) placed in an anti-cryostat inserted into the 300 mm bore of the superconductive coil. Two endoscopes are used for observation: one for the light and another for the camera.

The lower part of the experimental cell is made of a sapphire tube of 30 mm diameter and 100 mm length closed by two copper flanges equipped with heat exchangers, inside which the gaseous cold helium may flow. The temperature of the two flanges is controlled within 0.01K with heaters. The upper part of the cell is the stainless steel tube. The pressure in the cell is measured with an accuracy of 10 mb. Liquid oxygen is obtained by condensing pure O 2 gas (99.995%) in the cell. The cell is brought to the desired temperature by regulating the power of heaters and helium flow rate in the heat ex-changers that are situated at the top and at the bottom of the transparent part of the cell.

The disk is made of pure copper. Its diameter is 20.02 mm and its thickness is 2.98 mm. The weight of the disk is 7.68 g, its surface in contact with liquid is S = 8.07 cm 2 and its RMS roughness is 0.8 µm. A hole was drilled to introduce a platinum thermometer (supplier: Heinz Messtechnic, model W60-50, class A; diameter 1.8 mm, length 11 mm) in its heart. After each quench experiment, the disk may be risen into a hot zone that situates 0.6 m above the transparent part of the cell with a screw/nut system driven by a stepper motor. The hot zone of the cell is regulated at 270 K. Once the disk attains this temperature, the grip is loosened again and the disk drops into the liquid. Meanwhile, the disk temperature is acquired at a frequency of 100 Hz.

III. GRAVITY COMPENSATION

The principle of magnetic compensation of gravity in a pure non-magnetic substance is an application of the magnetic volume force opposite to the Earth gravity [START_REF] Beaugnon | [END_REF][7][8].

This force per unit volume is given by the expression

F m = χ m 2µ 0 ∇(B 2 ), ( 3
)
where B is the magnetic field, µ 0 is the magnetic permeability of vacuum and χ m the magnetic susceptibility of the material. Since χ m is proportional to the density ρ of the nonsubstance, one may introduce a quantity α = χ m /ρ independent of the density. It means that the buoyancy force is compensated if one applies the field with dB 2 /dz = 2µ 0 g/α, where g is the Earth gravity acceleration and z is the vertical coordinate. For O 2 at 90K, 2µ 0 g/α ∼ 8 T 2 /m. These conditions are achieved in OLGA facility near the bottom of the superconductive coil.

It was demonstrated [9] that the magnetic force (3) cannot be made constant in a volume. In our experimental cell, Earth gravity may be exactly compensated at a single point called complete compensation point. Such conditions are called hereafter 0g and the corresponding coil current is denoted I 0g . In the vicinity of this point, a residual acceleration acts on O 2 [10]. As the magnetic force is a way too small to compensate the gravity for copper, the disk falls with the Earth gravity acceleration. Its stop point situates at z = -242 mm with respect to the center of the coil. At this location, Fig. 5 gives the calculated iso-values of residual gravity acceleration modulus at 0g conditions. One can see that at the edge of the disk, the residual acceleration is smaller than 2% of g.

To achieve 0g at the required place inside the cell, a single bubble is generated by the boiling process with a heater situated at the bottom of the cell. Then the current in the coil is adjusted until the bubble center coincides with the optimal previously calculated point (Fig. 6).

The effective gravity level [8] may be varied by changing the coil current I. It is known [11] that the effective gravity level g * at the geometrical point of complete compensation is given by the equation

g * = 1 - I 2 I 2 0g
g.

FIG. 6. Image of O2 bubble at 0g inside the cell at equilibrium as seen in the periscope mirror used to adjust the position of the exact compensation point with respect to the cell. For this case, the temperature is 92.24K and the current in the coil is I0g = 241 A. The elongated shape of the bubble is due to the weak coupling of the magnetic field with the gas-liquid interface shape [8].

The magnetic susceptibility of oxygen varies with temperature, and I 0g depends weakly on the liquid temperature. For the isothermal conditions, the currents for the complete gravity compensation are I 0g = 239.7 A at 90.07K (pressure of 1 bar, at saturation) and I 0g = 241 A at 92.24K (subcooling 5K for the pressure of 2 bars). When the disk is dropped into liquid O 2 , the liquid is slightly heated up, resulting in a small decrease of the magnetic susceptibility. By observing the movement of the gas generated by boiling around the disk (Fig. 7), we could see that the current had to be slightly increased to obtain 0g in the cell during the disk cooling. In this case, I 0g = 242.1 A at 90.07K, and I 0g = 245.18 A at 92.24K. Quenching experiments were performed for two O 2 thermodynamic conditions: saturated O 2 and subcooled O 2 with respect to its saturation temperature for given pressure that was controlled independently since the cell was connected to a large oxygen-filled vessel.

A. Case of saturated liquid

For saturated O 2 experiments, the temperature of the cell was regulated at 90.07K and the vessel pressure was maintained at 1 bar (which is the saturation pressure corresponding to this temperature). 15 quenches were performed.

1. Cooling curves of the disk Fig. 8 shows the cooling curves of the disk for different levels of gravity. We can see that the lower the gravity, the larger is the cooling time. The 0g cooling time is very large with respect to other cases and the corresponding curve is truncated: it took more than 10 minutes to cool down the disk. The complete curve is shown in the inset to Fig. 13 below.

Heat transfer curves

Using the above described method, we have calculated the heat transfer curves for different levels of gravity (Fig. 9). The film boiling points were calculated by making a centered moving average on 50 points for the film boiling and from 10 to 20 points for the transitional/nucleate boiling. One can clearly see the influence of the gravity level on the transferred heat flux in film boiling regime and on the maximum and minimum heat flux values. [START_REF] Brentari | Boiling heat transfer for oxygen, nitrogen, hydrogen and helium[END_REF] with our experiment under 1g. One can see that the quench method gives the maximal heat flux ∼ 12 W/cm 2 ), which is two times lower than the CHF value 24.9 W/cm 2 obtained with the Kutateladze correlation. The latter value reasonably agrees with the experimental value (22-23 W/cm 2 ) obtained in pool boiling experiments [START_REF] Lyon | [END_REF]. The low value of the maximum heat flux reflects the transient nature of the quench cooling as discussed in Fig. 2(sec. I A).

B. Case of subcooled liquid

To measure the subcooling influence, experiments where performed for several subcooling values. The subcooling was obtained with the following method (cf. Fig. 11). The temperatures of two lower heat exchangers of the cell (cf. Fig. 4) are regulated at the same value T 1 . The temperature T 2 of the upper heat exchanger is regulated at a higher value equal to the saturation temperature corresponding to the gas pressure. Under these conditions, the equilibrium level of liquid forms in the region of cell maintained at the temperature T 2 . A thermal gradient exists in the liquid between two upper heat exchangers. The liquid is isothermal at T 1 in between two lower heat exchanges. The final drop point of disk is situated in this initially isothermal region.

To achieve the 5K subcooling, T 1 = 92.24K and T 2 = 97.24K. The latter temperature corresponds to the pressure of 2 bars. For the 7K subcooling (see below), T 1 = 90.24K.

Cooling curves

Fig. 12 shows the raw cooling curves of the disk for 6 levels of gravity. The temperature of cell was regulated at 92.24K and the pressure was adjusted at 2 bars which results in a subcooling of 5K. The temperature acquisition frequency was 100 Hz. Similarly to the saturated case, the lower the gravity level, the longer is the cooling time. Unlike the saturated case, the 0g cooling time is only twice larger than that of the 1g case.

Impact of subcooling

Fig. 13 shows that the subcooling reduces drastically (20 times) the general cooling time at 0g while the effect at any nonzero gravity is much weaker (cf. Figs. 8 and12). Such a reduction may be explained as follows. Fig. 7b shows the formation at 0g of a vapor bubble that surrounds the solid and completely insulates it from the liquid. The bubble can recondense if the liquid is subcooled and cannot recondense if the liquid is saturated. For the saturated case, the heat exchange is thus produced via the weak heat conduction through the vapor and through the almost non-conductive at cryogenic temperatures stainless steel rod that supports the disk (cf. Fig. 4). When the liquid is subcooled, the recondensation of the vapor at the external bubble interface is very efficient and the gas layer surrounding the disk is much thinner. This effect improves the heat exchange thus favoring cooling.

Heat transfer curves

We have plotted the heat transfer curves for 6 levels of gravity (Fig. 14) using the same method as the one used for the saturated case. One can see that unlike the saturated quench, the maximum and the minimum heat flux values are both nonzero for the 0g case.

C. Minimum heat flux

Fig. 15 shows the minimum heat flux as a function of the gravity level. The minimum value was obtained by fitting the raw data with a second degree polynomial. The MHF occurrence is usually associated with the Taylor instability of the vapor film [14,15] so that q MHF ∼ g 1/4 should be satisfied. Fig. 15 shows an exponent close to this value.

V. CONCLUSION

Quench cooling has been analyzed under different gravity levels. Earth gravity was compensated with magnetic forces. Measurements were performed by quenching a copper disk of 20 mm diameter and 3 mm thickness from 270K to 90K. The analysis of different boiling regimes has shown that the study of the critical heat flux (CHF) via quenching experiments is hardly possible because of the transient nature of the quench cooling: the maximum heat flux value is twice smaller than the CHF measured by other authors in the stationary boiling experiments.

It has been shown that cooling via quenching in microgravity takes usually exceedingly long time because the vapor generated during the boiling envelopes completely the cooled body (film boiling regime), thus thermally insulating it from the liquid. It is shown that the artificial subcooling (for instance, by transient pressurization) can drastically accelerate the cooling. This speeding up is explained by the recondensation of the vapor envelope at the vapor-liquid interface in case of subcooling.

The dependence of the heat transfer and, in particular, of the minimum heat flux of boiling (MHF) on gravity is analyzed. As expected, the heat transfer improves with the gravity level. The MHF gravity dependence shows an exponent close to that of the Zuber theory (exponent 1/4).

FIG. 1 .

 1 FIG. 1. Temperature evolution of a copper disk initially at about 270K dipped into O2 bath at 90K under Earth gravity.

FIG. 2 .

 2 FIG. 2. (Color online) Schematic heat transfer curves for stationary boiling experiments and transient quenching.

FIG. 3 .

 3 FIG. 3. (Color online) A sketch and a general view of the OLGA facility.

  FIG. 4. (Color online) A sketch of the cell and a photo of its lower part (the upper part of the cell is removed). The disk is shown in the upper position on the sketch and in the lower on the photo.

FIG. 5 .

 5 FIG. 5. (Color online) Iso lines of the residual acceleration modulus in % of g inside the cell of OLGA at the ending point of the disk fall at 0g.

FIG. 7 .FIG. 8 .

 78 FIG. 7. (Color online) Film boiling at different gravity (its direction is shown with vertical arrows) at initial liquid O2 temperature (before disk drop) 92.24K. a) I = 241 A: the residual gravity is directed downwards and the gas surrounding the disk slowly rises. b) I = 245.18 A: the gas recondenses without rising and surrounds the disk during all the cooling time (zero gravity). c) I = 252.7 A: the residual gravity is directed upwards and the gas goes downwards.

FIG. 9 .

 9 FIG. 9. (Color online) Heat transfer curves for saturated O2 quench cooling at 1 bar and different gravity levels.

Fig. 10

 10 Fig. 10 is a comparison of Kutateladze and Breen &Westwater correlations[START_REF] Brentari | Boiling heat transfer for oxygen, nitrogen, hydrogen and helium[END_REF] with our experiment under 1g. One can see that the quench method gives the maximal heat flux ∼ 12 W/cm 2 ), which is two times lower than the CHF value 24.9 W/cm 2 obtained with the Kutateladze correlation. The latter value reasonably agrees with the experimental value (22-23 W/cm 2 ) obtained in pool boiling experiments[START_REF] Lyon | [END_REF]. The low value of the maximum heat flux reflects the transient nature of the quench cooling as discussed in Fig.2(sec. I A).

FIG. 10 1 T 1 T 1 FIG. 11

 1011111 FIG. 10. (Color online) Comparison of the Kutateladze (broken line) and Breen&Westwater (circles) correlations with our data from Fig. 9 corresponding to 1g (solid line).

FIG. 12

 12 FIG. 12. (Color online) Evolution of the disk temperature at different gravity levels for the 5K subcooled case.

FIG. 13

 13 FIG. 13. (Color online) Evolution of the disk temperature at different subcooling levels at 0g. The same curves are shown in the inset in the semi logarithmic scale.

FIG. 14

 14 FIG. 14. (Color online) Heat transfer curves for the 5K subcooling case at different gravity levels.

FIG. 15 .

 15 FIG. 15. (Color online) Minimum heat flux versus gravity level for saturated and subcooled O2. Inset: same in log-log scale; the 1/4 slope is shown for comparison.

ACKNOWLEDGMENTS

The authors gratefully acknowledge Jerôme Chartier, Jean-Marc Mathonnet, Patrick Bonnay and Stephane Garcia (CEA/SBT) for their contribution throughout this work.

This work was partially supported by CNES.


Puechmorel Delahaye Daniel
email: puechmor@recherche.enac.fr

A New Genetic Algorithms Working on State Domain Order Statistics

Keywords: order statistics, state domains, genetic algorithm

This paper presents a new concept of Genetic Algorithm in which an individual is coded as a domain of the state space and is evaluated with the help of order statistics. For this first version only continuous criteria has been investigated. An hypercube domain of the state space is associated with each individual and is randomly sampled according to a distribution for which asymptotic extremes are known. Regular fitnesses are computed for all the samples in each domain and are combined to produce a prospectiveness criterion. A regular GA and this new GA are compared on classical N dimensional functions such as Sphere, Step, Ackley, Griewank for different values of N. A final comparison is given on the classical Lennard-Jones Molecular Conformation problem with 30 atoms. For both versions, a regular GA has been used; the first one works on state points and the other one on state domains. For all tests, and for the same number of criterion evaluations, this new algorithm performs much better than the classical one.

Introduction

At the beginning of GA, binary coding was mainly used to encode optimization problem (a good description of this kind of coding may be found in [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF]). This coding is easy to manipulate and is very adapted for discrete problems with binary decision variables. The schema theory can be applied to this kind of GA and produce some convergence theorems. Unfortunately, this coding is not adapted for real vector optimization problems and direct real coding has been developed [START_REF] Michalewicz | Genetic algorithms + data structures = evolution programs[END_REF]. Usually GAs work on state points and evaluate fitnesses on those points with succession of exploration and exploitation phases. The exploitation phase is guided by the selection and the exploration phase is guided by the crossover and mutation operators.

In order to enforce and speed up those two phases, the present algorithm works on state domains instead of state points. The principle is the same as the one used for the Branch and Probability Bound [START_REF] Zhigljavsky | Theory of global random search[END_REF]. For this B&B, the global state domain is split into smaller parts which are evaluated with the help of order statistics. The probability that the global optimum belongs to a domain is then used to guide the branching phase.

The present GA works the same way by coding each individual by an N dimensional hypercube and uses order statistics to produce the associated prospectiveness and then guide the exploration or the exploitation.

In a first part, a brief description of order statistics is given. The second part describes how order statistics may be used to strongly enhance the performance of a regular GA and how coding, fitness evaluation and operators have been implemented. Finally, the third part gives some comparisons of both algorithms on hard to optimize mathematical functions.

Order Statistics

Let (x 1 , . . . x N ) be a sample of N iid random variables. The orders statistics associated with that sample is ( 1 , . . . , N ) where i is the i-th largest value of (x 1 , . . . , x n ). It's well known from classical extreme statistics theory [START_REF] Gumbell | Statistics of extremes[END_REF] that under mild assumptions the distribution of the max of N iid random variable with common density function will converge point wise to one of the Frechet, Weibull or Gumbell distributions. More precisely, one may find two sequences (a n ) n∈N , a n > 0 and (b n ) n∈N such that :

∀x lim n→+∞ F n (a n x + b n ) = G(x)
where G is one of :

G 1α : x → exp(-x -α ) Frechet G 2α : x → exp(-(-x) -α ) Weibull G 3 : x → exp(-exp(-x)) Gumbell
where α is a positive real number named "Tail index". In most cases α has to be estimated, which requires a considerable amount of samples to be accurate, but for optimization purpose we will see that this value can be set a priori. Now, if we let M (F ) = inf {x|F (x) < 1} we have the following result [START_REF] Reiss | Approximate distributions of order statistics[END_REF] :

-If M (F ) = +∞ , lim t→+∞ [1 -F (tx)]/[1 -F (t)] = x -α for x > 0 -If M (F ) < +∞ , lim t→0 [1 -F (M (F ) -tx)]/[1 -F (M (F ) -t)] = (-x) α for x < 0 -As a limiting case , lim t→M (F ) [1 + F (T + g(t)x)]/[1 -F (t)] = e -x for x ∈ R
and :

g(t) = M (F ) t (1 -F (u))du/(1 -F (t))
which gives necessary and sufficient conditions for weak convergence of distributions.

If we now go back to the problem of maximizing a function f : R n → R and assuming that f can be approximated by a definite negative quadratic form in a sufficiently small ball centered at each point of its domain, then one can prove that conditions yielding to a Weibull limiting distribution are fulfilled and that α = 2/n. If the function f is exactly a definite negative quadratic form :

f : x → x, HX
and if samples are drawn according to a probability measure µ, then we have :

F : t → 1 -µ( t QΛ -1/2 B(0, t 1/2 ))
where Q and Λ are respectively orthogonal and diagonal matrices arising from the decomposition :

H = t QΛQ
and the notation B(0, t) stands for the open ball of center 0 and radius t. Now by linearity, we have directly the de Haan condition [START_REF] De Haan | On regular variation and its application to weak convergence of sample extremes[END_REF] :

(1 -F (ut))/(1 -F (t)) = u n/2
for u > 0. Of course, this result can be straightforwardly extended to the case of a non zero maximum M (F ), yielding again to a Weibull limiting distribution for the maximum :

lim N →+∞ F N (M (F ) + a N x) = exp(-(-x) 2/n ) with : a N = M (F ) -inf {t|1 -F (t) ≤ N -1 } + o(1)
If the function f has several equivalent maxima, each of them satisfying the de Haan condition locally, the Weibull distribution with 2/n as tail index is again the limiting distribution of the maximum.

In the following, we will assume that the previous quadratic model is valid, at least locally so that we may take α = 2/n. Note anyway that after a sufficient number of samples has been drawn, it's possible to have a sufficiently accurate estimator of α which may speed up the algorithm. This refinement has not been implemented in our test algorithm but will be in the final version. Several estimators exists for α and can be found in [START_REF] Zhigljavsky | Theory of global random search[END_REF].

Prospectiveness

In the following sampling will be done on an hyper-rectangle D defined by its principal diagonal (x 1 , x 2 ) with x 1 , x 2 ∈ R n . Let ( 1 , . . . , N ) be the order statistics associated with sample (f (x 1 ), . . . , f (x N )), the points x 1 , . . . , x N being drawn from D using uniform sampling (Markovian sampling may be used too and gives the same limiting distribution). The prospectiveness criterion of the domain D is an estimation of the probability that the true maximum of f occurs in D. It can be computed using the first order statistics by :

c(D) = 1 - M -1 M -k 2/n k
with k = min(5, N/10) (this value is based on numerical experiments. see [START_REF] Zhigljavsky | Theory of global random search[END_REF]) and M the maximum value of f observed so far (this value is updated after each generation of the GA).

A high prospectiveness indicates that the domain must be exploited while a low value shows that the domain is unpromising and must be either dropped or expanded.

GA and order statistics

Introduction

The general scheme of our GA is quite the same as a regular GA: it first generates an initial population, applies a selection to identify the best individuals and diversify the population by applying operators. The main difference come from the coding, and the fitness evaluation.

This GA works on state domains instead of state points and is able to identify the prospectiveness of a domain by computing order statistics. Then a low fitness will be given to a domains with a low prospectiveness. The power of order statistics enable to evaluate "large" domains with a "few" samples on it and produces a fitness which summarize the property of an entire zone.

It may be noticed that this improvement may be adapted to any GA which is working on continuous state space. It is independent of the selection, the scaling, the sharing etc ....

The main adaption have to be done on the coding, the operators and on the fitness evaluation.

Coding

The chromosomes used for our GA are hypercubes of a N-dimensional state space which are encoded with two points on a diagonal (see figure 1)

After the building of a chromosome, only the field P 1 and P 2 are updated. The other fields will be addressed by the evaluation of the associated fitness and will be used by the operators in order to enhance the exploitation and the exploration. During the evaluation, the domain is sampled and order statistics are computed.

State Domain

Fig. 1. Hypercube domain coding

Operators

Crossover Different crossover operators have been implemented but the one which produce the best results uses the geometrical properties of the parent domains (see figure 2).

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ P P C 1 C C C 1 P 1 P 1 2 2 2 2

Fig. 2. Crossover operator

Two parent domains P 1 and P 2 are randomly selected and 2 situations have to be investigated:

1. the two parents share a common sub domain (left part of figure 2). This shared sub domain (shaded area) becomes one of the children (C1) and the second one is given by the smaller domain which encompass the two parents (dashed line). 2. the two parents are independents (right part of figure 2). An "inter-subdomain" is then built (shaded area) to create the first child C 1 and, as in the first case, the second child is given by the smaller domain which encompass the two parents (dashed line).

Depending on the parents, this crossover enhance both, exploration and exploitation.

Mutation For 20% of mutation a full new drawing of the domain is applied and in the others cases the mutation operator is guided by the prospectiveness C which is computed during the fitness evaluation. C ∈ [0, 1], and is maximum when the confidence about the computed max is maximum. So, when C is close to "1" one has to enforce the exploitation (centering on P m and contraction). On the other side, when C is close to "0", the exploration has to be enforced (the operator draws a new domain). Finally when C is not on extrema (0 or 1) a compromise between exploration and exploitation is done (centering and extension).

An example of the three previous situations is given on figure(see figure 3).

Pm

Fitness Evaluation

In order to apply the selection, individuals have to be evaluated with a fitness function. When a chromosome has to be evaluated, its domain is randomly sampled (30 draws) with the original criterium function (the one we want to optimize). The order statistics and the associated confidence are then computed.

The fitness given to the chromosome is then computed with the help of the observed max and the confidence.

Test Functions

Different test functions have been used in order to compare our method with a classical GA :

-Sphere function :

f 1 (x) = i=N i=1 x 2 i , -50.0 ≤ x i ≤ 50.0 -Step function : f 2 (x) = i=N i=1 x i + 0.5 2 , -50.0 ≤ x i ≤ 50.0 -Ackley function : f 3 (x) = -c 1 . exp (S 1 (x)) -exp (S 2 (x)) + c 1 + e with S 1 = -c 2 1 N N i=1 x 2 i S 2 = 1 N N i=1 cos (c 3 .x i ) c 1 = 20 c 2 = 0.2 c 3 = 2π -30.0 ≤ x i ≤ 30.0 -Griewank function :f 4 (x) = 1 400 * N N i=1 x 2 i - N i=1 cos ( xi √ i ) + 1 -600.0 ≤ x i ≤ 600.0 -Rosenbrook function :f 5 (x) = N -1 i=0 100 * (x 2 i -x i+1 ) 2 + (1 -x i ) 2 -30.0 ≤ x i ≤ 30.0
-Lennard-Jones function : This function is coming from the famous Lennard-Jones Molecular Configuration problem: the problem of finding the structure or relative positions of a cluster of atoms that minimizes the potential energy of the structure.The Lennard-Jones problems assume that the potential energy of the molecule is given by the sum of the pairwise interaction between atoms. The position of atoms being given in the three dimensional space, a problem with K atoms has N = 3.K real variables to be optimized. The potential energy is then given by the following function:

f 6 (x) = N/3 i=0 i-1 j=0 1 d 6 ij -2. 1 d 3 ij d ij = (x i -x j ) 2 + (x i+1 -x j+1 ) 2 + (x i+2 -x j+2 ) 2 min f 5 (x) = -128.287 N = 90(30 atoms) -2.0 ≤ x i ≤ 2.0
All the function have to be minimized and have their minimum at 0 unless the Lenard-Jones function for which only an experimental min is used (best known min=-128.287 for 30 atoms). It must be noticed that both algorithm use the same selection scheme (stochastic remainder without replacement which is not the best) and do not use any scaling or sharing operators. From this point of view both algorithm may be still enhanced.

Our goal being to compare the influence of domain chromosome and order statistics we wanted them to work exactly the same way form the selection point of view.

The number of evaluations being different at each generation for those two algorithms, the number of generation has been adapted in order to maintain the same number of evaluations for all experiments.

It must be noticed that the following curves have been adjusted in order to represent both result on the same graph. Those adjustments have been done on both axis. The "x" axis address the number of evaluations for our GA and must be scaled for the standard GA (x 20). The "y" axis represent the fitness given by both algorithms. The given results given are so differents that a logarithm scale has been used to see both curves.

The parameters used for our GA are the following:

individuals 100 generations 500 probability of crossover 0.4 probability of mutation 0.3 For the Rosenbrook, Lennard-Jones the number of generation has been extended to 1500 and 2500 respectively. The experiments have been done on a PentiumII 300 MHz and last 7 minutes for N=200 and 14 minutes for N=2000 (N:dimension of the state space). It must be noticed that other experiments has been done for the same functions with the optimum moved in the state space (without symmetries) and the given results are quite the same.

Function f 1 f 2 f 3 f 4 Standard AG -N=200

Conclusion

This paper has given a good application of the gain which could be given by the mix of different methods. On one side, the main advantage of order statistics for optimization is their ability to summarize the properties of an entire domain with a "small" sample. On the other side, the evolution process of GA is able to build the most adapted chromosome to environment given by the fitness landscape. The mix of both methods really increase the performances of GA by guiding the exploration and exploitation phases. For all tests, the results given by this new GA, are much better than the ones given by a standard GA. It must be noticed that this algorithm may be still improved the following way:

-a better selection scheme may be used; -the order statistics may also be used in order adapted the drawing random law in the domains; -pools of samples may be used for different domains which could be randomly open on it. Those pools may be updated every K generations.

Fig. 3 .

 3 Fig. 3. Mutation operator

-5 are the five first statistics, P m is the position of the max in the domain and C is the relative confidence of the first statistic.


D Delahaye

S Puechmorel

Air Traffic Complexity Versus Control Workload

T his paper presents a new air traffic complexity metric based on dynamical systems. Based on a set of radar observations (position and speed) a vector field interpolating these data is constructed. Once the field has been obtained, the Lyapunov spectrum of the associated dynamical system is computed on points evenly spaced on a spatial grid. The results of the computations are summarized on complexity maps, with high values indicating areas to avoid or to carefully monitor. A first approach based on linear dynamical system enable to compute an aggregate complexity metric. In order to produce complexity maps, an extension of the previous approach have been developed. Based on such complexity metric, a control workload predictor has been developed validated on a real operational airspace.

Introduction

In a control sector, the higher the number of aircraft, the more the control workload increases (in a non-linear manner). A limit exists after which the controllers in charge of a control sector are unable to accept additional aircraft, obliging these new aircraft to travel around the sector, moving through less charged neighboring sectors. In this case, the sector is said to be saturated. This critical state should be avoided, as it provokes a cumulative overloading phenomenon in preceding sectors which can back up as far as the departure airport. The saturation threshold is very difficult to estimate, as it depends on the geometry of routes traversing a sector, D.Delahaye MAIAA , French Civil Aviation University, 7 Ave Edouard Belin 31055 Toulouse FRANCE, email: delahaye@recherche.enac.fr S.Puechmorel MAIAA , French Civil Aviation University, 7 Ave Edouard Belin 31055 Toulouse FRANCE, email: stephane.puechmorel@.enac.fr the geometry of the sector itself, the distribution of aircraft along routes, the performances of the control team, etc. One widely accepted threshold is fixed at 3 conflicts and 15 aircraft for a given sector. This maximum load should not last for more than ten minutes as it places the controllers under considerable stress, with the risk that they will no longer be able to manage traffic in optimal safety conditions.

The control workload measurement is critical in many domains of ATM as it is at the heart of optimization processes. Examples include the following applications:

• Airspace comparison (US/Europe).

• Validation of future concepts (SESAR,NEXTGEN,etc.).

• Analysis of traffic control action modes (situation before and after control).

• Optimization of sectorization.

• Optimization of sector grouping and de-grouping (pre-tactical alert: anticipation of an increase in congestion in a group of sectors in order to carry out degrouping in an optimal manner). • Optimization of traffic assignment.

• Determination of congestion pricing zones.

• Organic control assistance tools.

• Generation of 4D trajectories.

• Prediction of congested zones.

• etc.

The operational capacity of a control sector is currently measured by the maximum number of aircraft able to traverse the sector in a given time period. This measurement does not take account of the orientation of traffic and considers geometrically structured and disordered traffic in the same manner. Thus, in certain situations, a controller may continue to accept traffic even if operational capacity has been exceeded (structured traffic); in other situations, controllers may be obliged to refuse additional traffic even though operational capacity has not yet been reached (disordered traffic). Thus, a measurement in terms of the number of aircraft per unit of time constitutes an insufficient metric for the representation of the difficulty level associated with a particular traffic situation.

In the context of operational control, the ideal would be to find a metric which precisely measures the level of mental effort needed to manage a set of aircraft. Without going quite so far, it is possible to find complexity metrics which go beyond a simple measurement of the number of aircraft. We shall begin by clarifying two essential notions for use in the rest of this chapter:

• Control workload: measurement of the difficulty for the traffic control system of treating a situation. This system may be a human operator or an automatic process. In the context of operational control, this workload is linked to the cognitive process of traffic situation management (conflict prediction and resolution, trajectory monitoring, etc.).

• Traffic complexity: intrinsic measurement of the complexity associated with a traffic situation. This measurement is independent of the system in charge of the traffic and is solely dependent on the geometry of trajectories. It is linked to sensitivity to initial conditions and to the inter-dependency of conflicts. Incertitude with respect to positions and speeds increases the difficulty of predicting future trajectories. In certain situations, this incertitude regarding future positions can increase exponentially, making the system extremely complex in that it is virtually impossible to reliably extrapolate a future situation. When a future conflict is detected, a resolution process is launched which, in certain situations, may generate new conflicts. This inter-dependency between conflicts is linked to the level of mixing between trajectories. As an example, figure 1 presents three traffic situations classed according to increasing level of difficulty as a function of the level of predictability and of inter-dependency between trajectories.

Low sensitivity

No conflict

Easy situation

High sensitivity

Conflicts with no interaction between solutions

Average situation

High sensitivity

Potential conflicts with interaction between solutions

Difficult situation

Fig. 1 Three traffic situations classed by increasing order of complexity One way of interpreting these notions is to imagine oneself to be in charge of each situation in a context where the radar imaging equipment has ceased to work. Naturally, our attention is immediately focused on the situation on the right, as it is difficult to predict (in terms of the appearance of conflicts) and presents a high level of inter-dependency between trajectories. The middle solution, which presents a significant risk of conflict, is easy to manage as the same direction order must simply be given to all of the aircraft (+90 or -90 degrees) in order to place them into safe roundabout trajectories. Finally, in the situation on the left, the trajectories do not present any difficulties and the relative distance between the aircraft will be maintained, at least for the immediate future.

Research into air traffic complexity metrics has attracted considerable attention in recent years, particularly in the United States and in Europe. The first projects were launched in Germany in the 1970s, and since then the subject has continued to develop. Currently, NASA, MIT and Georgia Tech are involved in work on the subject within the framework of the NextGen project. In Europe, the DSNA, the DLR and the NLR are involved in similar activities linked to SESAR.

The objective of most of this work is to model the control workload associated with given traffic situations. The main approaches are as follows:

• Workload model based on traffic level [START_REF] Schmidt | A queueing analysis on the air traffic controller's workload[END_REF]. This approach defines the workload as the proportion of control time over an hour, taking account of the average duration of routine control tasks for an aircraft, the average time taken to resolve conflicts per aircraft, the average rate of arrivals in a sector per hour and the average rate of conflicts in a sector per hour. • Queue-based model [START_REF] Maugis | Techniques de détermination de la capacité des secteurs de contrôle de l'espace aérien: statistiques et simulations[END_REF]. In this case, a control sector is modeled as a system receiving an input (airplanes) and providing a service, allowing the aircraft to traverse the sector in safety. The sector may then be modeled as a service center including one or more servers and an airplane queue. By applying queuing theory, this approach allows us to determine a maximum acceptable arrival rate for a sector. • Model based on airspace structure [START_REF] Janic | En route sector capacity model[END_REF]. In this case, the capacity of a sector is based solely on its structure (flight levels, routes, route intersections, etc.). • Dynamic density [START_REF] Laudeman | Dynamic density : an air traffic management metric[END_REF][START_REF] Sridhar | Airspace complexity and its application in air traffic management[END_REF][START_REF] Chatterji | Measure for air traffic controller workload prediction[END_REF]. This model, developed by NASA, consists of measuring a set of traffic characteristics (number of changes in direction, changes in speed, changes in altitude, etc.) and the workload experienced by a controller, then carrying out linear regression in order to adjust the model to the experienced workload as precisely as possible. This model of control workload presents the two following drawbacks:

1. incapacity for generalization to new sectors 2. modeling is highly dependent on the controllers used to infer the model.

Other approaches [START_REF] Granger | A traffic complexity approach through cluster analysis[END_REF][START_REF] Lee | Air traffic complexity : An input-output approach[END_REF] model the complexity of a traffic situation using automatic conflict resolution algorithms, for which we measure the number of trajectory modifications required in processing a given situation. In the same way as before, these methods are highly dependent on the type of algorithm used to resolve conflicts. These considerations have led us to develop intrinsic traffic complexity metrics which are only linked to trajectory structure, and not to the system used to process them.

The next section presents the new trajectory complexity metric we have developed based on a non linear dynamical system modeling of the air traffic. An application the the real operation Reims airspace is then presented in the following section.

2 Trajectory Complexity Metric Based on Non Linear Dynamical Systems

Non linear Dynamical System modeling of Air traffic

A non linear dynamical system is summarized by the following equation :

Ẋ(t) = f(X) (1)
where X is the state vector of the system (X = [x, y, z] T ) and f : C 2 vector field, describe systems which integral curves may fit the observed trajectories. This equation associates a vector speed Ẋ to a position in the space coordinate X and then synthesis a particular vector field. Based on the observations of the aircraft (positions, speed vectors), the dynamical system has to be adjusted with the minimum error.This fitting is done with a Least Square Minimization (LMS) method for which the following criterion is used :

E 1 = i=N ∑ i=1 V i -f(X i ) 2 (2)
where N is number of observations. If we consider criterion E 1 only, it can be shown that there is an infinite number of vector fields f which can be adjusted to the observations. In order to keep the smoothest one, another criterion is added which has to be minimized, the so-called "div-curl" criterion :

E 2 = R 3 α ∇divf(X) 2 + β ∇curlf(X) 2 dX (3)
with α, β positive weights controlling the smoothness of the approximation by focusing on constant divergence or constant curl. In the following, we will consider α = β = 0.5; in such case :

E 2 = ∆ f(X) 2
where ∆ f(X) is the Laplacian of the vector field f. The joint minimization of E 1 and E 2 induces a unique optimum vector field [START_REF] Amodei | A vector spline approximation[END_REF] :

f(X) = N ∑ i=1 φ ( X -X i ).a i + A.X + B
with a i parameter vectors (one for each observation),

φ ( X -X i ) = Q( X -X i 3 )
and

Q =   γ 0 0 0 γ 0 0 0 γ   with γ = ∂ 2 xx + ∂ 2 yy + ∂ 2 zz
The resulting adjustment is done without error (⇒ minE 1 = 0). When α = β = 0.5, the vector spline function φ has the following structure :

φ ( X -X i ) = 12. X -X i (r i = X i -X ).
It must be noticed, that farthest observations has more weight in such calculation. Then, in order to compute the smoothest vector field which fit exactly the measures, all observations have to be taken into account in the computation.

Lyapunov Exponents

The metric chosen for complexity computation relies on a measure of sensitivity to initial conditions of the underlying dynamical system called Lyapunov exponents.

In order to figure out what Lyapunov exponents are, let consider a point and look at its evolution when transported by the dynamical system. Let x 0 be fixed (initial point) and let γ be a point trajectory of the dynamical system associated to the vector field f given by :

γ(t, x 0 ) = x 0 + t 0 f(u, γ(u, x))du (4)
Assume now that trajectory is disturbed by a small perturbation ε, we have :

γ(t, x 0 + ε) = γ(t, x 0 ) + ∇ x f(γ(t, x)).ε + o( ε )
where ∇ x f(t, γ(t, x)) is the differential of the vector field f at x. Divergence to nominal trajectory with respect to time is thus γ(t, x 0 )γ(t, x) = D(t, s) (see figure 2). γ(t, x) being defined as a flow :

∂ γ(t, x) ∂t = f(t, γ(t, x)) γ(0, x) = x
with f a smooth vector field, it is possible to show that D(t) satisfies a differential equation also. Given a nominal trajectory γ(t, x 0 ), then divergence of nearby trajectories can be found up to order one in xx 0 by solving :

∂ D(t, x) ∂t = ∇ x f(t, γ(t, x)).D(t, x) D(0, x) = x -x 0
If the three space dimensions are considered (x, y, z), and since the previous equation is linear, it can be extended to the matrix form :

0 V γ γ 0 x x (t,x ) (t,x) γ D(t, )
Fig. 2 Time evolution of a reference trajectory and a perturbed trajectory

dM(t) dt = ∇ x f(t, γ(t, x)).M(t) M(0) = Id
Where each column of the M matrix corresponds to the divergence associated to the principal coordinate axis. This equation is called the variational equation of the flow. The variational equation describes in some sense a linear dynamical system "tangent" to the original one. Let U t (t)Σ (t)V (t) = M(t) be the singular value decomposition of M(t). The Lyapunov exponents are mean values of the logarithms of the diagonal elements of Σ (t):

κ(s) = - 1 T T 0 log(Σ ii (t))dt ∀Σ ii (t) ≤ 1
When Lyapunov exponents are high, the trajectory of a point under the action of the dynamical system is very sensitive to initial conditions (or to the parameters on which the vector field may depend), so that, situation in the future is unpredictable. On the other hand, small values of the Lyapunov exponents mean that the future is highly predictable (expected to be comfortable for a controller).

So, the Lyapunov exponent map determines the area where the underlying dynamical system is organized. It identifies the places where the relative distances between aircraft do not change with time (low real value) and the ones where such distance change a lot (high real value).

Let us now describe the practical algorithm for computing complexity maps.

1. Dynamical regression based on radar observations. 2. The vector field and its associated gradient is then computed on each point of a cube of airspace. 3. Lyapunov exponents are computed at each point of the cube by the mean of a Runge-Kutta integration.

The figure 3 shows an example of Lyapunov exponents map for which full organized miles in trail trajectories (from south west to north east) cross two random traffic situations. This figure shows clearly a complexity valley on the miles in trail Fig. 3 Miles in trail traffic between disordered areas direction. This organization may have been detected even if the miles in trail trajectories would have been structured on a curve trajectory. That is the strong point of this metric: Lyapunov exponents are able to identify any kind of trajectory organization (when aircraft follow the same trajectory at the same speed).

This approach is able to produce a picture of the complexity at a given time. By summing the Lyapunov exponents in a given airspace on can extract a scalar metric which has been used in the Reims experiments presented in the next section.

Workload Prediction based on complexity metric

Based on the previous complexity algorithm, we propose to compare our metric to the workload measured on a real control sector. To do such experiments, we have used some real traffic sample of control sector belonging to the Reims airspace area (see figure 4). Then, we have considered the traffic crossing this airspace for the 14 February 2013 from 10:00 till 12:00AM. The controller in charge of such traffic have been invited to quantify the workload they feel by using a dedicated HMI. Such workload measurements have been done every five minutes. These measurements have been recorded and aggregated in vector a {X i } i=1..N . For the same traffic samples we have measured the associated traffic complexity at the same step and a vector of complexity measures has also been built : {Y i } i=1..N .

The objective of the study is to establish the level of correlation between those two time series. The first guess was to used the standard correlation measure :

C XY = N ∑ i=1 (X i -X)(Y i -Ȳ ) (5)
It is maximum when (X -X) = λ (Y -Ȳ ). However, since workload values given by the controllers are only relative, such a correlation may gives false results, we have decided to used a rank correlation approach which indicates if both series evolve in the same direction.

First, we consider any pair of samples in both series and we introduce the following variable :ε i j = (X i -X j )(Y i -Y j ) and four situations can be identified :

ε i j (X i -X j ) (Y i -Y j ) + + + + + + - + - - - +
In the first two cases the method classify the times series samples in the same way.

Then, we introduce the variable σ = 1 N(N-1) ∑ i> j ε i j which indicate how often both series vary the same direction. If σ is close to one or minus ones both series are classified in the same way but if σ is close to zero no conclusion can be given.

The confidence level was checked using a Spearman test [START_REF] Fieller | ests for rank correlation coefficients[END_REF].

H 0 : σ = 0 Pr { reject H 0 /H 0 true}
To quantify such hypothesis the P value is computed. If such P value is close to zero, then the correlation is not linked to hazard.

Based on the traffic samples the overall rank correlation value reaches 0.78 which is very good number based on the data we get from operation.

The associated P Value = 10 -6 meaning that such correlation meaningful. The evolution of both series is given on figure 5. Those measures are then aggregated on the blue curve (with a time window filtering in order to get a smooth curve).The green curve represents the complexity measures by our algorithm for the same time samples.

Conclusion

We have presented in this paper a new air traffic complexity metric based on vector field model of air traffic. Most of the work has been devoted to improvements in interpolating splines used to fit a vector field to the observation. Unlike linear models which produce mean complexity indicators, the non-linear one may give local information, thus providing a way of displaying maps of complexity. Based on this complexity metric, a new control workload predictor has been developed and validated on a real operational airspace.

Fig. 4

 4 Fig. 4 Airspace used for our experiments

Fig. 5

 5 Fig.5Times series evolution. Those curves represent the evolution of the control workload and the associated complexity with time. To built those curves we consider the Reims airspace and the associated real traffic. Every two minutes we request the controller to rank the workload he feels. Those measures are then aggregated on the blue curve (with a time window filtering in order to get a smooth curve).The green curve represents the complexity measures by our algorithm for the same time samples.