Sonic boom propagation in the Planetary Boundary Layer 1.A brief overview of sonic boom

Sonic boom is a complex phenomenon [START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF][START_REF] Maglieri | Sonic Boom[END_REF][START_REF] Plotkin | State of the art of sonic boom modeling[END_REF][START_REF] Coulouvrat | The challenges of defining an acceptable sonic boom overland[END_REF][START_REF] Maglieri | Sonic Boom Six Decades of Research[END_REF] that couples compressible fluid mechanics near the aircraft, acoustic propagation in the far-field and human perception. It occurs when a body is moving at a speed above the local sound speed. Since the body is moving faster than the perturbation it creates, conic shocks (Mach cones) which are attached to the geometry are created [START_REF] Whitham | The flow pattern of a supersonic projectile[END_REF]. These shocks then propagate in the atmosphere from the aircraft down to the ground. During this propagation, the waveform evolves from near-field to far-field as shown on Fig. 1.1: due to cumulative nonlinear effects, higher pressure parts of the signal travel at a slightly higher speed than the ambient sound speed and lower pressure parts travel slightly Figure 1.1: Schematics of waveform distortion. Extracted from [START_REF] Plotkin | State of the art of sonic boom modeling[END_REF].

slower. The shocks tend to coalesce leading to the classical N-wave shape on the ground.

In the near field, a linearised supersonic flow theory [START_REF] Whitham | The flow pattern of a supersonic projectile[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF] was developed by G.B. Whitham to compute the generation of the perturbation. Whitham's theory has been the foundation for sonic boom theory since the late 50s. It has been extended to take into account lift effects for non axisymmetric bodies [START_REF] Walkden | The Shock Pattern of a Wing-Body Combination, Far From the Flight Path[END_REF]. It remains today at the basis of matching between CFD and acoustics [START_REF] George | Reduction of sonic boom by azimuthal redistribution of overpressure[END_REF], and sonic boom minimization [START_REF] Seebass | Minimum sonic boom shock strengths and overpressures[END_REF][START_REF] George | Propagation of Sonic Booms and other Weak Nonlinear Waves through Turbulence[END_REF]. However, this method cannot handle the case of blunt or hypersonic bodies for which the dynamic of the flow is highly nonlinear, though it gives accurate far-field results otherwise. Nowadays, CFD allows computations of the flow around complex geometries including low-boom design [START_REF] Pawlowski | Origins and overview of the shaped sonic boom demonstration program[END_REF] or hypersonic speeds [START_REF] Loubeau | Effects of Meteorological Variability on Sonic Boom Propagation from Hypersonic Aircraft[END_REF][START_REF] Henneton | Numerical Simulation of Sonic Boom from Hypersonic Meteoroids[END_REF][START_REF] Luquet | Sonic Boom Assessment of a Hypersonic Transport Vehicle with Advanced Numerical Methods[END_REF]. In the far-field, where the perturbations are small and regular enough to use acoustic theory, the propagation is well described by the nonlinear geometrical acoustic theory [START_REF] Whitham | On the propagation of weak shock waves[END_REF][START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Plotkin | State of the art of sonic boom modeling[END_REF]. The coupling between mid-and far-fields can be done directly or by using the multipole matching [START_REF] Page | An efficient method for incorporating computational fluid dynamics into sonic boom prediction[END_REF] which extracts from CFD the part of the pressure field satisfying assumptions of Whitham theory. Recently, anisotropic mesh adaptation has been used to propagate the shocks far enough from the aircraft to reach sufficiently distant zones where the assumptions of geometric acoustics are satisfied [START_REF] Loseille | Anisotropic Mesh Generation Application to High-Fidelity Simulation in CFD[END_REF].

During the propagation, the atmosphere plays a major role on both the localization of the sonic boom footprint and on the shape and amplitude of ground waveforms. Its influence can be divided in two: (i) the macro effects due vertical gradients of atmospheric pressure, temperature and wind; (ii) and the micro effects due to the small scale turbulence. The macro effects result in refraction by bending the ray paths. It defines intensity, location and lateral extent of the ground boom. As shown on Fig. 1.2, upward refraction due to negative temperature gradients leads to a finite width carpet under the flight path surrounded by a shadow zone were no boom is heard. After this Figure 1.2: Schematics of sonic boom exposure. Extracted from [START_REF] Maglieri | Sonic Boom Six Decades of Research[END_REF].

shadow zone, a secondary boom exists resulting from rays that have been reflected by the upper atmosphere [START_REF] Le Pichon | Acoustic propagation and atmosphere characteristics derived from infrasonic waves generated by the Concorde[END_REF][START_REF] Gainville | Modélisation de la propagation atmosphérique des ondes infrasonores par une méthode de tracé de rayons non linéaire[END_REF]. Large scale variability due to meteorological variations of the atmosphere stratification has been quantified for various climates [START_REF] Blumrich | Variability of focused sonic booms from accelerating supersonic aircraft in consideration of meteorological effects[END_REF]. Other macro effects are the absorption and dispersion due to thermoviscosity and molecular relaxation of nitrogen and oxygen [START_REF] Sutherland | Atmospheric absorption in the atmosphere up to 160 km[END_REF] which modify the shock structure [START_REF] Cleveland | Propagation of sonic booms through a real, stratified atmosphere[END_REF][START_REF] Cleveland | Comparison of computer codes for the propagation of sonic boom waveforms through isothermal atmospheres[END_REF][START_REF] Cleveland | Time-domain modeling of finite-amplitude sound in relaxing fluids[END_REF]. As will be detailed later on, micro effects due to atmospheric turbulence can significantly alter the waveform especially in terms of peak overpressure and rise time. They have also an effect on the lateral boom by scattering noise in the shadow zone.

The classical approach of geometrical acoustics has several limitations in the sonic boom case because: (i) it cannot take into account diffraction and (ii) it is a high frequency approximation. Sonic boom focusing occurs for an accelerating or manoeuvring aircraft [START_REF] Wanner | Theoretical and Experimental Studies of the Focus of Sonic Booms[END_REF]. The ray tube area vanishes and ray tracing predicts infinite amplitude. Locally, around the caustic, diffraction has to be taken into account in addition to nonlinearities to remove the singularity [START_REF] Guiraud | Acoustique géométrique, bruit balistique des avions supersoniques, focalisations[END_REF]. A numerical solver for the resulting nonlinear Tricomi equation has recently been developed [START_REF] Auger | Numerical Simulation of Sonic Boom Focusing[END_REF][START_REF] Marchiano | Numerical simulation of shock wave focusing at fold caustics, with application to sonic boom[END_REF][START_REF] Salamone | Solution of the Lossy Nonlinear Tricomi Equation Applied to Sonic Boom Focusing[END_REF]. Another limitation of ray theory happens for shadow zones. Geometrical acoustic predicts no ray there and thus no signal at all. To alleviate this problem, creeping wave theory has been developed [START_REF] Onyeowu | Diffraction of sonic boom past the nominal edge of the corridor[END_REF][START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF][START_REF] Berry | Controlled experiments of the diffraction of sound by a curved[END_REF]. It has been applied to sonic boom [START_REF] Coulouvrat | Sonic boom in the shadow zone: A geometrical theory of diffraction[END_REF] along with proper matching between ray tracing and creeping waves. Finally, in the Planetary Boundary Layer (PBL), the wavelength of sonic boom is comparable to the size of the turbulent structure and the high frequency hypothesis is not valid [START_REF] Maglieri | Some effects of airplane operations and the atmosphere on sonic-boom signatures[END_REF]. This is the subject of the present study and a detailed review of this problem is now detailed in the following Sec. 1.1.2 Left: low-wind conditions (U 0 ≈ 5 m.s -1 ); right: high-wind conditions (U 0 ≈ 14 m.s -1 ). Extracted from [START_REF] Hilton | Sonic-boom measurements during bomber training operations in the chicago area[END_REF].

Influence of the turbulence on sonic boom

As stated in the previous section, atmospheric turbulence causes large variations in boom signatures altering shock overpressure and rise time. This high variability of the waveforms is illustrated by Fig. 1.3 which plots booms recorded by five microphones on two days, one with low-wind conditions (stable low atmosphere), and the other one with high-wind conditions (unstable low atmosphere).

First, this waveform variability has been investigated using flight tests [START_REF] Hilton | Sonic-boom measurements during bomber training operations in the chicago area[END_REF][START_REF] Garrick | Atmospheric effects on the sonic boom[END_REF][START_REF] Maglieri | Variability in sonic-boom signatures measured along an 8000-foot linear array[END_REF][START_REF] Lee | Sonic Booms Produced by United States Air Force and United States Navy Aircraft: Measured Data[END_REF][START_REF] Willshire | Preliminary results from the White Sands Missile Range sonic boom propagation experiment[END_REF]. They have shown that perturbations are random as expected from propagation through turbulence. Most perturbations are located after the shock. The peak shock overpressure can be doubled in some cases, but perturbations decrease rapidly after the shock, thus leading to spiked waveforms. Such signals are visible on the high wind case of Fig. 1.3. In other cases, on the contrary, rounded waveforms with lessened amplitude are recorded. Also, the leading and trailing shocks are affected similarly indicating that the turbulence can be considered as frozen. Recently, D-SEND (Drop test for Simplified Evaluation of Non-symmetrically Distributed sonic boom) flight tests [START_REF] Jaxa | [END_REF] have been realized, measuring the boom from a dropped out projectile. In particular, assuming thermal turbulence only, Kanamori et al. [START_REF] Kanamori | Development of nonlinear acoustic propagation analysis tool toward realization of loud noise environment prediction in aeronautics[END_REF] compared measured statistical distributions, to a numerical model combining ray tracing and diffraction approach based on the 2D HOWARD method [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF] at the basis of the present work. An example of simulated and measured pressure waveforms is shown on

Models for sonic boom propagation in the PBL

At the end of the 60s this subject was tackled analytically by many researchers. First, Pierce [START_REF] Pierce | Spikes on Sonic-Boom Pressure Waveforms[END_REF] proposed a model based on geometrical acoustics taking into account diffraction for the low frequencies. This model was able to explain the apparition of spiked and rounded waveform but the finer details of the perturbation could not be quantified. Crow [START_REF] Crow | Distortion of sonic bangs by atmospheric turbulence[END_REF] developed a shock scattering model. It is a first order scattering model inspired by linear random scattering theory [START_REF] Chernov | Wave Propagation in a Random Medium[END_REF][START_REF] Tatarski | The Effects of the Turbulent Atmosphere on Wave Propagation[END_REF][START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF] with each shock modelled as a step shock. It predicts that for a Kolmogorov spectrum of turbulence, the envelope of the perturbation varies as:

ψ 2 (t) = t c t 7/6 (1.1)
with ψ the ratio between the scattered pressure and the peak overpressure, t is the time and t c is a characteristic time in the range of 0.5 -1.5 m.s -1 for usual atmospheric conditions. The power -7/6 is directly related to the inertial range of turbulence modelled by Kolmogorov spectrum. Kamali and Pierce [START_REF] Kamali | Time Dependence of Variances of Sonic Boom Waveform[END_REF] were able to verify this -7/6 scaling law against a large set of sonic boom signatures. One major problem with Crow's theory is that it does not give bounds for the shock overpressure and does not explain the anomalously long rise time observed during flight tests [START_REF] Lilley | The structure of shock waves at large distances from bodies traveling at high speeds[END_REF]. To handle this, Plotkin and George [START_REF] Plotkin | Propagation of weak shock waves through turbulence[END_REF] extended the theory to finite time shock waves. This lead them to a Burgers' equation where viscosity is linked to the medium turbulence. This allowed them to compute an upper bound for rise time. Pierce raised some concerns about this model and proposed [START_REF] Pierce | Statistical Theory of Atmospheric Turbulence Effects on Sonic-Boom Rise Times[END_REF] a statistical model based on multiple wave front folding. Ffowcs Williams and Howe also tried to tackle the problem [START_REF] Williams | On the possibility of turbulent thickening of weak shock waves[END_REF].

