One-way approximation for the simulation of weak shock wave propagation in atmospheric flows.

One-way approximation for the simulation of weak shock wave propagation in atmospheric flows Louis-Jonardan Gallina) ^ de France, 91297 Arpajon, France Commissariat a l’ energie atomique et aux energies alternatives–DAM–Ile

Gaudard nier and Eric Mathieu Re Universit e Pierre et Marie Curie Paris 06, Institut Jean le Rond d’Alembert (Unit e Mixte de Recherche du Centre National de la Recherche Scientifique 7190), 4 Place Jussieu 75005 Paris, France

Thomas Farges ^ de France, 91297 Arpajon, France Commissariat a l’ energie atomique et aux energies alternatives–DAM–Ile

gis Marchiano and François Coulouvrat Re Universit e Pierre et Marie Curie Paris 06, Institut Jean le Rond d’Alembert (Unit e Mixte de Recherche du Centre National de la Recherche Scientifique 7190), 4 Place Jussieu 75005 Paris, France

(Received 25 July 2013; revised 12 January 2014; accepted 11 March 2014) A numerical scheme is developed to simulate the propagation of weak acoustic shock waves in the atmosphere with no absorption. It generalizes the method previously developed for a heterogeneous medium [Dagrau, Renier, Marchiano, and Coulouvrat, J. Acoust. Soc. Am. 130, 20–32 (2011)] to the case of a moving medium. It is based on an approximate scalar wave equation for potential, rewritten in a moving time frame, and separated into three parts: (i) the linear wave equation in a homogeneous and quiescent medium, (ii) the effects of atmospheric winds and of density and speed of sound heterogeneities, and (iii) nonlinearities. Each effect is then solved separately by an adapted method: angular spectrum for the wave equation, finite differences for the flow and heterogeneity corrections, and analytical method in time domain for nonlinearities. To keep a one-way formulation, only forward propagating waves are kept in the angular spectrum part, while a wide-angle parabolic approximation is performed on the correction terms. The numerical process is validated in the case of guided modal propagation with a shear flow. It is then applied to the case of blast wave propagation C 2014 Acoustical Society of America. within a boundary layer flow over a flat and rigid ground. V [http://dx.doi.org/10.1121/1.4869685] PACS number(s): 43.28.Fp, 43.28.Gq, 43.25.Cb, 43.25.Jh [PBB]

I. INTRODUCTION

The deployment of the International Monitoring System (IMS) of the comprehensive nuclear-test-ban treaty (CTBT) requires the capacity to efficiently simulate the propagation of acoustic shock waves inside the terrestrial boundary layer atmosphere over important distances. The essential influence of atmospheric conditions on the propagating waves has been known since the beginning of the twentieth century (Evers and Haak1 give a historical review of loud sound sources). The characteristics of the detected acoustic signatures either natural or human-made vary greatly depending on the impact of the atmospheric effects on the traveling waveforms: This renders analysis and diagnostic processes of the emitting sources very difficult. Several sources of aerial shock waves, such as explosions,2,3 meteors,4 volcanoes,5 lightning and sprites,6–9 and sonic booms10 are known to produce important nonlinear effects, while propagating in a windy atmosphere. A review of the major physical effects involved in the propagation is given by Norris et al.11 Ray

a)

Author to whom correspondence should be addressed. Also at: Universite Pierre et Marie Curie Paris 06, Institut Jean le Rond d’Alembert. Electronic mail: [email protected]

J. Acoust. Soc. Am. 135 (5), May 2014

Pages: 2559–2570

tracing techniques are often used to propagate infrasonic waveforms over important distances either with nonlinear effects or not. Their computational efficiency allows one to perform statistical analyses.12,13 However, the ray theory presents several major limitations. It is based on a high frequency approximation and it cannot describe any diffraction effect such as penetration of lower frequency waves inside shadow zones14 or local amplification near caustics.15 On the other hand, resolving all scales to solve the complete set of Navier-Stokes equations is numerically prohibitive.16–18 The construction of alternative numerical schemes, which are more precise than the ray tracing method but still remain computationally appealing, is a challenge. Parabolic approximations (like the KZK equation19) go beyond the ray theory limitations and have been extremely successful at describing finite amplitude narrow beams. Aver’yanov et al.20 proposed a generalization of the KZK equation including flow motion at first order. This approach has been used to study the scattering of weak shock waves in a turbulent flow.21,22 However the standard parabolic approximation limits the validity of the model to narrow angle propagation. This raises the need to go beyond parabolic approximation to simulate nonlinear acoustic waves diffracted at important angles, for instance in case of elongated sources, or of

0001-4966/2014/135(5)/2559/12/$30.00

C 2014 Acoustical Society of America V

2559

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 131.156.157.31 On: Thu, 29 Jan 2015 10:33:45

scattering by wind flow stratification and atmospheric turbulence. This is the objective of the present work. It presents the numerical implementation of the nonlinear wave model equation with wind established by Coulouvrat,23 by means of an original resolution method first proposed by Dagrau et al.24 and based on the work of Christopher and Parker.25 This method named HOWARD in Ref. 24 now becomes FLHOWARD (flow and heterogeneous one-way approximation for resolution of diffraction). The presented approach has been purposely developed to keep full angle validity (þ/90 ) at least in the homogeneous and nonlinear cases. A wide angle approximation is made only on heterogeneities and wind correction terms. As a consequence, the numerical resolution is performed in a one-way approach, which keeps the numerical efficiency of parabolic approximations but neglects backscattering. It is to be noted that no absorption effect, neither classical thermoviscosity nor molecular relaxation, is taken into account yet. Even though absorption is weak in the considered illustration case (propagation of a 50 Hz central frequency signal over a few kilometers) it could nevertheless influence especially the high frequency spectrum contained within the shock structure. However, the present work focuses on introducing flow motion in the numerical method. Adding an additional effect would only obscure this objective. The theoretical model will be recalled in Sec. II. It will then be reformulated in a way suitable for the one-way approach (Sec. III). The marching algorithm (Sec. IV) and the numerical discretization (Sec. V) will be presented. Section VI is dedicated to validation tests for the new wind flow terms. Finally, an example of application in the context of atmospheric sound propagation is described and interpreted in Sec. VII.

winds or turbulent fluctuations are generally much smaller. It is difficult to quantify the magnitude order of turbulence in general as it may be strongly dependent on the particular meteorological situation. However turbulent wind fluctuations are unlikely to be larger than a few meters per second. A high precision Doppler lidar experiment26 is used to measure the vertical wind component in a daytime convective boundary layer with a high sampling rate over a flat ground. Measured amplitudes are of order 3 m/s at the most, which is only about 1% of the sound speed. Consequently, nonlinear quadratic effects relative to turbulent fluctuations will be negligible in many similar cases. Two-dimensional density and sound speed fluctuations are assumed to be of relative order M2. Indeed, temperature fluctuations are on the order of 6 C–10 C over the atmospheric boundary layer, which induces sound speed fluctuations of about 1% ¼ O(M2). It is then possible to separate a mean component (noted with a bar) from a spatially fluctuating component (noted with a prime) of order M2,

II. MODEL EQUATION

Notation Ds/Dt ¼ @/@t þ V0 r is used for the convected derivative associated with the mean sheared flow. The first two terms of the left-hand side(lhs) of Eq. (4) constitute the usual wave equation in a heterogeneous medium, convected in the sheared flow V0 as if it were uniform. The second-order convected derivative introduces quadratic terms with respect to V0x (of order M2). The influence of the non-uniformity of the sheared flow is taken into account by the term proportional to dV0x/dz in the lhs of Eq. (4). Note these terms are only of first order M. Hence Eq. (4) is of mixed order (intermediate between orders 1 and 2) with respect to the shear flow Mach number, as it keeps the quadratic terms with respect to the shear flow associated with convection, but not those associated with the flow gradient. A fully consistent equation of second order with quadratic terms (proportional to V0xdV0x/dz) could be established [see Eq. (19) of Ref. 23], but it introduces third-order space derivatives and may get unstable. However, numerical comparisons23 for the mixed-order equation with solutions of the exact Lilley’s equation27 show that the error introduced by Eq. (4) on the phase velocity of guided modes is of order 104 only for Mach numbers of order 0.1. This is much smaller than the error (of order 102) for the equation strictly of order 1, and not much higher than the error (of order 105) for the full second-order equation. This is explained by the fact that convection terms play a dominant role in

We consider a two dimensional inviscid fluid, with x the horizontal coordinate and z the vertical one. For a given point x ¼ (x, z) at time t with ex the unit horizontal vector, we note q(x, t) the density, p(x, t) the pressure, v(x, t) the flow speed, and c0(x) the sound speed. Time scales for significant variations of atmospheric parameters are assumed to be much larger than the characteristic period of acoustic waves. Hence the propagation medium can be considered as “frozen” from an acoustical point of view, and the following decomposition can be introduced, with index 0 for quantities associated to the ambient medium that are independent from time, and index a for time varying acoustic perturbations, f ðx; tÞ ¼ f0 ðxÞ þ fa ðx; tÞ

with

f ¼ ðq; v; pÞ:

(1)

The flow motion is separated into a vertically sheared flow V0(z) ¼ V0x(z) ex with a Mach number M typically of order 0.1 at the most, and a smaller perturbation u0(x), including vertical winds and/or turbulent fluctuations, of smaller Mach number, v0 ðxÞ ¼ V0 ðzÞ þ u0 ðxÞ ¼ V0x ðzÞex þ u0 ðxÞ:

(2)

These hypotheses are reasonable for atmospheric studies since usual storm winds rarely exceed 40 m/s, and vertical 2560

J. Acoust. Soc. Am., Vol. 135, No. 5, May 2014

c0 ðxÞ ¼ c0 þ c00 ðxÞ and

0 þ q00 ðxÞ: q 0 ðx Þ ¼ q

(3)

With all these approximations, it is possible to establish a scalar wave equation [Eq. (20) of Ref. 23] describing nonlinear atmospheric acoustic propagation, ð dV0j t @ 2 pa ðx; t0 Þ 0 1 D2s pa rpa q r dt þ 2 0 q0 dz 1 @z@xj c20 Dt2 ð @u0j t @ 2 pa ðx; t0 Þ 0 2 @rpa 2 ¼ 2 u0 dt @t @xi 1 @xi @xj c0 þ

b @ 2 p2a : q0 c40 @t2

(4)

Gallin et al.: Propagation of shock waves in flows


determining the acoustical phase velocity and are fully taken into account in Eq. (4). On the right-hand side (rhs) of Eq. (4), one can see terms associated successively with the convection of the turbulent flow u0, with its gradient and finally nonlinear terms with respect to the acoustical field. Here b is the usual nonlinear parameter, equal to (c þ 1)/2 for a perfect gas with c the ratio of specific heats. Formally, the partial time derivative appearing in the nonlinear term should be the convected time derivative. However, the acoustical Mach number associated with nonlinear effects is most frequently much smaller than the flow Mach number. In air, for a very intense shock wave of amplitude 2000 Pa (or 180 dB) it is around 0.01, e.g., of order M2. Hence, nonlinear convected terms of order M are negligible because they are of order M3 at most. It is therefore consistent to identify in the last term of Eq. (4) Ds/Dt with @/@t. III. ONE-WAY EQUATION

The numerical resolution of Eq. (4) generalizes the method previously developed in the heterogeneous case24 to the moving case. That method can be summarized as follows: Eq. (4) is written in the form of a linear wave equation in a homogeneous medium, plus a perturbation P of order M including all other terms (associated to the flow motion, density and sound speed heterogeneities and nonlinearities), 1 @ 2 pa @ 2 pa @ 2 pa 2 2 ¼ P: @x @z c20 @t2

(5)

Second-order derivatives with respect to x prevent Eq. (5) to be solved numerically in a forward one-way marching method, from initial plane x ¼ 0 toward increasing values of x. To do this for a one-way numerical approach, Eq. (5) is rewritten in a time frame moving with the mean sound speed by introducing the retarded time s ¼ t x/ c 0 in the selected horizontal direction x, 2 @ 2 pa @ 2 pa @ 2 pa 2 ¼ P0: @z c0 @x@s @x2

(6)

Several one-way approximations can now be obtained. Standard parabolic approximation would amount to neglecting all second-order derivatives with respect to x in Eq. (6). The Claerbout28 approximation or a wide-angle parabolic approximation29 would amount to approximating the second-order derivatives with respect to x with the help of the standard parabolic approximation: ð @ 2 pa c0 s @ 3 pa 0 ds : (7) @x2 2 1 @x@z2 In the present method, the above wide-angle approximation Eq. (7) is used, but only in the perturbation term P 0 of order M and not in the lhs of Eq. (6) which is the dominant term and which is left unchanged. This yields a more precise equation, that is only partially one-way. The next step consists in replacing the acoustic pressure pa by the pseudopotential variable /(x): pa(x) ¼ ð@/=@tÞ (x). The advantage of using the pseudo-potential is that it is continuous through J. Acoust. Soc. Am., Vol. 135, No. 5, May 2014

shocks (while the acoustic pressure pa is discontinuous) and well-adapted to numerical the resolution of the nonlinear part of the equation.24,30 Performing all these steps is tedious but straightforward. This finally yields the partially one-way equation, where for simplicity the terms depending on the wind fluctuations u0(x) were ignored, c0 @ 2 / @ 2 / @2/ V0x ðzÞ @ 2 / V0x ðzÞ @ 2 / þ ¼ þ 2 @x2 @x@s @z2 c0 @x@s c20 @s2 ðs 2 2 dV0x @/ dV0x @ / 0 V0x ðzÞ @ 2 / c0 ds þ dz @z dz @z@x 2 c 30 @s2 0 V 2 ðzÞ @ 2 / c0 ðx; zÞ @ 2 / þ 0x þ 2 c 0 @z2 c2 @s2 0 1 @q0 @/ @/ @q0 @/ c0 c0 þ @z @z 2q0 ðx; zÞ @x @s @x " # 2 b @ @/ þ : (8) @s 2 q 0 c30 @s At last all terms involving an odd derivative of / being numerically unstable, they are transformed into a more conservative form by increasing the order of derivation of / to the next even power. For instance a term such as dV0x/dz @//@z in Eq. (8) is expressed in Eq. (10) under the following form: dV0x @/ @ @/ @2/ ¼ V0x (9) V0x 2 : dz @z @z @z @z Doing so leads to the numerically more stable equation Eq. (10), @2/ V0x ðzÞ @ 2 / V0x ðzÞ @ 2 / c0 @ 2 / @ 2 / þ ¼ þ 2 @x2 @x@s @z2 c0 @x@s c20 @s2 ðs 3 ^ 2 @ / @ / ds0 V0x ðzÞ 2 þ c0 V0x ðzÞ @z @x@z2 2 2 V0x @/ ðzÞ @ 2 / V0x ðzÞ @ 2 / @ V0x ðzÞ þ þ 2 c 0 @z2 @z @z 2 c 30 @s2 ðs 0 c ðx; zÞ @ 2 / @ @2/ V0x ðzÞ ds0 þ 0 2 c0 @z @z@x c0 @s2 c0 1 @q0 @/ @q0 @/ þ 2q0 ð x; zÞ @x @s 2q0 ð x; zÞ @x @x c0 c0 @ 2 / @ @/ q0 ðx; zÞ þ 2q0 ðx; zÞ @z 2 @z2 @z " # 2 b @ @/ þ : (10) @s 2 q 0 c30 @s Equation (10) can be written in a more compact form by splitting its rhs into physical effects, @2/ ¼ ½D/ þ ½H/ þ ½N /: @x@s