Using multiple scattering theory, they derived a diffusion equation for the distribution of energy in wavenumber space. However, this model fails to take into account spatial relocation of acoustic energy caused by phase scrambling. Several reviews of these theories can be found [START_REF] Pierce | Effects of Atmospheric Irregularities on Sonic-Boom Propagation[END_REF][START_REF] Plotkin | Review of sonic boom theory[END_REF][START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part III: Validation of sonic boom propagation models[END_REF].

Laboratory scale and numerical experiments

Validating these models is difficult since the only available data emanate from flight tests for which the atmospheric turbulence parameters cannot be controlled and measurement are scant. To overcome this difficulty, numerous laboratory scale experiments have been carried out. Bauer and Bagley [START_REF] Bauer | Sonic boom modeling investigation of topographical and atmospheric effects[END_REF] performed ballistic experiments with supersonic bullets fired in an open chamber. A turbulent flow was used to model the atmospheric turbulence. Large distortions of the waveforms and increases in rise time were observed. Davy and Blackstock [START_REF] Davy | and Diffraction Wave by a Gas-Filled Soap Bubble[END_REF], motivated by Pierce's theory [START_REF] Pierce | Spikes on Sonic-Boom Pressure Waveforms[END_REF] based on focusing and defocusing effect of turbulence, studied short N-waves produced by an electrical spark propagating through a soap bubble filled with different gases. Argonfilled bubbles acted as converging lenses leading to spiked waveforms, while heliumfilled bubbles (diverging lenses) caused rounded waveforms. Ganjehi et al. [START_REF] Ganjehi | Ondes de choc acoustiques en milieu hétérogène, des ultrasons au bang sonique[END_REF][START_REF] Ganjehi | Evidence of wave front folding of sonic booms by a laboratory-scale deterministic experiment of shock waves in a heterogeneous medium[END_REF] reproduced nonlinear wavefront folding due to deterministic heterogeneities in a water tank with a 1:100,000 ratio with sonic boom. Ribner [START_REF] Ribner | Laboratory simulation of development of superbooms by atmospheric turbulence[END_REF] investigated spiking and rounding of N-waves within a jet. He found that waves (generated in a shock tube) travelling upstream were focused while waves travelling downstream were defocused.

Tubb [START_REF] Tubb | Measured effects of turbulence on the rise time of a weak shock[END_REF] also used a shock tube to measure the effects of grid-generated turbulence on the rise time of a weak shock. Statistically, passage through turbulence causes a doubling of the rise time.

Still relying on a shock tube, Hesselink and Sturtevant [START_REF] Hesselink | Propagation of weak shocks through a random medium[END_REF] propagated a plane step shock through a random medium composed of helium and refrigerant 12. They have shown that waveform distortion is mainly due to focusing and defocusing by random heterogeneities. They supported this result with a ray tracing approach explaining most of their results.

Raspet et al. [START_REF] Raspet | Statistical and numerical study of the relationship between turbulence and sonic boom characteristics[END_REF] analysed flight test results and developed a scattering centerbased model (based on Born approximation) to propagate sonic boom through single scale gaussian turbules. They reproduced the spiked and rounded waveforms along with the increase in rise time. Boulanger et al. [START_REF] Boulanger | Sonic boom propagation through a realistic turbulent atmosphere[END_REF] extended this model to a more realistic atmosphere based on a von Kármán spectrum. This improved the comparison to measurement but the rise time was still underestimated. Finally, Kelly et al. [START_REF] Kelly | Scattering of sonic booms by anisotropic turbulence in the atmosphere[END_REF] applied this model to an anisotropic turbulence. This did not improve the comparison with the measurement, indicating that this single-scatter method cannot completely describe sonic boom propagation through turbulence.

Lipkens and Blackstock [START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part I. General Reults[END_REF][START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part II. Effect of turbulence intensity and propagation distance through turbulence[END_REF] scaled their experiments to sonic boom for both the N-waves and characteristics of the atmosphere. N-waves produced by sparks were propagated through kinematic turbulence within a jet. Plane and spherical waves were both studied. Waveform distortion of actual sonic boom was successfully reproduced. Also the average peak pressure was shown to slightly decrease through turbulence while the rise time largely increases. This dataset was also analysed by Lipkens [START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part III: Validation of sonic boom propagation models[END_REF] to validate the models of Crow [START_REF] Crow | Distortion of sonic bangs by atmospheric turbulence[END_REF], Plotkin and George [START_REF] Plotkin | Propagation of weak shock waves through turbulence[END_REF], Pierce [START_REF] Pierce | Statistical Theory of Atmospheric Turbulence Effects on Sonic-Boom Rise Times[END_REF] and Ffowcs-Williams and Howe [START_REF] Williams | On the possibility of turbulent thickening of weak shock waves[END_REF]. It was found that Crow's theory yielded a maximum for the Root Mean Square (RMS) perturbations with an overestimation of a factor 5 compared to the experiments. Its extension by Plotkin and George largely overpredicts rise time but no clear conclusion could be drawn as the model experiment was outside the model range of validity. Pierce's model of multiple wave front folding was successful in reproducing the rise times of the experiments with an error never exceeding 20%. Ffowcs-Williams and Howes multiple scattering model underestimated the rise time, which supports the fact that a phase scrambling mechanism is responsible for the large rise times.

Blanc-Benon et al. [START_REF] Blanc-Benon | Propagation of finite amplitude sound through turbulence: Modeling with geometrical acoustics and the parabolic approximation[END_REF] relied on numerical methods to investigate shock wave propagation through random turbulence. The turbulence was modelled as multiple realizations of a random temperature and velocity field. First, geometrical acoustic was used along with a nonlinear transport equation to compute the waveform. Then a one-way model was worked out under the form of a KZK equation modified to take into account the effect of a random temperature field. Both models confirmed the observations of Lipkens et al. [START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part I. General Reults[END_REF][START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part II. Effect of turbulence intensity and propagation distance through turbulence[END_REF]: average peak pressure decreases while rise time increases. Also simulations lead to waveform distortions similar to experimental ones. The main mechanism for waveform distortion and variability in peak pressure and rise time is shown to be the formation of caustics (wavefront focusing). Continuing this work, Averiyanov et al. [START_REF] Averiyanov | Parabolic equation for nonlinear acoustic wave propagation in inhomogeneous moving media[END_REF] developed a KZK equation extended to a slightly inhomogeneous and slowly moving medium. This equation was then solved numerically at two dimensions to simulate propagation of an N-wave through kinematic turbulence by Averiyanov et al. [START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF]. The same features as Blanc-Benon et al. [START_REF] Blanc-Benon | Propagation of finite amplitude sound through turbulence: Modeling with geometrical acoustics and the parabolic approximation[END_REF] were retrieved. The importance of taking into account diffraction and nonlinearity was shown. The influence of the turbulence model was also assessed. The mean peak pressure depends strongly on the the large scale heterogeneities while rise time depends on both large and small scale heterogeneities. The presence of a turbulent transverse flow is negligible but it becomes important when increasing the size of heterogeneities, such as for a uniform transverse flow. In parallel, Averiyanov et al. [START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF] performed a series of laboratory scale experiments using high amplitude and short duration spherical N-waves generated by spark pulses. The turbulent medium was produced by a jet and was fully characterized. This work confirmed the influence of the large scale heterogeneities on the average peak pressure. Also similarly to linear acoustic [START_REF] Blanc-Benon | Intensity fluctuations of spherical acoustic waves propagating through thermal turbulence[END_REF], the probability density of the peak pressure is shown to follow a generalised gamma distribution which smoothly varies from a log-normal to an exponential distribution as the turbulence level is increased. Salze et al. [START_REF] Salze | Laboratory-scale experiment to study nonlinear N-wave distortion by thermal turbulence[END_REF] performed analogous experiments for thermal turbulence. Results corroborate the previous experiments and the generalised gamma distribution was again retrieved for the probability density of the peak overpressure.

Another experiment of shock waves propagation through kinematic turbulence was carried by Sasoh et al. [START_REF] Sasoh | Statistical behavior of post-shock overpressure past grid turbulence[END_REF]. The decrease of average peak pressure was recovered along with a large standard deviation almost proportional with the RMS of the turbulent velocity.

All these studies provide a good understanding of the phenomenon. However, one of the difficulties is that the laboratory scale experiments cannot scale perfectly all the sonic boom characteristics. For instance, as can be seen on Table . 1.1 for the work of Averiyanov et al. [START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF], the ratio between wavelength and integral length scale is not conserved. For Lipkens and Blackstock set-up [START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part I. General Reults[END_REF], nonlinear effects are about ten times larger than sonic boom ones. In all aerial small scale experiments, thermoviscous effects are strongly enhanced. This can be compensated in water [START_REF] Ganjehi | Evidence of wave front folding of sonic booms by a laboratory-scale deterministic experiment of shock waves in a heterogeneous medium[END_REF] but in this case, only deterministic heterogeneities can be investigated.

Numerical methods in nonlinear acoustics

Several numerical methods are proposed in the literature for simulating propagation of acoustic shock waves in both homogeneous, heterogeneous or moving media. Most of them were developed either for atmospheric propagation (sonic boom, infrasound) or more frequently for biomedical applications (HIFU, lithotripsy).

Geometrical acoustics

One of the first methods is the ray tracing [START_REF] Blokhintzev | The Propagation of Sound in an Inhomogeneous and Moving Medium I[END_REF][START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Candel | Numerical solution of conservation equations arising in linear wave theory: application to aeroacoustics[END_REF][START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF]. It solves the equations of linear geometrical acoustics to obtain the ray path. In the nonlinear regime, a nonlinear transport equation similar to Burgers' equation is solved along each ray to get the waveform evolution. This has been developed first purposely for sonic boom [START_REF] Whitham | On the propagation of weak shock waves[END_REF][START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Thomas | Extrapolation of sonic boom pressure signatures by the waveform parameter method[END_REF][START_REF] Taylor | The TRAPS sonic boom program[END_REF]. Nowadays, infrasound investigations also routinely rely on it [START_REF]Infrasound Monitoring for Atmospheric Studies[END_REF]. However this method has important limitations as (i) it is a high frequency approximation and (ii) it does not take into account diffraction effects. This makes it impossible to predict phenomena such as shadow zone, focusing or scattering by turbulence. Diffraction effects therefore need to be introduced by solving a nonlinear wave equation using various methods which are now presented.

One-way approaches

One important class of methods is the one way approximation. It assumes that there is a main direction of propagation and that the backscattered field is negligible. As a consequence it can be applied only for weakly heterogeneous media. The original nonlinear one way wave equation is the so called Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation [START_REF] Zabolotskaya | Quasi-plane waves in the nonlinear acoustics of confined beams[END_REF][START_REF] Kuznetsov | Equations of nonlinear acoustics[END_REF]. It is a nonlinear version of the standard parabolic equation valid only for small angles of propagation (±15 • ). For a review about linear parabolic approximation see [START_REF] Lee | Parabolic Equation Development in the Twentieth Century[END_REF]. The first numerical algorithms were developed for axisymmetric narrow beams in former USSR during the seventies, alternating between frequency domain for absorption and diffraction and time domain for nonlinear effects [START_REF] Bakhvalov | Nonlinear Theory of sounds Beams[END_REF]. On the contrary, Aanonsen et al. [START_REF] Aanonsen | Distortion and harmonic generation in the nearfield of a finite amplitude sound beam[END_REF] worked entirely in the frequency domain. Coordinate transformations [START_REF] Hamilton | Nonlinear effects in the farfield of a directive sound source[END_REF][START_REF] Hart | Nonlinear effects in focused sound beams[END_REF] allow to follow more easily beam divergence or focusing. First extensions to 3D were proposed by Kamakura et al. [START_REF] Kamakura | Harmonic generation in finite amplitude sound beams from a rectangular aperture source[END_REF]. A fully time domain algorithm was achieved by Lee and Hamilton [START_REF] Lee | Time-domain modeling of pulsed finite-amplitude sound beams[END_REF]. Nonlinear effects are solved by the exact Poisson's solution. As a consequence, the advancement step must be small enough so as to allow absorption effects to prevent shock formation. This can be very penalising in case of a weakly absorbing medium or a low-frequency case such as sonic boom or infrasound. Nevertheless this algorithm was applied to simulate a high frequency, laboratory scale experiment [START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part I. General Reults[END_REF] of sonic boom in a turbulent temperature field [START_REF] Blanc-Benon | Propagation of finite amplitude sound through turbulence: Modeling with geometrical acoustics and the parabolic approximation[END_REF]. Three dimensional extension of that algorithm is proposed in 2007 by Jing and Cleveland [START_REF] Jing | Modeling the propagation of nonlinear threedimensional acoustic beams in inhomogeneous media[END_REF]. Many HIFU applications rely on the KZK equation, an extensive literature review being out of the scope of the present work. An alternative of the KZK equation is the Nonlinear Progressive-wave Equation (NPE) [START_REF] Mcdonald | Time-domain solution of the parabolic equation including nonlinearity[END_REF] which still relies on a narrow angle parabolic approximation but with advancement in time instead of space. One application to sonic boom wavefront folding has been performed by Piacsek [START_REF] Piacsek | Atmospheric turbulence conditions leading to focused and folded sonic boom wave fronts[END_REF].

One of the main drawbacks of these methods is the small angle limitation. It can be removed using for the linear part the angular spectrum method. The first approach in this way was the work of Christopher and Parker [START_REF] Christopher | New approaches to nonlinear diffractive field propagation[END_REF] which a Discrete Hankel Transform (DHT) implemented to propagate axisymmetric beams with nonlinear effects handled in the frequency domain. DHT was replaced by Rayleigh integral [START_REF] Tavakkoli | Modeling of pulsed finite-amplitude focused sound beams in time domain[END_REF] and then by the more efficient Fast Fourier Transform (FFT) at 3D by Zemp et al. [START_REF] Zemp | Modeling of nonlinear ultrasound propagation in tissue from array transducers[END_REF]. Both of them introduced a second order split-step between diffraction and nonlinearities. This approach has been extended in 2D to the heterogeneous case by Varslot and Taraldsen [START_REF] Varslot | Computer simulation of forward wave propagation in soft tissue[END_REF] and more recently through the so-called HOWARD (Heterogeneous One-Way Approximation for the Resolution of Diffraction) method [START_REF] Dagrau | Simulation de la propagation du bang sonique : de la CFD à l'acoustique non linéaire[END_REF][START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF]. Contrarily to previous works, these two studies solve the nonlinear effects in the time domain by means of Poisson's solution. Varslot and Taraldsen solve the heterogeneities in sound speed along with nonlinearities while Dagrau et al. handle them separately. In this last work nonlinear effects are solved for the potential through the Burgers-Hayes method [START_REF] Burgers | Further statistical problems connected with the solution of a simple nonlinear partial differential equation[END_REF][START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Coulouvrat | The challenges of defining an acceptable sonic boom overland[END_REF] which enables to manage inviscid shock waves which are important for atmospheric applications. Yuldashev et al. [START_REF] Yuldashev | Simulation of three-dimensional nonlinear fields of ultrasound therapeutic arrays[END_REF] use the same angular spectrum approach for 3D homogeneous simulations. The heterogeneities are introduced only by phase screens [START_REF] Yuldashev | Propagation d'ondes de choc dans les ilieux aléatoires avec des inhomogénéités dans l'espace ou dans une couche mince[END_REF].

Full wave approaches

For cases where the assumption of a main direction of propagation is not valid such as for strongly heterogeneous media, a full nonlinear wave equation has to be solved. This approach has first been tackled by Sparrow and Raspet [START_REF] Sparrow | A numerical method for general finite amplitude wave propagation in two dimensions and its application to spark pulses[END_REF] at 2D for investigating Mach reflection based on a Finite Differences Time Domain (FDTD) method. Since, several authors have adapted the FDTD approach to various problems. These algorithms are listed in Table . 1.2 along with their main characteristics: problem dimensionality (2D or 3D); scheme order; applications involving cases either weakly (for instance harmonic generation) or strongly (shock waves) nonlinear.

Marsden et al. [START_REF] Marsden | A study of infrasound propagation based on high-order finite difference solutions of the Navier-Stokes equations[END_REF] solve the 2D full Euler equations in a stratified atmosphere with no flow using high order numerical methods developed for aeroacoustics with a shock sensor. A Finite Volume Method (FVM) [START_REF] Velasco-Segura | A finite volume approach for the simulation of nonlinear dissipative acoustic wave propagation[END_REF] has been proposed to easily handle complex geometries. Recently, a Discontinuous Galerkin (DG) method [START_REF] Tripathi | Numerical Method for Propagation of Acoustical Shock Waves in Complex Geometry Using Discontinuous Galerkin Method[END_REF] has been developed so as to handle shock waves while remaining high order thanks to artificial viscosity and shock sensor. DG method also allows simulations in complex

Reference

Dimension Weakly nonlinear Order or Shock Sparrow and Raspet [START_REF] Sparrow | A numerical method for general finite amplitude wave propagation in two dimensions and its application to spark pulses[END_REF] 2D Shock 4 th Hallaj and Cleveland [START_REF] Hallaj | FDTD simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound[END_REF] 2D Weakly 4 th Ginter et al. [START_REF] Ginter | Full-wave modeling of therapeutic ultrasound: nonlinear ultrasound propagation in ideal fluids[END_REF] 2D Shock 4 th DRP Pinton et al. [START_REF] Pinton | A heterogeneous nonlinear attenuating full-wave model of ultrasound[END_REF] 3D Weakly 4 th Del Pino et al. [START_REF] Del Pino | 3D Finite Volume simulation of acoustic waves in the earth atmosphere[END_REF] 2D Shock up to 17 th Okita et al. [START_REF] Okita | The role of numerical simulation for the development of an advanced HIFU system[END_REF] 3D Weakly 6 th Table 1.2: Main characteristics of FDTD schemes for nonlinear wave equation.

geometries. The k-wave pseudo-spectral method [START_REF] Treeby | Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method[END_REF][START_REF] Treeby | Modeling nonlinear wave propagation on nonuniform grids using a mapped k-space pseudospectral method[END_REF][START_REF] Jaros | Use of multiple GPUs on shared memory multiprocessors for ultrasound propagation simulations[END_REF][START_REF] Jaros | Full-wave nonlinear ultrasound simulation on distributed clusters with applications in high-intensity focused ultrasound[END_REF] is an alternative way to relax the constrains on the spatial mesh by computing spatial gradients with a Fourier collocation method. A promising method is the Fourier continuation [START_REF] Albin | Fourier continuation methods for high-fidelity simulation of nonlinear acoustic beams[END_REF] which seems to successfully get high order schemes and handling the presence of shock wave while allowing for complex geometries. These methods are attractive as physics is well represented. However, it is difficult to combine simultaneously high order schemes necessary for the acoustic propagation and the presence of shock waves.

Moving medium

Note that in all the above cited works (except for the ray-tracing method), the considered medium can be heterogeneous but is always at rest. There have been relatively few investigations on weak shock waves propagation in moving media. With this in view, Averiyanov et al. generalized the KZK equation [START_REF] Averiyanov | Parabolic equation for nonlinear acoustic wave propagation in inhomogeneous moving media[END_REF] to a slowly moving medium. A numerical scheme extends the 2D time domain algorithm of Lee and Hamilton [START_REF] Lee | Time-domain modeling of pulsed finite-amplitude sound beams[END_REF] with axial flow convection and absorption handled in the frequency domain [START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF]. Also the Poisson's solution for the resolution of inviscid Burgers' equation is replaced by a spatially centred Godunov scheme. Very recently, Sabatini et al. [START_REF] Sabatini | Numerical simulation of infrasound propagation in the Earth's atmosphere: Study of a stratospherical arrival pair[END_REF] extended Marsden et al. work to incorporate a 2D stratified flow for infrasound application. Finally, to our knowledge, the only other reference about shock waves in moving media is the so-called FLHOWARD method [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF][START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF] that takes into account influence of a mean shear flow at 2D. This additional effect is handled similarly to heterogeneities in the HOWARD method.

Objective

The previous analysis shows there is currently no numerical method able to fully simulate shock waves (including sonic booms) propagation in a realistic three dimensional turbulent atmosphere. Most of the existing algorithms are two dimensional. Only very few include wind effects and none at 3D. Concerning sonic boom applications, existing software have been applied mostly to compare with laboratory scale experiments. Therefore, they were carried out with different parameters scaling: more intense sources of higher frequency so more nonlinear problems with more absorption; and a smaller ratio of acoustic wavelength to integral length scale of turbulence.

The first main objective of this thesis is to develop such a method. Given the high numerical cost of simulating 3D shock waves propagation, a one-way approach has been selected as a trade-off between accuracy and efficiency. Developing a 3D propagation code requires an efficient parallelization methodology. The so-called 2D FLHOWARD method with stratified flow has been extended and implemented with the following new features:

• fully three-dimensional;

• parallel programming using the Message Passing Interface (MPI) library;

• 3D flow effects;

• atmospheric absorption (thermoviscosity and molecular relaxation);

• Dirichlet or von Neuman boundary conditions instead of periodicity, in particular for rigid ground reflection.

The resulting software has been named FLHOWARD3D.

The second main objective is to apply this newly developed software to quantify influence of turbulence on sonic boom propagation in the Planetary Boundary Layer. This has been achieved within the frame of the European Commission funded project ATLLAS II (Aero-Thermodynamic Loads on Lightweight Advanced Stuctures II, 2011-2015, contract ACP0-GA-2010-263913, 7th Framework Programme). ATL-LAS II project aims at building tools to design a new aircraft for high speed flight with a cruise Mach number around 5. The overall project is split up along the vehicle design process and the development of advanced light-weight, high-temperature material. The globally optimized vehicle targets aerodynamic, propulsive, structural and thermal objectives but nevertheless takes into account restrictions imposed by emissions regulations and sonic boom mitigation. With this last objective, sonic boom of the ATLLAS II configuration is fully assessed in the present work. In particular, this assessment is achieved in a statistical way taking into account variability resulting from atmospheric turbulence. To our knowledge, such approach has never been used at the early design stage.

Content

The first part of the present manuscript details the FLHOWARD3D algorithm (Chapter 2). The model equation describing nonlinear shock wave propagation through a moving heterogeneous atmosphere [START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF] is first recalled. It is transformed into a partially one-way equation adapted to the numerical resolution. Resulting dispersion relation is compared to the exact one. The numerical algorithm is then detailed. It is based on split-step approach separating diffraction, flow, absorption and heterogeneities effects, and nonlinearities. Handling of the boundary conditions is then addressed. The resulting dispersion relation is established and discussed. Finally, parallel implementation of the algorithm is outlined.

In a second part, the FLHOWARD3D code is validated (Chapter 3). This is achieved by comparing the results of the code either to known analytical solutions, or to simulations performed with other codes. Various test cases are chosen to validate the different effects implemented in FLHOWARD3D: linear homogeneous diffraction, medium heterogeneities, flow motion, nonlinearities, absorption and boundary conditions. When possible, cases coupling several effects are investigated.

Sonic boom of the ATLLAS II hypersonic aircraft is assessed in Chapter 4. A model for the Planetary Boundary Layer is first proposed mixing Monin-Obukhov Similarity Theory (MOST) for the mean stratified atmosphere (wind and temperature), along with homogeneous isotropic fluctuations satisfying a von Kármán energy spectrum. Sonic boom is then evaluated first in the standard atmosphere with no turbulence as the reference case. This is achieved by matching near-field CFD simulations based on automatic mesh refinement to nonlinear ray tracing describing long range propagation. Sonic boom of ATLLAS II configuration is quantified both ground track and laterally at the edges of the geometrical carpet, both in terms pressure signature, peak overpressure and annoyance metrics. Statistical effect of turbulent wind fluctuations is then investigated for ground track boom by means of FLHOWARD3D software for two intensities of turbulence. Typical distorted waveforms and spectra are presented and peak overpressure variability is analysed. The sonic boom penetration in the shadow zone is also reproduced numerically comparing the cases with and without turbulence.

According to catastrophe theory, and in the linear regime, caustics are structurally stable, thus should persist in an heterogeneous medium. However, for shock waves, linear catastrophe theory is not sufficient and nonlinear effects have to be taken into account to predict the caustic field. Using the tools previously developed (FL-HOWARD3D software and turbulence synthesis), we are able to investigate numerically structural stability of cusped caustics in the unexplored case of nonlinear shock waves as detailed in Chapter 5. In this chapter, a model equation describing nonlinear shock wave propagation through a moving heterogeneous atmosphere [START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF] is first recalled. Then it is transformed into a partially one-way equation adapted to the numerical resolution. Resulting dispersion relation are computed and compared to the exact one and to the parabolic equation. The numerical algorithm is then detailed. It is based on splitstep approach separating diffraction, flow, absorption and heterogeneities effects, and nonlinearities. Handling of the boundary conditions is then addressed. The resulting numerical dispersion relation is established and discussed. Finally, parallel implementation of the algorithm is outlined.

Model

Equation for nonlinear propagation through a moving heterogeneous medium

A viscous fluid in a Cartesian domain is considered; x is the main propagation direction, z the vertical coordinate, and y the transverse one. At a position x and time t, we denote ρ(x, t) the density, c 0 (x, t) the sound speed, v(x, t) the flow velocity, and p(x, t) the pressure. Since the atmospheric characteristic time scale is much larger than the acoustical one (a few hours for atmospheric state, a few minutes for turbulent fluctuations and a fraction of seconds for a sonic boom signal), the medium can be considered as frozen which means that the quantities associated with the ambient medium (• 0 ) can be separated from the time varying acoustical fluctuations (• a ):

f (x, t) = f 0 (x) + f a (x, t) (2.1)
with f = (ρ, v, p). The ambient flow is then separated between a mean vertically stratified flow V 0 (z) with only horizontal components V 0 • e z = 0, which is of order of the ambient flow Mach number M , and turbulent fluctuations u 0 (x) of smaller order, approximately

M 2 : v 0 (x) = V 0 (z) + u 0 (x). (2.2)
This assumption is consistent with the situation in the atmospheric boundary layer where the mean wind velocity rarely exceeds 20 m/s, corresponding to a Mach number M = V 0 /c 0 of 0.06. Its fluctuations are much lower: a few meters per second [START_REF] Stull | An introduction to boundary layer meteorology[END_REF][START_REF] Wilson | Acoustic propagation through anisotropic, surface-layer turbulence[END_REF]. Hence, the Mach number of the fluctuation is approximately of the order of M 2 = 0.004. Finally, the temperature fluctuations are of order 5 °C to 10 °C which implies that the sound speed and density fluctuations are also of the order of one percent and thus O(M 2 ). Therefore they are also separated between a mean component (•) and a spatially varying one (• ′ ):

c 0 (x) = c0 + c ′ 0 (x) (2.3) ρ 0 (x) = ρ0 + ρ ′ 0 (x). (2.4)
Using these assumptions it is possible to model the nonlinear propagation of acoustic waves in a moving inhomogeneous medium with a nonlinear scalar equation [START_REF] Coulouvrat | New equations for nonlinear acoustics in a low Mach number and weakly heterogeneous atmosphere[END_REF]:

1 c 2 0 D 2 s p a Dt 2 -ρ 0 ∇ • ∇p a ρ 0 + 2 dV 0j dz t -∞ ∂ 2 p a (x, t ′ ) ∂z∂x j dt ′ = - 2 c 2 0 u 0 • ∂∇p a ∂t + β ρ 0 c 4 0 ∂ 2 p 2 a ∂t 2 + δ c 4 0 ∂ 3 p a ∂t 3 (2.5)
where the operator D s /Dt = ∂/∂t + V 0 • ∇ is the convective derivative associated to the ambient mean flow. Eq. 2.5 takes into account quadratic nonlinearities and atmospheric absorption; β = (γ + 1)/2 and δ = 4

3 µ + µ B + κ(c -1 v -c -1
p /ρ 0 are respectively the nonlinear parameter and diffusivity of the medium. Here µ is the shear viscosity, µ b the bulk viscosity and κ the thermal conductivity. c p and c v are the specific heats at constant pressure and constant volume, respectively. γ = c p /c v is the ratio of the specific heats and is equal to 1.4 in air so that β = 1.2. Eq. 2.5 is of mixed order between 1 and 2 in Mach number as it contains O(M 2 ) terms for V 0j in the convective derivative while neglecting the O(M 2 ) terms corresponding with the effect of shear associated with dV 0j dz . The first ones are involved in phase effects which are dominant as soon as frequency is high enough while the other ones are involved only in amplitude effects and therefore play a smaller role for the considered frequencies. These approximations are quantified numerically in [START_REF] Coulouvrat | New equations for nonlinear acoustics in a low Mach number and weakly heterogeneous atmosphere[END_REF][START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF][START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF]. The same work is done on the turbulent fluctuations where terms with du 0j dx i are neglected. Moreover, formally, the time derivatives associated with nonlinearities and absorption should be convective derivatives but the effect of convection can be neglected as explained in [START_REF] Coulouvrat | New equations for nonlinear acoustics in a low Mach number and weakly heterogeneous atmosphere[END_REF].

Partially one-way equation

Eq. 2.5 is not easily handled numerically so some transformations have to be performed. First, it is written in the form of a homogeneous wave equation with a perturbation term on the right hand side:

1 c0 2 ∂ 2 p a ∂t 2 - ∂ 2 p a ∂x 2 - ∂ 2 p a ∂y 2 - ∂ 2 p a ∂z 2 = P.
(2.6)

The right hand side P contains all terms linked to flow motions, medium heterogeneities, sound absorption and nonlinearities. They all are of order M at most. In a second stage, in order to use a one-way propagation approach, the retarded time τ = t -x/ c0 is introduced. Eq. 2.6 is rewritten in a time frame moving with the mean sound speed in the x-direction.

2 c0 2 ∂ 2 p a ∂x∂τ - ∂ 2 p a ∂x 2 - ∂ 2 p a ∂y 2 - ∂ 2 p a ∂z 2 = P ′ . (2.7)
As a third step, a wide-angle parabolic approximation is applied only on the perturbation terms P ′ rather than on the full wave equation 2.7. To do this, the second order derivative in x is replaced using the linear homogeneous equation written in retarded time:

∂ 2 p a ∂x 2 = 2 c0 2 ∂ 2 p a ∂x∂τ - ∂ 2 p a ∂y 2 - ∂ 2 p a ∂z 2 + O(M ). (2.8)
This wide-angle approximation is of order M 2 since it is applied only on terms that are already of order M at most. It is therefore of a higher order than applying a wideangle approximation to the full equation as P ′ itself is of order M . As a counterpart, a second order derivative in x remains on the homogeneous part of the equation. Thus it is only a partially one-way equation. This second order derivative will be dealt with during the numerical resolution in 2.2.2. Finally, in order to handle shock waves, the acoustic pressure p a is replaced by a pseudo-potential:

p a (x) = ∂φ ∂t (x).
(2.9)

The pseudo-potential has the advantage to remain continuous through shocks and is well adapted for the resolution of the nonlinear part of the equation [START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF]. It results finally in the following equation:

∂ 2 φ ∂x∂τ (x, τ ) = Dφ(x, τ ) + Hφ(x, τ ) + N φ(x, τ ) + Aφ(x, τ ), (2.10)
where D is an operator representing diffraction effects:

Dφ(x, τ ) = c0 2
∂ 2 φ ∂x 2 + ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 , ( 2.11)
The heterogeneities and wind effects are described by operator H decomposed as:

Hφ(x, τ ) = F LH (1) (x, τ ) + F LH (2) (x, τ ) + T H(x, τ ).
(2.12)

F LH (1) is the operator for the effect of the mean vertically stratified flow on propagation:

F LH (1) φ(x, τ ) = V 0x c0 2 ∂ 2 φ ∂τ 2 - V 0x c0 ∂ 2 φ ∂x∂τ - V 0y c0 ∂ 2 φ ∂y∂τ + V 2 0x 2 c 0 ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 - V 2 0x 2 c0 3 ∂ 2 φ ∂τ 2 + V 0x V 0y c0 1 c0 ∂ 2 φ ∂τ ∂y - ∂ 2 φ ∂x∂y - V 2 0y 2 c 0 ∂ 2 φ ∂y 2 + c0 V 0x τ -∞ ∂ 3 φ ∂x∂z 2 dτ ′ -c0 τ -∞ ∂ ∂z V 0x ∂ 2 φ ∂x∂z dτ ′ + d dz V 0x ∂φ ∂z -V 0x ∂ 2 φ ∂z 2 + c0 V 0y τ -∞ ∂ 3 φ ∂y∂z 2 dτ ′ -c0 τ -∞ ∂ ∂z V 0y ∂ 2 φ ∂y∂z dτ ′ , ( 2.13)
First line of Eq. 2.13 describes linear convection terms by the mean stratified flow.

The second and third lines are for the nonlinear quadratic convection terms. The last three lines emanate from the gradient of the ambient flow. F LH (2) is the operator for the effect of the turbulent fluctuations on propagation:

F LH (2) φ(x, τ ) = u 0x c0 2 ∂ 2 φ ∂τ 2 - u 0x c0 ∂ 2 φ ∂x∂τ - u 0y c0 ∂ 2 φ ∂y∂τ - u 0z c0 ∂ 2 φ ∂z∂τ . (2.14)
Eq. 2.14 takes only into account the linear convection effects due to the flow fluctuations. T H is the operator describing the influence of the sound speed and density of the medium inhomogeneities on the propagation:

T Hφ(x, τ ) = 1 2ρ 0 ∂ρ 0 ∂x ∂φ ∂τ -c0 ∂ρ 0 ∂x ∂φ ∂x - c0 2ρ 0 ∂ ∂y ρ 0 ∂φ ∂y + ∂ ∂z ρ 0 ∂φ ∂z + c0 2
∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 + ǫ 2 c0 3 ∂ 2 φ ∂τ 2 , (2.15) with ǫ = 2 c0 c ′ 0 + c ′2 0 . (2.16)
The first line of Eq. 2.15 stands for the density heterogeneities in the main propagation direction x, the second line for those in the transverse directions y and z. The last term takes into account all effects due to sound speed heterogeneities. Note that in the definitions of F LH (1) (Eq. 2.13), F LH (2) (Eq. 2.14) and T H (Eq. 2.15) the components V 0x , u 0x and ∂ρ 0 /∂x appear differently from the other ones because of the introduction of retarded time in the main propagation direction x. A simplified FLHOWARD equation can be written by considering only the dominant linear convection flow effects for both the mean wind and its fluctuations. In this case, the H operator has to be replaced by H (s) defined as:

H (s) φ(x, τ ) = F LH (s) (x, τ ) + T H(x, τ ), (2.17)
with

F LH (s) φ(x, τ ) = v 0x c0 2 ∂ 2 φ ∂τ 2 - v 0x c0 ∂ 2 φ ∂x∂τ - v 0y c0 ∂ 2 φ ∂y∂τ - v 0z c0 ∂ 2 φ ∂z∂τ . (2.18)
Here, only the dominant linear convection effects are taken into account for the total ambient velocity field v 0 defined in Eq. 2.2. Note that for the velocity fluctuations u 0 (x), operators F LH (s) and F LH (2) are identical. Operators H and H (s) differ only for the mean stratified flow V 0 (z). Finally, N is the operator associated to the nonlinearities:

N φ(x, τ ) = β 2 ρ0 c0 3 ∂ ∂τ ∂φ ∂τ 2 , ( 2.19)
and A to the absorption:

Aφ(x, τ ) = ∂ ∂τ δ 2 c0 3 ∂ 2 φ ∂τ 2 .
(2.20)

More details on the establishment of Eq. 2.10 are given in Appendix A. The main physical effects are clearly separated in equation 2.10. We will take advantage of this for the numerical resolution of the problem in 2.2.1.

Dispersion relation

The validity of the model equation 2.10 in linear propagation (operators N φ(x) and Aφ(x) are set to zero) can be examined precisely in two cases by writing its dispersion relation for a bi-dimensional plane wave:

φ(x, y, z, τ ) = A exp ik 0 ( kx -1)x + kz z -c0 τ (2.21)
where k 0 = ω 0 /c 0 is the wave number, and kx = k x /k 0 , kz = k z /k 0 are the dimensionless components of the wave vector. Two cases will be investigated.

Sound speed heterogeneities:

The first case was studied by Dagrau et al. [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF], it is a no-flow case with ρ ′ 0 = 0 and c ′ 0 is constant. In this case, for Eq. 2.10 we have:

kx = 1 -ǫ/ c0 2 -kz 2 (2.22)
while the exact dispersion relation is:

kx = c0 2 /c 2 0 -kz 2 .
(2.23)

To retrieve the exact dispersion relation and get a higher precision, ǫ can be redefined as:

ǫ(x) = c0 2 1 - c0 2 c 2 0 (x) = 2 c0 c ′ 0 + O(c ′2 0 ) = 2 c0 c ′ 0 + O(M 4 ). (2.24)
Velocity heterogeneities: The second case investigated is the propagation of a plane wave in a uniform flow parallel to the direction of propagation (V 0x = cte, V 0y = 0). The exact equation is given by:

D 2 s p a Dt 2 -c 2 0 ∆p a = 0 (2.25)
which corresponds to the dispersion relation:

kx 2 + 2M 1 -M 2 kx + kz 2 -1 1 -M 2 = 0. (2.26)
Flhoward equation for F LH (1) (in this case F LH (2) = 0) is written:

∂ 2 φ ∂x∂τ = c0 2
∂ 2 φ ∂x 2 + ∂ 2 φ ∂z 2 + V 0x c0 2 ∂ 2 φ ∂τ 2 - V 0x c0 ∂ 2 φ ∂τ ∂x + V 2 0x 2 c0 ∂ 2 φ ∂y 2 - 1 c0 2 ∂ 2 φ ∂τ 2 (2.27)
and its associated dispersion relation is:

kx 2 + 2M kx + kz 2 -1 1 + M 2 = 0. (2.28)
Equations 2.26 and 2.28 differ only by a term of order O(M 3 ) which is consistent with the assumption made in establishing Eq. 2.10. To obtain a higher order of precision, the Mach number M can be replaced by M/ 1 -M 2 , in a way similar as was done previously for the sound speed velocity. The difference between the two dispersion relations then becomes a term of order O(M 4 ) [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF][START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF]. The lower order FLHOWARD equation for F LH (s) :

∂ 2 φ ∂x∂τ = c0 2
∂ 2 φ ∂x 2 + ∂ 2 φ ∂z 2 + V 0x c0 2 ∂ 2 φ ∂τ 2 - V 0x c0 ∂ 2 φ ∂τ ∂x (2.29)
has the dispersion relation:

kx 2 + 2M kx + kz 2 -1 = 0, (2.30)
which clearly differs from Eq. 2.26 by terms of orderM 2 . Finally, the last equation to compare is the parabolic equation of Averyianov et al. [START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF][START_REF] Averiyanov | Parabolic equation for nonlinear acoustic wave propagation in inhomogeneous moving media[END_REF][START_REF] Averiyanov | Propagation des ondes acoustiques à travers un milieu turbulent: Etudes théorique et expériementale des effets de diffraction et des effets non linéaire[END_REF]:

∂ 2 φ ∂x∂τ = c0 2
∂ 2 φ ∂z 2 + V 0x c0 2 ∂ 2 φ ∂τ 2 (2.31)
which leads to the dispersion relation:

kx -1 + M + kz 2 2 = 0. (2.32)
On Fig. 2.1, the various dispersion relations are compared for M = [0.0, 0.5]. It can be seen that the full FLHOWARD equation with F LH (1) is highly accurate as expected.

Simplified FLHOWARD equation with F LH (s) is also more accurate than the classical parabolic equation. The improvement is due mostly to the diffraction term ∂ 2 φ ∂x 2 and not to the higher order flow terms especially for the low Mach numbers, since there are only small differences between F LH (1) and F LH (s) . Also, it is to be noted that both approximations behave really well even with relatively high Mach numbers, even if the Mach numbers were supposed to be small for the derivation of the equations.

Numerical resolution 2.2.1 Split-step method

As explained in the previous section, Eq. 2.10 is set in a time frame moving at the mean sound speed in the propagation direction. Thus a spatially advancing scheme in the x-direction is needed. A split-step method is chosen [START_REF] Yanenko | The method of fractional steps[END_REF]. It allows to solve each part of Eq. 2.10 separately and then to couple the solution by assembling the results. The L.H.S. is the coupling term.

φ(x + ∆x, y, z, τ ) = φ D ∆x • φ H+A ∆x • φ N ∆x (x, y, z, τ ) + O(∆x). (2.33) with • is the operator composition: g • f (x) = g(f (x)). φ X ∆x is solution of the formal equation ∂ 2 φ ∂x∂τ (x, τ ) = Xφ(x, τ
) on the ∆x step with operator X = D, H + A or N . Equation (2.33) is a first order split-step. For numerical reasons explained in 2.2.5, absorption and heterogeneities are solved during the same sub-step. Alternatively, the following second order split-step [START_REF] Strang | On the construction and comparison of difference schemes[END_REF][START_REF] Yanenko | The method of fractional steps[END_REF] can be used:

φ(x + ∆x, y, z, τ ) = φ N ∆x 2 • φ D ∆x 2 • φ H+A ∆x • φ D ∆x 2 • φ N ∆x 2 (x, y, z, τ ) + O(∆x 2 ). (2.34)
The main advantage of the split-step scheme is that each part of the equation is solved using an efficient numerical method adapted for each physical effect:

• Diffraction using the angular spectrum method [START_REF] Goodman | Introduction to Fourier Optics[END_REF] in the 3D spectral space (ω, k y , k z ),

• Heterogeneities, wind and absorption in the frequency domain (ω, y, z): an analytic solution is used for phase effects, the remaining terms are solved using finite differences, • Nonlinear effects in the physical domain (τ, y, z) with the semi-analytic Burgers-Hayes method [START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF].

Between these steps, Fourier transforms are used to change the space in which the resolution is done.

φ i (t, y, z) φ1 (ω, ky, kz) φ2 (ω, ky, kz) φ 3 (ω, y, z) φ 4 (ω, y, z) φ 5 (t, y, z) φ i+1 (t, y, z) 3 D F F T 2 D F F T 1 D F F T Diffraction
Heterogeneities and Flows Nonlinearities Figure 2.2: First order split-step scheme.

Fig. 2.2 shows the algorithm of the first order split-step used for the resolution of Eq. 2.10.

Resolution of diffraction effects

Diffraction effects are contained in the operator D defined in Eq. 2.11. It results into the equation:

∂ 2 φ ∂x∂τ = c0 2
∂ 2 φ ∂x 2 + ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 . (2.35)
Eq. 2.35 is solved using the angular spectrum method in the spectral domain (k y , k z , ω).

The time Fourier transform along the τ axis is defined:

φ(y, z, ω) = ∞ -∞ φ(y, z, τ )e -iωτ dτ. (2.36)
The 2D spatial Fourier transform in the plane orthogonal to the main propagation direction x is:

φ(k y , k z , τ ) = ∞ -∞ ∞ -∞
φ(y, z, τ )e -ikyy e -ikzz dydz.

(2.37)

Using these transforms, equation 2.35 becomes:

d 2 φ dx 2 -2ik d φ dx -k 2 y + k 2 z φ = 0 (2.38)
with k = ω/ c0 . This is a second order ordinary differential equation which can be solved analytically. Its discriminant is:

∆ = 4 -k 2 + k 2 y + k 2 z . (2.39)
Eq. 2.38 has two solutions: one for propagation in the positive x direction and one for the propagation in the negative x direction. Only the solution propagating in the positive x direction is selected in the algorithm: the backscattered field is cancelled. Thus the numerical solution is fully one-way. If -k 2 + k 2 y + k 2 z > 0, only the evanescent wave propagating and decaying in the positive x-direction is selected. The solution is:

φ(x + ∆x) = φ(x)exp ∆x ik --k 2 + k 2 y + k 2 z . (2.40) If -k 2 + k 2 y + k 2 z < 0, only the wave propagating in the positive x-direction is selected. The solution is: φ(x + ∆x) = φ(x)exp ∆x ik -i k 2 -(k 2 y + k 2 z ) . (2.41)
The solution is then retrieved in the physical space using the inverse Fourier transforms. Practically, the Fourier transforms are implemented using Fast Fourier Transform (FFT) algorithm provided by the FFTW library [START_REF]FFTW official[END_REF]. This method has the advantages to be fast due to the speed of the FFT algorithm while keeping a spectral accuracy.

Resolution of flows and heterogeneities effects

Heterogeneities and flows effects are contained in operator H which is rewritten as

H = H 1 + H 2 with: H 1 = V 0x c0 2 ∂ 2 φ ∂τ 2 - V 0x c0 ∂ 2 φ ∂x∂τ - V 2 0x 2 c0 3 ∂ 2 φ ∂τ 2 + u 0x c0 2 ∂ 2 φ ∂τ 2 - u 0x c0 ∂ 2 φ ∂x∂τ + ǫ 2 c0 3 ∂ 2 φ ∂τ 2 + 1 2ρ 0 ∂ρ 0 ∂x ∂φ ∂τ -c0 ∂ρ 0 ∂x ∂φ ∂x , ( 2.42)
and

H 2 = - V 0y c0 ∂ 2 φ ∂y∂τ - u 0y c0 ∂ 2 φ ∂y∂τ - u 0z c0 ∂ 2 φ ∂z∂τ + V 2 0x 2 c 0 ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 + V 0x V 0y c0 1 c0 ∂ 2 φ ∂τ ∂y - ∂ 2 φ ∂x∂y - V 2 0y 2 c 0 ∂ 2 φ ∂y 2 + c0 V 0x τ -∞ ∂ 3 φ ∂x∂z 2 dτ ′ -c0 τ -∞ ∂ ∂z V 0x ∂ 2 φ ∂x∂z dτ ′ + c0 V 0y τ -∞ ∂ 3 φ ∂y∂z 2 dτ ′ -c0 τ -∞ ∂ ∂z V 0y ∂ 2 φ ∂y∂z dτ ′ + d dz V 0x ∂φ ∂z -V 0x ∂ 2 φ ∂z 2 + c0 2
∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 - c0 2ρ 0 ∂ ∂y ρ 0 ∂φ ∂y + ∂ ∂z ρ 0 ∂φ ∂z . (2.43)
For a plane wave, H 1 handles the phase effects which are preponderant as long as the frequency is high enough, which is the case for the considered applications. It corresponds to the terms of Eq. 2.12 containing only x and τ derivatives. It is worth noting that phase effects in operator H 1 involve only: (i) convection in the main propagation direction x (terms linear or nonlinear with respect to V 0x and linear with respect to u 0 , (ii) sound speed heterogeneities, (iii) gradient of density in the main propagation direction. The remaining terms of Eq. 2.12 are contained in the operator H 2 which represents the coupling between diffraction, heterogeneities and wind for a non planar wave. Amplitude effects in operator H 2 involve: (i) linear convection due to the transverse component of the velocity field (first line), (ii) some quadratic convection terms (second line), (iii) gradient of ambient flow (lines 3 to 5), (iv) density gradients in the transverse directions (last line). The effects contained in operator H are handled in the frequency domain (ω, y, z). For stability and accuracy reasons [START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF], all phase effects are solved analytically. The coupling of the multiple effects is achieved through a second order finite differences scheme with semi-implicit second order Crank-Nicolson advancement scheme. For the simplified FLHOWARD equation case, operator H 1 is replaced by operator

H (s)
1 defined as:

H (s) 1 = v 0x c0 2 ∂ 2 φ ∂τ 2 - v 0x c0 ∂ 2 φ ∂x∂τ + ǫ 2 c0 3 ∂ 2 φ ∂τ 2 + 1 2ρ 0 ∂ρ 0 ∂x ∂φ ∂τ -c0 ∂ρ 0 ∂x ∂φ ∂x , ( 2.44)
and operator H 2 is replaced by operator H (s) 2 defined as:

H (s) 2 = - v 0y c0 ∂ 2 φ ∂y∂τ - v 0z c0 ∂ 2 φ ∂z∂τ + c0 2
∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 - c0 2ρ 0 ∂ ∂y ρ 0 ∂φ ∂y + ∂ ∂z ρ 0 ∂φ ∂z . ( 2

.45)

In this simplified case, phase effects in operator H (s) 1 involve only: (i) linear convection in the main propagation direction x, (ii) sound speed heterogeneities, (iii) gradient of density in the main propagation direction. Amplitude effects in operator H 2 involve only: (i) linear convection due to the transverse components of the velocity field (first line), (ii) density gradients in the transverse directions (last line).

Phase effects

Phase effects are defined through operator H 1 in Eq. 2.42 and are solved analytically. The equation to solve is:

∂ 2 φ ∂x∂τ = V 0x c0 2 ∂ 2 φ ∂τ 2 - V 0x c0 ∂ 2 φ ∂x∂τ - V 2 0x 2 c0 3 ∂ 2 φ ∂τ 2 + u 0x c0 2 ∂ 2 φ ∂τ 2 - u 0x c0 ∂ 2 φ ∂x∂τ + ǫ 2 c0 3 ∂ 2 φ ∂τ 2 + 1 2ρ 0 ∂ρ 0 ∂x ∂φ ∂τ -c0 ∂ρ 0 ∂x ∂φ ∂x . (2.46)
In the frequency domain, it can be cast as a first order linear ordinary differential equation:

iω 1 + V 0x c0 + u 0x c0 + c0 2ρ 0 ∂ρ 0 ∂x ∂ φ ∂x = ω 2 V 2 0x 2 c0 3 - V 0x c0 2 - u 0x c0 2 - ǫ 2 c0 3 + iω 2ρ 0 ∂ρ 0 ∂x φ.
(2.47)

Solution for this equation is given explicitly by:

φ(x + ∆x, y, z, ω) = φ(x, y, z, ω)exp x+∆x x Θ(ζ)dζ , ( 2.48)
with:

Θ(ζ) = ω 2 V 2 0x 2 c0 3 -V 0x c0 2 -u 0x c0 2 -ǫ 2 c0 3 + iω 2ρ 0 ∂ρ 0 ∂x iω 1 + V 0x c0 + u 0x c0 + c0 2ρ 0 ∂ρ 0 ∂x . (2.49)
To compute the integral, the second order trapezoidal rule (Simpson's rule) is used:

φ(x + ∆x, y, z, ω) = φ(x, y, z, ω) exp ∆x 2 [Θ(x) + Θ(x + ∆x)] + O(∆x 2 ). (2.50)
For the simplified FLHOWARD equation case, Θ(ζ) is replaced by Θ (s) (ζ):

Θ (s) (ζ) = -ω 2 v 0x c0 2 + ǫ 2 c0 3 + iω 2ρ 0 ∂ρ 0 ∂x iω 1 + v 0x c0 + c0 2ρ 0 ∂ρ 0 ∂x . (2.51)

Coupling effects

Operator H 2 which represents the coupling between diffraction, heterogeneities and wind for a non planar wave (Eq. 2.43) cannot be solved analytically. The finite difference method in frequency domain is employed. It is of second order accuracy in the y and z directions. In the x-direction, a semi-implicit Crank-Nicolson scheme is used. This scheme is unconditionally stable and of second order accuracy. We resort to an Alternate Direction Implicit (ADI) method to separate the y and z directions. The two resultant matrices are tridiagonal, so that they can be solved using the Thomas' algorithm given in [START_REF] Quarteroni | Numerical mathematics[END_REF]. Details are given in Appendix B. The method is similar for the simplified FLHOWARD equation for operator

H (s) 2 .

Resolution of nonlinear effects

The nonlinear effects are taken into account by the operator N (Eq. 2.19) which results into an inviscid Burgers' equation:

∂φ ∂x = β 2 ρ0 c0 3 ∂φ ∂τ 2 .
(2.52)

To solve it, the quasi-analytical Burgers-Hayes methods is chosen. It relies on the implicit Poisson solution. This solution is exact while it is single-valued (no shock).

In case of a multi-valued solution (after the shock formation), the weak shock theory is employed. The classical equal area rule of Landau [START_REF] Landau | On Shock Waves at Large Distances from the Place of Their Origin[END_REF] is computationally too expensive to be practical. Instead, Burgers-Hayes [START_REF] Burgers | Further statistical problems connected with the solution of a simple nonlinear partial differential equation[END_REF][START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF] method handles it using the fact that the physical solution for the potential corresponds to the maximum value of the multivalued Poisson's solution as can be seen on Fig. 2.3. This condition is necessary to satisfy the second principle of thermodynamics.

φ(x, y, z, τ ) = max φ(x, y, z, θ) - β∆x ρ0 c0 3 ∂φ(x, y, z, θ) ∂θ τ = θ - β∆x ρ0 c0 3 ∂φ(x, y, z, θ) ∂θ . ( 2.53)
A numerical interpolation is needed to retrieve the solution Eq. 2.53 on the initial retarded time τ grid. Practically, the order of this interpolation makes no differences thus only a first order interpolation is used.

Absorption and relaxation

The operator A deals with the absorption of the medium. In the atmosphere, at audible frequencies or below, the thermoviscous absorption included in Eq. 2.20 is not the main attenuation mechanism. Dominant absorption terms are due to relaxation of molecular nitrogen and oxygen and is controlled by the water vapour content [START_REF] Hamilton | Nonlinear acoustics[END_REF][START_REF] Sutherland | Atmospheric absorption in the atmosphere up to 160 km[END_REF][START_REF] Cleveland | Time-domain modeling of finite-amplitude sound in relaxing fluids[END_REF][START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF]. Thus an ad-hoc expression is added to include attenuation and dispersion due to relaxation:

∂φ ∂x = L τ (φ). (2.54)
L τ is a linear operator containing the effects of attenuation and dispersion in a thermoviscous fluids containing multiple relaxation mechanisms. It is given by:

L τ (φ) = δ 2 c0 3 ∂ 2 φ ∂τ 2 + j c ′ j c0 2 τ -∞ ∂ 2 φ ∂τ ′2 e -(τ -τ ′ )/t j dτ ′ (2.55)
with δ the diffusivity of sound defined in 2.1, t j the relaxing time of the j th relaxation process and c ′ j the increase in phase velocity associated with this mechanism when the frequency goes from 0 to infinity. Eq. 2.54 is handled in the frequency domain (ω, y, z):

∂ φ ∂x =   -ω 2 δ 2 c0 3 - ω c0 2 j c ′ j ωt j 1 + iωt j   φ.
(2.56)

As for the heterogeneities, we solve this ordinary differential equation:

φ(x + ∆x) = φ(x)exp   - ω 2 δ 2 c0 3 - ω c0 2 j c ′ j ωt j 1 + iωt j   (2.57)
Practically, since the method is the same as for the phase effects of the heterogeneities H 1 , the resolution is integrated in the same step to limit the computation time. Details on the absorption coefficients in the atmosphere [START_REF] Iso | Acoustics -Attenuation of sound during propagation outdoors -Part 1: Calculation of the absorption of sound by the atmosphere[END_REF] are given in Appendix C.

Boundary conditions

The code is designed to perform simulations of acoustic propagation in the atmosphere near the ground (Fig. 2.4). Hence, at least two kinds of boundary conditions are required. The perfectly rigid ground is modelled using a perfectly reflecting boundary condition. The other boundary conditions are handled as artificial absorbing layers mimicking infinite medium in the vertical (z > 0) and lateral (y) directions.

Reflecting boundary conditions

The perfectly rigid condition is given by:

∂p a ∂z z=0 = 0. (2.58)
The difficulty to implement this boundary condition comes from the use of Fourier transform to solve the diffraction operator. The classical Fourier transform enforces periodicity on the boundaries. The simplest way to take into account this boundary conditions is the image method [START_REF] Morse | Theoretical acoustics[END_REF] which consists in solving the symmetric problem of the one we are solving. Unfortunately it is numerically inefficient both in terms of computation time and of memory usage since one has to solve the problem on a domain twice as big as needed. This solution was previously implemented in FLHOWARD 2D [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF] but is not tractable in 3D due to the memory and computation time increase. An other common solution is to solve the problem on a grid generated by using Chebyshev or Legendre polynomials [START_REF] Fornberg | A practical guide to pseudospectral methods[END_REF]. With this method, the grid is no longer Cartesian, which is more complicated to implement. To overcome this problem, FLHOWARD 3D uses the cosine transform to enforce von Neumann boundary conditions (Eq. 2.58):

φ(y, k z , τ ) = ∞ -∞ φ(y, z, τ )cos(ik z z)dz. (2.59)
Discrete Cosine Transform (DCT) naturally enforces even parity at both of its extremities [START_REF]FFTW official[END_REF] and can be computed efficiently using FFTW. An other advantage is that Dirichlet boundary condition p a (z = 0) = 0 can also be modelled using the Sine transform. Also, combination of Cosine and Sine transform should allow us to extend this work to the propagation over a finite impedance ground as has been done in electromagnetism [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF].

Absorbing boundary layer

The current state of the art for enforcing non reflecting boundary conditions is the Perfectly Matched Layers (PML) [START_REF] Collino | Perfectly matched absorbing layers for the paraxial equations[END_REF][START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF]. With PML, waves arriving with any incidence angle will not reflect on the boundary. Instead of PML, the simpler Absorbing Boundary Layer (ABL) reduces the reflection on the border of the domain. It consists in introducing an artificial absorption term in an upper layer in order to absorb the incident wave as shown on Fig. 2.4. This artificial absorption term is implemented in addition to the physical absorption. The absorption coefficient is quadratically increasing from the beginning of the layer to the domain limit. It is the same as implemented in FLHOWARD 2D [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF][START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF]:

φ(x + ∆x, y, z, ω) = φ(x, y, z, ω)exp (-α(z)) (2.60) with α(z) ==      21 (z -L) 2 h 3 if z > L, 0 else, (2.61)
where L is the height of the boundary and h is the width of the ABL.

Dispersion relation

On Sec. 2.1.3, the dispersion relation of Eq. 2.10 was examined. This did not take into account the effect of the numerical resolution. Here, the effect of the split-step method on the dispersion relation is investigated. As was done in 2.1.3, a plane wave solution under the form of Eq. 2.21 is injected in the studied equation. First, for the angular spectrum method associated to operator D:

∂ 2 φ ∂x∂τ = c0 2
∂ 2 φ ∂x 2 + ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 , ( 2.62)
one gets:

k x1 2 = 1 -kz 2 .
(2.63)

Then equation associated to operator F LH (1) reduces to:

∂ 2 φ ∂x∂τ = V 0x c0 2 ∂ 2 φ ∂τ 2 - V 0x c0 ∂ 2 φ ∂x∂τ + V 2 0x 2 c0 ∂ 2 φ ∂z 2 - 1 c0 2 ∂ 2 φ ∂τ 2 . (2.64)
Its dispersion relation is:

k x2 = 1 1 + M 0 1 - M 2 0 2 kz 2 + 1 . (2.65)
The total dispersion relation for the split-step method is given by:

kx = k x1 + k x2 -1, (2.66)
and so for the combination of operators D and F LH (1) the dispersion relation is:

kx = 1 -kz 2 + 1 1 + M 0 1 - M 2 0 2 kz 2 + 1 -1. (2.67)
Similarly, for the operators D and F LH (s) , we get:

kx = 1 -kz 2 + 1 1 + M 0 -1. (2.68)
As was done in 2.1.3 a comparison is made with the parabolic approximation of Averyianov Eq. 2.31 whose dispersion relation is unaffected by the split-step method since it involves first order effects which are completely uncoupled:

kx = 1 - kz 2 2 -M 0 . (2.69)
The first thing that can be seen on Fig. 2.5, is that the dispersion relation for the FLHOWARD method is indeed one way. For M < 0.3 both equations for F LH (1) and F LH (s) agree really well with the exact dispersion relation. Starting with M = 0.3, some discrepancies begin to appear for all wave numbers for F LH (1) while, surprisingly, only propagation at large angles induces an error for F LH (s) . So with this numerical method, adding the higher order terms leads to a decrease in accuracy for the case of a uniform flow. For comparison, the parabolic equation has a good accuracy only for the same range of Mach number, as F LH (s) but, as expected, also only for small angles of propagation.

To better understand why F LH (1) induces a larger error at large Mach numbers, a Taylor expansion with respect to Mach number M 0 up to order O(M 3 0 ) is performed on the above numerical dispersion relations. They become:

• for the exact convected wave equation:

kx = 1 -kz 2 -M 0 + M 2 0   1 -kz 2 + 1 2 kz 2 1 -kz 2   ; (2.70)
• for FLHOWARD with operator F LH (1) :

kx F LH (1) = 1 -kz 2 -M 0 + M 2 0 2 1 -kz 2 ;
(2.71)

• for FLHOWARD with operator F LH (s) :

kx F LH (s) = 1 -kz 2 -M 0 + M 2 0 ; (2.72)
• for Averiyanov parabolic equation:

kxP E = 1 - kz 2 2 -M 0 . (2.73)
The first three ones are strictly identical up to order M 0 thanks to the angular spectrum method. As expected, parabolic equation is valid only for small angles of propagation.

At second order M 2 0 , due to the combined effects of equation approximations and numerical split step, the F LH (1) operator induces a systematic error of order M 2 0 even for small angles of propagation when kz << 1. On the contrary, the F LH (s) turns out to be precise up to order M 2 0 kz 4 . Therefore it remains much more accurate as long kz is not to large.

Consequently, for practical applications this analysis proves the use of FLHOWARD equation with simplified operator F LH (s) is recommended. This is due to the numerical split-step method, combined with the underlying approximations of the FL-HOWARD equation.

Parallel computing

Even if the one-way methods are computationally cheaper than DNS calculations, for three dimensional computations, it is necessary to use distributed memory architectures. For an acoustic wave at a frequency f 0 = 10 Hz in a domain (x, y, z) = (3 × 3 × 3 km 3 ) using ten points per wavelength we have δx = 3.4 m. It means that around one thousand points in each directions are required. Around the same number of points is needed for the temporal discretization. Thus the problem has around 10 12 degrees of freedom. Distributed memory computers have many processors, each with its own memory. They allow to tackle bigger computations with the restitution time remaining acceptable. Currently, this kind of computation is classical and a lot of finite-differences, finite-elements or finite-volumes code implement it. The domain can be decomposed so as to allocate one sub domain to one processor which performs its computations on this small part of the numerical domain. It is also known as the Single Program Multiple Data (SPMD) paradigm. However, here, this strategy is complicated due to the presence of Fourier transforms in the algorithm. A Fourier transform is a non-local operation incompatible with the SPMD paradigm. To handle this, the 1D decomposition (also known as slice decomposition) is implemented. The domain is only decomposed in one direction, the two others remaining local to the processor as shown in Fig. 2.6. With this method the Fourier transform can be applied on two directions. For the last direction, a parallel transpose is done so that the decomposed direction is changed. This parallel transpose requires a lot of data to be transferred between the processors: it is a MPI_ALL_TO_ALL instruction which means that each processor broadcasts the data to the others. Communication between processors is a time consuming task. This heavy exchange of data is the bottleneck for the proper scaling of the code. However as it will be shown in Sec. 3.5, it turns out that its effect on the scaling is lower than expected. Practically, the 1D decomposition is implemented thanks again to the FFTW library which provides FFT and DCT operators, and parallel transposition [START_REF]FFTW official[END_REF].

Figure 2.6: Scheme for 1D domain decomposition on two processors. In red, the direction which is local for the processor.

Conclusion

An original three-dimensional one-way method to model and simulate the nonlinear propagation of acoustic shock waves in the atmosphere has been developed. This method strives to be as high order as possible while remaining computationally tractable.

To do so, it relies on the fractional step method which permits to solve different simpler problems using methods that are as efficient as possible. When possible, spectral or analytical solutions are employed. For the other terms, finite differences method is used. To remain consistent, both the finite differences and split-step schemes are of second order accuracy. A lot of work has been devoted to the implementation of boundary conditions: rigid ground and free field. Contrary to most of the current nonlinear codes which are implemented on single processor or shared-memory computers [START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF][START_REF] Jing | Modeling the propagation of nonlinear threedimensional acoustic beams in inhomogeneous media[END_REF][START_REF] Yuldashev | Simulation of three-dimensional nonlinear fields of ultrasound therapeutic arrays[END_REF], the method has been implemented for high performance computing on the distributed memory architecture. The complete algorithm is illustrated Fig. End In this chapter, the validation of the FLHOWARD3D code is presented. This is achieved by comparing the results of the code either to known analytical solutions, or to simulations performed with other codes. The validation process is similar to the one used to validate HOWARD [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF][START_REF] Dagrau | Simulation de la propagation du bang sonique : de la CFD à l'acoustique non linéaire[END_REF] and FLHOWARD [START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF][START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF] but with new complementary test cases. These various test cases are chosen to validate the different effects implemented in FLHOWARD3D: linear homogeneous diffraction, medium heterogeneities, flow motion, nonlinearities, absorption and boundary conditions. When possible, cases coupling several effects are investigated. The list of selected tests and the part of the algorithm that they assess is summarized in Table 3.1.

φ i (y, τ, z) φ 0 (y, τ, z) φ 1 (y, τ, z) φ1 (ω, k y , k z ) φ2 (ω, k y , k z ) φ2 (ω, y, z) φ3 (ω, y, z) φ3 (ω, k y , k z ) φ4 (ω, k y , k z ) φ 4 (y, τ, z) φ i+1 (y, τ, z) φ n (

Three-dimensional circular acoustic piston in a linear homogeneous medium

The first part of the validation process is the 3D propagation of a pure tone wave in a linear homogeneous medium. It is chosen to validate at 3D the diffraction part (operator D) of the algorithm and the Absorbing Boundary Layer (ABL). The case of a circular piston in three-dimensions (Fig. 3.1) is considered. The analytical solution of this problem [START_REF] Morse | Theoretical acoustics[END_REF][START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF] is given in Appendix D.1. Note that the name piston is somewhat inappropriate in this case since pressure is imposed on the surface whereas Section Operators Dimension Reference solution 3.1 (1) or D + F LH (s) ) + Rigid 2D Semi-analytical 3.4.2 D + F LH (2) 2D Numerical for a true piston the velocity would be imposed. For the classical parabolic equation it would be meaningless to make the difference but for higher-order methods such as FLHOWARD, it is not equivalent. The parameters are chosen to correspond to highintensity focused ultrasound (HIFU) treatment: c 0 = 1500 m/s, ρ 0 = 1000 kg/m 3 . The incident wave is sinusoidal with a frequency f = 1 MHz and a wavelength λ = 0.0015 m. The radius of the piston is r = 0.006 m which corresponds to 4λ. A uniform pressure of 5 Pa is applied on this surface and the pressure is set to zero outside. The computational domain extends from 0 m to 0.25m (0 to 250λ) in the propagation direction (x), -0.07 m to 0.07m (-33λ to 33λ) in the transverse directions (y, z) and from 0 s to 10 -6 s (one period) in time. Numerically, each direction (x, y, z, τ ) is discretized using 1024 points. Therefore the problem has 10 12 degrees of freedom. Absorbing boundary layers of size 0.02 m are used so that the physical domain is reduced to -0.05 m to 0.05 m in the transverse directions. Fig. 3.2 shows the pressure amplitude radiated from the piston in the (x, y) plane where the characteristic diffraction pattern can be seen. An interesting thing is that there is no reflected wave at the top and bottom boundaries. This shows that the ABLs are working properly. This can also be seen on the axis pressure Fig. analytical solution in Fig. 3.3 shows a good agreement even in the near-field which is usually not reproduced correctly by the parabolic approach.

D + ABL 3D Analytical 3.2 D + T H 3D Analytical 3.3 N + A 1D Analytical 3.4.1 (D + F LH

Scattering of a plane wave by a a spherical heterogeneity

The scattering of a plane wave by a discontinuous spherical heterogeneity in sound speed shown Fig. m to 1020m (0 to 15λ) in the propagation direction (x), -2040 m to 2040m (-30λ to 30λ) in the transverse directions (y and z) and from 0 s to 0.2 s (one period) in time. Numerically, each direction (x, y, z, τ ) is discretized using 1024 points. Reflecting boundary conditions are used. They are far enough from the heterogeneity so as not to interfere with the transmitted field in the considered domain. This case is much more difficult than the atmospheric application the model was designed for. First the heterogeneity is discontinuous whereas the parameters of the atmosphere are varying continuously. Second, the intensity of the heterogeneity, 5% in sound speed, corresponds to a variation of temperature of 34 °C, extremely large for the atmosphere over such a short distance. Fig. 3.5 shows the pressure amplitude scattered by the heterogeneity in the (x, y) plane. We can see the characteristic diffraction pattern. Comparison with the analytical results in Fig. 3.6 shows good agreement. The main differences are: (i) the amplitude value of the low pressure just after the heterogeneity in the defocusing zone, (ii) the lack of oscillations before the heterogeneity in our result. The difference of amplitude in the defocusing area is dependant on the number of points used and full convergence is hard to achieve. This is probably due to the difficulty of meshing properly a sphere with a cartesian grid. Note that the position of the minimum is well reproduced. The absence of oscillations is due to the one-way approach, as the backscattered wave is not taken into account. This hypothesis is acceptable if the heterogeneity is small (less than 5 %). An other result (not shown here) is that taking into account or omitting the coupling term of the algorithm (operator H 2 in 2.2.3, solved by finite differences) has no influence in this case. This indicates that phase effects are dominant to explain the diffraction pattern. However, the difference in computational time between the cases with or without operator H 2 is quite significant. FLHOWARD, Analytical. The grey zone marks the heterogeneity location.

Nonlinear propagation in a thermoviscous medium

To validate the nonlinear and absorption part, we use the 1D case of propagation of a pure tone acoustic wave of angular frequency ω 0 in a nonlinear thermoviscous medium. This problem has an analytical solution given by Mendousse [START_REF] Hamilton | Nonlinear acoustics[END_REF]:

P = 4 Γ ∞ n=0 Γ(-1) n+1 nI n (0.5Γ)exp(-n 2 αX)sin(nω 0 τ ) I n (0.5Γ) + 2(-1) n I n (0.5Γ)exp(-n 2 αX)cos(nω 0 τ ) (3.1)
with P = p a /P 0 the acoustical pressure normalized by the pressure amplitude at the source P 0 , X = σρ 0 c 0 βP 0 ω 0 the propagation distance, α = βP 0 ω 0 Γρ 0 c 3 0 the absorption coefficient and Γ = βP 0 ω 0 αρ 0 c 3 0 the Gold'berg number which measures the ratio between nonlinearity and absorption. Notation σ is for the dimensionless distance and σ = 1 corresponds to the shock distance in the inviscid case. Here, the parameters are chosen to be Γ = 50 and σ = 3 so that the final point is well beyond the shock distance and the shock wave is well formed with a characteristic saw-tooth shape. The large Gold'berg number implies that the primary effect is the nonlinearity which strongly dominates over the absorption. Numerically 512 points are used to discretize the waveform in retarded time. least-square error metric defined by:

error % = 100 (P F LHOW ARD (τ ) -P (τ )) 2 P (τ ) 2 (3.2)
is plotted against the number of discretization points in the axial direction N x used for propagation in Fig. 3.9. The error remains small, even if a small number of points is used for the propagation. As expected [START_REF] Loubeau | Effects of Meteorological Variability on Sonic Boom Propagation from Hypersonic Aircraft[END_REF], the second order split-step induces a great improvement compared to the first order split step: the error is nearly always under 1%, even for a very small number of points (as soon as N x is larger than 20, which is equivalent to say that an error less than 1% is achieved with 7 points per shock formation distance).

Flow Validation

The last part to validate is the propagation in a moving medium. Two cases are considered: the first one to validate the stratified shear flow represented by the operator F LH (1) and the second one for the velocity fluctuations of the operator F HL (2) . Assumptions used for these tests are in agreement with those used to derive the model. Only linear propagation in an inviscid medium is considered here, so that operators A and N are omitted.

Linear acoustic modal propagation in a two-dimensional waveguide with shear flow

The first case of validation of the flow effect deals only with a stratified shear flow. It is the linear case of propagation of an acoustic mode in a waveguide in presence of a shear flow as illustrated on Fig. 3.10. This validation as already been used for the 2D version of the code [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF]. There are two main advantages for this validation case. First, the modes can be compared to the solution of Lilley's equation [START_REF] Lilley | The Generation and Radiation of Supersonic Jet Noise[END_REF], which is an exact equation for the propagation through sheared flows. Second, modes are only characterized by their transverse profile and their phase speed. In this case the numerical error will cumulate along the propagation and thus this is a very demanding test case. The main problem of [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF] was the difficulty to adapt the angular spectrum to this test case. Only FFT were used in the transverse direction, and so only periodic boundary conditions could be imposed. Here, the use of the DCT overcomes the problem and allows us to model a realistic duct. It was shown [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF] that the profile of the sheared flow does not change much the numerical results, so we will restrict to a Poiseuille flow with M = +0.1. A reference solution is obtained by searching for a solution in a modal form p a (x, z, τ ) = Π(z)exp [(k x -k 0 )x -ω 0 τ ] for FLHOWARD equation neglecting nonlinearities, absorption and heterogeneities. For the flow effects, the operator F LH (1) is used. It results into:

           d 2 Π(z) dz 2 + 2k x dM/dz (1 + M 2 ) dΠ(z) dz + k 2 0 1 - k 2 x + 2M k x 1 + M 2 Π(z) = 0 dΠ dz z=-L = dΠ dz z=L = 0. (3.3)
For comparison, the modal form of Lilley's equation is given by:

         d 2 Π(z) dz 2 + 2k x dM/dz (1 -k x M ) dΠ(z) dz + k 2 0 (1 -k x M ) 2 -k 2 x Π(z) = 0 dΠ dz z=-L = dΠ dz z=L = 0. (3.4)
Eq. 3.3 is then solved using a shooting method to obtain the profiles of the modes Π(z), and the associated wavenumber k x at a given frequency ω 0 . This also gives the phase speed: c φ = ω 0 /k x . The duct half-height is chosen L = 1 m, the sound speed is c 0 = 340 m/s. We can remark that modes Π(z) are computed numerically in the line x = 0 but since they are modal solution, their propagation is given analytically for x > 0. Consequently, the reference solution here is semi-analytical. Fig. 3.11 displays the two selected test modes. Mode 1 is the first symmetric mode and mode 3 is the second antisymmetric one. Comparison with Lilley's exact modes can be found in [START_REF] Coulouvrat | New equations for nonlinear acoustics in a low Mach number and weakly heterogeneous atmosphere[END_REF][START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF][START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF] and is not reproduced here.

Each mode profile is then used as input condition in the plane x = 0 for the code. Modes are numerically propagated over a long distance: 1000λ in the x-direction. Numerical dissipation and dispersion are assessed using respectively the amplitude of the mode and its phase speed. The studied frequencies are f = [400, 600, 800, 1000] Hz. The mesh is chosen with eight points per wavelength in the propagation direction (N x = 8000), N z = 512 points in the transverse direction and N τ = 256 for the temporal signal. Fig. 3.12 shows the relative error on the amplitude after the propagation over one thousand wavelengths for the selected two modes. First, for high frequencies, the error tends to decrease. Near the cut-off, at low frequencies, the model gives erroneous results in amplitude but the profile is conserved. These results are consistent with the fact that a wide-angle approximation has been performed on the flow part of the equation. Therefore, low frequency modes near cut-off propagate quite perpendicularly to the guide axis and are poorly described due to the numerical errors as shown in the analysis of the dispersion relation performed in Sec. 2.2.7. The other interesting result to investigate is the importance of taking into account the operator H 2 and the corresponding finite differences part of the algorithm that is used to solve it. Without operator H 2 , error is really small (nearly always under 10%) which indicates that this operator is negligible as it was already shown for the heterogeneities in Sec. 3.2. Moreover, even using the conservative form recommended by Gallin et al. [START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF], the scheme is unstable when taking into account H 2 as a cumulative increase in amplitude can be seen. Concerning the dispersion error, Fig. 3.13 shows that FLHOWARD gives a very accurate phase speed with a relative error smaller than the percent even near the cut-off. Since the FLHOWARD algorithm will be used primarily for atmospheric propagation, this validation test case shows that it can be used with confidence. In the atmosphere there will be no hard walled duct, the ducting effect will be due to variation in the meteorological conditions which are continuous. Also, due to the variability of the atmosphere, the distances of ducting will be shorter and there will be no such cumulating of error effects. It is to be noted that this validation test case has been carried out before the numerical analysis of Sec. 2.2.7 was achieved. Therefore no validation test based on the simplified but more precise F LH (s) operator could be done.

Scattering of a plane wave by a finite-circulation vortex

The second case of flow validation is the propagation of a pure tone plane wave through a vortex shown by Fig. 3.14. This problem has been tackled using both analytical [START_REF] Shea | Sound scattering by a potential vortex[END_REF][START_REF] Ford | Scattering of acoustic waves by a vortex[END_REF][START_REF] Howe | On the scattering of sound by a rectilinear vortex[END_REF] or numerical methods such as parabolic equations [START_REF] Candel | Numerical solution of wave scattering problems in the parabolic approximation[END_REF] or DNS [START_REF] Colonius | The scattering of sound waves by a vortex: numerical simulations and analytical solutions[END_REF]. It has been chosen as validation case for the Wide-Angle Parabolic Equation of Dallois et al. [START_REF] Dallois | Wide angle parabolic equations in moving media: Sound diffraction by a core vortex[END_REF] and the FDTD model of Cheinet et al. [START_REF] Cheinet | Unified modeling of turbulence effects on sound propagation[END_REF][START_REF] Ehrhardt | Modélisation en domaine temporel de la propagation acoustique[END_REF]. The vortex has a finite-circulation, its tangential speed is given by:

v θ (r) = Γ 2πr 1 -exp -α r 2 L 2 (3.5)
where r is the distance to the center of the vortex, Γ the circulation and L the size of the vortex. α = 1.256431 is a constant chosen so that the maximum velocity occurs at r = L. The radial velocity is set to zero. This vortex is an exact solution of the incompressible Navier-Stokes equation, it is called the Lamb-Oseen vortex. In this work we will present only one configuration: Γ = 1511, M = 0.25 and L = 2 m. The wave is generated 40 m before the vortex, its frequency is f = 43 Hz. These parameters are chosen to be the same as Colonius [START_REF] Colonius | The scattering of sound waves by a vortex: numerical simulations and analytical solutions[END_REF] and Cheinet [START_REF] Cheinet | Unified modeling of turbulence effects on sound propagation[END_REF]. It is a difficult test case since the Mach number is relatively high while the derivation of the FLHOWARD equation assumed a small Mach number as is expected in outdoor acoustics. The computational domain extends from -40 m to 40m in the propagation direction (x), -90 m to 90m in the transverse direction (y) and 0 s to 0.2 s (one period) in time. The boundary conditions used in y are reflections since it was found that they do not interfere with the validation process. The domain is discretized using 2048 points in both x and y directions. 512 points are used for the temporal signal. The mesh is voluntarily oversampled to evaluate only the model. Two computations were performed. The first one tested the FLHOWARD model with operator F LH (2) including the effect of transverse flow H = H 1 + H2. Note that in this case F LH (2) = F LH (s) . The second omitted the coupling operator H 2 of the algorithm, which is the only one that includes the effect of the transverse flow relative to the main propagation direction. In this case we have H = H 1 . The results are presented using what Colonius called the root-meansquare (RMS) of the scattered field. It consists in subtracting the plane incident field from the results of the computations and then take its RMS value (for a sinusoidal signal corresponds to dividing its amplitude by √ 2). Fig. 3.15 presents the RMS scattered field resulting from FLHOWARD simulation. Two well known characteristics of this case can be observed. First, the effect of refraction due to the finite circulation of the vortex is clearly seen. The scattered field is not zero before entering the vortex.

Second, an asymmetry can be observed with an interference pattern due to the fact that propagation is in the flow direction for y < 0 and opposite for y > 0. This validates qualitatively the use of FLHOWARD model. Quantitatively, Figure 3.16 compares the two FLHOWARD computations (with or without the operator H 2 ) with Colonius DNS results [START_REF] Colonius | The scattering of sound waves by a vortex: numerical simulations and analytical solutions[END_REF]. The value of the field is extracted on a circle of radius 20 m centred on the vortex. Before going through the vortex, the refraction is perfectly taken into account.

The main maxima and minima are well reproduced even though there are differences mostly on the amplitude. The FLHOWARD simulations lead to slightly smoother fluctuations of the pressure field in a way similar to what was observed for scattering by a speed of sound heterogeneity. The main differences between FLHOWARD and DNS are mostly due to the wide-angle approximation on the correction terms that take into account the flow effects. As expected, the transverse flow described by H 2 has only a small effect, only slightly shifting the curves slightly closer to the DNS simulations. The levels of maxima and minima are nearly the same with or without transverse flow.

Performance

Two metrics are used to quantify the performance of the code. The first is the strong scaling. It is defined as how the restitution time varies with the number of processors for a fixed overall problem size. Ideally, the restitution time should be halved when doubling the number of cores. Strong scaling is a measure needed when the goal is to tackle a problem quicker when using more computational cores i.e. for cpu-bound computation. The second is weak scaling. It is defined as how the restitution time varies with the number of processors for a fixed problem size per processor. Ideally, the restitution time should remain constant when both the mesh number of points and the number of cores are doubled. It is of interest when a bigger mesh is needed, i.e. for memory-bound computations.

These performances tests are performed on the cluster of Institut Jean le Rond d'Alembert (UMR7190 UPMC CNRS) which has 29 nodes of four 6-cores (AMD Opteron 2435 processors) cadenced at 2.6 GHz for a total of 696 cores. Its theoretical peak processing power is 3 Tflops. Each node has 64 GB of memory resulting in a total of 1,8 TB of RAM for the entire cluster. The nodes are interconnected with Infiniband connections. It runs on a Linux CentOS operating system. The procedure used is the following. The simulation is the propagation of an N-wave of amplitude 100 Pa and duration 0.2 s through an acoustical lens resulting from a 3D sound speed heterogeneity:

c 0 = c0 1 -0.3 exp (x -x 0 ) 2 σ 2 x + (y -y 0 ) 2 σ 2 y + (z -z 0 ) 2 σ 2 z , (3.6)
with σ x = 50 m and σ y = σ z = 100 m. Such a test case has already been selected [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF] for mesh convergence studies at 2D. It is here chosen because it involves all the main operators D, N and H 1 of the algorithm. The computational domain extends from 0 m to 680 m in the propagation direction (x), and from -340 m to 340 m in the transverse directions (y, z). The time window is 1.5 s long. All the writing and reading operations on the hard-disk are deactivated so that only the efficiency of the computation is quantified in the test. The results are presented in times spent for performing one advancement step from x to x + ∆x. To compute this, hundred steps are performed and the computational time is averaged to reduce the errors of measurement.

Strong scaling: For the strong scaling, two configurations were tested. The first one has a domain (N τ × N y × N z ) of (1024 × 1024 × 1024) resulting in 1 billion points while the second one has 4 billions points with a domain of (2048 × 2048 × 1024). These two domain sizes were chosen because they are the targeted domain sizes for the application of the code. Computations were performed on up to 512 cores. Results presented on Fig. 3.17 show a good strong scaling with a 1.7 decrease in computation time when doubling the number of cores compared to the ideal linear case of factor 2. For these sizes of mesh, the computation cost is much higher the communications between cores even using MPI_ALL_TO_ALL calls. Let us recall that such a collective communication instruction is requested here only for transposition operations as explained in section 2. is not as good (factor 2.4) as for the domains of bigger size (factor 1.8). This is because, for a small number of points, communications between cores become more important in the restitution time.

For both weak and strong scaling, the results do not depend on the specific value of N τ , N y or N z but only on the overall number of points (N τ × N y × N z ). The parallelization strategy therefore seems to show good results for the considered applications which target domain sizes between 1 and 4 billions points distributed on around 256 cores. : 2 20 , 2 22 , : 2 24 points per core.

Conclusion

This chapter detailed the validation of the FLHOWARD3D code. In addition to the already performed test cases at 2D [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF][START_REF] Dagrau | Simulation de la propagation du bang sonique : de la CFD à l'acoustique non linéaire[END_REF][START_REF] Gallin | Caractérisation acoustique des éclairs d'orage[END_REF][START_REF] Gallin | One-way approximation for the simulation of weak shock wave propagation in atmospheric flows[END_REF] new complementary and quantitative test cases have been investigated covering all involved operators and boundary conditions. All validations turned out to be satisfactory even in cases were the FL-HOWARD equation is beyond its theoretical range of validity. The parallelization performance tests showed also satisfying behaviour for both weak and strong scaling. This allows us to target numerical domains of the order of a few billion points necessary for contemplated sonic boom applications. Moreover, the theoretical analysis of numerical dispersion summarised by Fig. 2.5 proves that the split-step method leads to less error for the low order equation with operator H (s) than for the full equation with operator H. It is more precise numerically to take into account only the linear effect of flow motion and neglect quadratic convection terms. The numerical test cases with the waveguide show that there is little benefit to take into account the coupling terms described by operators H 2 or H (s) 2 which especially include the effect of flow gradient. In some cases, omitting this coupling operator even significantly reduces the numerical error. The case of wave scattering by a vortex confirms these observations and also shows that the transverse components of the flow motion have little effect on the wave diffraction pattern. As a conclusion, we can say that the lowest order terms for the flow motion is the most precise from a numerical point of view. In the following, only operator H (s) 1 will be considered for heterogeneous and flow effects.

This chapter begins by proposing a model for the Planetary Boundary Layer mixing

Monin-Obukhov Similarity Theory (MOST) for the mean stratified atmosphere (wind and temperature), along with homogeneous isotropic turbulent fluctuations satisfying a von Kármán energy spectrum. Sonic boom is then evaluated first in the standard atmosphere with no turbulence as the reference case. This is achieved by matching near-field CFD simulations based on automatic mesh refinement to nonlinear ray tracing describing long range propagation. Sonic boom of ATLLAS II Mach 5 configuration is quantified both ground track and laterally at the edges of the geometrical carpet in terms pressure signature, peak overpressure and annoyance metrics. Statistical effects 55 of turbulent wind fluctuations are then investigated for ground track boom by means of FLHOWARD3D software for two intensities of turbulence. Typical distorted waveforms and spectra are presented and peak overpressure variability is analysed. Finally, the sonic boom penetration in the shadow zone is reproduced numerically comparing the cases with and without turbulence. The Planetary Boundary Layer (PBL) or atmospheric boundary layer is the lower part of the atmosphere. Its structure depends on the hour during the day, and on the meteorological conditions. It is also heavily influenced by its contact with the ground. There exist two major states for the PBL as can be seen on Fig. 4.1 from Kaimal [START_REF] Kaimal | Atmospheric boundary layer flows: their structure and measurement[END_REF]. During the night, the buoyancy flux at the surface is negative and damps the turbulence so that the stratification is stable. It is solely driven by the shear turbulence produced by the geostrophic wind of the free atmosphere. During the day, the atmospheric boundary layer is frequently unstable due to the positive buoyancy flux at the surface which generates thermal instabilities. It is called the convective boundary layer and can be divided in two zones [START_REF] Stull | An introduction to boundary layer meteorology[END_REF]. Near the ground, the surface layer extends vertically to about 100-200m. There, turbulence is mostly due to shear stress. Above lies the mixed layer which extends from the surface layer to 0.5-2 km depending on the meteorological conditions. In the mixed layer, the turbulent diffusivity tends to be largest and mean gradients of wind and conserved scalars smallest [START_REF] Wyngaard | Turbulence in the Atmosphere[END_REF]. There, the temperature decreases almost linearly with the altitude.

The Planetary Boundary Layer

In this study, we deal only with the day case for which the turbulent level is the highest. For the propagation of sonic boom, the wind and temperature fields have to be known. Measured data cannot be used because they are usually provided only at points, and not over the entire field. So one has to model the PBL. The most correct way would be to simulate the PBL directly using the Navier-Stokes equations (DNS). However, this is not feasible due to the computational cost. A DNS computation needs a mesh of the order of Re 3/4 in each direction; Re = U L ν is the Reynolds number which is used to characterize the flow. For the PBL, we have a typical flow velocity U = 15 m.s -1 , the characteristic length is its height L = 1000m and the kinematic viscosity is ν = 1.4510 -5 m 2 .s -1 . This corresponds to a Reynolds number of 10 9 . This would lead to a mesh of the order of 10 20 points. An other method would be to perform large eddy simulation (LES) [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF], which consists in solving an approximation of Navier-Stokes equations where the small scales of turbulence are modelled instead of computed, thus reducing the size of the computational grid. This method is currently used in meteorology [START_REF] Moeng | A Large-Eddy-Simulation Model for the Study of Planetary Boundary-Layer Turbulence[END_REF][START_REF] Sullivan | A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows[END_REF][START_REF] Wyngaard | Turbulence in the Atmosphere[END_REF] and has recently been applied to electromagnetic wave propagation in a turbulent atmosphere [START_REF] Dipankar | A new phase-screen method for electromagnetic wave propagation in turbulent flows using large-eddy simulation[END_REF][START_REF] Dipankar | Trajectory of an optical vortex in atmospheric turbulence[END_REF]. However it is computationally expensive to have a fine enough mesh for the eddy sizes needed for the acoustic propagation. Indeed, from the acoustical point of view, the required numerical resolution is very demanding for shock waves. Several numerical tests of convergence for shock waves simulation can be found in literature. For shock wave propagation in turbulent medium, [START_REF] Averiyanov | Propagation des ondes acoustiques à travers un milieu turbulent: Etudes théorique et expériementale des effets de diffraction et des effets non linéaire[END_REF] recommends 25 points per wavelength in the transverse direction. For a typical Concorde sonic boom, this would lead to a resolution around 3 meters. Sonic boom focusing [START_REF] Marchiano | Numerical simulation of shock wave focusing at fold caustics, with application to sonic boom[END_REF] induces a characteristic scale of diffraction h Dif f = (2/(c 2 0 T 2 sig R caus )) -1/3 with T sig the signal duration (0.3s for Concorde) and R caus the radius of curvature of the caustic (in the range of 10 to 100 kilometres for a stratified atmosphere, probably less in a turbulent medium). This gives a value of around 300 meters, and with a required mesh of about a few 100 points over this distance, this leads again to a resolution of the order of the meter. Finally, from the perception point of view, the frequency of maximum sensitivity of human ear is about 3 kHz which corresponds to a wavelength about 10 cm. So for our contemplated applications, resolution of the order of the meter or less is targeted. Coupling between LES of the PBL and acoustic propagation has been performed in two-dimensions for linear acoustics [START_REF] Wilson | Characterization of uncertainty in outdoor sound propagation predictions[END_REF][START_REF] Wilson | Dependence of predictive skill for outdoor narrowband and broadband sound levels on the atmospheric representation[END_REF] to assess the efficiency of simpler atmospheric models with a resolution of the order of 1 meter again. However it cannot be used routinely. The simpler method practiced in this study is to model the mean atmosphere (mean temperature and shear wind in the propagation direction) by the well established Monin-Obukhov Similarity Theory (MOST) [START_REF] Stull | An introduction to boundary layer meteorology[END_REF][START_REF] Kaimal | Atmospheric boundary layer flows: their structure and measurement[END_REF][START_REF] Wyngaard | Turbulence in the Atmosphere[END_REF] as summarized in Section 4.2. Then a synthetic turbulence approach [START_REF] Frehlich | Simulation of Three-Dimensional Turbulent Velocity Fields[END_REF] will be detailed in Section 4.3 for generating the turbulent fluctuations. Evaluation of the sonic boom from ATLLAS II configuration in the non turbulent case will be the object of Section 4.4. It will be used as an input for quantifying the impact of turbulence on the undertrack boom in Section 4.5 and on the lateral boom in Section 4.6.

The mean atmosphere

Monin-Obukhov Similarity Theory is used to determine the vertical profiles of mean temperature and wind velocity within the atmospheric surface layer whose thickness is roughly 100-200 m in an unstable stratification. However, since the wind in the mixed layer is nearly constant and the temperature adiabatic [START_REF] Stull | An introduction to boundary layer meteorology[END_REF], it is reasonable to use MOST even for higher altitudes. MOST assumes that characterization of the surface layer needs only two parameters that can be determined from ground measurements: the friction velocity u * and the surface sensible heat flux Q s . Thus, the wind profile V 0 (z) and temperature profile T (z) can be determined:

T (z) = T r -(z -z r )Γ d + P t T * κ v [ln(z/z r ) -Ψ h (z/L 0 ) + Ψ h (z r /L 0 )] (4.1) V 0 (z) = u * κ v [ln(z/z 0 ) -Ψ m (z/L 0 ) + Ψ h (z 0 /L 0 )] . ( 4.2)
Here, z r is a reference height, T r = T (z r ) is the value of temperature at this height, Γ d = 0.0098 K / m is the dry adiabatic lapse rate, P t = 0.95 is the turbulent Prandtl number, T * is a temperature scale, κ v = 0.40 is the von Kármán constant. The Obukhov length is L 0 = -u 3 * T r ρ 0 /(κ v Γ d Q s ), and Ψ h,m (ξ) are universal function profiles given by [START_REF] Wilson | An alternative function for the wind and temperature gradients in unstable surface layers[END_REF]:

Ψ h,m (ξ) =    2 ln 1+ √ 1+a h,m |ξ| 2/3 2 , if ξ < 0 -b h,m ξ, if ξ ≥ 0 (4.3)
where a h = 7.9, a m = 3. The value of z 0 was chosen to represent a flat ground with short grass. This example was chosen according to Ostashev et al. [START_REF] Ostashev | Source localization from an elevated acoustic sensor array in a refractive atmosphere[END_REF].

Turbulent fluctuations

Characteristics of homogeneous and isotropic turbulence

As already explained in section 4.1, the best way to describe the turbulence would be by Direct Numerical Simulation or Large Eddy Simulation but it is not possible due to the computational resources needed as explained in the introduction. However, turbulence can be described statistically at least under the assumption of homogeneity and isotropy. These assumptions are not realistic from a meteorological point of view because of the atmosphere vertical stratification and shear wind. Some works have been devoted to estimate the influence of anisotropy or intermittency on linear acoustics [START_REF] Wilson | Calculated coherence and extinction of sound waves propagating through anisotropic, shear-induced turbulent velocity fluctuations[END_REF][START_REF] Wilson | Sound-wave coherence in atmospheric turbulence with intrinsic and global intermittency[END_REF] but none on nonlinear acoustics.

A vector velocity fluctuation field whose ensemble average is null: u 0i (r) = 0 is considered where u 0i (r) is the ith component of the fluctuation velocity vector u 0 , is the ensemble average and r is the position vector. The statistics of this field are described by the covariance tensor which is determined by the spatial correlations:

B ij (r 1 , r 2 ) = u 0i (r 1 )u 0j (r 2 ) (4.4)
where r 1 and r 2 are two points of the space. In case the turbulence is homogeneous, it depends only of the separation between the points s = r 2r 1 :

B ij (s) = u 0i (r)u 0j (r + s) . (4.5)
If u 0 is an incompressible homogeneous isotropic vector field, B ii (0) = σ 2 u is the variance of the velocity component u 0i . Its total variance is u 2 0 = 3 u 2 01 = 3σ 2 u . The velocity spectral density is given by the Fourier transform of the correlations:

φ ij (k) = 1 8π 3 R 3 B ij (s) exp (-ik • s)ds (4.6) B ij (s) = R 3 φ ij (k) exp (ik • s)dk. (4.7)
Here notation k is used only in this section for the turbulent wave vector. It should not be mistaken with the acoustical wave vector used in the rest of the document. For a scalar field n, the same definitions apply and we have the correlation:

B(r 1 , r 2 ) = n(r 1 )n(r 2 ) , ( 4.8)
which reduces to: B(s) = n(r)n(r + s) , (4.9)

if n is homogeneous and to

B(s) = n(r)n(r + s) , (4.10)
if n is also isotropic (only the distance between the points matter). The variance is obtained by B(0) = σ 2 n . Its spectral density is given by the Fourier transform of the correlations:

φ(k) = 1 8π 3 R 3 B(s) exp (-ik • s)ds (4.11) B(s) = R 3 φ(k) exp (ik • s)dk. (4.12)
One way to handle some of the differences between a vector and a scalar field is to introduce their energy spectra E(k) which is given by:

φ ij (k) = E u (k) 4πk 4 (δ ij k 2 -k i k j ) (4.13)
for a non-divergent vector velocity field and:

φ(k) = E n (k) 2πk 2 (4.14)
for the scalar field. The energy spectrum is also easier to interpret than the spectral density since it relates the intensity of the vortex to the wavenumber k and thus to its size 2π/k. It is to be noted that in this case, the vector field is isotropic but each of its component is anisotropic. The turbulent field energy depends only on the wave number k (assumption of isotropy) whereas its spectral densities depend on the wave vector k. The scalar and vector fields can follow the same energy spectrum E(k) [START_REF] Wilson | Turbulence models and the synthesis of random fields for acoustic wave propagation calculations[END_REF]. But since σ 2 is taken as the variance of each component of the vector field we have

E u (k) = 3E(k) and E n (k) = E(k).
Different models for the energy spectrum are proposed in the litterature [START_REF] Hinze | Turbulence[END_REF][START_REF] Ostashev | Propagation of sound in a turbulent medium. I. Plane waves[END_REF][START_REF] Tennekes | A First Course in Turbulence[END_REF]. Here the three most used ones are presented. First the Gaussian spectrum:

E(k) = σ 2 k 4 L 5 0 24 √ π exp - k 2 L 2 0 4 ; (4.15)
then the classical Kolmogorov spectrum:

E(k) = 4Γ(17/6) 3 √ πΓ(1/3) σL 0 (kL 0 ) -5/3 ; (4.16)
and the more complete von Kármán one:

E(k) = 4Γ(17/6) 3 √ πΓ(1/3) σ 2 k 4 L 5 0 (1 + k 2 L 2 0 ) 17/6 .
(4.17)

In the preceding equations, Γ is Euler's gamma function. The parameters of these spectra are the variance σ 2 of the field and the characteristic length L 0 . Kolmogorov and Gaussian spectra have often be used in several acoustic studies [START_REF] Wilson | Acoustic scattering and the spectrum of atmospheric turbulence[END_REF][START_REF] Ostashev | Acoustics in Moving Inhomogeneous Media[END_REF][START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF] but they are not fully consistent with some of the features of turbulence. First, it is known that in the inertial range of the spectrum, the energy of the vortices is decreasing as k -5/3 [START_REF] Hinze | Turbulence[END_REF][START_REF] Ostashev | Propagation of sound in a turbulent medium. I. Plane waves[END_REF][START_REF] Tennekes | A First Course in Turbulence[END_REF]. This decay is well reproduced by both von Kármán and Kolmogorov spectra but not by the Gaussian one for which the decay is exponential. Also, for atmospheric turbulence, there exits a length scale L 0 which is the limit of the biggest vortices. This length scale is present in both the Gaussian and von Kármán spectra which can therefore represent the large scale inhomogeneities. The Kolmogorov spectrum does not reproduce this energy-containing regions and thus is not accurate for the large scale inhomogeneities. All these features can be seen on Fig. 4.4. Thus the von Kármán spectrum can be seen as more representative than the two others and will be used in this work.

Generation of turbulence

Taylor's hypothesis of frozen turbulence will be used since the characteristic time of acoustic propagation is much smaller than the characteristic time of turbulent changes. Acoustical waves are propagated through multiple realizations of turbulence.

To obtain a random field of turbulence (either temperature or velocity) using an energy spectrum, several methods exist: the random Fourier modes of Kraichnan [START_REF] Kraichnan | Diffusion by a Random Velocity Field[END_REF][START_REF] Chevret | A numerical model for sound propagation through a turbulent atmosphere near the ground[END_REF][START_REF] Karweit | Simulation of the propagation of an acoustic wave through a turbulent velocity field: A study of phase variance[END_REF][START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF], the random phase generation [START_REF] Gilbert | Calculation of turbulence effects in an upward-refracting atmosphere[END_REF][START_REF] Wilson | Acoustic scattering and the spectrum of atmospheric turbulence[END_REF], the turbules method [START_REF] De Wolf | A random-motion model of fluctuations in a nearly transparent medium[END_REF][START_REF] Mcbride | Scattering of sound by atmospheric turbulence: Predictions in a refractive shadow zone[END_REF] and the recent quasi-wavelet method [START_REF] Goedecke | Quasi-wavelet model of von Kármán spectrum of turbulent velocity fluctuations[END_REF][START_REF] Goedecke | Quasi-Wavelet Models of Turbulent Temperature Fluctuations[END_REF][START_REF] Wilson | Quasi-wavelet formulations of turbulence with intermittency and correlated field properties[END_REF][START_REF] Wilson | Sound-wave coherence in atmospheric turbulence with intrinsic and global intermittency[END_REF][START_REF] Wilson | Quasi-wavelet formulations of turbulence and other random fields with correlated properties[END_REF]. In this work, the random field generation method [START_REF] Frehlich | Simulation of Three-Dimensional Turbulent Velocity Fields[END_REF] is used to generate realizations of a homogeneous and isotropic turbulent field. It was chosen for its ease of implementation and its computational efficiency. According to Frehlich, [START_REF] Frehlich | Simulation of Three-Dimensional Turbulent Velocity Fields[END_REF], a divergence free, vectorial field can be expressed as:

u 0j (x, y, z) = Nx mx=0 Ny my=0 Nz mz=0 w j (m x , m y , m z ) exp 2πi xm x N x + ym y N y + zm z N z .
(4.18) Here, k j = 2πm j /N j are the components of the discrete wave vector, and N j are the number of discretization points of the turbulent field in each direction j = x, y, z. The complex random vector w is generated by:

   w 1 w 2 w 3    =    H 11 0 0 H 12 H 22 0 H 13 H 23 H 33    •    N 1 N 2 N 3    (4.19)
with N a vector of arrays given by:

N j (m x , m y , m z ) = a j (m x , m y , m z ) + ib j (m x , m y , m z ) a j (m x , m y , m z ) 2 = b j (m x , m y , m z ) 2 = 1 a j (m x , m y , m z )b j (m x , m y , m z ) = 0. (4.20)
Here, a j (m x , m y , m z ) and b j (m x , m y , m z ) are uncorrelated arrays of zero mean Gaussian random numbers. The coefficients of the matrix in 4.19 are:

H 11 (m x , m y , m z ) = φ 11 (k x , k y , k z )
H 13 (m x , m y , m z ) = φ 13 (k x , k y , k z ) ∆k x ∆k y ∆k z φ 11 (k x , k y , k z ) (4.24) H 13 (m x , m y , m z ) = φ 23 (k x , k y , k z )∆k x ∆k y ∆k z -H 12 (m x , m y , m z )H 13 (m x , m y , m z ) φ 22 (k x , k y , k z ) (4.25) H 33 (m x , m y , m z ) = φ 33 (k x , k y , k z )∆k x ∆k y ∆k z -H 2 13 (m x , m y , m z )H 2 23 (m x , m y , m z ) (4.
26) where ∆k j = 2π/(δ j N j ) with δ j the grid spacing in the j-direction and j = x, y, z. We recall here, that quantities φ ij are the spectral densities of the turbulent field define by Eq. 4.6. In this method, variance of the velocity amplitude σ u and characteristic length of the vortices L 0 are the only two input parameters needed to generate a homogeneous isotropic field of turbulence. The method was presented for the generation of a random vector field but can also be applied to a scalar field by using only one component in the above equations.

An example of realization of both a von Kármán and a Gaussian spectrum is shown on Fig. 4.5. The smooth aspect of the Gaussian spectrum realization is due to the exponential decay of small vortices. For acoustic applications, this difference is important because smaller vortices will influence the higher frequency part of the wave spectrum whose acoustical wavelengths are comparable to the vortices size.

Sonic boom in the standard atmosphere 4.4.1 Geometry and CFD computation

The vehicle studied is the ATLLAS2 hypersonic configuration shown on Fig. 4.6. It is made of a main fuselage, a delta wing. Three high bypass turbofans are integrated either in the fuselage or in two nacelles. They are equipped with variable geometry mixed intakes. A CFD computation of this configuration was performed by ONERA [START_REF] Luquet | Sonic Boom Assessment of a Hypersonic Transport Vehicle with Advanced Numerical Methods[END_REF] using their in-house multi-physics polyhedral platform CEDRE [START_REF] Chevalier | CEDRE : Development and validation of a multiphysic computational software[END_REF] and its CFD solver CHARME. The plane is cruising at a Mach number M = 5 and an altitude of 26 km. The flow is assumed inviscid and the Euler equations are solved using an AUSM+ second order flux scheme. In order to compute accurately the shock in the mid-field, an automatic mesh adaptation was used. The Feflo.a software developed by INRIA [START_REF] Loseille | Anisotropic Mesh Generation Application to High-Fidelity Simulation in CFD[END_REF] was used to determine the mesh requirement. The computed pressure signal is then extracted on a cylinder around the vehicle as presented on Fig. 4.7 along with the final adapted mesh. Cylinders of extraction are of radius comparable to the aircraft length. On the nearest extraction from the body, the flow is shown to exhibit a strong 3D behavior resulting from the complexity of the configuration: azimuthal variations of the pressure field can be very abrupt at 30 meters. This justifies the necessity to use the CFD and automatic mesh adaptation. On Fig. 4.8, the impact of the mesh adaptation can be observed on the dimensionless field √ r(p -p 0 )/p 0 taking into account the expected cylindrical decay of sonic boom as 1/ √ r with r the distance from the flight path. Without it, the signal is completely damped before reaching the mid-field and cannot be used as an input for the acoustic propagation phase down to the ground. For the present study, only the most distant extraction performed at 100 meters will be used.

Nonlinear ray tracing method

ATLLAS II sonic boom is evaluated at the ground level in the geometrical carpet for the standard atmosphere [START_REF]ICAO. Manual of the ICAO Standard Atmosphere (extended to 80 kilometres (262500 feet)[END_REF] with no wind and a perfectly reflecting ground. Attenuation relies on the example of humidity vertical profile proposed by ISO [START_REF] Iso | Acoustics -Attenuation of sound during propagation outdoors -Part 1: Calculation of the absorption of sound by the atmosphere[END_REF]. As already explained, the matching between CFD and ray tracing is done at R = 100 m. No multipole matching process [START_REF] Page | An efficient method for incorporating computational fluid dynamics into sonic boom prediction[END_REF] is applied. Due to the temperature decrease with the altitude, the sonic boom is refracted upwardly and the carpet has a well defined lateral edge. The sonic boom prediction method is based on a standard nonlinear ray tracing method [START_REF] Hayes | Sonic boom propagation in a stratified atmosphere with computer program[END_REF][START_REF] Thomas | Extrapolation of sonic boom pressure signatures by the waveform parameter method[END_REF][START_REF] Plotkin | State of the art of sonic boom modeling[END_REF] in a stratified, moving atmosphere [START_REF] Blokhintzev | The Propagation of Sound in an Inhomogeneous and Moving Medium I[END_REF][START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF]. The algorithm is described in Blumrich et al. [START_REF] Blumrich | Meteorologically induced variability of sonic-boom characteristics of supersonic aircraft in cruising flight[END_REF]. To avoid singularities near cut-off when rays are grazing, the position of a point along a given ray is parametrized by the eikonal function rather than by its altitude. To calculate the geometrical ray tube area On top, without mesh adaptation, on bottom with mesh adaptation (with courtesy of ONERA).

along the rays in a stratified medium, the implementation follows the method of [START_REF] Candel | Numerical solution of conservation equations arising in linear wave theory: application to aeroacoustics[END_REF] and requires the numerical integration of 13 differential equations: 4 for ray tracing, 8 for determination of the ray-tube cross-section (4 for the derivative of ray equations with respect to each of the two ray coordinates parametrizing a single ray, e.g., time of emission and azimuthal angle) and 1 for the nonlinear age variable. Integration is numerically achieved by using standard Runge-Kutta algorithm. Along a given ray, the pressure field satisfies a Generalized Burgers' Equation (GBE) [START_REF] Cleveland | Time-domain modeling of finite-amplitude sound in relaxing fluids[END_REF][START_REF] Cleveland | Comparison of computer codes for the propagation of sonic boom waveforms through isothermal atmospheres[END_REF] including nonlinearities, atmospheric absorption and geometrical spreading. Nonlinear effects are deemed essential in the long-range propagation of finite-amplitude waves such as sonic booms [START_REF] Whitham | On the propagation of weak shock waves[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. They are responsible for the slow evolution of the waveform until the typical "N" shape is achieved which is frequently recorded at ground level. After an algebraic transformation through the introduction of the age variable, geometrical spreading is eliminated. The resulting equation is solved with the same method as used here (see section 2), combining Burgers-Hayes method for nonlinear effects (section 2.2.4), analytic solution in the frequency domain for absorption (section 2.2.5) and second order split-step method (section 2.2.1). The algorithm is identical to the one developed in the first ATLLAS project [START_REF] Loubeau | Effects of Meteorological Variability on Sonic Boom Propagation from Hypersonic Aircraft[END_REF] and has been validated in section 3.3.

Sonic boom at ground level

To compute the ground sonic boom, the pressure signature is discretized with 2 14 points over a time window of duration 0.95 s. This leads to a maximum computed frequency equal to 8600 Hz, well beyond the frequency of maximum sensitivity of human ear.

Fig. 4.9 presents the lateral extension of the carpet which is 136 km wide. Due to the size of the vehicle and the high altitude, the duration of the signal shown on Fig. 4.10 is longer than for existing aircraft (Concorde, military fighters, Shaped Sonic Boom Demonstrator) and thus has a lower frequency content around 1.5 Hz. The maximum peak overpressure for under track boom is 70.7 Pa which is slightly smaller than sonic boom recorded for Concorde (around 100 Pa). The signal waveform is very similar to a classical N-wave shape for the compression phase. However the expansion phase is not symmetric due to lift effect which is known to induce a signal with non-zero mean value. The near-field signal results from the multiple shocks emanating from the intakes, the wing edges, the exhausts and the tail. These different shocks progressively coalesce away from the aircraft. At ground level, the shocks are slightly smeared out by the atmospheric absorption as can be seen on zoom of Fig. 4.10. It displays the characteristic shape of a shock wave in a relaxing medium. The frequency spectrum shows the typical boom shape [START_REF] Shepherd | A loudness calculation procedure applied to shaped sonic booms[END_REF] with a 6 dB decay per octave until 2000 Hz and then a 12 dB decay per octave. At the geometrical cut-off, the sonic boom is of much smaller amplitude 18 Pa because of the reduced influence of lift at the source, as can be seen on Fig. 4.8. Also due to this smaller lift effect, there remain two shocks in the compression phase, the second one of very small amplitude. As for the undertrack boom the shocks are slightly smeared out by atmospheric absorption which is not fully negligible even at such low frequency because of the very long propagation path and weaker nonlinearities. The frequency spectrum shows the fast decay of 12 dB per octave occurs at much lower frequency around 100 Hz. Beyond the cut-off, the creeping wave theory leads to an exponential decay and a fast smoothing of the shock waves as explained in [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF][START_REF] Coulouvrat | Sonic boom in the shadow zone: A geometrical theory of diffraction[END_REF].

Sonic boom metrics

One important feature when analyzing sonic boom, is the quantification of the annoyance caused to the people whether they perceive the boom outdoor or indoor. This annoyance is produce by the audible perception of sonic boom itself. However most of reported annoyance from sonic boom is due to startle, indirect induced rattle noise, vibrations of the structures and possible damages to buildings [START_REF] Maglieri | Sonic Boom[END_REF][START_REF] Coulouvrat | The challenges of defining an acceptable sonic boom overland[END_REF][START_REF] Maglieri | Sonic Boom Six Decades of Research[END_REF]. In this section, two metrics which are loudness based will be investigated: the A-weighted Sound Exposure Level (ASEL) and the C-weighted Sound Exposure Level (CSEL). They are classical metrics used for sonic boom impact estimation. An other common metric is the Perceived Loudness deciBel (PLdB). It was not chosen because it has shown a strong correlation with ASEL while being much more complex to compute. At the present time, there is no agreement of the scientific community on which metric is the most suitable and on any level of acceptability. A-weighting is the most common weighting. It simulates the human ear frequency response to noise. A good correlation with annoyance is indicated from community surveys for military fighters [START_REF] Leatherwood | Summary of recent NASA studies of human response to sonic booms[END_REF] and from sonic boom simulators for low-boom designed business jets [START_REF] Loubeau | A New Evaluation of Noise Metrics for Sonic Booms Using Existing Data[END_REF]. It puts emphasis on the boom high frequencies (above 200 Hz), which correspond to the audible part of the signal due to the shock wave and the associated startle effect. The C-weighting has been recommended by the US Committee on Hearing, Bioacoustics and Biomechanics (CHABA) [START_REF] Chaba | Community Response to High-Energy Impulsive Sounds: An Assessment of the Field Since[END_REF] study for high energy impulsive sound. It is putting weight on lower frequencies than the A-weighting in order to better take into account annoyance related to vibrations. These metrics are based both on a frequency weighting of the power spectrum. Computation of the Sound Exposure Level (SEL) relies on Parseval theorem for the pressure field energy:

+∞ -∞ p a (t) 2 dt = 2 +∞ 0 | pa (f )| 2 df (4.27)
where pa (f ) is the Fourier transform of the signal and f is the frequency. To take into account the frequency dependant response, a weighting H(f ) is applied to the Fourier transform of the signal. The weighted power spectrum is then transformed into a dB unit so as to get:

HSEL = 10 log 2
T ref +∞ 0 | pa (f )H(f )| 2 df p 2 ref (4.28)
where T ref = 1 s and p ref = 2 10 -5 Pa. A and C weighting are given respectively by [START_REF]IEC. Electroacoustics -sound level meters[END_REF]:

H ASEL (f ) = 12200 2 f 4 10 2/20 (f 2 + 20.6 2 )(f 2 + 12200 2 ) (f 2 + 107.7 2 )(f 2 + 737.9 2 ) (4.29)
and: According to the present knowledge [START_REF] Maglieri | Sonic Boom Six Decades of Research[END_REF] based on existing supersonic aircraft, Nlike booms with amplitudes in the range 50 Pa to 150 Pa and with a duration from 100 ms to 300 ms are deemed highly annoying by a percentage from 13% to 33% of the population. Even though ATLLAS II configuration leads to a boom slightly longer (around 400 ms), so with a lower frequency content, we expect that the percentage of population highly annoyed by ATLLAS II configuration would be in a similar range given the similar peak overpressure at least over the central part of the carpet.

H CSEL (f ) = 12200 2 f 2 10 0.06/20 (f 2 + 20.6 2 )(f 2 + 12200 2 ) . ( 4

Effect of a turbulent atmosphere on the undertrack boom 4.5.1 Computational process

In this section, the influence of the turbulence encountered in the planetary boundary layer as modeled in Sec. 4.1 is investigated. The computation process is sketched on Fig. 4.11. The ground signal from the ray tracing computation (see previous section) is used as an input at the top of the planetary boundary layer. This waveform is then numerically propagated with the FLHOWARD3D software in the vertical direction down to the ground. Several assumptions are made in this process to allow for an easy matching with ray tracing. First, geometrical spreading is neglected as the input is taken as a plane wave. Second, refraction is not taken into account, a vertical straight line main path of propagation being chosen. These assumptions are of little importance in the present study as they are applied for the undertrack boom only over a small part of the overall propagation path from the aircraft to the ground. For propagation over the PBL thickness, geometrical spreading and refraction are small compared to turbulence effects. Moreover geometrical spreading and refraction have already been taken into account in the ray tracing phase down to the ground but in a way uncoupled to the turbulence influence. However, as explained in Appendix F, even for this low amplitude signals, nonlinearity as to be taken into account. The meteorological conditions are: a mean sound speed c0 = 340 m.s -1 , a density ρ 0 = 1.2 kg.m -3 and a nonlinear coefficient β = 1.2. The turbulent velocity fluctuations are considered frozen and are generated as explained in Sec. 4.3. The outer scale of the vortices L 0 is 100 meters and two kinematic turbulence intensities following a von Kármán spectrum are chosen: σ v = 0.5m.s -1 and σ v = 2m.s -1 . Thermal turbulence is not considered here: c ′ 0 = 0. The main direction of propagation is along the x vertical direction. The propagation distance is chosen as 1500 meters which is characteristic for a fully developed convective boundary layer. Both the transverse directions y and z are 1500 meters long in order to contain large structures of turbulence. Each spatial direction uses 1024 discretization points (see Appendix E for a discussion on the mesh convergence). Periodic boundary conditions are assumed laterally in both y and z directions. This is coherent with the periodicity requested to synthesize turbulence. For the time variable τ , the numerical window is [-0.6, 1.6] seconds. It is chosen (i) in order to be long enough to capture all the scattered signals arriving after the ballistic wave and (ii) so that these scattered waves do not interfere with the main signal through the artificial periodicity imposed by the numerical method. Due to limited numerical resources, the temporal window involves only a maximum number of 4096 points. This leads to an overall 4 billions degrees of freedom per propagation step in the vertical direction. As a consequence, the time signals have a frequency cut-off around 900 Hz when considering the Nyquist limit. This is much lower than the 8600 Hz frequency cut-off of the ray tracing output. This one is extended over the time window by a zero padding and is then linearly interpolated over the new numerical grid. Consequently, a significant part of the audible frequency content cannot be simulated. In particular, the frequency cut-off is too low to compute precisely the loudness metrics. Computing the ASEL value for this too limited frequency spectrum for the input signal would give 81.94 dBA, 9 dB less than the true value. Achieving a sufficiently high frequency resolution would require ten times more points in the time direction which is intractable in terms of computational power. FLHOWARD3D is used with the simplified version

H = H (s)
1 which was found to be the best compromise between precision and efficiency in Sec. 3. This implies in particular that only the vertical component of velocity fluctuations is considered. The computations were carried out on 128 cores on the cluster described in Sec. 3.5. A full computation needs around 30 hours.

Pressure variability

Fig. 4.12 displays the maximum peak pressure in a vertical plane along the propagation direction from the top of the PBL (x = 0 m) to the ground (x = 1500 m). For the two turbulence intensities, a strong variability is obvious. For the strongest turbulence, overpressure doubling can be observed at some points. Variability increases the peak overpressure much more frequently and much more efficiently than it decreases it. The structure of the field is clearly anisotropic as strongly influenced by the propagation direction. This is not due to neglecting the horizontal components of the turbulent fluctuations [START_REF] Averiyanov | Propagation des ondes acoustiques à travers un milieu turbulent: Etudes théorique et expériementale des effets de diffraction et des effets non linéaire[END_REF]. This anisotropic variability is due to the wave focusing in regions of low effective sound speed (c 0 + v 0x ), and defocusing in regions of high effective sound speed [START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF]. Note however that the the notion of random caustic formation [START_REF] Blanc-Benon | On the appearance of caustics for plane sound-wave propagation in moving random media[END_REF] is here applicable mostly to the high frequency part of the spectrum as the acoustical main wavelength is comparable to the turbulence outer scale. This might reduce the probability of occurrence of very strong amplifications. In the horizontal directions, the scattering is clearly isotropic as illustrated by Fig. 4. [START_REF] Berry | Singularities in waves and rays[END_REF] showing the peak overpressure at altitude x = 750 m. All these features are also visible on Fig. 4.14 showing two transverse extractions in the two horizontal directions y = 0 and z = 0 at mid-altitude x = 750 m. The two extractions show similar variability, no statistical difference is visible between the two directions. The trend according to which turbulence mostly amplifies the peak overpressure is also undeniable. Large pressure fluctuations occur over a distance of 100 m. We hypothesize this scale is directly related to the turbulence outer scale L 0 . However, as the boom main wavelength is of the same order, it could also influence this variability scale.

Due to the effects of diffraction, nonlinearities and scattering by random inhomogeneities, waveforms change during the propagation. On Fig. 4.15 some signals are presented. They were chosen because they are typical for observations of sonic boom propagation through atmospheric turbulence [START_REF] Lee | Sonic Booms Produced by United States Air Force and United States Navy Aircraft: Measured Data[END_REF]. From top to bottom we have: (i) the classic N-wave which was used as input for the computation (in black); (ii) a spiked waveform characteristics for the passage of the wavefront through a caustic; (iii) a multiple spiked signal which has traveled through several caustics; and (iv) finally a rounded one with a longer rise time. The energy spectrum of these waves is also displayed on the same figure (in black for the input signal, in red for the other ones). The input spectrum is a classical N-wave spectrum with a cut-off due to the mesh and interpolation at around 600 Hz below the theoretical value at 900 Hz. At 600 Hz, it is about 65 dB less than at 1 Hz. Due to the longer time window the spectrum is smoother at low frequencies (below 10 Hz). On the spiked signal it can be seen that the mid-range frequencies, roughly between 10 Hz and 100 Hz corresponding to the shocks spikes are enhanced up to 5 dB corresponding to the strong peak pressure spike. These spikes are of a duration 0.01 s in agreement with the 100 Hz observed frequency. As a counterpart some energy dips are observed around 270 Hz and 570 Hz probably resulting from random destructive interference. For the multiple spiked waveform, some similar trends are observed with some frequency ranges enhanced (here between 40 Hz and 300 Hz) and some frequency ranges lessened (here between 10 Hz and 30 Hz). Finally, the rounded waveform has a significant decrease of energy above 30 Hz linked to the long rise time of the shocks around 0.015 s. These spectra are characteristic of the scattering of the shocks. Their energy is redistributed spatially increasing or decreasing the amplitude. The variance of the effective sound speed is σ v = 2 m.s -1 , less than 1% of the mean sound speed value and the size of the largest vortices is about L 0 = 100 m. As a consequence the main acoustical wavelength is of the order of 2L 0 and is almost unaffected by turbulence. For visible effects of the turbulence, it is necessary that the acoustical wavelength be small enough compared to the size L 0 (here about 30 Hz). This is due to the fact that the ATLLAS II configuration leads to a longer signal and thus to a lower spectral content (below 5 Hz) than usual sonic booms.

Statistical analysis

Statistics on the maximum peak pressure are now investigated. To obtain statistics, an ergodic hypothesis is used. Instead of performing many simulations for several turbulent field realizations, statistics are computed spatially across the horizontal (transverse) directions for a single realization. Fig. 4.16 and Fig. 4.17 presents the four first statistical moments for peak pressure along the vertical propagation axis for the turbulence levels. Average peak pressure shows an unexpected behaviour, first increasing before reaching a plateau and then beginning to decrease. This is different from observations of laboratory scale experiments for kinematic turbulence [START_REF] Lipkens | Model experiment to study sonic boom propagation through turbulence. Part III: Validation of sonic boom propagation models[END_REF][START_REF] Averiyanov | Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media[END_REF][START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF]. Note that similar experiments have recently been performed for thermal turbulence [START_REF] Salze | Laboratory-scale experiment to study nonlinear N-wave distortion by thermal turbulence[END_REF]. The main parameters of these ones are summarized in Table . 4.2. These experiments used spark sources to generate a high frequency shock wave propagating through a turbulent jet. The ratio of the turbulence outer scale L 0 to the thickness of the turbulent layer (H turb ), the turbulence intensity and the signal peak overpressure are comparable to typical sonic boom values. However, the experiment was constrained by the acoustical source so that the main acoustic wavelength λ = c 0 /f is relatively small compared to L 0 and H 0 . Moreover, high frequency and intensity implied a shock formation distance L shock = ρ 0 c 3 0 /(2πβf P 0 ) and a viscous absorption distance L viscous = ρ 0 c 3 0 /(2π 2 f 2 µ) that are around the meter so not so large compared to H turb . Consequently attenuation due to nonlinear and absorption effects is several orders of magnitude larger than for ATLLAS II sonic boom.

These different scaling ratios might explain the diverging pressure growth rates. Damping effects are negligible in the present computation case which explains to our opinion the initial growth of the average peak pressure. This one is likely to be due to the superposition of high frequencies scattered waves to the unaffected low frequency boom. For the low level of turbulence, the mean peak pressure initial growth is much slower than for the intense case. For both turbulence levels, the standard deviation shows a quick increase at the beginning of the propagation, before decreasing slowly. It is significant even for the low level of turbulence which supports the high variability already seen on the pressure map on Fig. 4.12. The peak of standard deviation seems to be proportional to the turbulent intensity: about 3 Pa for σ u = 0.5 m.s -1 and 12 Pa for σ u = 2 m.s -1 . The next two statistical moments are skewness and kurtosis. The 4.2: Comparison between characteristic scales for classical sonic boom, ATLLAS II configuration and laboratory scale experiments of Averiyanov et al. [START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF].

skewness measures the asymmetry of the distributions about its mean. In this case, it grows rapidly in the positive which means that the right tail will be longer than the left one indicating the higher probability of high peak pressures than low ones. This effect is shown to increase with the level of turbulence. The kurtosis represents the peakedness of the distribution. It will be positive if most of the standard deviation is due to infrequent extreme events compared to a Gaussian distribution. It has a similar trend as the skewness, increasing quickly and amplified by a high level of turbulence. The rapid increase in standard deviation, skewness and kurtosis is due to the rapid formation of caustics which widen the distribution and lead to an asymmetry with an excess in the high pressure tail. As refocusing of a wavefront that has already gone through a caustic is difficult due to phase scrambling [START_REF] Pierce | Statistical Theory of Atmospheric Turbulence Effects on Sonic-Boom Rise Times[END_REF], peak pressure begin to decrease slowly after achieving a maximum at a propagation distance which depends on the intensity of the turbulence. All these results can be seen clearly on the histograms of the probability density function (PDF) which are plotted on Fig. 4.18 for five distances of propagation. Moreover, an important result is that, as was already shown [START_REF] Blanc-Benon | Intensity fluctuations of spherical acoustic waves propagating through thermal turbulence[END_REF][START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF][START_REF] Salze | Laboratory-scale experiment to study nonlinear N-wave distortion by thermal turbulence[END_REF], the PDF can be fitted by a gamma distribution:

P DF (p a ) = ξ a p (a-1) a Γ(a) exp(-ξp a ) (4.31)
where ξ and a are two fit parameters that depend on propagation distance and turbulent intensity. The fact that even with completely different conditions, numerical simulations (present work along with [START_REF] Ehrhardt | Evaluating a linearized Euler equations model for strong turbulence effects on sound propagation[END_REF] for linear acoustics) and experimental data [START_REF] Blanc-Benon | Intensity fluctuations of spherical acoustic waves propagating through thermal turbulence[END_REF][START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF][START_REF] Salze | Laboratory-scale experiment to study nonlinear N-wave distortion by thermal turbulence[END_REF] follow the gamma distribution strongly supports its universality for the propagation of acoustic waves (linear or nonlinear) in random media.

Rise time

This section investigates the behaviour of rise time θ 0.1-0.9 for pressure waveforms. For a simulated waveform p a (τ ) at a given position (x, y, z), rise time is defined as the time the pressure needs to increase from 10 % to 90% of the peak overpressure. Note that with the chosen time resolution the rise time is not fully converged (in particular rise times equal or shorter than the mesh size are set to 0) but the results are nonetheless interesting. The rise time in a vertical plane (x, y, z 0 = 0m) is presented on top of Fig. 4.19 for the large turbulent intensity σ v = 2 m.s -1 . A high variability is observed. As for the peak pressure it is clearly anisotropic and inhomogeneous along the propagation direction. The main difference with the behaviour of the peak overpressure is that the rise time spatial variations are less smooth. It can vary strongly from one point to the next. The bottom of Fig. 4.19 representing rise time in a horizontal plane (x 0 = 750m, y, z) reinforces this observation of large variations over short distances. In these transverse directions y, z, the rise time variations appear isotropic and homogeneous, similarly to peak overpressure variations. Average and standard deviation (statistics are computed over the 2 transverse y, z directions) of the rise time along the propagation are plotted on Fig. 4.20. As expected, the average rise time increases as the signal propagate so does its variance. The comparison of the results for 2 turbulent intensities shows that the increase in average and variance of rise time is proportional to the turbulence intensity σ v . The size of the numerical domain does not enable us to say whether mean rise time reaches a plateau (as modelled by Plotkin and George [START_REF] Plotkin | Propagation of weak shock waves through turbulence[END_REF]) or not (as suggested by Pierce [START_REF] Pierce | Statistical Theory of Atmospheric Turbulence Effects on Sonic-Boom Rise Times[END_REF]). Histograms of the probability density functions of the rise time are presented on Fig. 4.21 at different propagation distances for the 2 turbulence intensities. They show that the distribution has an asymmetry with a large excess for the long rise time tail especially for the strong turbulence intensity. Along the propagation, there is an enlargement of the statistical distribution along with a decrease in the probability amplitude of a given value.

The relation between peak overpressure and rise time is also investigated. To do this, a scatter plot of the peak pressure in function of the rise time is displayed on Fig. 4.19 for 2 propagation distances and for largest turbulence intensity. First, the characteristics observed previously on histograms of rise time and peak overpressure are retrieved. Second, for a high peak overpressure, the rise time is low. This is likely to correspond to spiked waveforms. Large rise times occur only for small peak overpressure, which we expect to correspond to rounded waveforms. Nevertheless, short rise times can occur even for low amplitude signals. In case boom annoyance is linked to rise time, this could be an important feature as it might induce large variability in annoyance even for low amplitude N-waves. Also it is interesting to note that in his thesis, Salze [START_REF] Salze | Propagation acoustique non-linéaire en milieu inhomogène avec effets de sol: expériences à l'échelle du laboratoire[END_REF] presented similar scatter plots for the laboratory scale experimental propagation of shock waves in thermal turbulence.

Sonic boom in the shadow zone

The case of sonic boom penetration into the shadow zone at the edge of the primary carpet is now studied as illustrated by Fig. 4.23. In case of negative vertical temperature gradient or adverse wind, the sonic boom is refracted upwardly and the carpet has a finite width. In the shadow zone, the ray theory cannot predict the signal which is due to diffraction and scattering by turbulence. The case of an idealized shadow zone with no turbulence was studied by Coulouvrat [START_REF] Coulouvrat | Sonic boom in the shadow zone: A geometrical theory of diffraction[END_REF]. In the present work, the additional influence of turbulence is investigated. A plane wave is propagated in the upwardly refracting atmosphere presented in 4.1. This atmosphere is modeled by using Monin-Obukhov similarity presented in Sec. 4.2. The case of a mostly sunny Turbulence is assumed homogeneous and isotropic and obeys to a von Kármán energy spectrum. The characteristic scale is L 0 = 100 m and three levels of variance σ u = 1.2, 2.4 and 3.6 m.s -1 were tested. The temporal signal of the incoming wave is extracted from the ray tracing computation performed in Sec. 4.4 at the cut-off. Note that this waveform results from a simulation in a standard atmosphere. The matching is therefore not fully consistent as the mean vertical wind profile has not been taken into account for ray-tracing. This would require to match MOST profile with atmospheric data above the PBL which is beyond the scope of the present study. The main propagation axis is now horizontal. Propagation is simulated over 4 km, the height and width of the simulation domain are 4 km. In the vertical direction, an Absorbing Boundary Layer (ABL) of thickness 400 m is imposed on the top of domain. Rigid boundary conditions are imposed on all other boundaries. The duration of the time window is 2.2 s. The temporal signal, the altitude and the transverse direction are sampled with 1024 points each. As for undertrack boom, the input signal results from zero padding and interpolation of the ray-tracing. The maximum frequency of time mesh is therefore 233 Hz. Higher frequency sampling could not be achieved because: (i) the ground rigid boundary condition makes simulation about twice as long as the periodic case;

(ii) simulation of multiple realizations is required to achieve statistical convergence. The computational domain has more than 1 billions degrees of freedom; 1024 points are used in the propagation direction. On Fig. 4.24 we can observe the progressive decay of the ground positive peak pressure in the central vertical plane (y = 0). On the top figure, the atmosphere is stratified but not turbulent and the shadow zone can be clearly seen. Just above the shadow zone, a caustic is formed due to the curvature of the initially plane wavefront induced by atmospheric refraction. On the bottom figure, for σ v = 3.6m.s -1 , the turbulence is shown to create multiple random variability including focusing and defocusing areas in a way similar to the undertrack boom. This process of random scattering redirects energy into the shadow zone which turns out to be less pronounced. As a counterpart, the geometrical caustic is strongly perturbed.

As can be seen on Fig. 4.25, the positive peak pressure on the ground surface is still globally decaying with the distance but due to turbulence, a high variability with areas of enhanced amplitudes is nevertheless observed relatively far in the shadow zone. Fig. 4.26 presents two simulated temporal waveforms at x = 2 km, either in the shadow zone (z = 0) or near the caustic (z ≈ 250 m) in the cases with or without turbulence. The presence of oscillations before the first shock is due to the periodicity of the time window: the perturbations leaving the time window at the right side are re-injected on the left side. In the shadow zone, without turbulence, the initial N-wave profile progressively smears out along the ground due to exponential attenuation of creeping waves [START_REF] Berry | Controlled experiments of the diffraction of sound by a curved[END_REF][START_REF] Coulouvrat | Sonic boom in the shadow zone: A geometrical theory of diffraction[END_REF]. In the turbulent case, the positive peak pressure is increased and high frequencies associated to the shocks are more strongly scattered from the illuminated to the shadow zones than low frequencies [START_REF] Crow | Distortion of sonic bangs by atmospheric turbulence[END_REF]. This results into a very large increase of the wave spectrum with turbulence for frequencies above 5 Hz. We can therefore expect a large increase of the ASEL level of the boom in the shadow zone due to turbulence although we are not able to quantify this because of the too low frequency limit of the present numerical simulations. Near the caustic the typical U-wave observed without turbulence is strongly affected in the turbulent realization because of the random phase scrambling [START_REF] Pierce | Statistical Theory of Atmospheric Turbulence Effects on Sonic-Boom Rise Times[END_REF]. As a consequence, the positive peak overpressure is reduced. The spectrum modification is similar to what was observed for undertrack boom with changes mostly above 20 Hz.

Finally, to quantify the decrease of positive peak pressure qualitatively seen in Fig. 4.25, a statistical study is conducted. The positive peak pressure is averaged on the ground along the transverse direction. To get statistical convergence several realizations of turbulence are necessary, here up to 10 computations for each level of turbulence. Fig. 4.27 shows the average peak overpressure and its standard deviation compared to the case without turbulence. Both show the classical exponential decay. However, the average positive peak pressure increases with the turbulent intensity as expected because of random scattering from the illuminated zone to the shadow zone. But even for high turbulent intensities, this increase is relatively small: no more than 5 Pa. The standard deviation first increases over distances in the range from 500 m to 1000 m before slowly decaying. This behaviour is similar to what happens undertrack but the amplitude of the deviation is much weaker. Again, the more intense the turbulence, the more standard deviation increases.

Conclusion

An extensive sonic boom evaluation including the effects of turbulence in the atmospheric boundary layer has been performed for ATLLAS II hypersonic configuration. Using aerodynamic data provided by ONERA, a standard nonlinear ray tracing method enabled us to quantify the overall sonic boom carpet. For the considered configuration, it leads to a sonic boom slightly lower than Concorde's one in terms of peak overpressure, of longer duration and with a much larger carpet because of highest speed and flight altitude. Finally, FLHOWARD3D has been applied to evaluate the effect of atmospheric turbulence on sonic boom propagation. Two cases were investigated: (i) ground track sonic boom in cruise conditions, (ii) sonic boom penetration into the lateral shadow zone. The first case provides information on the global annoyance caused by sonic boom, while the second precise the geographical extent of the boom.

In the first case of undertrack boom, a statistical study is performed for two intensities of turbulence. For a high intensity of turbulence, the mean peak pressure is shown to increase before reaching a plateau and then decreases. This behaviour differs from the one observed for laboratory scale experiments. This difference is explained by the fact that at boom scale, absorption and nonlinear effects are much weaker and the wavelength to integral length scale is larger by one order of magnitude. For a low turbulence intensity, pressure nearly remains unchanged. For both turbulence intensities, the standard deviation is relatively high, indicating that there is a high probability of enhanced peak pressure. For the second case of propagation in the shadow zone, in the non turbulent case, the exponential wave decay associated to creeping waves penetration is recovered along with formation of a caustic in altitude. The turbulent case quantifies the effects of random scattering. A statistical analysis for various turbulence intensities shows that the mean peak pressure slightly increases, with a variability augmenting with turbulence intensity.

As a conclusion, turbulence is not decreasing the sonic boom impact as was previously thought. This is due to the specific low-frequency waveform generated by the ATLLAS II configuration. It generally leads to a higher mean peak pressure on the ground with a probability of occurrence of amplitudes up to twice the unperturbed signal. Nevertheless quantification in terms of loudness metrics adapted to sonic boom remains to be done and would deserve future studies with much more computational resources. Nonetheless, relying on the estimations performed in this study, overland flight of the investigated configuration is likely to be deemed highly annoying by a significant percentage of the impacted population.

Introduction

Focusing of weak shock waves at a caustic is a fundamental phenomenon. A caustic is an amplitude singularity which is structurally stable. Different kinds of caustics exist: fold, cusp, swallow-tail etc. They are described and classified by the theory of catastrophes [START_REF] Thom | Structural stability and morphogenesis[END_REF][START_REF] Berry | Waves and Thom's theorem[END_REF][START_REF] Berry | Singularities in waves and rays[END_REF]. This theory provides the generic form for the high-frequency limit of the pressure field around the caustics in the linear regime, where diffraction cannot be neglected and is the main physical mechanism. A review of the acoustical caustics in linear regime can be found in [START_REF] Marston | Geometrical and catastrophe optics methods in scattering[END_REF]. If the incoming wave is a shock wave, the solution given by the linear theory is invalid since it predicts infinite amplitude. To regularize the field at the caustic, nonlinearities needs to be added to diffraction [START_REF] Guiraud | Acoustique géométrique, bruit balistique des avions supersoniques, focalisations[END_REF]. For weak shock waves, two kinds of caustics have been studied: the fold ans the cusp. For the fold caustic, Guiraud [START_REF] Guiraud | Acoustique géométrique, bruit balistique des avions supersoniques, focalisations[END_REF] established the nonlinear Tricomi equation satisfied by the pressure field. Numerical resolution of this equation indicates that taking into account nonlinear effect along with diffraction is sufficient to get a bounded field [START_REF] Auger | Numerical Simulation of Sonic Boom Focusing[END_REF][START_REF] Marchiano | Numerical simulation of shock wave focusing at fold caustics, with application to sonic boom[END_REF][START_REF] Salamone | Solution of the Lossy Nonlinear Tricomi Equation Applied to Sonic Boom Focusing[END_REF]. Similarly, for the cusped caustics, Cramer ans Seebass [START_REF] Cramer | Focusing of weak shock waves at an arete[END_REF] showed that for 93 a step shock, the acoustic field satisfies locally the KZK equation. Coulouvrat [START_REF] Coulouvrat | Focusing of weak acoustic shock waves at a caustic cusp[END_REF] generalized the previous study by showing that the KZK equation is valid for any incident shock wave. He also formulated the boundary conditions of the problem for its numerical resolution. Marchiano et al. [START_REF] Marchiano | Nonlinear focusing of acoustic shock waves at a caustic cusp[END_REF] used these results to investigate the pressure field resulting from an N-wave focusing at a caustic cusp; again it was found that nonlinearity allows the field to remain bounded on the caustic.

Along with the numerical work of Marchiano et al. [START_REF] Marchiano | Nonlinear focusing of acoustic shock waves at a caustic cusp[END_REF], shock wave focusing on a caustic was explored to determine if wavefront folding could occur for weak shock waves. First, experimentally, Sturtevant and Kulkarny and Beasley et al. [START_REF] Sturtevant | The focusing of weak shock wave[END_REF][START_REF] Beasley | A laboratory investigation of n-wave focusing[END_REF] produced focusing around caustic cusp. Then, Piacsek [START_REF] Piacsek | Atmospheric turbulence conditions leading to focused and folded sonic boom wave fronts[END_REF] used the Nonlinear Progressive wave Equation (NPE) to study the issue numerically. All theses studies showed evidence of wavefront folding.

The present work aims to study structural stability of a cusped caustic for weak shock wave under the influence of kinematic turbulence. To do so, a focused wavefront similar to the one proposed Piacsek is propagated through a turbulent medium as shown on Fig. 5.1. First, a 2D computation of an N-wave focusing in a homogeneous medium is performed to assess the validity of the FLHOWARD model to handle this case. Then, the same focusing is simulated in a 3D turbulent medium. Results are compared to the homogeneous case. The third direction, along which the field is invariant in the homogeneous non turbulent case, allows us to analyse variability in the turbulent case 5.2.

Focusing in a homogeneous medium

Computational process

A 2D curved wavefront is propagated in a homogeneous atmosphere. To simplify the boundary conditions and mesh convergence, the wavefront with a parabolic ripple given by Piacseck is used:

τ c (y) = L c c 0 1 + L c 2R c y L c 2 -1 . (5.1)
The parabolic ripple is given in terms of arrival time τ c and is function of the maximal depth L c = 25, 5 m and minimum radius of curvature R c = 510 m. The temporal signal is an ideal N-wave of peak overpressure 400 Pa and duration 0.2 s. The considered medium of propagation is inviscid air, with a mean sound speed c0 = 340 m.s -1 , a density ρ 0 = 1.2 kg.m -3 and a nonlinear coefficient β = 1.2. The shock formation distance computed for a sine wave of same amplitude and 0.2 s period is 920 m. Given the fact that the initial wave has already shocks and that its amplitude tends to increase because of focusing, this value shows that nonlinear effects play a significant role all along the propagation. The main direction of propagation is along the x direction. The propagation distance is chosen as 2720 m so that the whole focusing process occurs in the computational domain. 2048 points are used to discretize it. The transverse direction y is 680 m long. It is discretized with 512 points. Periodic boundary conditions are assumed laterally in the y directions. Because of focusing, periodicity does not influence the focal region. For the time variable τ , the numerical window is [-0.2, 0.3] s. As a result of the convergence study carried out in the next section, the temporal window involves 4096 points. As a consequence, the time signals have a frequency cut-off around 4000 Hz when considering the Nyquist frequency.

Nonlinear caustic

The peak overpressure (max τ (p(x, y, τ ))) in the homogeneous medium is presented on Fig. 5.3 for linear (β = 0) and nonlinear cases. Note that simulation for linear case for which convergence is theoretically not attainable, is made for N τ = 16384 points for the time window. The focusing occurs on a line at the center of the domain y = 0. Its maximum occurs after around 575 m of propagation in the linear case, but after 950 m in the nonlinear case. It is around 3 times the initial peak overpressure in the nonlinear case compared to 6.5 times in the linear case. Comparison can be made with a geometrical ray tracing computation shown on Fig. 5.4 for which the rays focus at around 500 m which is the beginning of the focal spot for nonlinear case. As expected, in the linear case the focus is slightly shifted away from the geometrical caustic due to diffraction. On the contrary, we observe a strong shift of the focus position and a reduced amplification in the nonlinear case. All this indicates a strongly nonlinear configuration. The characteristic cusped shape is clearly visible for the linear peak overpressure field, while this one is strongly perturbed in the nonlinear case.

The transverse profile of maximum peak overpressure are presented on Fig. 5.5 for 5 distances of propagation (540 m, 1080 m, 1630 m, 2170 m and 2720 m) in the linear (540 m, 1080 m only) and nonlinear cases. As expected in the linear case, no side lobe is visible before the focus while 2 local maxima appear beyond, one for each branch of the cusp. In the nonlinear case, this pattern is considerably modified. The width of the main lobe is enlarged, its amplitude is reduced so that the 2 side lobes disappear. This explains the differences between linear and nonlinear cases of Fig. 5.3. In Fig. 5.6, the nonlinear wavefronts in the (τ, y) plane are displayed at 6 different propagation distances (0m, 540 m, 1080 m, 1630 m, 2170 m and 2720 m). The first one is the initial N-wave with the parabolic ripple in retarded time (Eq. 5.1). Then as the propagation distance increases, the wavefront is progressively folding until a swallow tail shape is reached for both the head and tail shock, in agreement with the catastrophe theory. This folding goes along with the amplification of the pressure field. Note the swallow tail shape already appears at 540 m, much beyond the focal point of maximum amplitude but around the geometrical focus given by ray theory. This confirms that the caustic keeps stable in the nonlinear regime even though the cusp shape is no more visible when looking only at the pressure amplitude. Note that on the last plot at x = 2720 m, parasite reflections occur on the lateral boundaries due to the periodic boundary conditions. However, they do not influence the focal spot.

Grid convergence

The homogeneous computational case is also used to achieve a grid convergence study. It is a highly demanding case because the amplitude of the shock is governed by a balance between diffraction and nonlinear effects. Since results are mostly sensitive to time discretization [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF], it is chosen to only vary the discretization of the time window, In linear regime and nonlinear regime .

N τ in the range 64 to 16384. On Fig 5 .7 the maximum of the peak pressure in the domain is plotted against the number of temporal points. It shows the maximum pressure increases along with mesh refinement. However, this increase is bounded in the nonlinear case according to Guiraud's assumption, but is unbounded in the linear case. This confirms that FLHOWARD nonlinear code can properly capture a cusped caustic. In the nonlinear case, the slope begins to decrease after 1024 points in time discretization with a convergence (error smaller than 5%) after 4096 points. Therefore, this number of points is chosen as a trade-off between accuracy and efficiency on the rest of this Chapter. Since the chosen configuration is inviscid (no absorption nor relaxation), the rise time is not bounded because it is controlled by the numerical dissipation only.

Focusing in a turbulent medium

3D propagation medium

The computation uses the same numerical parameters as the previous homogeneous case. However the propagation medium is now a 3D turbulent one. First a third direction z is introduced with the same parameters as the y direction: 780 m long discretized with 512 points. Periodicity is also enforced on its boundaries as was done for the y direction. This is coherent with the periodicity requested to synthesize turbulence. Also, to reduce the computational cost of the simulation, the wave is propagated only over 1360 m as the maximum peak pressure has been reached before that distance in the homogeneous case. It is therefore discretized using only 1024 points in this direction thus keeping the same axial resolution as for the homogeneous case. The turbulent conditions chosen for this study are similar to the ones selected in the Sec. 4.5 of the previous chapter. The turbulent velocity fluctuations are considered frozen and are generated as explained in Sec. 4.3. The outer scale of the vortices L 0 is 100 meters and the turbulent intensity following a von Kármán spectrum is set to a high value σ v = 2m.s -1 . Thermal turbulence is not considered here: c ′ 0 = 0.

Caustic stability

Fig. 5.8 presents the peak overpressure along the propagation. On the top, the homogeneous case is reproduced while on the bottom, the peak overpressure in the central plane max τ (P (x, y, z 0 = 0, τ )) is plotted. Similarly to the plane wave and shadow zone analysed in Chapter 4, the turbulence creates multiple random focusing and defocusing areas. Nevertheless, the main focal spot remains clearly identifiable in the vicinity of the unperturbed one. Turbulence shifts the position of the focal area randomly, changes its amplitude and its width. Fig. 5.9, shows the positive peak overpressure in the transverse plane (y, z) at the propagation distance x 0 = 1360 m (end of the numerical domain): max τ (P (x 0 , y, z, τ )). The position of the wavefront focusing is clearly shown to stay localized around the line y = 0. Turbulence is therefore inducing significant random variations of the peak overpressure in terms of both amplitude and axial position but it is only slightly shifting the position of the focus in the y transverse direction and preserves the existence of the caustic. Wavefront folding in the turbulent case is illustrated by Fig. 5.10 for 6 different distances x in the central plane z = 0. Compared to the homogeneous case, one observes that the overall shape of the wavefront is conserved. Folding still begins to occur mostly beyond the geometrical focus x = 540 m and further develops. This therefore numerically confirms the structural stability of the caustic. Nevertheless, due to turbulence, additional random focusing are visible at various positions. So one can say that the process of weaker random focusing due to turbulence superimposes to the main geometrical focusing process. Fig. 5.11 presents some temporal signals and their spectra. First the input Nwave signal is presented at the top. Its spectrum is classical with a cut-off at around 3000 Hz (not shown here) below the theoretical cut-off of the mesh. Then, a signal outside the focal area is presented in the turbulent case. It appears as a rounded Nwave with small perturbations localized just after the head and tail shocks. Regarding the signal spectrum, some frequencies are enhanced while others are lessened due to constructive or destructive interferences resulting from random scattering. These features are identical to those detailed in Chapter 4. The third signal extracted is at the focus for the homogeneous medium. It exhibits the characteristic U shape associated to shock wave focusing. Its spectrum is similar to the N-wave one, only being globally increased for all frequencies because of amplification. Finally the signal at the focal point in a turbulent medium is shown (point of maximum amplitude in the central plane z = 0). The U-wave shape is again retrieved, turbulence having little effects on this signal: only small high frequencies ripples are visible. Let us nevertheless recall that the position of this focus (mainly the axial one) and the peak overpressure are strongly influenced by turbulence. Accordingly, the spectra is weakly modified compared to the homogeneous focus.

Thus by observing the peak overpressure maps, the wavefronts and some temporal waveforms, one can conclude that the characteristics of the caustic in an homogeneous medium persist in the turbulent case. In particular, the focusing still takes place along a line near the axis y = 0, main wavefront folding persists and the typical U-shaped waveforms are retrieved at the focus. Consequently, at least for the parameters chosen for this computation, the caustic remains stable. Nevertheless turbulence strongly influences the axial position of the focus and its amplitude, and adds local random focusing processes to the geometrical one.

Statistical analysis

In the previous section, the caustic has been shown to keep stable in a turbulent medium even though it is not unaffected. Now the effects of the turbulence on some of the caustic characteristic are synthesised statistically. With this in view, the third z direction of the simulation is used to performs statistics. This is possible because the homogeneous problem is invariant in this direction.

Peak pressure along the propagation

The main characteristic of the caustic investigated here is the peak pressure amplitude as caustics are defined as amplitude singularities for geometrical acoustics. The average peak pressure is computed as the mean along the transverse direction z of the maximum peak overpressure for each propagation distance x: P z (x) = mean z (max y,τ (P (x, y, z, τ ))). This quantity is plotted on Fig. 5.12 as the black dotted line. It is observed that this average peak pressure has a shape globally similar to the homogeneous case (black solid line). It increases up to the geometrical focus around 550m. Then it keeps more or less constant over the numerical domain. It is likely to decrease at larger distances which are outside the numerical domain. Compared to the homogeneous case the pre-focus growth is larger due to the multiple turbulent focusing. Also the corresponding standard deviation (red dotted line) has the same behaviour as the average: increasing up to geometrical focus before reaching a plateau. Note that we are not absolutely sure at present stage that the beginning of the plateau around the geometrical focus is not linked to the probability of occurrence of the first caustic. It is interesting to note that along the propagation, the maximum peak pressure P max (x) = max y,z,τ (P (x, y, z, τ )) in a turbulent field (green solid line) is always much higher than for the homogeneous case and can reach really high amplitude: more than twice the homogeneous case value. Again, the overall maximum is max x (P max (x)) is reached near the geometrical focus.

Focal width

Another interesting characteristic of the caustic is the width in the transverse direction y of the focal area along the propagation distance x. In this study we define it as the width W (x, z) of the area over which the peak overpressure max τ (P (x, y, z, τ )) is higher than the half of the maximum peak overpressure max y,τ (P (x, y, z, τ )). It is then averaged over the z direction to obtain W z (x). The evolution of this average width along the propagation is observable on Fig. 5.13. In the homogeneous medium (solid black line), it is equal to the size of the domain as long as there is no focusing. Then it reaches its minimum value near the focal point before slowly increasing. The behaviour is similar in the turbulent case (dashed black line) but occurs sooner because of random focusing. Again the minimum value is close to the geometrical focus. Beyond it, its value is very similar to the homogeneous case indicating again a weak sensitivity to turbulence in the transverse direction. The standard deviation increases rapidly in the early stage when focusing is only due to the turbulence and then decreases rapidly before reaching a plateau of about 30 m.

Focus position

The last characteristic of the caustic examined in the present work is the focus position. Here, the focus is defined as the position (x f (z), y f (z)) of the maximum peak overpressure in a given plane z = cte. An example for plane z = 0 is given on Fig. 5.8 with the corresponding focused signal shown on Fig. 5.11 (last line). For each plane defined by a given z value this position is reported on Fig. 5.14. The variations of focal position in the transverse direction y turn out small: always less than ±20 m away from its position in the homogeneous medium. This confirms statistically our previous observations from Figs. 5.8 and 5.9 according to which caustics slightly move away in the transverse direction y. On the contrary, in the propagation direction x, the variation is quite large with values comprised in a range of 1000 m. The mean axial position is 888 m with a standard deviation 258 m. Note that this mean value is shifted backward compared to the homogeneous case but is nevertheless not equal to the geometrical focus position. Ellipses representing the average distance in both directions from the focal point in a homogeneous medium (elongated dark pink ellipse) and 3 times its variance (light pink ellipse) are also shown on this figure.

Conclusion

A numerical study about the structural stability of a nonlinear cusped caustic has been performed by propagating a rippled wavefront in a turbulent medium. To achieve this goal, tools developed in the thesis: FLHOWARD 3D and synthesis of turbulent field, were used. First, the caustic formation in a homogeneous medium was studied in the linear and nonlinear regimes. It enabled us to qualitatively validate the use of the FLHOWARD method for this study. Then the caustic stability relative to a random perturbation is established for the chosen parameters: the cusped focusing process is shown to persist in the turbulent medium. Finally, a statistical study was performed on three parameters of the caustic: amplitude, width and position. In terms of amplitude, pressure amplification occurs sooner compared to the unperturbed case up to the geometrical focus. Mean value of the peak overpressure keeps comparable to the homogeneous case but with a high variability. Focus position turns out to vary mostly in the axial direction.

Chapter 6

Conclusion and Perspective

Conclusion

In this thesis, an original three-dimensional one-way method to model and simulate the nonlinear propagation of acoustic shock waves in the atmosphere has been developed. It relies on a scalar wave equation which includes diffraction, flow and heterogeneities effects, nonlinearity, thermoviscous absorption and molecular relaxation. The model strives to be as high order as possible while remaining computationally tractable. To do so, its numerical resolution relies on the fractional step method which permits to solve different simpler problems using algorithms that are as efficient as possible. When possible, spectral or analytical solutions are employed. For the other terms, finite differences method is used. To remain consistent, both the finite differences and split-step schemes are of second order accuracy. To asses the validity and accuracy of the resulting FLHOWARD3D software, its dispersion relation, including the split-step scheme, was established and shown to be exact up to order M 2 0 kz 4 so of second order relative to the flow Mach number M 0 and fourth order relative to the propagation angle measured by the dimensionless wavenumber kz . Rigid ground and free field boundary conditions have been implemented. The algorithm is implemented for high performance computing on distributed memory architecture.

In Chapter 3, FLHOWARD3D has been validated using quantitative test cases covering all involved operators and boundary conditions. All validations turned out to be satisfactory even in cases were the model equation is beyond its theoretical range of validity. The parallelization performance tests showed also satisfying behaviour for both weak and strong scaling. This allows us to target numerical domains of the order of a few billion points necessary for contemplated sonic boom applications. Along to the establishment of dispersion relations, the waveguide test case helps us to choose the model level of accuracy. It is more precise numerically to take into account only the linear effect of flow motion and neglect quadratic convection terms. Also, there is little benefit to include some coupling terms (described by operators H 2 or H (s)

2 ) which especially include the effect of flow gradients. In some cases, omitting this coupling operator even significantly reduces the numerical error. The case of wave scattering by a vortex confirms these observations and also shows that the transverse components of the flow motion have little effect on the wave diffraction pattern. Thus, only the simplest operator (H (s) 1 ) was considered for heterogeneous and flow effects in the rest of the thesis.

Chapter 4 is dedicated to quantify the effects of turbulence in the atmospheric boundary layer on sonic boom. This analysis has been performed for ATLLAS II hypersonic configuration. Using aerodynamic data provided by ONERA, the overall sonic boom carpet was quantified using a standard nonlinear ray tracing method. For the considered configuration, it leads to a sonic boom slightly lower than Concorde's one in terms of peak overpressure, of longer duration and with a much larger carpet because of highest speed and flight altitude. The tool developed in the first two chapter, FLHOWARD3D software, has been applied to evaluate the effect of atmospheric turbulence on sonic boom propagation. Two cases were investigated: (i) ground track sonic boom in cruise conditions, (ii) sonic boom penetration into the lateral shadow zone. The first case provides information on the global annoyance caused by sonic boom, while the second one precises the geographical extent of the boom. For undertrack boom, a statistical study is performed for two intensities of turbulence. For a high intensity of turbulence, the mean peak pressure is shown to increase before reaching a plateau and then decreases. This behaviour differs from the one observed for laboratory scale experiments. This difference is explained by the fact that, at boom scale, absorption and nonlinear effects are much weaker and the wavelength to integral length scale ratio is larger by one order of magnitude. For a low turbulence intensity, average peak overpressure nearly remains unchanged. For both turbulence intensities, the standard deviation is relatively high, indicating that there is a high probability of enhanced peak pressure. Peak overpressure statistical distribution follow a gamma distribution thus confirming observations from laboratory scale experiments. For the second case of propagation in the shadow zone, in the non turbulent case, the exponential wave decay associated to creeping waves penetration is recovered along with formation of a caustic in altitude. The turbulent case quantifies the effects of random scattering. A statistical analysis for various turbulence intensities shows that the mean peak pressure slightly increases, with a variability augmenting with turbulence intensity. As a conclusion for ATLLAS II configuration, turbulence is not decreasing the sonic boom impact as could have been thought a priori. This is due to the specific low-frequency waveform generated by the ATLLAS II configuration. Turbulence statistically leads to an enhanced mean peak pressure on the ground with a probability of occurrence of amplitudes up to twice the unperturbed signal.

Finally, the structural stability of a nonlinear cusped caustic has been investigated numerically by propagating a rippled wavefront in a turbulent medium. First, the caustic formation in a homogeneous medium was studied in the linear and nonlinear regimes. It enabled us to qualitatively validate the use of the FLHOWARD3D method for this study. Then the caustic stability relative to a random kinematic perturbation is established for the chosen parameters: the cusped focusing process is shown to persist in the turbulent medium. Finally, a statistical study was performed on three parameters of the caustic: amplitude, width and position. In terms of amplitude, pressure amplification occurs sooner compared to the unperturbed case up to the geometrical focus. Mean value of the peak overpressure keeps comparable to the homogeneous case but with a high variability. Focus position turns out to vary mostly along the axial direction.

Perspective

New development could be added to the FLHOWARD3D software. One of the major assumptions in the method is that the ground is plane and perfectly reflecting. Relying on previous works performed for the electromagnetic split-step Fourier parabolic equation [START_REF] Dockery | An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation[END_REF], a finite impedance model could be implemented by mixing cosine and sine transforms. Also topography could be implemented first by knife-edge diffraction, which is a zeroth order approximation but has shown satisfying results [START_REF] Donohue | Propagation modeling over terrain using the parabolic wave equation[END_REF][START_REF] Dockery | An overview of recent advances for the TEMPER radar propagation model[END_REF]. Performances of the code could be increased by adapting the parallelization so as to run FLHOWARD3D on several thousands of cores. With this in view, a 2D domain decomposition could be used instead of the present 1D one. This has been shown to be a viable way to improve scaling [START_REF] Laizet | Incompact3d: A powerful tool to tackle turbulence problems with up to O (105) computational cores[END_REF] for FFT dominated algorithms. An other possible development could be to modify the method by replacing plane wave retarded time by cylindrical or spherical retarded times so as to better follow wavefronts for line or point sources.

During the study of ATLLAS II sonic boom propagation in the Planetary Boundary Layer, an increase of the mean peak overpressure was outlined. This was assumed to be due to the ratio between acoustic wavelength and integral scale of turbulence, which was quite high here compared to previous studies. An investigation of the effect of the size of the heterogeneities would be needed to better understand this phenomenon. Also, the temporal mesh in this study was not fine enough to obtain fully converged results for rise time and to allow the computation of sonic boom metrics. Thus performing the same study on finer meshes especially for time variable, would be interesting. For sonic boom, it is needed to investigate other aircraft configurations: for low boom design, the waveform is shaped so as to be less annoying. Effects of atmospheric turbulence on this kind of pressure signals is currently unknown [START_REF] Maglieri | Sonic Boom Six Decades of Research[END_REF] and would deserve further investigations. Comparison with flight test campaigns is an indispensable step to fully qualify the present method. Sensitivity of the boom variability to the turbulence model remains an unexplored issue, in particular effects of PBL inhomogeneity and anisotropy. The outlined gamma distribution for pressure characterised by only two parameters, could be used as a starting point for model reduction, quantifying the dependence of these parameters relative to turbulence and waveform characteristics. The present study about caustic stability could be applied to sonic boom focusing resulting from transonic acceleration, which remains one of the main obstacles to overland supersonic flight. Also, FLHOWARD3D capabilities are planned to be applied to other acoustical shock waves studies, in particular resulting from geophysical sources.

Appendix A

Flhoward equation

Details of the development of the Eq. 2.5 is presented. It is recalled that the goal is to obtain an equation that can be numerically solved using a one-way method. First Eq. 2.5 is recalled:

1 c 2 0 D 2 s p a Dt 2 -ρ 0 ∇ • ∇p a ρ 0 + 2 dV 0j dz t -∞ ∂ 2 p a (x, t ′ ) ∂z∂x j dt ′ = (A.1) - 2 c 2 0 u 0 • ∂∇p a ∂t -2 ∂u 0j ∂x i t -∞ ∂ 2 p a (x, t ′ ) ∂x i ∂x j dt ′ + β ρ 0 c 4 0 ∂ 2 p 2 a ∂t 2 + δ c 4 0 ∂ 3 p a ∂t 3 .
It is written in the form of an homogeneous wave equation with a perturbation term to become Eq. 2.6.

To do this, each term is developed. For the convective derivative:

A = D 2 s p a Dt 2 = ∂ 2 p a ∂t 2 + 2V 0 • ∇ ∂p a ∂t + V 0 • ∇ (V 0 • ∇p a ) . (A.2)
Since V 0z = 0 et V 0 = V 0 (z), then:

A = ∂ 2 p a ∂t 2 + 2 V 0x ∂ 2 p a ∂x∂t + V 0y ∂ 2 p a ∂y∂t + V 0x ∂ ∂x

A.1 Retarded time

To allow for a one way numerical resolution, retarded time τ = t -x/ c0 is introduced:

P ′ = - 1 c0 2 2 V 0x ∂ 2 p a ∂x∂τ - 1

A.2 Wide angle parabolic approximation

Since second order derivative is not adapted to one-way method, a wide angle parabolic approximation is applied on the perturbation term: So the perturbation term is now:

∂ 2 p a ∂x 2 = 2 c0
P ′ = - 1 c0 2 2 V 0x ∂ 2 p a ∂x∂τ - 1

A.3 Potential

To handle shocks, a pseudo-potential is introduced:

p a (x) = ∂φ ∂t (x). (A.15)
It results in: Appendix B

P ′ = - 1 c0 2 2 V 0x ∂ 2 φ ∂x∂τ - 1 c0
∂ 2 φ ∂τ 2 + V 0y ∂ 2 φ ∂y∂τ +V 2 0x - ∂ 2 φ ∂y 2 - ∂ 2 φ ∂z 2 + 1 c0 2 ∂ 2 φ ∂τ 2 +2V 0x V 0y ∂ 2 φ ∂x∂y - 1

Finite differences and tridiagonnal system

In this appendix the term that corresponds to the coupling between diffraction, heterogeneities and wind for a non plane wave are numerically solved thank to finite differences in the frequency domain. It corresponds to the following equation: First, we go to the frequency space. For the sake of simplicity, the index corresponding to the frequency will be ignored in the following.

∂ 2 φ ∂x∂τ = - V 0y c0 ∂ 2 φ ∂y∂τ + V 2 0x 2 c 0 ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 + V 0x V 0y c0 1 c0 ∂ 2 φ ∂τ ∂y - ∂ 2 φ ∂x∂y - V 2 0y 2 c 0 ∂ 2 φ ∂y 2 + c0 V 0x τ -∞ ∂ 3 φ ∂x∂z 2 dτ ′ -c0 τ -∞ ∂ ∂z V 0x ∂ 2 φ ∂x∂z dτ ′ + c0 V 0y τ -∞ ∂ 3 φ ∂y∂z 2 dτ ′ -c0 τ -∞ ∂ ∂z V
iω ∂ φ ∂x = -iω V 0y c0 ∂ φ ∂y + V 2 0x 2 c 0 ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 + V 0x V 0y c0 iω c0 ∂ φ ∂y - ∂ 2 φ ∂x∂y - V 2 0y 2 c 0 ∂ 2 φ ∂y 2 + c0 V 0x iω ∂ 3 φ ∂x∂z 2 - c0 iω ∂ ∂z V 0x ∂ 2 φ ∂x∂z + c0 V 0y iω ∂ 3 φ ∂y∂z 2 - c0 iω ∂ ∂z V 0y ∂ 2 φ ∂y∂z + d dz V 0x ∂ φ ∂z -V 0x ∂ 2 φ
∂ 3 φ ∂x∂z 2 = φ i j,k+1 -2 φ i j,k + φ i j,k-1 -φ i-1 j,k+1 + 2 φ i-1 j,k -φ i-1 j,k-1
∆x∆z 2 (B.9)

∂ 3 φ ∂y∂z 2 = θ φ i j+1,k+1 -2 φ i j+1,k + φ i j+1,k-1 -φ i j-1,k+1 + 2 φ i j-1,k -φ i j-1,k-1 2∆y∆z 2
+ (1 -θ) φ i-1 j+1,k+1 -2 φ i-1 j+1,k + φ i-1 j+1,k-1 -φ i-1 j-1,k+1 + 2 φ i-1 j-1,k -φ i-1 j-1,k-1 2∆y∆z 2 . (B.10) Due to the conservative form, the following terms cannot be discretized using classical schemes while allowing the resolution of a tridiagonal sytem of equation. Instead, half-point schemes [START_REF] Fornberg | Generation of finite difference formulas on arbitrarily spaced grids[END_REF] are used.

∂ ∂z V 0x ∂ 2 φ ∂x∂z = V i 0x k+1/2
∆x∆z 2 φ i j,k+1 -φ i j,k -φ i-1 j,k+1 + φ i-1 j,k

- V i 0x k-1/2 ∆x∆z 2 φ i j,k -φ i j,k-1 -φ i-1 j,k + φ i-1 j,k-1 (B.11) ∂ ∂z V 0y ∂ 2 φ ∂y∂z = θ V i 0y k+1/2 2∆y∆z 2 φ i j+1,k+1 -φ i j+1,k -φ i j-1,k+1 + φ i j-1,k -θ V i 0y k-1/2 2∆y∆z 2 φ i j+1,k -φ i j+1,k-1 -φ i j-1,k + φ i j-1,k-1 + (1 -θ) V i-1 0y k+1/2 2∆y∆z 2 φ i-1 j+1,k+1 -φ i-1 j+1,k -φ i-1 j-1,k+1 + φ i-1 j-1,k -(1 -θ) V i-1 0y k-1/2 2∆y∆z 2 φ i-1 j+1,k -φ i-1 j+1,k-1 -φ i-1 j-1,k + φ i-1 j-1,k-1 (B.12) 123 ∂ ∂z V 0x ∂ φ ∂z = θ V 0 k+1/2 ∆z 2 φ i j,k+1 -φ i j,k - V 0 k-1/2 ∆z 2 φ i j,k -φ i j,k-1 + (1 -θ) V 0 k+1/2 ∆z 2 φ i-1 j,k+1 -φ i-1 j,k - V 0 k-1/2 ∆z 2 φ i-1 j,k -φ i-
+ (1 -θ) ρ 0 j,k+1/2 ∆z 2 φ i-1 j,k+1 -φ i-1 j,k - ρ 0 j,k-1/2 ∆z 2 φ i-1 j,k -φ i-1 j,k-1 . (B.15)
All this scheme are of second order accuracy in y and z. In the propagation direction x, a Crank-Nicolson scheme (θ = 1/2) is used. It is of second order while remaining unconditionally stable which allows to have the same size of step as the others part of the algorithm. This scheme results in the solving of a tridiagonal linear system. This is done using the Thomas algorithm [START_REF] Quarteroni | Numerical mathematics[END_REF] which is really efficient. Note: if V 0x or V 0y has no index, it is implicitly the k index. In the following, the cross derivative terms (of higher orders) are neglected. It allows the use of the Alternating Direction Implicit (ADI) method [START_REF] Leveque | Finite Difference Methods for Ordinary and Partial Differential Equations[END_REF]. ADI is part of the splitting methods. It consist in splitting the equation in the two directions so that two "small" tridiagonal linear system are solved instead of one "big" sparse system.

Beginning with the z-direction: This equation take the form of the following linear system:

φ i-1 j,k-1 c0 V 0x ω 2 ∆z 2 - c0 V 0x j,k-1/2 ω 2 ∆z 2 -φ i-1 j,k-1 (1 -θ) - ∆xV 2 0x 2i c0 ω∆z 2 + ∆xV 0x iω∆z 2 + c0 ρ 0k-1/2 ∆x 2iρ 0 ω∆z 2 - ∆x c0 2iω∆z 2 - V 0xj,k-1/2 ∆x iω∆z 2 + φ i-1 j,k 1 - 2 c0 V 0x ω 2 ∆z 2 + c0 V 0x j,k-1/2 ω 2 ∆z 2 + c0 V 0x j,k+1/2 ω 2 ∆z 2 -φ i-1 j,k (1 -θ) ∆xV 2 0x i c0 ω∆z 2 - 2∆xV 0x iω∆z 2 + V 0xk-1/2 ∆x iω∆z 2 + V 0xk+1/2 ∆x iω∆z 2 -φ i-1 j,k (1 -θ) - c0 ρ 0j,k+1/2 ∆x 2iρ 0 ω∆z 2 - c0 ρ 0j,k-1/2 ∆x 2iρ 0 ω∆z 2 + ∆x c0 iω∆z 2 + φ i-1 j,k+1 c0 V 0x ω 2 ∆z 2 - c0 V 0x j,k+1/2 ω 2 ∆z 2 -φ i-1 j,k+1 (1 -θ) - ∆xV 2 0x 2i c0 ω∆z 2 + ∆xV 0x iω∆z 2 + c0 ρ 0j,k+1/2 ∆x 2iρ 0 ω∆z 2 - ∆x c0 2iω∆z 2 - V 0xk+1/2 ∆x iω∆z 2 = φ i j,k-1 c0 V 0x ω 2 ∆z 2 - c0 V i 0x k-1/
A • u = q (B. 17
)
where A is the tridiagonal matrix: Its coefficients are given by: q coefficients are:

A =                 a j,
a j,k = 1 + c0 ω 2 ∆z 2 -2V 0x + V i 0x k+1/2 + V i 0x k-1/2 + iθ∆x ω∆z 2 - V 2 0x c0 + 2V 0x + c0 ρ 0j,k+1/2 2ρ 0 + c0 ρ 0j,k-1/2 2ρ 0 -c0 -V 0xk+1/2 -V 0xk-1/2 (B.19) b j,k = c0 ω 2 ∆z 2 V 0x -V 0x k-1/2 + iθ∆x ω∆z 2 V 2 0x 2 c0 -V 0x - c0 ρ 0j,k-1/2 2ρ 0 + c0 2 + V 0xk-1/2 , (B.20) c j,k = c0 ω 2 ∆z 2 V 0x -V 0x k+1/2 + iθ∆x ω∆z 2 V 2 0x 2 c0 -V 0x - c0 ρ 0j,
q j,k = φ i-1 j,k-1 c0 ω 2 ∆z 2 V 0x -V 0x j,k-1/2 -φ i-1 j,k-1 (1 -θ)i∆x ω∆z 2 V 2 0x 2 c0 -V 0x - c0 ρ 0j,k-1/2 2ρ 0 + c0 2 + V 0xk-1/2 + φ i-1 j,k 1 + c0 ω 2 ∆z 2 -2V 0x + V 0x j,k-1/2 + V 0x j,k+1/2 -φ i-1 j,k (1 -θ)i∆x ω∆z 2 - V 2 0x c0 + 2V 0x -V 0xk+1/2 -V 0xk-1/2 -φ i-1 j,k (1 -θ)i∆x ω∆z 2 + c0 ρ 0j,k+1/2 2ρ 0 + c0 ρ 0j,k-1/2 2ρ 0 -c0 + φ i-1 j,k+1 c0 ω 2 ∆z 2 V 0x -V 0x j,k+1/2 -φ i-1 j,k+1
(1 -θ)i∆x

ω∆z 2 + V 2 0x 2 c0 -V 0x - c0 ρ 0j,k+1/2 2ρ 0 + c0 2 + V 0xk+1/2 . (B.22)
Now on the y directions:

φ i-1 j-1,k - V 0x V 0y 2iω c0 ∆y -φ i-1 j-1,

Fig. 1 . 1 Fig. 2 . 1

 1121 Fig. 1.1 Schematics of waveform distortion. Extracted from [148]. . . . . . . Fig. 1.2 Schematics of sonic boom exposure. Extracted from [122]. . . . . . Fig. 1.3 Example of variability of boom waveform for a B-58 flight test in 1965. Left: low-wind conditions (U 0 ≈ 5 m.s -1 ); right: high-wind conditions (U 0 ≈ 14 m.s -1 ). Extracted from [83]. . . . . . . . . . . . . . Fig. 1.4 Comparison of computed and experimental signals for the D-SEND flight tests. Extracted from [97]. . . . . . . . . . . . . . . . . . . . . . . Fig. 2.1 Comparison of dispersion relations. . . . . . . . . . . . . . . . . . . Fig. 2.2 First order split-step scheme. . . . . . . . . . . . . . . . . . . . . . Fig. 2.3 Illustration of Burgers-Hayes method. Top : pressure -Down : potential. Solid black lines : initial values. Dotted red lines : multivalued Poisson solutions. Solid red lines : maximum of Poisson solution (potential) and corresponding shock pressure waveform. . . . . . . . . . Fig. 2.4 Example a computation case for outdoor propagation over the ground. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 2.5 Comparison of dispersion relations including the split-step effect. . Fig. 2.6 Scheme for 1D domain decomposition on two processors. In red, the direction which is local for the processor. . . . . . . . . . . . . . . . Fig. 2.7 Implemented FLHOWARD3D algorithm. The first variable is the one decomposed among the processors. It is also the one that is contiguous in memory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 3.1 3D pure tone piston in an homogeneous medium. . . . . . . . . . . Fig. 3.2 Pressure amplitude (in Pa -Color level) radiated by a pure tone pressure piston. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 3.3 Pressure amplitude on the axis of a 3D pressure piston. . . . . . . . Fig. 3.4 Scattering of a plane wave by a heterogeneous sphere. . . . . . . . . Fig. 3.5 Pressure amplitude (in Pa -Color level) resulting from scattering of a plane wave by a heterogeneous sphere. . . . . . . . . . . . . . . . . Fig. 3.6 Axial pressure amplitude resulting from scattering of a plane wave by a heterogeneous sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 3.7 Pressure amplitude for the propagation of a plane wave in a nonlinear thermoviscous medium. . . . . . . . . . . . . . . . . . . . . . . . . Fig. 3.8 Pressure amplitude of the first 100 harmonics for the propagation of a plane wave in a nonlinear thermoviscous medium. . . . . . . . . . . Fig. 3.9 Effect of the split-step and discretization on the error for the propagation of a plane wave in a nonlinear thermoviscous medium. . . . . . .

Fig. 3 .Fig. 4 . 1

 341 Fig. 3.10 Linear acoustic modes propagation in a two-dimensional waveguide with shear flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 3.11 Selected acoustic modes of FLHOWARD equation for the propagation in a two-dimensional waveguide with Poiseuille shear flow. . . . . Fig. 3.12 Amplitude error for a propagation over 1000λ. . . . . . . . . . . . . Fig. 3.13 Dispersion error for a propagation over 1000λ. . . . . . . . . . . . . Fig. 3.14 Scattering of a plane wave by a finite-circulation vortex. . . . . . . Fig. 3.15 Normalized RMS pressure amplitude radiated by the scattering of a plane wave by a finite-circulation vortex. . . . . . . . . . . . . . . . . . Fig. 3.16 Normalized RMS pressure amplitude radiated by the scattering of a plane wave by a finite-circulation vortex on a circle placed at 2.5λ. . . Fig. 3.17 Strong scaling: averaged restitution time for one step versus number of cores for a given overall domain size. . . . . . . . . . . . . . . . . Fig. 3.18 Weak scaling: averaged restitution time for one step versus number of cores for a given domain size per core. . . . . . . . . . . . . . . . . . . Fig. 4.1 Representation of the planetary boundary layer by Kaimal [94]. . . Fig. 4.2 Classical profiles in the convective boundary layer by Stull [163]. . Fig. 4.3 Example of vertical profiles of mean sound speed (left) and flow velocity (right) for Monin-Obukhov Similarity Theory. . . . . . . . . . . Fig. 4.4 von Kármán, gaussian and Kolmogorov spectra. . . . . . . . . . . . Fig. 4.5 Example of axial velocity u 0x fluctuation for a random realization of turbulence in the (x, y) plane obtained using a von Kármán spectrum (top) and a Gaussian spectrum (bottom). σ u = 2 m.s -1 , L 0 = 100 m. . . Fig. 4.6 Skin mesh of ATLLAS II configuration (with courtesy of ONERA). Fig. 4.7 CFD extractions of the pressure signal on cylinders (with courtesy of ONERA). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 4.8 Normalized pressure signal √ r(p -p 0 )/p 0 extracted from CFD versus distance along the aircraft (x in meters). From left to right: azimuth 0°(vertical), 60°and 90°(horizontal) for 3 radii of extraction: 30 m (black), 50 m (red), 100 m (blue). On top, without mesh adaptation, on bottom with mesh adaptation (with courtesy of ONERA). . . . . . . Fig. 4.9 Sonic boom carpet with pressure amplitude (Pa) as color level computed using ray tracing for ATLLAS II configuration in the standard atmosphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 4.10 Left: pressure signals on the ground computed using ray tracing for ATLLAS II configuration in the standard atmosphere with zoom on the first shock. Right: corresponding frequency spectrum. . . . . . . . . Fig. 4.11 Sketch of the numerical process for evaluating sonic boom propagation through turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 4.12 Maximum positive peak pressure (color level in Pa) for ATLLAS II configuration along vertical propagation in a synthetic PBL. . . . . . Fig. 4.13 Maximum positive peak pressure (color level in Pa) for ATLLAS II configuration in a horizontal plane (x = 750 m). . . . . . . . . . . . . Fig. 4.14 Maximum positive peak pressure along 2 horizontal lines. . . . . . .

Fig. 4 . 90 Fig. 5 . 1 98 Fig. 5 . 5 1 Shock waves and sonic boom in inhomogeneous media Contents 1 . 1 1 . 2

 49051985511112 Fig. 4.15 Examples of signals (left) and their energy spectrum (right). Form top to bottom: initial N-wave, spiked waveform, multiple spiked waveform and rounded waveform. In black the initial signal, in red signals after propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Fig. 4.16 Average (top) and standard deviation (bottom) of the positive peak pressure along the propagation. . . . . . . . . . . . . . . . . . . . . . . . 78 Fig. 4.17 Skewness (top) and kurtosis (bottom) of the positive peak pressure along the propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Fig. 4.18 Histograms of the positive peak pressure at increasing propagation distances. Left: σ = 0.5 in red and right σ = 2.0 in green. The blue line correspond to the fitted gamma distribution. . . . . . . . . . . . . . . . 80 Fig. 4.19 Rise time (color level in s) for ATLLAS II configuration in a synthetic PBL for σ v = 2 m/s. Along vertical propagation in plane z = 0 (top) and in a horizontal plane (x = 750 m) (bottom). . . . . . . . . . . 82 Fig. 4.20 Average (top) and standard deviation (bottom) of the rise time along the propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Fig. 4.21 Histograms of the rise time for increasing propagation distances. Left: σ = 0.5 m.s -1 in red and right σ = 2.0 m.s -1 in green. . . . . . . . . 84 Fig. 4.22 Peak overpressure (Pa) against rise time (s) at x = 750 m (top) and x = 1500 (bottom) for σ = 2.0 m.s -1 . . . . . . . . . . . . . . . . . . 85 Fig. 4.23 Shadow zone of a sonic boom . . . . . . . . . . . . . . . . . . . . . 86 Fig. 4.24 Sonic boom penetration into the shadow zone. Vertical map of peak overpressure (color level in Pa): (above) no turbulence, (below) one realization of turbulence with σ v = 3.6m.s -1 . . . . . . . . . . . . . . 87 Fig. 4.25 Sonic boom penetration into the shadow zone. Ground map of peak overpressure (color level in Pa) for one realization of turbulence with σ v = 3.6m.s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Fig. 4.26 Temporal signatures at 2 km (left) and their spectrum (right). The signal is taken on the ground (top) and near the caustic (bottom). The red curve indicates a propagation through a realization of turbulence σ v = 3.6m.s -1 compared to the non-turbulent case in black. . . . . . . . 89 Fig. 4.27 Statistics of the positive peak pressure along the propagation in the shadow zone for different turbulent levels. Average(top) and standard deviation (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Fig. 5.1 Sketch of the wavefront focusing in a turbulent medium. . . . . . . 95 Fig. 5.2 3D view of the caustic. In blue: incident wavefront, in red: geometrical caustic in the homogeneous case, in black: sketch of turbulent structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Fig. 5.3 Maximum positive peak pressure (color level in Pa) in the propagation plane in an homogenous medium. Top: linear simulation, bottom: nonlinear simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Fig. 5.4 Same as bottom of Fig.5.3 superposed with distribution of acoustic rays. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

Figure 1 . 3 :

 13 Figure 1.3: Example of variability of boom waveform for a B-58 flight test in 1965. Left: low-wind conditions (U 0 ≈ 5 m.s -1 ); right: high-wind conditions (U 0 ≈ 14 m.s -1 ).Extracted from[START_REF] Hilton | Sonic-boom measurements during bomber training operations in the chicago area[END_REF].

Figure 1 . 4 :

 14 Figure 1.4: Comparison of computed and experimental signals for the D-SEND flight tests. Extracted from [97].

Fig. 1

 1 Fig. 1.4.

13

 13

Figure 2 . 1 :

 21 Figure 2.1: Comparison of dispersion relations.

Figure 2 . 3 :

 23 Figure 2.3: Illustration of Burgers-Hayes method. Top : pressure -Down : potential. Solid black lines : initial values. Dotted red lines : multivalued Poisson solutions. Solid red lines : maximum of Poisson solution (potential) and corresponding shock pressure waveform.

Figure 2 . 4 :

 24 Figure 2.4: Example a computation case for outdoor propagation over the ground.

Figure 2 . 5 :

 25 Figure 2.5: Comparison of dispersion relations including the split-step effect.

Figure 3 . 1 :

 31 Figure 3.1: 3D pure tone piston in an homogeneous medium.

Figure 3 . 2 :

 32 Fig.3.2 shows the pressure amplitude radiated from the piston in the (x, y) plane where the characteristic diffraction pattern can be seen. An interesting thing is that there is no reflected wave at the top and bottom boundaries. This shows that the ABLs are working properly. This can also be seen on the axis pressure Fig.3.3 where no oscillations due to reflections are visible in the far-field. The comparison with the

Figure 3 . 3 :

 33 Figure 3.3: Pressure amplitude on the axis of a 3D pressure piston. FLHOWARD computation, Analytical solution and Analytical solution for the parabolic equation. The near-field is zoomed.

Figure 3 . 4 :

 34 Figure 3.4: Scattering of a plane wave by a heterogeneous sphere.

Figure 3 . 5 :

 35 Figure 3.5: Pressure amplitude (in Pa -Color level) resulting from scattering of a plane wave by a heterogeneous sphere. The grey zone marks the heterogeneity location.

Figure 3 . 6 :

 36 Figure 3.6: Axial pressure amplitude resulting from scattering of a plane wave by a heterogeneous sphere. FLHOWARD, Analytical. The grey zone marks the heterogeneity location.

Fig. 3 .Figure 3 . 7 :Figure 3 . 8 :

 33738 Figure 3.7: Pressure amplitude for the propagation of a plane wave in a nonlinear thermoviscous medium. FLHOWARD, Analytical.

Figure 3 . 9 :

 39 Figure3.9: Effect of the split-step and discretization (N x number of discretization points in the propagation direction) on the error for the propagation of a plane wave in a nonlinear thermoviscous medium.Second order, First order.

Figure 3 . 10 :

 310 Figure 3.10: Linear acoustic modes propagation in a two-dimensional waveguide with shear flow.

Figure 3 . 11 :

 311 Figure 3.11: Selected acoustic modes of FLHOWARD equation for the propagation in a two-dimensional waveguide with Poiseuille shear flow.

Figure 3 .Figure 3 . 16 :

 3316 Figure 3.15: Normalized RMS pressure amplitude (P RM S -Color level) radiated by the scattering of a plane wave by a finite-circulation vortex. The black circle has a radius of 2.5λ and shows where the pressure is extracted for Fig. 3.16 .

Figure 3 . 17 :

 317 Figure 3.17: Strong scaling: averaged restitution time for one step versus number of cores for a given overall domain size. domain (1024 × 1024 × 1024), : domain (2048 × 2048 × 1024), dashed lines: linear scaling.

Figure 3 . 18 :

 318 Figure 3.18: Weak scaling: averaged restitution time for one step versus number of cores for a given domain size per core.: 2 20 , 222 , : 2 24 points per core.

77 4. 6

 776 Planetary Boundary Layer . . . . . . . . . . . . . . . . . . . . . 4.2 The mean atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Turbulent fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Sonic boom in the standard atmosphere . . . . . . . . . . . . . . . . 4.4.1 Geometry and CFD computation . . . . . . . . . . . . . . . 64 4.4.2 Nonlinear ray tracing method . . . . . . . . . . . . . . . . . 65 4.4.3 Sonic boom at ground level . . . . . . . . . . . . . . . . . . 66 4.4.4 Sonic boom metrics . . . . . . . . . . . . . . . . . . . . . . . 67 4.5 Effect of a turbulent atmosphere on the undertrack boom . . . . . . 4.5.1 Computational process . . . . . . . . . . . . . . . . . . . . . 70 4.5.2 Pressure variability . . . . . . . . . . . . . . . . . . . . . . . 71 4.5.3 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . 75 4.5.4 Rise time . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sonic boom in the shadow zone . . . . . . . . . . . . . . . . . . . . 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4 . 1 :

 41 Figure 4.1: Representation of the planetary boundary layer by Kaimal [94].

Figure 4 . 2 :

 42 Figure 4.2: Classical profiles in the convective boundary layer by Stull [163].

  6, b h = 8.4, and b m = 5.3 are numerical constants. Thanks to these universal functions, MOST can be used regardless of the state of the PBL: stable, unstable or neutral. It can be seen as a generalization of the logarithmic wind profile derived for neutral conditions. A profile of temperature and one of flow velocity are presented on Fig. 4.3 for a mostly sunny day with strong wind: Q s = 200 W.m 2 and u * = 0.7 m.s -1 . The parameters were chosen z r = 1 m, T r = 288.15 K and z 0 = 0.1 m.

Figure 4 . 3 :

 43 Figure 4.3: Example of vertical profiles of mean sound speed (left) and flow velocity (right) for Monin-Obukhov Similarity Theory.

Figure 4 . 4 :

 44 Figure 4.4: Three turbulent spectra: von Kármán, Gaussian, and Kolmogorov.

Figure 4 . 5 :

 45 Figure 4.5: Example of axial velocity u 0x fluctuation for a random realization of turbulence in the (x, y) plane obtained using a von Kármán spectrum (top) and a Gaussian spectrum (bottom). σ u = 2 m.s -1 , L 0 = 100 m.

Figure 4 . 6 :

 46 Figure 4.6: Skin mesh of ATLLAS II configuration (with courtesy of ONERA).

Figure 4 . 7 :

 47 Figure 4.7: CFD extractions of the pressure signal on cylinders (with courtesy of ONERA).

Figure 4 . 8 :

 48 Figure 4.8: Normalized pressure signal √ r(p -p 0 )/p 0 extracted from CFD versus distance along the aircraft (x in meters). From left to right: azimuth 0°(vertical), 60°a nd 90°(horizontal) for 3 radii of extraction: 30 m (black), 50 m (red), 100 m (blue).On top, without mesh adaptation, on bottom with mesh adaptation (with courtesy of ONERA).

  at -68 km and 68 km Cut-off at -68 km and 68

Figure 4 . 9 :

 49 Figure 4.9: Sonic boom carpet with pressure amplitude (Pa) as color level computed using ray tracing for ATLLAS II configuration in the standard atmosphere.

Figure

  Figure 4.10: Left: pressure signals on the ground computed using ray tracing for ATLLAS II configuration in the standard atmosphere with zoom on the first shock. Right: corresponding frequency spectrum.undertrack, at the cut-off.

. 30 )

 30 ASEL and CSEL values of the ATLLAS II configuration are reported in Table.

4 . 1 .

 41 They are compared to an ideal N-wave of amplitude 100 Pa, duration 0.27 s and rise time 1 ms typical for Concorde boom undertrack.

Figure 4 . 11 :

 411 Figure 4.11: Sketch of the numerical process for evaluating sonic boom propagation through turbulence.

Figure 4 . 12 :

 412 Figure 4.12: Maximum positive peak pressure (color level in Pa) for ATLLAS II configuration along vertical propagation in a synthetic PBL for 2 turbulence intensities: top σ v = 0.5 m/s, bottom σ v = 2 m/s.

Figure 4 . 13 :

 413 Figure 4.13: Maximum positive peak pressure (color level in Pa) for ATLLAS II configuration in a horizontal plane (x = 750 m) for 2 turbulence intensities: top σ v = 0.5 m/s, bottom σ v = 2 m/s.

Figure 4 . 14 :

 414 Figure 4.14: Maximum positive peak pressure along 2 horizontal lines( y = 0, z = 0) for 2 turbulence intensities: left σ v = 0.5 m/s, right σ v = 2 m/s.

Figure 4 . 15 :

 415 Figure 4.15: Examples of signals (left) and their energy spectrum (right). Form top to bottom: initial N-wave, spiked waveform, multiple spiked waveform and rounded waveform. In black the initial signal, in red signals after propagation.

Figure 4 . 16 :

 416 Figure 4.16: Average (top) and standard deviation (bottom) of the positive peak pressure along the propagation. σ = 0.5 and σ = 2.0

Figure 4 . 17 :

 417 Figure 4.17: Skewness (top) and kurtosis (bottom) of the positive peak pressure along the propagation. σ = 0.5 and σ = 2.0

Figure 4 . 18 :

 418 Figure 4.18: Histograms of the positive peak pressure at increasing propagation distances. Left: σ = 0.5 in red and right σ = 2.0 in green. The blue line correspond to the fitted gamma distribution.

Figure 4 . 19 :

 419 Figure 4.19: Rise time (color level in s) for ATLLAS II configuration in a synthetic PBL for σ v = 2 m/s. Along vertical propagation in plane z = 0 (top) and in a horizontal plane (x = 750 m) (bottom).

Figure 4 . 20 :Figure 4 . 21 :

 420421 Figure 4.20: Average (top) and standard deviation (bottom) of the rise time along the propagation.σ v = 0.5 m.s -1 and σ v = 2.0m.s -1 .

Figure 4 . 22 :

 422 Figure 4.22: Peak overpressure (Pa) against rise time (s) at x = 750 m (top) and x = 1500 (bottom) for σ = 2.0 m.s -1 .

Figure 4 . 23 :

 423 Figure 4.23: Shadow zone of a sonic boom

Figure 4 . 24 :

 424 Figure 4.24: Sonic boom penetration into the shadow zone. Vertical map of peak overpressure (color level in Pa): (above) no turbulence, (below) one realization of turbulence with σ v = 3.6m.s -1 .

Figure 4 . 25 :Figure 4 . 26 :

 425426 Figure 4.25: Sonic boom penetration into the shadow zone. Ground map of peak overpressure (color level in Pa) for one realization of turbulence with σ v = 3.6m.s -1 .

Figure 4 . 27 :

 427 Figure 4.27: Statistics of the positive peak pressure along the propagation in the shadow zone for different turbulent levels. Average(top) and standard deviation (right).

93 5. 2

 932 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Focusing in a homogeneous medium . . . . . . . . . . . . . . . . . . 94 5.2.1 Computational process . . . . . . . . . . . . . . . . . . . . . 94 5.2.2 Nonlinear caustic . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2.3 Grid convergence . . . . . . . . . . . . . . . . . . . . . . . . 96 5.3 Focusing in a turbulent medium . . . . . . . . . . . . . . . . . . . . 101 5.3.1 3D propagation medium . . . . . . . . . . . . . . . . . . . . 101 5.3.2 Caustic stability . . . . . . . . . . . . . . . . . . . . . . . . 102 5.3.3 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . 106 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: Sketch of the wavefront focusing in a turbulent medium.

Figure 5 . 3 :

 53 Figure 5.3: Maximum positive peak pressure (color level in Pa) in the propagation plane in an homogenous medium. Top: linear simulation, bottom: nonlinear simulation.

Figure 5 . 4 :

 54 Figure 5.4: Same as bottom of Fig.5.3 superposed with distribution of acoustic rays.

Figure 5 . 5 :

 55 Figure 5.5: Transverse slice of the peak overpressure at various propagation distances. Top figure: nonlinear simulation at 540 m, 1080 m, 1630 m, 2170 m and 2720 m (from black to light blue); bottom figure: linear simulation at 540 m and 1080 m.

Figure 5 . 6 :

 56 Figure 5.6: Pressure field in the (τ, y) plane. Extracted at a distance of propagation (from left to right and top to bottom): 0 m, 540 m, 1080 m, 1630 m, 2170 m and 2720 m.

Figure 5 . 7 :

 57 Figure 5.7: Maximum pressure as a function of the number of points in retarded time. In linear regime and nonlinear regime .

Figure 5 . 8 :

 58 Figure 5.8: Maximum positive peak pressure (color level in Pa) along the propagation. Top: no turbulence, bottom, one realization of turbulence.

Figure 5 . 9 :Figure 5 . 10 :

 59510 Figure 5.9: Maximum positive peak pressure (color level in Pa) in a transverse plane.

Figure 5 . 11 :Figure 5 . 12 :Figure 5 . 13 :

 511512513 Figure 5.11: Some signals (left) and their spectra (right). From top to bottom (i) input N-wave, (ii) rounded N-wave outside focal spot, (iii) nonlinear focused signal (maximum amplitude) for the nonturbulent case, (iv) nonlinear focused signal (maximum amplitude) for the turbulent case in the central plane.

Figure 5 . 14 :

 514 Figure 5.14: Distribution of the focal position (x f (z), y f (z)). Each point corresponds to one z value. See text for legend of ellipses.

Figure F. 1 :

 1 Figure F.1: Average (top) and standard deviation (bottom) of the pressure along propagation.: Nonlinear (β = 1.2), : linear (β = 0).

Table 1 .

 1 1: Comparison between characteristic scales for a sonic boom and laboratory scale experiments of Averiyanov et al.[START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF].

		Atmosphere Laboratory scale experiments
	Turbulent layer width	1-2 km	0.3-1.8 m
	Integral length scale	100-200 m	16-18 cm
	RMS fluctuation velocity	0-2.5 m/s	0-3.4 m/s
	Pressure amplitude	10-100 Pa	10-100 Pa
	Signal duration	0.1-0.2 s	30-80 µs
	Wavelength	34-68 m	1-2.7 cm

Chapter 2 Flhoward method Contents 2

  .1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

	2.1.1 Equation for nonlinear propagation through a moving het-
	erogeneous medium . . . . . . . . . . . . . . . . . . . . . . . 14
	2.1.2 Partially one-way equation . . . . . . . . . . . . . . . . . . . 15
	2.1.3 Dispersion relation . . . . . . . . . . . . . . . . . . . . . . . 17
	2.2 Numerical resolution . . . . . . . . . . . . . . . . . . . . . . . . . . 19
	2.2.1 Split-step method . . . . . . . . . . . . . . . . . . . . . . . . 19
	2.2.2 Resolution of diffraction effects . . . . . . . . . . . . . . . . 21
	2.2.3 Resolution of flows and heterogeneities effects . . . . . . . . 22
	2.2.4 Resolution of nonlinear effects . . . . . . . . . . . . . . . . . 25
	2.2.5 Absorption and relaxation . . . . . . . . . . . . . . . . . . . 26
	2.2.6 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . 26
	2.2.7 Dispersion relation . . . . . . . . . . . . . . . . . . . . . . . 28
	2.2.8 Parallel computing . . . . . . . . . . . . . . . . . . . . . . . 30
	2.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Table 3 .

 3 1: List of the validation test cases with tested operators and boundary conditions (ABL for Absorbing Boundary Layer -Rigid for the perfectly reflecting boundary).

  2.8. The change in slope between 16 (where the slope is almost 2) and 32 cores for the smallest domain occurs because for 16 cores only one node is necessary, whereas for more than 24 cores several nodes are used. Then communications are no longer local to the node and use the Infiniband network. During the weak scaling test, three different numbers of points per core were tested: 2 20 , 2 22 and 224 . It corresponds respectively to a domain size (N τ × N y × N z ) of (64 × 128 × 128), (256 × 128 × 128) or (256 × 256 × 256) points per core. Computations were performed using from 8 up to 512 cores. Good overall weak scaling is achieved as it can be seen on Fig.3.18. The restitution time increases by a factor around 2 only while the overall domain sizes and number of cores increase by a factor 64. As expected, the performance for the domain of smallest size (2 20 points per core)

	Weak scaling:

Table 4 .

 4 4.10: Left: pressure signals on the ground computed using ray tracing for ATLLAS II configuration in the standard atmosphere with zoom on the first shock. 1: ASEL and CSEL values simulated undertrack and at cut-off for ATLLAS II configuration compared to Concorde-like N-wave.

	Right: corresponding frequency spectrum.	undertrack,	at the cut-off.

  0 c j,0 0 . . .

	b j,1 0 . . . . . . . . .	. . . . . . c j,k . . . . . . b j,nz-2 a j,nz-2 c j,nz-2 . . . 0 . . . . . . . . . . . . . . . b j,k a j,k . . . . . . . . . . . . . . . . . . . . . . . . . . . 0	               	.	(B.18)
	0 . . . . . . . . .	0	b j,nz-1 a j,nz-1

  k (1 -θ) -

						∆xV 2 0x 2i c0 ω∆y 2 -	∆xV 0y 2 c0 ∆y	+	V 0x V 0y ∆x 2∆y c0 2 +	V 2 0y ∆x 2i c0 ω∆y 2
	-φ i-1 j-1,k (1 -θ) +	c0 ρ 0j-1/2,k ∆x 2iρ 0 ω∆y 2 -	∆x c0 2iω∆y 2
	+ φ i-1 j,k
	-φ i-1 j,k (1 -θ)	∆xV 2 0x i c0 ω∆y 2 -	V 2 0y ∆x i c0 ω∆y 2 -	c0 ρ 0j+1/2,k ∆x 2iρ 0 ω∆y 2 -	c0 ρ 0j-1/2,k ∆x 2iρ 0 ω∆y 2 +	∆x c0 iω∆y 2
	+ φ i-1 j+1,k	V 0x V 0y 2iω c0 ∆y
	-φ i-1 j+1,k (1 -θ) -	∆xV 2 0x 2i c0 ω∆y 2 +	∆xV 0y 2 c0 ∆y	-	V 0x V 0y ∆x 2∆y c0 2 +	V 2 0y ∆x 2i c0 ω∆y 2
	-φ i-1 j+1,k (1 -θ) +	c0 ρ 0j+1/2,k ∆x 2iρ 0 ω∆y 2 -	∆x c0 2iω∆y 2
	= φ i j-1,k -	V 0x V 0y 2iω c0 ∆y
	+ φ i j-1,k θ -	∆xV 2 0x 2i c0 ω∆y 2 -	∆xV 0y 2 c0 ∆y	+	V 0x V 0y ∆x 2∆y c0 2 +	V 2 0y ∆x 2i c0 ω∆y 2
	+ φ i j-1,k θ +	c0 ρ 0j-1/2,k ∆x 2iρ 0 ω∆y 2 -	∆x c0 2iω∆y 2
	+ φ i j,k
	+ φ i j,k θ	∆xV 2 0x i c0 ω∆y 2 -	V 2 0y ∆x i c0 ω∆y 2 -	c0 ρ 0j+1/2,k ∆x 2iρ 0 ω∆y 2 -	c0 ρ 0j-1/2,k ∆x 2iρ 0 ω∆y 2 +	∆x c0 iω∆y 2
	+ φ i j+1,k	V 0x V 0y 2iω c0 ∆y
	+ φ i j+1,k θ -	∆xV 2 0x 2i c0 ω∆y 2 +	∆xV 0y 2 c0 ∆y	-	V 0x V 0y ∆x 2∆y c0 2 +	V 2 0y ∆x 2i c0 ω∆y 2
	+ φ i j+1,k θ +	c0 ρ 0j+1/2,k ∆x 2iρ 0 ω∆y 2 -	∆x c0 2iω∆y 2 .	(B.23)
	A coefficients are:
	a j,k = 1 +	iθ∆x ω∆y 2 -	V 2 0x c0	+	V 2 0y c0	+	c0 ρ 0j+1/2,k 2ρ 0	+	c0 ρ 0j-1/2,k 2ρ 0	-c0 , (B.24)

136Appendix F. Effect of nonlinearity for sonic boom propagation through turbulence

List of Tables

Tab. [START_REF] Aanonsen | Distortion and harmonic generation in the nearfield of a finite amplitude sound beam[END_REF].1 Comparison between characteristic scales for a sonic boom and laboratory scale experiments of Averiyanov et al. [START_REF] Averiyanov | Random focusing of nonlinear acoustic N-waves in fully developed turbulence: laboratory scale experiment[END_REF]. . . . . . . . . . . . Tab. 1.2 Main characteristics of FDTD schemes for nonlinear wave equation.

q coefficients are:

Appendix C

Relaxation and absorption parameters

First we introduce the new parameters:

. The ISO 9613-1 norm (1993) give the following values for dioxygen :

For diazote : :

Frequency for thermoviscous absorption is given by:

At the reference pressure and temperature, with h m = 1%, we obtain:

Analytical solution D.1 Linear pure tone acoustic piston in an homogeneous medium

This case is called a piston source but stricly speaking, it is not. For a piston is the velocity is imposed on a surface whereas we use the pressure due to a simpler implementation in the FLHOWARD method. This leads to some differences in the analytical solution [131] [141].

We start with the Kirchhoff theorem in the frequency domain:

where p a (x) is the pressure at the observation point, p a (x 0 ) is the pressure on the piston, Ĝ(x, x 0 ) is the Green function of the problem, S 0 the piston surface and n the normal to the piston. To impose the pressure condition, the Green function Ĝ(x, x 0 ) has to be zero over the piston surface S 0 . The natural choice is: Ĝ

which is null if x 1 = x 2 = x 0 . Integral D.1 can be written:

The Green function derivative is:

Thus pressure is given by:

Appendix D. Analytical solution

The integral can be solved numerically to obtain the pressure field. Here, since the solution is needed only on the axis we will use the cylindrical coordinates and select only the axis to obtain:

D.2 Scattering of a plane wave by a spherical heterogeneity

The analytical solution of Anderson [START_REF] Anderson | Sound scattering from a fluid sphere[END_REF] is briefly recalled. A spherical heterogeneity of radius a, sound speed c h , density ρ h and wavenumber k h is placed in a medium of sound speed c 0 , density ρ 0 and wavenumber k 0 . The problem is decomposed in three wave: the incident plane wave (r > a) p 0 , the scattered one (r > a) p sc and the transmitted (r < a) p t with p = p 0 + p sc in the infinite medium. The solution is obtained by representing these fields using spherical harmonics.

with s = c h /c 0 and t = ρ h /ρ 0 .

Appendix E

Grid convergence for sonic boom propagation through turbulence

The aim of this appendix is to show that the grid used for the undertrack boom computation in Sec. 4.5 is fine enough. To demonstrate this point, a second simulation with coarser grid is performed and results are compared. The largest turbulent intensity σ u = 2 m.s -1 is chosen. In each direction(x, y, z), the number of points is divided by 2. Thus 512 points are used in each spatial direction. Statistics for the positive peak pressure are compared for the two meshes. On : fine mesh, : rough mesh.

Appendix F

Effect of nonlinearity for sonic boom propagation through turbulence

This Appendix deals with the importance of taking into account nonlinearity for the computation of the effect of turbulence on undertrack boom in Sec. 4.5. A computation for the largest turbulent intensity σ u = 2 m.s -1 is performed with nonlinear parameter β set to 0 (so equivalent to linear propagation). Average and standard deviation of the positive peak pressure are presented on Fig. F.1. For the linear computation an increase up to 2 Pa more than the nonlinear computation is seen. One interesting result is that most of the difference for the standard deviation is around the distance where most of the caustics appear. In linear regime, the pressure on the caustics generated by turbulence is theoretically infinite (see Chapter 5 for details). Thus, in linear regime, the value of the peak overpressure strongly depends on the mesh refinement: the finer the mesh, the higher the amplitude. This behaviour is not acceptable and highlights why nonlinearity is needed to accurately simulate sonic boom propagation in the planetary boundary layer.