(11)

In Eq. (11), operator [D] describes the diffraction effect associated to the linear wave equation in a homogeneous Gallin et al.: Propagation of shock waves in flows

2561


medium without flow [first two terms of the rhs of Eq. (10)]. That dominant part of the equation does not include any oneway approximation yet, as it keeps second-order derivatives with respect to x. Then, operator [H] stands for the “heterogeneities and flow corrections” operator, including the third to fifteenth terms of the rhs of Eq. (10). The wideangle approximation Eq. (7) has been applied on this term, so that no second-order derivatives @ 2//@x2 appear there. Finally, the nonlinear operator [N ] corresponds to the last term of Eq. (10). In order to obtain a numerically suitable expression, [H] is decomposed respectfully to the space derivatives of /, @/ @2/ þ ½H2zz 2 @x @z @3/ @/ @2/ ½ ½ þ V : þ½H3xzz þ V 1z 2xz @x@z2 @x @z@x

½H/ ½H0 / þ ½H1x

(12)

Expressions of the various operators are detailed in the Appendix. It is worth noting that operators [H0 ], [H1x ], and [V 1z ] already appear in the no-flow case.24 The flow terms modify the coefficients of these operators and introduce three new ones: [H2zz ], [H3xzz ], and [V 2xz ]. When adapted to our hypotheses and notations (no sound diffusivity, no transverse flow along the direction Oz, use of a potential variable /) the nonlinear standard parabolic equation of Aver’yanov et al. with flow20 becomes c0 @ 2 / V0x ðzÞ @ 2 / @2/ b @/ @ 2 / ¼ þ þ : 2 @z2 @x@s 0 c30 @s @s2 c20 @s2 q

(13)

The three terms on the rhs of Eq. (13) correspond respectively to terms two, three, and 16 of the rhs of Eq. (10). In Eq. (13) no quadratic term relative to the flow V0x is present, and several linear ones are neglected. The validity of Eq. (10) can be quantified precisely in the linear case of a uniform flow V0x (the case of a stratified flow will be investigated later on by means of numerical simulation), by comparing its dispersion equation with the one of the convected wave equation. A plane wave solution is sought as /ðx; z; sÞ ¼ A exp½ik0 ððkx 1Þx þ kz z c0 sÞ;

ðkx 1Þ þ k2z þ M ¼ 0:

(17)

The two dispersion relations Eqs. (15) and (16) differ by a term of order O(M3) in agreement with the precision of the model Eq. (4), that has not been downgraded by the wide angle approximation. An even higher precision can also be recovered simply by replacing the Mach number M by M/(1 M2). The dispersion relation of Eq. (10) will then differ from the exact one by a quantity of order M4 only. Note that a similar modification of the definition of the sound speed fluctuation has been proposed24 by replacing c 0 )2 1 by 1 ( c 0 /c0)2 to recover the exact dispersion (c0/ relation in the heterogeneous case. Both transformations will be operated. Figure 1 shows the comparison of the exact and approximate dispersion relations in the cases M ¼ 0, M ¼ 0.1, M ¼ 0.3 and M ¼ 0.5. Only for the highest Mach number, with values much larger than the contemplated atmospheric cases, can the two curves be distinguished from one another. One of the advantages of the present Eq. (10) compared to parabolic approximations is that the O(M2) precision (at least when the flow is uniform) is kept for any propagation angle. Moreover, the equation also admits evanescent waves for large values of the transverse wave numbers, while all transverse wave numbers will be artificially propagated in the parabolic case. Equation (17) is simpler but less precise, as (i) it is of order one O(M) only with respect to the ambient flow, (ii) the parabolic approximation restricts its application to a propagation direction close to the axial one (standard parabolic approximation), (iii) all the wave numbers are propagating (parabolic approximation). IV. MARCHING ALGORITHM AND SPLIT-STEP

Equation (11) is solved numerically by a marching algorithm in the x direction. Horizontal axis Ox, vertical axis Oz and retarded time axis Os are regularly discretized respectively [x1, xNx] with horizontal step Dx, [z1, zNz] with vertical

(14)

c 0 the wave number and (kx ¼ kx/k0, kz ¼ kz/k0) with k0 ¼ x0/ the dimensionless components of the wavevector. Writing c 0 , the dispersion relation for the exact convected V0x ¼ M wave equation in a uniform flow is 2 kx þ

2 2M kz 1 þ ¼ 0: k x 1 M2 1 M2

(15)

Equation (10) leads to the approximate dispersion relation, 2 2 (16) kx þ 2Mkx þ kz 1 ð1 þ M2 Þ ¼ 0: At last parabolic equation Eq. (13) yields 2562


FIG. 1. Dispersion relation in a uniform flow for various Mach numbers indicated by figures. Light dashed-dottted line: standard parabolic approximation Eq. (17). Solid lines: exact relation Eq. (15). Thick dashed lines: approximate relation Eq. (16). Gray dashed lines: backward propagating waves correctly described by Eq. (16), but ignored in the numerical framework. Gallin et al.: Propagation of shock waves in flows


step Dz and [s1, sNs] with time step Ds. Following the natural separation of the equation into a homogeneous part, a heterogeneous one and a nonlinear one, a split-step algorithm is used, the coupling between each effect being ensured by the common lhs of Eq. (11). The first-order split-step algorithm from distance xi to distance xiþ1 is written formally as D N ð Þ /ðxiþ1 ; z; sÞ ¼ /H Dx 䊊 /Dx 䊊 /Dx ðxi ; z; sÞ þ O Dx ;

(18)

where notation /O Dx (with operator O being any of the operator N , D, or H) stands for the numerical solution of the “sub-equation” @ 2//@x@s ¼ [O]/ over the step Dx. The ordering of equations to be solved is symbolized by 䊊 from right to left in the rhs of Eq. (18). The split-step as introduced in Eq. (18) is of first order with respect to the marching step Dx. In the present implementation, we use a more precise symmetric second-order split-step method (Fig. 2) first proposed by Strang,31 D H D N /ðxiþ1 ; z; sÞ ¼ /N Dx=2 䊊 /Dx=2 䊊 /Dx 䊊 /Dx=2 䊊 /Dx=2 ðxi ; z; sÞ

þ OðDx2 Þ:

(19)

The choice on the order of the half-substeps is justified in the next paragraph. One can denote that the last halfnonlinear substep at the i ! i þ 1 iteration can be merged with the first half-nonlinear substep at the following i þ 1 ! i þ 2 iteration. Thus only diffraction operator D is effectively split into two half-substeps.

Each of the equations @ 2//@x@s ¼ [O]/ is discretized numerically in an appropriate way. The nonlinear sub-equation (O ¼ N ) is solved in the time domain for a better handling of shock waves, using an analytical solution given by the Burgers-Hayes method,30,32,33 " 2 # bDx @/ðxi ; z;hÞ N ; /Dx ðxi ; z; sÞ ¼ max /ðxi ; z; hÞ @h 0 c30 q bDx @/ðxi ; z; hÞ : @h 0 c30 q

FIG. 2. Numerical scheme. J. Acoust. Soc. Am., Vol. 135, No. 5, May 2014

The homogeneous linear wave equation (O ¼ D) is also solved quasi-analytically by a one-way angular spectrum method34 in Fourier space (kz, x). Defining the spatial Fourier transform in the vertical direction, ð þ1 ð x; kz ; sÞ ¼ p1ffiffiffiffiffiffi /ðx; z; sÞeikz z dz; (22) / 2p 1 one gets ^ ^ 2ix d/ d2 / ^ ¼ 0; kz2 / þ 2 c0 dx dx

(23)

whose solution is ^ ðx þ Dx; x; k Þ / i z

V. NUMERICAL DISCRETIZATION

s ¼ h

Indeed, Eq. (20) is the implicit Poisson solution written for the potential. In the case of shock waves, that solution is multivalued, and the physically admissible solution for the potential is the maximum value, in order to satisfy the second principle of thermodynamics through shocks. The linear sub-equations (O ¼ D and O ¼ H) are solved in the frequency domain. The ordering of the operators in the second-order Strang split-step Eq. (19) is chosen so as to minimize the number of time Fourier tranform. The latter is defined as ð þ1 ^ ðx; z; xÞ ¼ p1ffiffiffiffiffiffi / /ðx; z; sÞeixs ds: (21) 2p 1

(20)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ 2 ðxi ; x; kz Þexp iDx k0 kz2 k0 : ¼/

(24)

Solution Eq. (24) is exact but assumes there is no backscattered field, induced for instance by the inhomogeneous part of the equation. Hence a one-way approximation is indeed performed also for the homogeneous part of the equation. Contrary to the heterogeneous part, that one-way approximation is not achieved by parabolic approximations, but numerically by selecting only forward propagating waves in the angular spectrum solution. However, the neglecting of the backscattered field is the major approximation (one-way approximation) of the present method. The previous two steps are strictly identical to the case of nonlinear wave propagation in an inhomogeneous but non-moving medium.24 Sub-equation @ 2//@x@s ¼ [H]/ is solved numerically over step Dx in the frequency domain. In order to keep the overall second-order numerical precision of the Strang splitting, for each frequency x, the equation is solved by an implicit Crank-Nicolson numerical scheme. Transverse first and second-order derivatives with respect to the vertical variable z are discretized by standard secondorder centered finite differences. The use of half points for the discretization of ½V 1z ð@/=@zÞ and ½V 2xz ð@ 2 /=@z@xÞ allows the implicit linear system to remain tridiagonal. The values of the ambient medium (c0, q0, V0x) at half grid points are interpolated at second order. For instance, one has V0x(zjþ1/2) ¼ [V0x(zj) þ V0x(zjþ1)]/2 þ O(Dz2). Gallin et al.: Propagation of shock waves in flows

2563


Free space boundaries are modeled by the same nonphysical absorbing layers as presented by Dragrau et al.24 These are written as an attenuation term in the heterogeneous part of Eq. (11). Note that physical bulk absorption, though not taken into account yet, could be introduced here in a similar way. The presence of a plane and rigid ground is also implemented in a different way depending on the physical effect. In the diffraction step the ground is replaced by a source image: this implies a doubling of the transverse axis Oz before applying the Fourier transform toward the space (kz, x). In the heterogeneous terms, the ground is modeled by imposing a normal velocity equal to zero in the finite differences scheme. No boundary conditions are necessary in the nonlinear part of the split-step algorithm since it is reduced to the inviscid Burgers’ equation. VI. VALIDATION TESTS

The proposed procedure to solve Eq. (11) is validated in the linear case of acoustic modes in a two-dimensional waveguide, for various types of shear flows with Mach number M(z)ex. The heterogeneous part has already been validated in a previous work.24 Several reasons justify the validation strategy. First, the numerical modes can be compared to the ones obtained from the exact Lilley equation.27 Second, because a mode is characterized only by its transverse profile and its phase speed that are kept unchanged during propagation, a numerical simulation over a very long propagation distance is a challenging test to check the properties of the numerical algorithm in terms of energy conservation and numerical dispersion. In the modal case, any numerical error generated by the scheme is always accumulated with no possibility of compensation. Such configurations can thus be interpreted as “worst cases” and provide a good estimation of the limits of our algorithm. One difficulty however in this process is to adapt the angular spectrum method for the diffraction step, to the guided case. Indeed, the method implies a spatial periodicity in the z direction. To make this compatible with the waveguide, we will consider a guide with pressure release boundary conditions /(z ¼ L) ¼ /(z ¼ L) ¼ 0 on both lateral boundaries of the waveguide of total width 2L. Moreover, we will consider only antisymmetric modes (modes of even order 2, 4, 6,…). A property of antisymetric modes is that their periodic prolongation with period 2L is C1 , while the periodic prolongation of symmetric modes (modes of odd order 1, 3, 5,…) have first-order derivatives discontinuous at z ¼ 6(2 n þ 1)L, hence are not adequate for the present numerical method. For a given frequency x0, a modal solution is sought as pa(x, z, s) ¼ P(z) exp [ik0((kx 1) x c0 s)]. Inserting it into Eq. (10) where the nonlinear effects are dropped, leads to the following eigenvalue problem: 8 " # 2 2 > 0 > < d P þ 2kx M dP þ k2 1 kx þ 2Mkx PðzÞ ¼ 0 0 ½1 þ M 2 (25) dz2 ½1 þ M2 dz > > : Pðz ¼ LÞ ¼ Pðz ¼ LÞ ¼ 0; with M0 ¼ dM/dz. By comparison, the exact modal equation for Lilley’s equation is 2564


h i d2 P 2kx M0 dP 2 ð Þ2 k2 PðzÞ ¼ 0: þ k þ 1 kM 0 dz dz2 ½1 kM (26) Equation (25) is solved numerically by a shooting method to determine the frequency x0, the different wave numbers kx, and the corresponding profiles P(z) of the propagating modes. The theoretical phase speed is then defined by c/ ¼ x0/kx. For each tested mode, the profile found by the shooting mode applied to Eq. (25) is used as a boundary condition at x ¼ 0 for the propagation software. Then the mode is numerically propagated over a distance of a thousand wavelengths 1000kxex (with kx ¼ 2p/kx a characteristic wavelength associated with the component along axis Ox of the wavevector). The huge distance of a thousand wavelengths is chosen purposely to test the numerical dispersion and dissipation. The discretization mesh is chosen with eight points per wavelength in the axial direction ex. An oversampled mesh is used along the transverse axis Oz (256 points per wavelength) and the time axis Os (128 points per period) in order to reach the numerical convergence of the mesh. The validation tests investigate the numerical conservation of the mode profiles P(z) during the propagation, and the good agreement of the numerical phase speed c/ is compared to the theoretical one. The presented results focus on two representative flows: a Poiseuille flow of Mach number M(z) ¼ M0(1.0 z2/L2) and a more sharply sheared flow, chosen to be uniform at the center of the wave guide: M(z) ¼ M0 for z僆[( 1 þ a)L, (1 a)L] and with shear gradients localized near the boundaries: M(z) ¼ M0 sin[p(1 z/L)/(2a)] for z 僆 [L(1 a), L] (and symmetrically for the other side). The relative thickness of the flow “boundary layer” is chosen equal to a ¼ 0.1. We consider both cases of positive and negative characteristic Mach number M0 ¼ 60.1. The numerical algorithm is tested for each of the four flows with frequencies x covering the range c 0 ¼ 20. Around the lower limit one from xL/ c 0 ¼ 5 to xL/ gets close to the cutoff frequency, and the numerical scheme shows some numerical discrepancies. For larger frequencies, the corresponding wave numbers kx get close to k0, the modes propagate mostly in the axial direction and algorithm does not show large changes when increasing the frequency. As an example, the profiles of the second mode (first antisymetric mode) for a Poiseuille flow of characteristic Mach numbers M0 ¼ 0.1 and M0 ¼ 0.1 are presented on Fig. 3 for xL/ c 0 ¼ 18.5. The agreement between the modes from the Lilley equation and equation Eq. (25) is excellent. Concerning the numerical algorithm, clearly the shapes of the mode profiles are conserved during the numerical propagation over 1000kxex. However, in the case of the positive Mach number, a small dissipative effect is also present: after a propagation of 1000kxex, the relative error of the maximum of the profile compared to the one at the initial step is 10%. The evolution of this error is illustrated in Fig. 4. One can observe a small but cumulative dissipative effect. For the negative Mach number, no dissipation effect is observed. Figure 5 shows in a synthetic way the relative errors of the amplitudes of the modes at the final propagation stage (a thousand wavelengths in the axial direction) for the Gallin et al.: Propagation of shock waves in flows


FIG. 3. Pressure transverse profile of mode 2 for the frequency xL/ c 0 ¼ 18.5 for a Poiseuille flow of positive Mach number M0 ¼ 0.1 (upper plot) and of negative Mach number M0 ¼ 0.1 (lower plot). Solid lines with dots: modes for exact Lilley equation. Solid lines with cross: modes of Eq. (25). Dashed lines: numerical modes after a propagation over 1000kxex. The profiles are normalized by the maximum amplitude of the Lilley mode.

different test cases (Poiseuille or sharply sheared flow, positive or negative Mach number, full frequency range). It thus appears that the numerical scheme does not conserve energy perfectly, and that some significant amplitude deviations may appear for low frequency modes propagating with a relatively large angle relative to the axis. In the Poiseuille flow case modes with relative error greater than 10% correspond to propagation angles greater than 30 . The case of the boundary layer flow is more sensible with amplification appearing for angles greater than 20 . This error is directly associated with the V 1z and V 2xz terms in Eq. (12). If these terms are dropped, energy conservation is ensured. This is linked to the fact that these terms are associated with uneven derivatives of the potential [see Eq. (9)], hence they may introduce dissipation or anti-dissipation. However, this effect is strongly reduced by using the conservative form Eq. (10) rather than the non-conservative form Eq. (8). For the modal propagation, a very small numerical dissipation or antidissipation error (typically here at worst 0.1% per wavelength) is systematically reported. However, for realistic outdoor propagation, because of energy leakage out of the possible waveguides, and because of atmospheric variability over long distances, such cumulative numerical errors are unlikely to occur.

FIG. 4. Normalized amplitude versus propagation distance for mode 2 of frequency xL/ c 0 ¼ 18.5 in the case of a Poiseuille flow with positive Mach number M0 ¼ 0.1. The relative error is 10% at the end of the propagation (1000kxex). J. Acoust. Soc. Am., Vol. 135, No. 5, May 2014

FIG. 5. Relative error of the amplitude (relative to the initial amplitude) for the Poiseuille flow (upper plot) and the boundary layer flow (lower plot) after a numerical propagation of 1000kxex. Solid lines: mode 2. Dashed lines: mode 4. (þ) signs: flow of positive Mach number M0 ¼ 0.1. () signs: flow of negative Mach number M0 ¼ 0.1.

Similarly, Fig. 6 shows in a synthetic way the relative errors of the phase speed of the modes at the final propagation stage. Here, in all cases, the error appears to remain always very small, less than 1%. Hence the numerical algorithm turns out to introduce very little numerical dispersion, even for modes propagating with a relatively large angle relative to the axis. VII. BLAST WAVE PROPAGATION WITHIN A BOUNDARY LAYER FLOW

An example of application is now provided in the context of atmospheric sound propagation. A line source is considered, of total height 170 m (e.g., 25k with k ¼ 6.8 m for a speed of sound equal to 340 m/s and a main signal frequency equal to 50 Hz). The line source, which extends down to the ground level, is tilted with an angle of 4 deg relative to the normal to the ground. The upper part of the segment (one quarter of the total length) is tapered with a decreasing cosinusoidal shape to smooth geometrical discontinuities. Each point along the source emits a blast wave of peak overpressure P0 ¼ 400 Pa and peak frequency f0 ¼ 50 Hz, with a pressure waveform described by35 P0 (1 t/td) et/td. The value td is chosen equal to 0.003 s: this corresponds to a peak frequency of 50 Hz in the spectrum of the waveform. This configuration aims at describing in a very primitive way the Gallin et al.: Propagation of shock waves in flows

2565


Three cases are investigated numerically: the case without wind, the case of a positive wind (“downward” propagation of sound in the same direction as wind), and at last the case of a negative wind (“upward” propagation of sound opposite to wind). The blast wave is numerically propagated away from the source over a distance of x ¼ 200k, which corresponds to a distance of 1360 m. The propagation axis Ox is meshed with eight points per wavelengths k, this corresponds to a total of 1600 points. The transverse domain Oz is 40k high and is meshed with 4096 points in total: all of them are used in the diffraction step to represent the Fourier transform of the physical space and its image. Only 2049 points are used in the finite differences scheme and the nonlinearstep, this corresponds to the physical space and the ground. 1024 points are used to discretize the moving time window. The thickness of the absorbing boundary layer is 10k with a damping factor of 1.4.24 Figures 7 and 8 present for each of the three cases the pressure fields as functions of altitude and retarded time, at distances respectively 50k (340 m) and 200k (1360 m). Figures 9–11 display the corresponding overpressure profiles at ground level z ¼ 0 at distances 50k (340 m), 150k (1020 m), and 200k (1360 m). Each wave profile at ground

FIG. 6. Relative error of the phase speed (relative to the theoretical one) for the Poiseuille flow (upper plot) and the boundary layer flow (lower plot) after a numerical propagation of 1000kxex. Solid lines: mode 2. Dashed lines: mode 4. (þ) signs: flow of positive Mach number M0 ¼ 0.1. () signs: flow of negative Mach number M0 ¼ 0.1.

lower part of a lightning channel which can develop during a return stroke, according to the line source model by Few.36 When no wind is present, the tilt angle is chosen so that nonlinear effects lead to a von Neumann reflection of the shock wave over the ground.37 The ground is considered to be perfectly plane and rigid. A wind flow is introduced, using the model proposed by Ostashev et al.38 based on the Monin-Obukhov similarity theory for the atmospheric boundary layer. Parameters correspond to one of the six benchmark meteorological regimes proposed by Ostashev and Wilson,39 the so-called “mostly sunny day with strong wind” case, with a friction velocity of 0.7 m/s, a surface roughness length of 1 cm (flat ground with short grass), a reference temperature near the ground of 25 C, and a surface heat flux of 200 W/m2, hence a negative (unstable stratification during day) Obukhov length of 151 m. The wind amplitude (see Fig. 4 in Ref. 38) varies from 0 m/s on the ground to a maximum of 16 m/s in altitude at 500 m, with a sharp gradient within the first 10 m above the ground. The maximum value corresponds to a flow Mach number M0 0.05, thus satisfying the low Mach assumption used in the present work. Since the present work focuses on the new flow terms of Eq. (10), the model proposed by Ostashev et al. for the temperature gradient is not used here. 2566


FIG. 7. Left column: overpressure field in pascals (gray color bar) in the domain (Os, Oz) after a nonlinear propagation over a distance of x ¼ 50k (340 m). Right column: pressure profiles at the altitudes 10, 50, 100, and 150 m. Vertical segments on the right represent a reference amplitude of 200 Pa. First line: no wind. Second line: positive strong wind. Last line: negative strong wind. Gallin et al.: Propagation of shock waves in flows


FIG. 8. Same plots as Fig. 7, but after a propagation over a distance of x ¼ 200k (1360 m).

level is computed either in the linear and nonlinear cases, and the pressure spectra are also systematically shown. Figures 7 and 8 show the strong effect of the boundary layer flow. Close to the source at 340 m (Fig. 7), the wavefield is dominated by the incident wavefront, the reflected wavefront having time to build up only near the ground. The differences in arrival times in Fig. 7 are due to the flow convection. In the no wind case, the nonlinear effect leads to a von Neumann reflection on the ground,37 with a very tiny Mach stem. Hence, at a 10 m altitude, one clearly sees only one shock going from Fig. 7 to Fig. 8. The ground signal shows the waveform has little evolved from the initial Kinney profile, even though the nonlinear process of energy transfers toward high and low frequencies has already begun in Fig. 9. In the case of positive wind, there is no nonlinear reflection, and at a 10 m altitude, two shocks (one incident, one reflected) of same amplitude are clearly visible. At higher altitudes, the wind is almost uniform, and its only effect is to advance the signal, but the waveforms with or without wind are almost identical (comparing the top and middle of Fig. 8). At ground level z ¼ 0, one observes a broadening of the time waveform (comparing Figs. 9 and 10). Indeed, because of the sharp wind gradient near the ground, the high frequencies (typically those above 50 Hz) tend to keep near the ground through a wave guide effect, which is frequency dependent. This dispersion is also visible in the oscillations of the spectrum above 80 Hz. In the case J. Acoust. Soc. Am., Vol. 135, No. 5, May 2014

FIG. 9. Case without wind at ground level z ¼ 0. Left column: overpressure on the ground (in pascals). Right column: corresponding spectrum in dB (reference: maximum of the spectrum from the initial signal at x ¼ 0). Propagation over a distance of x ¼ 50k (first line), x ¼ 150k (second line), and x ¼ 200k (third line). Solid lines: nonlinear case. Dashed lines: linear case. Light lines with dots: source spectrum at x ¼ 0.

of negative wind, the signals tend to turn up by upward refraction in Figs. 7 and 8. Near the ground in Fig. 11, the highest frequencies decay exponentially through the process of creeping wave propagation14,40 and the shock begins to smear out. The high frequency decay is very visible on the spectrum of the ground signal. The incident and reflected shocks merge not at ground level but in altitude, here at about 80 m in Fig. 8. Far from the source at 1360 m (Fig. 8), the reflected wavefront has fully built up. Once again, flow convection explains the strong differences between the arrival times. In the no wind case, the wavefront curvatures at the highest altitudes are due to the finite size of the source and diffraction at its upper edge. Incident and reflected signals are well observable and are almost similar (geometrical reflection). Near the ground, both of them begin to interfere and it is hard to tell whether nonlinear von Neumann reflection takes place or not. The ground signal in Fig. 9 clearly shows a difference between the linear and nonlinear case, this one showing (i) a lower peak amplitude (dissipation through shock waves), (ii) an earlier shock arrival (nonlinear speed of sound is greater than the linear one), hence a lengthening of the signal with a shift of the spectrum peak toward low frequencies (now around 30 Hz instead of 50 Hz), and (iii) oscillations in the high frequency part of the spectrum (the waveform progressively evolves from the blast-wave profile Gallin et al.: Propagation of shock waves in flows

2567


FIG. 10. Same as Fig. 9 in the case of a positive strong wind.

FIG. 11. Same as Fig. 9 in the case of a negative strong wind.

to an N-wave profile). In the case of positive wind, the part of the field at high altitudes, well above the waveguide effect, looks quite similar in Fig. 8 to the no wind case, except with earlier arrival times. On the contrary, between the ground and around 70 m, the waveguide effect due to downward propagation in a strongly sheared flow leads to strong dispersion effects. At a 50 m altitude, one still has clearly an incident and a reflected signal, but this reflected one is about twice as large. Because it is of larger amplitude, nonlinear effects are more important and it therefore looks closer to an N-wave. Compared to the no wind case, the signal at 0 and 10 m is about three times longer. The head part of the ground signal looks quite similar to an N-wave with two shocks, showing that nonlinear effects are quite strong. This is confirmed by comparing in Fig. 10 the ground signals computed with or without nonlinear effects, which display large differences (also in the long dispersive tail of the signal). Indeed, because of the wave guide effect due to wind shear, the acoustic energy tends to be more confined near the ground, and hence nonlinear effects are more important. To be noticed also Fig. 10 is the strong shift of the first trough in the spectrum: 80 Hz at 340 m (50k), 45 Hz at 1020 m (150k), and 20 Hz at 1360 m (200k). In the case of a negative wind in Fig. 11, the signal attenuation at ground level in the shadow zone is strong, and few of the frequencies larger than 50 Hz remain. The signal intensity is higher in Fig. 8 at altitudes between 50 and 100 m because of upward refraction. Because of the higher amplitude at these altitudes, the nonlinear effects are stronger, and the waveform looks quite

FIG. 12. Amplitude of the acoustic pressure field normalized to [0,1] for a given frequency (first line: 10 Hz, second line: 50 Hz, third line: 100 Hz). Cases of the positive (left column) and negative (right column) strong wind.

2568




close to an N-wave. Only well above the shadow zone limit can one clearly see the incident and reflected signals. Figure 12 presents maps of the pressure fields |^ p a (x, z, x)| at different frequencies f ¼ x/2p (10, 50, and 100 Hz) in the computational plane (Ox, Oz). In all cases, one can observe the interference fringes between the direct field radiated by the source, and the field reflected by the ground. The diffraction of the tapered upper part of the line source is also visible for all frequencies. The lowest frequencies fill all the space. For them there is no shadow zone, and there are few differences between cases with positive wind and those with negative wind. Such frequency dependent results could not be obtained by a ray tracing approach. In the case of a positive strong wind, at the higher frequencies 50 and 100 Hz, a wave guide effect is very visible near the ground (mostly within the first 10 m). There, the acoustic rays bounce between the rigid ground and the shear part of the wind flow when traveling away from the source. For the case of a negative wind flow, shadow zones are clearly defined and the upward refraction of the field away from the rigid ground is visible.

described by a ray tracing approach. Future research will aim at parallelizing the algorithm for the development of a three-dimensional version. Thermoviscous absorption and molecular relaxation will also have to be included. The numerical algorithm will be used to simulate atmospheric shock wave propagation in more realistic configurations, in terms of both meteorological data and source description. APPENDIX: OPERATORS DESCRIBING HETEROGENEITIES AND FLOWS

The operators [H0 ], [H1x ], [H2zz ], [H3xzz ], [V 1z ], [V 2xz ] introduced in Eq. (12) are ½H0 / ¼

þ ½H1x


2 V0x ðzÞ @ 2 / V0x ðzÞ @ 2 / ; 2 2 c0 @s 2 c 30 @s2

(A1a)

@/ @q0 @/ V0x ðzÞ @ 2 / c0 ¼ ; 2q0 ð x; zÞ @x @x @x c0 @x@s

(A1b)

@2/ @ 2 / V 2 ðzÞ @ 2 / c0 @ 2 / ¼ V0x ðzÞ 2 þ 0x þ ; (A1c) 2 @z @z 2 c 0 @z2 2 @z2 ðs 3 ^ @3/ @ / 0 ½H3xzz ¼ c V z ds ; (A1d) ð Þ 0 0x 2 @x@z @x@z2 c0 @/ @ @/ @ @/ ½V 1z ¼ V0x ðzÞ q0 ðx; zÞ ; 2q0 ðx; zÞ @z @z @z @z @z

VIII. CONCLUSION

The objective of the present work was to develop a oneway numerical approach to simulate nonlinear shock wave propagation in conditions valid for atmospheric flows. In particular, to be able to simulate wide angle propagation effects due for instance to extended sources, sources at altitude, or wave scattering effects by atmospheric inhomogeneities, a formulation was developed that goes beyond the parabolic approximation but keeps a one-way marching approach which is numerically efficient. It generalizes a method previously developed and validated in cases without flow. Flow motion is taken into account at an intermediate order between M and M2, where M is the ambient flow Mach number. In particular, convection effects are described at the highest M2 order. The resulting dispersion relation of this new approximate formulation turns out to be very precise for all angles (þ/90 ) and for flow Mach numbers up to 60.3. For sheared flows, the numerical method was validated by comparison with solutions of the exact Lilley equation in the case of guided modes. The comparison confirms that the proposed FLHOWARD numerical algorithm simulates the dispersion relations with good precision. As a counterpart, the algorithm does not turn out to be perfectly energy preserving. This comes out from two identified terms in the model equation. Further developments to improve the algorithm for these two terms would be useful. An example of blast wave propagation within a boundary layer flow shows that the numerical algorithm reproduces well-known features of atmospheric propagation, such as strong wave decays in shadow zones, upward refraction or wave guide propagation near the ground in the sheared layer. Moreover significant nonlinear effects are outlined, especially in the wave guide case with energy transfers toward the highest and lowest parts of the frequency spectrum, signal broadening, and formation of complex waveforms. The low frequency part of the signal is also shown to behave in a way that could not be

c00 ðx; zÞ @ 2 / 1 @q0 @/ þ 2 2 2q0 ðx; zÞ @x @s c0 @s

½H2zz

(A1e) ½V 2xz

@2/ ¼ c0 @z@x

ðs

@ @2/ V0x ðzÞ @z @z@x

ds0 :

(A1f)

1

L. Evers and H. W. Haak, “The characteristics of infrasound, its propagation and some early history,” in Infrasound Monitoring for Atmospheric Studies, edited by A. Le Pichon, E. Blanc, and A. Hauchecorne (Springer, Dordrecht, 2010), pp. 3–27. 2 K. L. Gee, V. W. Sparrow, M. M. James, J. M. Downing, C. M. Hobbs, T. B. Gabrielson, and A. A. Atchley, “The role of nonlinear effects in the propagation of noise from high power jet aircraft,” J. Acoust. Soc. Am. 123, 4082–4093 (2008). 3 O. Gainville, P. Blanc-Benon, E. Blanc, R. Roche, C. Millet, F. L. Piver, B. Despres, and P. F. Piserchia, “Misty picture: A unique experiment for the interpretation of the infrasound propagation from large explosive sources,” in Infrasound Monitoring for Atmospheric Studies, edited by A. Le Pichon, E. Blanc, and A. Hauchecorne (Springer, Dordrecht, 2010), pp. 575–598. 4 D. O. Revelle, “On meteor generated infrasound,” J. Geophys. Res. 81, 1217–1229, doi:10.1029/JA081i007p01217 (1976). 5 A. LePichon, E. Blanc, and D. P. Drob, “Probing high-altitude winds with infrasound,” J. Geophys. Res. 110, D20104, doi:10.1029/2005JD006020 (2005). 6 H. E. Bass, “The propagation of thunder through the atmosphere,” J. Acoust. Soc. Am. 67, 1959–1966 (1980). 7 T. Farges and E. Blanc, “Characteristics of infrasound from lightning and sprites near thunderstorm areas,” J. Geophys. Res. 115, A00E31, doi:10.1029/2009JA014700 (2010). 8 J. D. Assink, L. G. Evers, I. Holleman, and H. Paulssen, “Characterization of infrasound from lightning,” Geophys. Res. Lett. 35, L15802, doi:10.1029/2008GL034193 (2008). Gallin et al.: Propagation of shock waves in flows

2569


9

S. de Larquier and V. P. Pasko, “Mechanism of inverted-chirp infrasonic radiation from sprites,” Geophys. Res. Lett. 37, L24803, doi:10.1029/2010GL045304 (2010). 10 K. J. Plotkin, “State of the art of sonic boom modeling,” J. Acoust. Soc. Am. 111, 530–536 (2002). 11 D. Norris, R. Gibson, and K. Bongiovanni, “Numerical methods to model infrasonic propagation through realistic specifications of the atmosphere,” in Infrasound Monitoring for Atmospheric Studies, edited by A. Le Pichon, E. Blanc, and A. Hauchecorne (Springer, Dordrecht, 2010), pp. 541–573. 12 R. Blumrich, F. Coulouvrat, and D. Heimann, “Meteorologically induced variability of sonic-boom characteristics of supersonic aircraft in cruising flight,” J. Acoust. Soc. Am. 118, 707–722 (2005). 13 A. Le Pichon, M. Garces, E. Blanc, M. Barthelemy, and D. P. Drob, “Acoustic propagation and atmosphere characteristics derived from infrasonic waves generated by Concorde,” J. Acoust. Soc. Am. 111, 629–641 (2002). 14 F. Coulouvrat, “Sonic boom in the shadow zone: A geometrical theory of diffraction,” J. Acoust. Soc. Am. 111, 499–508 (2002). 15 R. Marchiano, F. Coulouvrat, and R. Grenon, “Numerical simulation of shock wave focusing at fold caustics, with application to sonic boom,” J. Acoust. Soc. Am. 114, 1758–1771 (2003). 16 V. W. Sparrow and R. Raspet, “A numerical method for general finiteamplitude wave propagation in 2 dimensions and its application to spark pulses,” J. Acoust. Soc. Am. 90, 2683–2691 (1991). 17 S. D. Pino, B. Despres, P. Have, H. Jourdren, and P.-F. Piserchia, “3D finite volume simulation of acoustic waves in the Earth atmopshere,” Comput. Fluids 38, 765–777 (2009). 18 G. Hanique-Cockenpot, “Etude numerique de la propagation non lineaire des infrasons dans l’atmosphe`re” (“Numerical study of nonlinear propagation of infrasounds in the atmosphere”), Ph.D. thesis, Ecole Centrale de Lyon, 2011. 19 V. P. Kuznetsov, “Equations of nonlinear acoustics,” Sov. Phys. Acoust. 16, 467–470 (1970). 20 M. V. Aver’yanov, V. A. Khokhlova, O. Sapozhnikov, P. Blanc-Benon, and R. Cleveland, “Parabolic equation for nonlinear acoustic wave propagation in homogeneous moving media,” Acoust. Phys. 52, 623–632 (2006). 21 M. Averiyanov, P. Blanc-Benon, R. Cleveland, and V. Khokhlova, “Nonlinear and diffraction effects in propagation of N-waves in randomly inhomogeneous moving media,” J. Acoust. Soc. Am. 129, 1760–1772 (2011). 22 M. Averiyanov, S. Ollivier, V. Khokhlova, and P. Blanc-Benon, “Random focusing of nonlinear acoustic N-waves in fully developed turbulence: Laboratory scale experiment,” J. Acoust. Soc. Am. 130, 3595–3607 (2011). 23 F. Coulouvrat, “New equations for nonlinear acoustics in a low Mach number and weakly heterogeneous atmosphere,” Wave Motion 49, 50–63 (2012).

2570


24

F. Dagrau, M. Renier, R. Marchiano, and F. Coulouvrat, “Acoustical shock wave propagation in a weakly heterogeneous medium: A numerical simulation beyond the parabolic approximation,” J. Acoust. Soc. Am. 130, 20–32 (2011). 25 T. Christopher and K. Parker, “New approaches to nonlinear diffractive field propagation,” J. Acoust. Soc. Am. 90, 488–499 (1991). 26 M. Lothon, D. Lenschow, and S. Mayor, “Coherence and scale of vertical velocity in the convective boundary layer from a Doppler lidar,” Boundary-Layer Meteorol. 121, 521–536 (2006). 27 G. M. Lilley, “The generation and radiation of supersonic jet noise,” Technical Report, Air Force Aero-Propulsion Laboratory, AFAPL-TR-7253 (1972). 28 J. F. Claerbout, Fundamentals of Geophysical Data Processing (McGraw-Hill, New York, 1976), pp. 194–207. 29 V. E. Ostashev, D. Juve, and P. Blanc-Benon, “Derivation of a wide-angle parabolic equation for sound waves in inhomogeneous moving media,” Acust. Acta Acust. 83, 455–460 (1997). 30 F. Coulouvrat, “A quasi exact shock fitting algorithm for general nonlinear progressive waves,” Wave Motion 46, 97–107 (2009). 31 G. Strang, “On the construction and comparison of difference schemes,” SIAM J. Numer. Anal. 5, 506–517 (1968). 32 J. M. Burgers, “Further statistical problems connected with the solution of a simple nonlinear partial differential equation,” Proc. Kon. Ned. Akad. Wetens., Ser. B. 57, 159–169 (1954). 33 W. D. Hayes, R. C. Haefeli, and H. E. Kulsrud, “Sonic boom propagation in a stratified atmosphere with computer program,” Technical Report, NASA CR-1299 (1969). 34 J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, New York, 1968), pp. 45–50. 35 F. G. Friedlander, “The Diffraction of sound pulses. 1. Diffraction by a semi-infinite plane” Proc. R. Soc. London, Ser. A 186, 322–344 (1946). 36 A. A. Few, “Power spectrum of thunder,” J. Geophys. Res. 74, 6926–6934, doi:10.1029/JC074i028p06926 (1969). 37 S. Baskar, F. Coulouvrat, and R. Marchiano, “Nonlinear reflection of grazing acoustical shock waves: Unsteady transition from von Neumann to Mach to Snell-Descartes reflections,” J. Fluid Mech. 575, 27–55 (2007). 38 V. E. Ostashev, M. V. Scanlon, D. K. Wilson, and S. N. Vecherin, “Source localization from an elevated acoustic sensor array in a refractive atmosphere,” J. Acoust. Soc. Am. 124, 3413–3420 (2008). 39 V. E. Ostashev and D. K. Wilson, “Relative contributions from temperature and wind velocity fluctuations to the statistical moments of a sound field in a turbulent atmosphere,” Acust. Acta Acust. 86, 260–268 (2000). 40 A. D. Pierce, Acoustics: An Introduction to its Physical Principles and Applications (Acoustical Society of America, New York. 1989), pp. 469–478.



Two-dimensional simulation of optical wave propagation through atmospheric turbulence.

Shock wave propagation along constant sloped ocean bottoms.

Rescaled Local Interaction Simulation Approach for Shear Wave Propagation Modelling in Magnetic Resonance Elastography.

Assessing viscoelasticity of shear wave propagation in cervical tissue by multiscale computational simulation.

finite element method.

Virtual cardiac monolayers for electrical wave propagation.

The robustness of the weak selection approximation for the evolution of altruism against strong selection.

DebtRank: A Microscopic Foundation for Shock Propagation.

Shock-induced wave propagation over porous and fractured borehole zones: theory and experiments.

Evanescent wave-based particle tracking velocimetry for nanochannel flows.

MEANS: python package for Moment Expansion Approximation, iNference and Simulation.

Self-controlled wave propagation in hyperelastic media.

Anisotropic wave propagation in nematic liquid crystals.

Acoustic precursor wave propagation in viscoelastic media.

Numerical simulation of the nonlinear ultrasonic pressure wave propagation in a cavitating bubbly liquid inside a sonochemical reactor.

Magnetohydrodynamic channel flows with weak transverse magnetic fields.

Nonlinear guided wave propagation in prestressed plates.

Immersive experimentation in a wave propagation laboratory.

Extracorporeal shock wave lithotripsy for biliary stones.

Topical anaesthesia for extracorporeal shock wave lithotripsy.

Extracorporeal shock wave lithotripsy for radiolucent stones.

Extracorporeal shock wave lithotripsy for gallstones.

Ring discretization of the wave spectrum for sea surface simulation.

An experimental study on the ultrasonic wave propagation in cancellous bone: waveform changes during propagation.