A full-field simulation methodology for sonic boom modeling on adaptive Cartesian cut-cell meshes

Journal Pre-proof A full-field simulation methodology for sonic boom modeling on adaptive Cartesian cut-cell meshes

R. Yamashita, L. Wutschitz, N. Nikiforakis

PII:

S0021-9991(20)30045-0

DOI:

https://doi.org/10.1016/j.jcp.2020.109271

Reference:

YJCPH 109271

To appear in:

Journal of Computational Physics

Received date:

17 September 2019

Revised date:

14 January 2020

Accepted date:

16 January 2020

Please cite this article as: R. Yamashita et al., A full-field simulation methodology for sonic boom modeling on adaptive Cartesian cut-cell meshes, J. Comput. Phys. (2020), 109271, doi: https://doi.org/10.1016/j.jcp.2020.109271.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.

Highlights • • • • •

A full-field simulation methodology for sonic boom modeling is developed. Adaptive mesh refinement and cut cell grids are used over a stratified atmosphere. A well-balanced MUSCL-Hancock scheme in a stratified atmosphere is formulated. Full-field simulation successfully reproduces the flight test for sonic boom. Simulation results agree with results of a waveform parameter method and flight test.

A full-field simulation methodology for sonic boom modeling on adaptive Cartesian cut-cell meshes R. Yamashita∗, L. Wutschitz, N. Nikiforakis Laboratory for Scientific Computing, Cavendish Laboratory, Department of Physics, University of Cambridge, UK

Abstract This paper develops a full-field direct numerical simulation methodology of sonic boom in a stratified atmosphere. The entire flow field, ranging from the near field around a supersonic body to the far field extending to the ground, is modeled by the three-dimensional Euler equations with a gravitational source term. Thus far, it has been solved using a structured grid, and an application of previous simulation to complex geometries has been limited. In this study, we realize a full-field simulation by employing the following four numerical approaches: (i) a hierarchical structured adaptive mesh refinement (AMR) method, (ii) a Cartesian cut cell method, (iii) a well-balanced finite volume method, and (iv) a segmentation method of the computational domain. A new well-balanced, MUSCL-Hancock scheme applied over Cartesian AMR and cut cell grids for a stratified atmosphere is formulated. The computational results of an oblique shock wave in a stratified atmosphere agree well with the exact solution. A full-field simulation successfully reproduces the Drop test for Simplified Evaluation of Non-symmetrically Distributed sonic boom (D-SEND) #1, conducted by the Japan Aerospace Exploration Agency (JAXA). The results of this simulation are in good agreement with those of the previous computational study, the waveform parameter method, and flight test measurements. The grid convergence study shows that the mesh size is fine enough to assess pressure signatures over the entire flow field. These results demonstrate that a full-field simulation with AMR and cut cell grids is a powerful tool for extensively analyzing three-dimensional shock wave propagation in a stratified atmosphere. Keywords: Sonic boom, Full-field simulation, Adaptive mesh refinement, Cut cell, Stratified atmosphere, Well-balanced finite volume method

1. Introduction Sonic boom arising from an aircraft flying at a supersonic speed is one of the primary problems in the development of commercial supersonic aircraft [1]. Generally, shock waves emanating from a supersonic ∗ Corresponding

author. Email addresses: [email protected] (R. Yamashita), [email protected] (L. Wutschitz), [email protected] (N. Nikiforakis)

aircraft coalesce into an N-shaped waveform (N-wave) because of the nonlinear effects over a fairly long 5

propagation distance of the waves, and an explosive-like sound occurs twice on the ground. Therefore, overland supersonic flight is currently restricted in many countries due to the sonic boom noise [2]. Low sonic boom design methods [3-5] have been developed for addressing the sonic boom problem. Flight tests conducted by NASA [6] and the Japan Aerospace Exploration Agency (JAXA) [7] demonstrated that the loudness level of sonic boom generated from low sonic boom aircraft can be significantly reduced, relative

10

to that generated from Concorde. The International Civil Aviation Organization (ICAO) has proposed an international standard to determine an acceptable level of sonic boom in overland supersonic flights. However, due to the lack of scientific data, it has been difficult to establish a clear standard over a variety of flight and atmospheric conditions from takeoff to landing. Therefore, flight tests for sonic boom were recently conducted by NASA [8-10], and the sonic boom prediction workshops [11-13] were held to evaluate accuracy

15

of the sonic boom prediction. Moreover, accurate prediction methods have been developed for investigating complex sonic boom phenomena in unsteady flights [14] and those in real atmospheric environments [15-17]. The accurate prediction of sonic boom is a challenging task due to the multi-scale nature of the phenomenon: The scale varies from millimeters around an aircraft to kilometers in a real stratified atmosphere [18]. In addition, sonic boom intensity significantly changes according to the flight pattern, such as accelera-

20

tion and maneuver, and atmospheric variability, such as atmospheric turbulence and wind [19]. Traditionally, the sonic boom prediction has been conducted by splitting an entire field into three sub-fields: the near field around a supersonic body, the far field extending to the ground, and the caustic-vicinity field where sonic boom focusing occurs. Near-field pressure waveforms have been obtained by wind tunnel experiments or computational fluid dynamics (CFD) analysis [2, 11, 12, 20, 21]. Far-field waveforms have been predicted

25

by one-dimensional analysis along a ray path, such as the waveform parameter method (WPM) [22], and the method of solving the Burgers’ equation [23, 24]. Sonic boom focusing in the caustic-vicinity field due to acceleration or maneuver has been predicted by two-dimensional analysis that solves the nonlinear progressive-wave equation [25] or the nonlinear Tricomi equation [14, 26]. The one-dimensional or twodimensional analysis described above is suitable for parametric studies that are used to design a low sonic

30

boom aircraft due to the low computational cost.

To improve state of the art techniques, a direct sim-

ulation methodology based on three-dimensional CFD analysis over the entire flow field, ranging from the near field around a supersonic aircraft to the far field extending to the ground, has been developed. An initial direct simulation under the assumption of uniform atmosphere has been performed by solving the three-dimensional Euler equations with unstructured grids [27]. Thereafter, a direct simulation method con35

sidering atmospheric stratification, known as a full-field simulation method, has been constructed to analyze sonic boom propagation through a real stratified atmosphere [28-30]. Accuracy of the full-field simulation was validated by comparison with flight test data and the WPM, which is a representative prediction method for sonic boom [29, 30]. Full-field simulation, extended to thermal nonequilibrium flow analysis, has been 2

performed to predict the rise time of sonic boom [30]. Furthermore, sonic boom cutoffs occurring due to 40

variation in atmospheric temperature with altitude have been successfully analyzed by full-field simulation [31, 32]. These results demonstrate that full-field simulation is a powerful tool for investigating complex phenomena, such as sonic boom cutoffs, that have been difficult to predict using conventional one-dimensional or two-dimensional analysis. To date, full-field simulation [28-32] has been performed using a solution-adapted structured grid in

45

boundary fitted coordinates (BFC), and thus it has been difficult to precisely capture all waves generated from complex geometries, such as three-dimensional aircraft configurations. To avoid this problem, we focus on the hierarchical structured adaptive mesh refinement (AMR) method [33] that can provide the advantage of a structured grid and enhance grid resolution only near shock waves. CFD analysis with AMR [2, 34, 35] has been conducted for predicting near-field waveforms under the assumption of uniform atmosphere.

50

However, it cannot be directly applied to sonic boom analysis in the far-field. This is because atmospheric stratification has a significant effect on sonic boom propagation in the far field, and the well-balanced finite volume method [36-38] for maintaining hydrostatic balance between the pressure divergence and gravitational source terms within discrete grid cells must be applied even in a regular Cartesian grid. In addition, fine-coarse boundaries in an AMR approach must be specially treated by considering fluctuation due to

55

waves and change due to atmospheric stratification with altitude, and the large computational cost in this method must be reduced for analyzing sonic boom propagation over distances of tens of kilometers. With these constraints in mind, the objectives of this study are to develop a full-field simulation method with AMR for three-dimensional sonic boom modeling in a stratified atmosphere and to achieve simulations in a reasonable computational time. The notable feature of this simulation method is to integrate the following

60

four numerical approaches: (i) a hierarchical structured AMR method for enhancing grid resolution around shock waves, (ii) Cartesian cut cell methods [2,34,35,39] for mesh generation around complex geometries, (iii) well-balanced finite volume methods in combination with a numerical correction method [29, 31] for maintaining hydrostatic balance in a stratified atmosphere, and (iv) a segmentation method of the computational domain [29, 31] for reducing the computational cost. This integrated simulation method includes a

65

novel, well-balanced finite volume method applied in Cartesian AMR grids in a stratified atmosphere. The simulations are validated by comparison with exact solutions, previous computational results, WPM, and flight test data. The rest of this paper is organized as follows: The computational methods, including the well-balanced finite volume method in a stratified atmosphere, the numerical correction method for preserving a freestream

70

condition in a stratified atmosphere, the AMR and cut cell methods, and the segmentation method of the computational domain, are described in section 2. An oblique shock wave generated from a simple wedge is analyzed to validate computational accuracy in a stratified atmosphere by comparison to the exact solution in section 3. In section 4, axisymmetric simulation in the near field around a supersonic 3

body is performed to demonstrate that all waves emanating from the supersonic body can be captured by 75

the AMR and cutcell grids, and two numerical schemes are compared. Thereafter, a full-field simulation is performed to demonstrate that it can reproduce the Drop test for Simplified Evaluation of Non-symmetrically Distributed sonic boom (D-SEND) #1 [7, 40] conducted by JAXA. The computational results are validated by comparison with existing computational results [29], WPM, and flight test data.

2. Computational method 80

2.1. Sonic boom modeling Sonic boom propagation for very long distances must be modeled by considering the effects of geometrical spreading, nonlinearity, and atmospheric stratification. Thus, the governing equations in Cartesian coordinates (x, y, z) are the three-dimensional Euler equations with a gravitational source term for considering a horizontally stratified atmosphere, ∂F ∂G ∂Q ∂E + + + = Sx + Sy , ∂t ∂x ∂y ∂z



ρ





ρu





ρv

           ρuv  ρu2 + p   ρu            Q =  ρv  , E =  ρuv  , F =  ρv 2 + p            ρvw  ρuw   ρw       (e + p)v (e + p)u e     0 0          ρgx   0          Sx =  0  , Sy =  ρgy  ,          0   0      ρgx u ρgy v e=



(1)



ρw

     ρuw       , G =  ρvw      ρw2 + p    (e + p)w

ρRT ρ + (u2 + v 2 + w2 ), γ−1 2

      ,     (2)

(3)

where u, v, and w are the velocity in the x-, y-, and z-directions, respectively. The gravity vector is vertical to the z-direction, and its components in the x- and y-directions are gx = −g sin ϕ and gy = g cos ϕ , respectively, where the acceleration of gravity is g = 9.80665 m/s2 , and ϕ is the angle between the gravity vector and the y-direction. Assuming an ideal gas, the gas constant is R = 287 J/(kgK) and the ratio of 85

specific heats is γ = 1.4.

4

2.2. Atmospheric modeling The atmospheric model is based on the standard atmosphere [41]. The atmospheric temperature is set to T∞ = T0 − βh,

(4)

where h is the altitude, and β is the temperature lapse rate. Subscripts ∞ and 0 denote the atmospheric and ground states, respectively. The atmospheric pressure and density are calculated using the equation of state for an ideal gas, p∞ = ρ∞ RT∞ ,

(5)

dp∞ = −ρ∞ g. dh

(6)

and the hydrostatic equation as

Consequently, the atmospheric pressure in the troposphere at β 6= 0 is calculated by p∞ = p0



T∞ T0

g  βR

.

(7)

2.3. Flow solver In a full-field simulation of sonic boom, shock waves propagating over distances of several tens of kilometers must be accurately resolved at low computational cost. Therefore, several numerical schemes for the 90

convective terms are investigated by an exact Riemann solver [42], which is extended to 2nd order accuracy by the MUSCL-Hancock scheme [42, 43] with the minmod and superbee limiters [42], and 3rd and 5th order monotonicity preserving weighted essentially non-oscillatory (MPWENO) schemes [44]. As a result of the axisymmetric analysis described in Sec. 4.1, the MUSCL-Hancock scheme with the superbee limiter is most suitable for analyzing sonic boom propagation in a reasonable computational time; thus, the further analyses

95

are conducted with only the MUSCL-Hancock scheme with the superbee limiter. The numerical schemes assuming a uniform atmosphere used in CFD analysis cannot be directly applied to analysis in a stratified atmosphere because hydrostatic balance between the pressure divergence and gravitational source terms given by Eq.(6) cannot be maintained through discrete grid cells. Therefore, well-balanced finite volume schemes [36-38] in a stratified atmosphere have been developed for avoiding

100

hydrostatic imbalance between discrete grid cells and for performing stable simulation in the fields of meteorology and astrophysics. The novel well-balanced finite volume scheme proposed in this study is described in Sec.2.3.1, and the numerical correction method for preserving a freestream condition in a stratified atmosphere is described in Sec.2.3.2.

5

2.3.1. Well-balanced finite volume scheme in a stratified atmosphere The principle of well-balanced finite volume schemes is to use a local time-dependent hydrostatic reconstruction within each cell [36]. In the case of a standard atmosphere given by Eqs.(4) and (7), the   reconstructed values q˜ = u ˜, v˜, w, ˜ T˜, p˜ in the x-direction are calculated as u ˜ = ui ,

(8)

v˜ = vi ,

(9)

w ˜ = wi ,

(10)

T˜ = Ti − β∆h,

(11)

p˜ = pi 105

T˜ Ti

g ! βR

,

(12)

where ∆h is a difference in altitude between cell i and the height at which the reconstructed values are calculated. The reconstructed values in the y- and z-directions are calculated in the same way. (i) Evaluation of the numerical flux A well-balanced first-order scheme is designed by replacing piecewise constant cell averages with reconstructed values given by Eqs. (8) to (12). The extension of a well-balanced scheme to higher-order is not straightforward, because it must be done by merging a piecewise higher-order reconstruction and a hydrostatic reconstruction for well-balancing in a stratified atmosphere. For this reason, a second order scheme is formulated by decomposing a given primitive variable into an equilibrium and a perturbation. The equilibrium reconstruction is given by Eqs. (8) to (12), and the perturbation is evaluated using piecewise linear reconstruction. Second-order, well-balanced schemes have been used in combination with second-order strong stability preserving (SSP) Runge-Kutta methods [38, 45]. However, the computational cost in the well-balanced schemes, in combination with the second-order Runge-Kutta method, is too high to perform a full-field simulation over distances of tens of kilometers in a reasonable computational time. This is because the computational cost for a stratified atmosphere is much higher than for a uniform atmosphere due to the hydrostatic reconstruction procedure, and because the Riemann problem is solved twice for a secondorder Runge-Kutta method. Therefore, to reduce the computational cost, we propose an alternative new well-balanced MUSCL-Hancock scheme, which is second-order in space and time. The Riemann problem

6

in this MUSCL-Hancock scheme is solved only once in each time step; thus, the computational cost can be reduced. The original MUSCL-Hancock scheme in the x-direction is described as

qiL∗ = qiL +

1 1 qiL = qi − ∆i , qiR = qi + ∆i , 2 2

(13)

∆i = ψ (qi − qi−1 , qi+1 − qi ) ,

(14)





 1 ∆t  1 ∆t  f qiL − f qiR , qiR∗ = qiR + f qiL − f qiR , 2 ∆x 2 ∆x

(15)

where q = (u, v, w, p, T ) is the primitive variable vector. Subscripts L and R denote the left and right L∗ cell interfaces, respectively, and ψ is the limiter function. The Riemann problem is then solved using qi+1

and qiR∗ . The well-balanced MUSCL-Hancock scheme in a stratified atmosphere is designed by modifying the original scheme with the hydrostatic reconstruction q˜ given by Eqs.(8) to (12), and the perturbation ∆q = q − q∞ , giving 1 ˜ R 1 ˜ q˜iL = q˜i−1/2 − ∆ ˜i = q˜i+1/2 + ∆ i, q i, 2 2

(16)

˜i = ψ (∆qi − ∆qi−1 , ∆qi+1 − ∆qi ) , ∆

(17)

     1 ∆t 1 ˜ 1 ˜ ∆ f qi − ∆ − f q + , i i i 2 ∆x 2 2      1 ˜ 1 ˜ 1 ∆t ∆ f qi − ∆ − f q + . = q˜iR + i i i 2 ∆x 2 2

q˜iL∗ = q˜iL + q˜iR∗

(18)

If uniform properties of the atmosphere are assumed in the entire computational domain, i.e., q∞ = const., the well-balanced MUSCL-Hancock scheme is exactly the same as the original scheme. Here, the Riemann 110

L∗ problem is solved using q˜i+1 and q˜iR∗ . The numerical fluxes in the y- and z-directions are calculated in the

same way. (ii) Evaluation of the gravitational source term The gravitational source term is evaluated to maintain hydrostatic balance within each cell through a

7

hydrostatic reconstruction of pressure given by Eq.(12) on the cell interfaces [36-38],   0       p˜i+1/2,j,k − p˜i−1/2,j,k   1   ˜ Sx,i,j,k = ,  0  ∆x      0  
p˜i+1/2,j,k − p˜i−1/2,j,k ui,j,k   0       0   1   S˜y,i,j,k = .  p˜i,j+1/2,k − p˜i,j−1/2,k  ∆y      0  
p˜i,j+1/2,k − p˜i,j−1/2,k vi,j,k

(19)

2.3.2. Numerical correction method To accurately capture shock waves in a reasonable computational time, the computational grid should be aligned with the shock waves. However, the freestream condition in a stratified atmosphere, i.e., u = u∞ , v = v∞ , w = 0, p = p∞ , T = T∞ , cannot be maintained when there is a non-zero, non-horizontal velocity component in the x- or y-direction (e.g. (ρ∞ u∞ )i+1/2,j,k − (ρ∞ u∞ )i−1/2,j,k 6= 0 when u∞ 6= 0 and ϕ 6= 0 deg, even though the hydrostatic balance between the pressure divergence and gravitational source terms is maintained by the well-balanced finite volume scheme. Therefore, to preserve a freestream condition in a stratified atmosphere, a numerical correction method, used in [29, 31], is applied. A slight modification to this method is made; the correction term SC is added to cancel the discretization error as   ∂Q ∂E ∂F ∂G = (Sx + Sy + SC )i,j,k , + + + ∂t ∂x ∂y ∂z i,j,k

˜∞,i+1/2,j,k − E ˜∞,i−1/2,j,k E F˜∞,i,j+1/2,k − F˜∞,i,j−1/2,k + ∆x ∆y  ˜ ∞,i,j,k+1/2 − G ˜ ∞,i,j,k−1/2  G + − S˜x∞ + S˜y∞ , ∆z i,j,k

(20)

SC,i,j,k =

where 115

∂Q ∂t

(21)

= 0. The numerical flux on a cell interface and the gravitational source term are calculated

by the hydrostatic reconstruction given by Eqs. (8) to (12) and (19), under the freestream condition in a stratified atmosphere. As the grid spacing approaches zero, the correction term also approaches zero. The Riemann problem is not solved to calculate the correction term in this study, although it is solved in the L∗ original numerical correction method without hydrostatic reconstruction. This is because qi+1 = qiR∗ when

hydrostatic reconstruction is applied. Finally, the correction term is evaluated as constant within each cell, 120

and the perturbation is computed in the simulation. 8

The well-balanced MUSCL-Hancock scheme given by Eqs.(16)-(18), and the numerical correction method given by Eqs.(20) and (21) are constructed to avoid the influence of discretization error in a stratified atmosphere. As the grid spacing approaches zero, the well-balanced MUSCL-Hancock scheme approaches the original scheme, and the correction term approaches zero. Therefore, if the grid convergence is confirmed as 125

needed in uniform properties of atmosphere, the well-balanced MUSCL-Hancock scheme in combination with the numerical correction method remains conservative, and the Rankine-Hugoniot relation is guaranteed. Computational accuracy is validated to analyze an oblique shock wave over a stratified atmosphere in section. 3. 2.4. Computational grid

130

The computational grid is composed of hierarchical structured AMR and cut cell grids in Cartesian coordinates as follows: 2.4.1. Adaptive mesh refinement The hierarchical, patch-based AMR approach [33] is used to adaptively refine Cartesian grids in regions where all waves generated from a supersonic body might need to be clarified. The flagging criterion to refine

135

the grids is set as dp/dx > 0, where shock or compression waves propagate. When new, finer-grid patches are created, according to the flagging criterion, they need to be populated with data. Therefore, to avoid generation of nonphysical waves due to hydrostatic imbalance on all hierarchical cells, the finer-grid data are set by interpolating a perturbation from the freestream condition on the coarse-grid cell. The ghost cells associated with each grid patch are similarly set by interpolating a perturbation from data on the parent

140

patch. Since the fluxes on a fine grid are calculated separately from those on a coarse grid, the flux into a coarse-grid cell at a fine-coarse boundary is replaced by the sum of the fluxes from adjacent fine-grid cells. However, the atmospheric states in a coarse-grid cell differ from those in fine-grid cells, and as a result, a numerical error occurs. Thus, the threshold value for avoiding amplification of the error and for stable analysis is set as that used in [31]; when a density perturbation is less than 10−7 kg/m3 , all perturbations

145

within the cell are treated as zero. 2.4.2. Cut cell method The localized proportional flux stabilization (LPFS) cut cell method [39] is used for mesh generation around geometries in a Cartesian grid. The maximum AMR level is always used around geometries, i.e., no fine-coarse boundaries on cut cell grids are used in this study. Since atmospheric stratification effects on the

150

flow field around a supersonic body are negligible (see Sec.4.2 for details) [29], the original LPFS cut cell method is directly applied to analysis in a stratified atmosphere, ignoring the gravitational source term on cut cell grids. 9

Figure 1: Overview of the segmentation method for the computational domain inclined at a Mach angle, in which the xy−plane is a symmetry plane of the supersonic body.

2.5. Segmentation of the computational domain The segmentation of the computational domain is useful for improving the efficiency of the computation in which shock waves propagate over very long distances. Therefore, the segmentation method used in [29, 31] is applied as follows. Figure 1 shows an overview of the segmentation method, for the computational domain used in a full-field simulation (see Sec.4.2 for details). In this figure, the xy-plane is a symmetry plane of the supersonic body, and the computational domain is inclined at a Mach angle [46] of   1 θ = arcsin , M∞

(22)

where M∞ is the freestream Mach number. The computational domain in the y-direction is split into 155

multiple subdomains with extended regions for avoiding the influence of the outer boundary at y = ymax from affecting the subdomains, and those in the x- and z-directions are enlarged to enclose all waves emanating from the body and to make their waves pass through the boundary at x = xmax = const. Since there are no waves ahead of the front shock wave, xmin and zmax are determined according to the location of the front shock wave at z = 0 and x = xmax , respectively. In the case of a stratified atmosphere, the

160

freestream Mach number decreases as the atmospheric temperature increases towards the ground. Therefore, the angle of each subdomain is changed according to the freestream Mach number at x = xmin and y = ymin , and the boundary between the subdomain (n) and subdomain (n-1) is set by extrapolating the perturbation obtained in the subdomain (n-1). Note that the difference in grid angle between the subdomain (n) and the subdomain (n-1) is set to a small value (less than 0.016 deg in this study), and thus the influence of

165

the extrapolation is much less than the grid resolution. The computation is performed from the near field around a supersonic body to the far field, extending to the ground, until the time-dependent variation in pressure at the boundary between the subdomain (n) and the subdomain (n-1) converges to less than 0.001 10

∆pmax at the boundary. Finally, the computational results obtained in all subdomains are analyzed to evaluate the sonic boom intensity.

170

3. Numerical simulation of an oblique shock wave The computational accuracy of shock wave propagation through a stratified atmosphere was validated by a comparison with the exact solution of an oblique shock wave. The two-dimensional Euler equations with a gravitational source term were solved to analyze the flow field around a wedge with an angle of θOSW = 5 deg at h = 2 km. The freestream Mach number at h = 2 km was set to M∞ = 1.4. The atmospheric condition was set to the ISO standard atmosphere given by Eqs.(4) to (6), where p0 = 101.325 kPa, T0 = 288.15 K, and β = 6.5 K/km; thus, the magnitude of the atmospheric variables in this simulation increases as altitude decreases. Under the aforementioned condition, the angle of the oblique shock wave at h = 2 km was βOSW = 52.78 deg, given by the oblique shock wave relation [46]: tan θOSW = 2 cot βOSW

2 sin2 βOSW − 1 M∞ . (γ + cos 2βOSW ) + 2

2 M∞

(23)

The non-inclined and inclined computational domains were set to confirm that hydrostatic balance in a stratified atmosphere can be maintained in an inclined grid, aligned with the oblique shock wave, as shown in Fig. 2. Each computational domain was solved using the regular Cartesian and AMR grids. The minimum grid spacing was set to 4 m in all grids, and the maximum AMR level was set to lmax = 2; i.e., the maximum 175

grid spacing in the AMR computation was 16 m. For the non-inclined grid, the computational domain has a dimensional domain size of 1,600 m × 2,240 m, and the regular Cartesian mesh is 400 × 560 cells. The surface of the body was evaluated by a slip wall boundary condition. The bottom and left boundaries were set to the freestream condition in a standard atmosphere, and the right boundary was evaluated by a 0th order extrapolation. For the inclined grid, the computational domain has a dimensional domain size of 512

180

m × 3,040 m, and the regular Cartesian mesh is 128 × 760 cells. The surface of the body was treated as a slip wall boundary condition. The bottom left boundary was set to the freestream condition in a standard atmosphere, and the top-right and bottom-right boundaries were evaluated by a 0th order extrapolation of the perturbation from the freestream condition. All boundaries in the simulation without well-balancing were treated as transmissive, to avoid divergence of the computation. However, in the computation without

185

well-balancing in an inclined grid, an oblique flow was generated in the entire computational domain, due to hydrostatic imbalance, and as a result, the simulation failed to properly predict an oblique shock wave. Shock-wave reflection at the ground was ignored. The CFL number was set to 0.8, and the computation was conducted until t = 30 s. Figure 2 shows the pressure distribution, and the AMR grids. Note that the magnitude of the atmospheric

190

variables in this simulation increases as altitude decreases, and the same flow field was obtained in the regular 11

(a) non-inclined grid with θ = 0 deg

(b) inclined grid with θ = 52.78 deg

Figure 2: Pressure distribution obtained by two-dimensional analysis using the AMR grid for the oblique shock wave generated from the wedge with θosw = 5 deg at h = 2 km in a standard atmosphere.

Cartesian and AMR grids. As shown in Fig. 2, an oblique shock wave is generated from the wedge at h = 2 km and propagates through the stratified atmosphere towards the ground. Since there are no other wave interactions in this configuration, atmospheric stratification is the only effect on the oblique shock wave. The shock-wave angle depends on the atmospheric temperature, which as expressed by Eq. (23), increases 195

towards the ground, and the pressure rise of the shock wave increases with increasing atmospheric pressure. Fine grids are generated only near the body and along the oblique shock wave, and the wave is accurately captured over the entire computational domain in both the non-inclined and inclined grids. Figure 3 shows the pressure waveform at the ground obtained by the non-inclined and inclined grids. Note that the pressure waveform is comparable in the regular Cartesian and AMR grids. The pressure increases due to the oblique

200

shock wave, and the configuration of pressure waveform is in good agreement in all cases. The oblique shock wave in the inclined grid is captured using only a few grid cells, whereas that in the non-inclined grid is captured with a larger number of cells; thus, the inclined grid is more suitable for high-resolution of the shock wave than the non-inclined grid. Tables 1 and 2 show, respectively, the pressure rise of an oblique shock wave and the pressure oscillation

205

scaled by atmospheric pressure at x = xmin at the ground. The pressure oscillation is caused by the 12

(a) non-inclined grid with θ = 0 deg

(b) inclined grid with θ = 52.78 deg

Figure 3: Pressure waveform at the ground obtained by two-dimensional analysis using the regular Cartesian grid for an oblique shock wave generated from the wedge in a standard atmosphere (WB: well-balancing, NCM: numerical correction method). Shock-wave reflection at the ground is ignored. Table 1: Pressure rise of an oblique shock wave at the ground obtained by two-dimensional analysis without considering shock-wave reflection at the ground in a standard atmosphere [kPa].

θ, deg

w/o WB∗

w/ WB

w/ WB and NCM#

w/ WB and NCM

0

25.31

25.32

25.32

25.34

52.78

-

25.33

25.33

25.33

maximum AMR level ∗

WB–well balancing

0 #

2

NCM–numerical correction method

hydrostatic imbalance between the pressure divergence and gravitational source terms, or by the generation of nonphysical waves at the inclined grid in a stratified atmosphere (see Sec.2.3 for details). Since the threshold value for avoiding numerical error at fine-coarse boundaries was set, as described in Sec.2.4.1, the pressure oscillation in the AMR computation was exactly zero. As shown in Table 1, the pressure rise of an 210

oblique shock wave was comparable in all cases. However, the pressure oscillation at x = xmin , unrelated to the pressure perturbation caused by the oblique shock wave, occurred when the well-balancing approach or the numerical correction method was not applied, as shown in Table 2. For the non-inclined grid, the pressure oscillation was of the order of machine precision, i.e. 10−16 , when the well-balanced approach was applied. For the inclined grid, a large pressure oscillation was generated even when the well-balanced

215

approach was applied, because the freestream condition in the stratified atmosphere was not preserved at the inclined grid (see Sec.2.3.2 for details). On the other hand, the pressure oscillation was at the order of machine precision when the well-balanced approach in combination with the numerical correction method 13

Table 2: Pressure oscillation scaled by atmospheric pressure at x = xmin at the ground, obtained by two-dimensional analysis of the oblique shock wave in a standard atmosphere. Shock-wave reflection at the ground is ignored.

θ, deg

w/o WB∗

w/ WB

w/ WB and NCM#

w/ WB and NCM

0

1.4E-05

2.9E-16

1.4E-16

0

52.78

-

3.9E-5

4.0E-16

0

maximum AMR level ∗

WB–well balancing

0 #

2

NCM–numerical correction method

was applied. Therefore, it can be seen that the well-balanced approach with the numerical correction method should be used for performing stable simulation using the inclined grid in the stratified atmosphere. Figure 4 shows the comparison between the computational results and the exact solution of oblique shock wave. Since the pressure rise of oblique shock wave is comparable in all cases, as shown in Table 1, the results obtained using the well-balanced approach with the numerical correction method in the regular Cartesian grid is depicted in Fig. 4. The exact solution of the oblique shock wave [46] can be calculated by
∆pOSW 2γ 2 M∞ sin2 βOSW − 1 , = p∞ γ+1 220

(24)

where the shock wave angle βOSW is obtained directly using the computational result, or by Eq.(23). In this study, θOSW was calculated from the velocity vector just behind the oblique shock wave in the computational results, and βOSW was calculated by Eq.(23). Then, the exact solution was obtained using Eq.(24). As shown in Fig. 4, the pressure rise of the oblique shock wave increases toward the ground. This is because the pressure rise of the oblique shock wave, given by Eq.(24), depends on three parameters: the atmospheric

225

pressure, the freestream Mach number, and the shock wave angle. The change in atmospheric pressure with altitude is more significant than those in the freestream Mach number and the shock wave angle; i.e., the denominator of the LHS of Eq.(24), i.e. the atmospheric pressure, increases towards the ground, although the RHS of Eq.(24) hardly changes. Consequently, the pressure rise of the oblique shock wave strongly depends on the atmospheric pressure and increases with decreasing altitude. The pressure rise in the inclined grid

230

changes smoothly as shown in Fig. 4(b), whereas with the non-inclined grid, it is slightly oscillating. This is because the grid lines of the non-inclined grid are not aligned with the oblique shock wave and thus the resolution of shock wave at the non-inclined grid is lower than that at the inclined grid as shown in Fig. 3. The relation between ∆pOSW and βOSW for both the non-inclined and inclined grids agrees well with the exact solution. The difference in pressure rise between the computational result and the exact solution at

235

the ground was 0.2 % in all cases when the well-balanced finite volume approach in combination with the numerical correction method was applied. These results show that the shock wave propagation through the stratified atmosphere was accurately captured by the well-balanced finite volume method in combination with the numerical correction method. 14

(a) non-inclined grid with θ = 0 deg

(b) inclined grid with θ = 52.78 deg

Figure 4: Comparison between the computational results obtained by two-dimensional analysis using the regular Cartesian grid for an oblique shock wave generated from the wedge in a standard atmosphere and the exact solution.The well-balanced finite volume method in combination with the numerical correction method is applied

4. Numerical simulation of sonic boom 240

The accuracy of numerical simulations for sonic boom was validated by comparison with previous computational results [29], the results of WPM, and flight test data [40]. The computation was used to reproduce D-SEND#1, conducted by JAXA, under the same conditions used in the previous computational study [29] which are summarized here. The flight test model was set to an axisymmetric N-wave model (NWM) with a length of L = 5.6 m. The tail fins attached to the NWM were neglected because they were assumed to have

245

a small effect on the sonic boom waveform in the far field. A steady level flight of M∞ = 1.43 at h = 6.039 km was assumed. The atmospheric conditions based on the meteorological data in the flight test was set by Eqs.(4) to (6), where p0 = 101.25 kPa, T0 = 275.95 K, β = 7.21 K/km. Thus, the atmospheric pressure and temperature at the flight altitude were p∞ = 44.83 kPa and T∞ = 232.41 K, respectively. An axisymmetric simulation in the near field around the NWM was performed to demonstrate that all waves emanating from

250

the NWM can be captured by AMR and cut cell grids, with two numerical schemes compared in Sec.4.1. Thereafter, a full-field simulation in the entire flow field, ranging from the near field, around the NWM, to the far-field, extending to the ground, was performed to demonstrate that it can reproduce real flight test results for sonic boom in Sec.4.2. The grid convergence was extensively studied to assess the sonic boom propagation over the entire flow field. The computation in this section was parallelized using Open MPI

255

and performed on an Intel Xeon E5-2640 v3 @ 2.60 GHz (16 cores).

15

(b) AMR grid (lmax = 6)

(a) regular Cartesian grid (1280×1280 cells)

Figure 5: Distribution of pressure perturbation from the atmospheric pressure in the near field around the NWM obtained by axisymmetric Euler analysis using the MUSCL-Hancock scheme with the superbee limiter.

Figure 6: Pressure waveforms at r/L = 1 obtained by axisymmetric Euler analysis using a MUSCL-Hancock scheme with the superbee limiter in the near field around the NWM. The previous computational result, obtained by the BFC grid [29], is shown (left: overall waveform, right: closeup of front shock wave).

16

(a) 3rd order

(b) 5th order Figure 7: Pressure waveforms at r/L = 1 obtained by axisymmetric Euler analysis using an MPWENO scheme in the near field around the NWM. The previous computational result, obtained by the BFC grid [29], is shown (left: overall waveform, right: closeup of front shock wave).

4.1. Near-field simulation A near-field numerical simulation for investigating computational accuracy of the AMR and cut cell grids, and an appropriate numerical scheme, was conducted by solving the axisymmetric Euler equations. The convective terms were evaluated by a MUSCL-Hancock scheme with the minmod and superbee limiters, 260

and 3rd and 5th order MPWENO schemes. The computational domain in axisymmetric coordinates (x, r) has a dimensional domain size of 20 m × 20 m. The axisymmetric simulation was performed using the regular Cartesian and AMR grids. The regular Cartesian grids are 320 × 320, 640 × 640, and 1,280 × 1,280 cells, and the grid spacing in each grid is ∆x = ∆r = 62.5, 31.25, and 15.625 mm, respectively. In the AMR computation, the maximum refinement level was varied between lmax = 1 and lmax = 6, where the minimum

265

grid spacing is ∆xmin = ∆rmin = 15.625 mm, regardless of the maximum refinement level. The axisymmetric 17

Figure 8: Comparison of pressure waveforms at r/L = 1 obtained by axisymmetric Euler analysis using the regular Cartesian and AMR grids in the near field around the NWM (left: overall waveform, right: closeup of front shock wave).

NWM was located at the top boundary, corresponding to the axisymmetric axis. The uniform properties of atmosphere at the flight altitude, i.e. p∞ = 44.83 kPa and T∞ = 232.41 K, were assumed through the entire computational domain. The surface of the NWM was evaluated using a cut cell grid and was treated as a slip wall boundary condition. The top boundary, except for the NWM, was evaluated with an axisymmetric 270

boundary condition. The left and bottom boundaries were set to a freestream condition in the uniform atmosphere at the flight altitude. The right boundary was evaluated by a 0th order extrapolation. The CFL number was set to 0.8, and the computation was conducted until t = 0.15 s. Figure 5 shows the distribution of the pressure perturbation from atmospheric pressure obtained by a MUSCL-Hancock scheme with the superbee limiter. Note that the configuration of the flow field is the same

275

regardless of the numerical scheme. As shown in Fig. 5, physical phenomena around the NWM consist of six waves as described in the previous computational study [29]: the front and rear shock waves emanating from the leading and trailing edges, respectively, the expansion waves from the front and rear ends of the cylindrical geometry, and the compression waves from the front cone and those for pressure restoration between the expansion waves. The distribution of the pressure perturbation obtained by the AMR grid is

280

similar to that from the regular Cartesian grid. As shown in Fig. 5(b), fine grids in the AMR computation are generated only near the shock waves and compression waves, because the flagging criterion to refine the grids was set as dp/dx > 0 as described in Sec. 2.4.1. Figures 6 and 7 show, respectively, the pressure waveform at r/L = 1 obtained by MUSUCL-Hancock and MPWENO schemes in the regular Cartesian grid. The horizontal axis is the relative time calculated

285

by the streamwise distance from the front shock wave and the freestream velocity. The relative time is set to zero when the pressure perturbation becomes one-tenth of the maximum pressure rise. The previous computational result obtained using the BFC grid [29] is depicted in Figs. 6 and 7 to validate computational 18

(b) computational time

(a) maximum pressure rise

Figure 9: Maximum pressure rise at r/L = 1 and computational time with respect to the minimum grid spacing obtained by axisymmetric Euler analysis in the near field around the NWM.

accuracy in the near field around the NWM. Note that, in the MUSCL-Hancock scheme, the pressure waveform is comparable regardless of the limiter function. As the grid resolution is enhanced, the pressure 290

rise caused by the shock waves and compression waves steepens. The pressure waveform obtained by the finest grid of 1280 × 1280 cells agrees well with that of the BFC grid, regardless of the numerical scheme. The results show that all waves generated from the NWM were accurately captured by the cut cell grid. For an MPWENO scheme, the pressure waveform is comparable at 3rd and 5th order, as shown in Fig. 7, although the pressure perturbation for 5th order is slightly higher than that for 3rd order. Figure 8 shows

295

the comparison of pressure waveforms at r/L = 1 between the regular Cartesian and AMR grids obtained by a MUSCL-Hancock scheme with a superbee limiter. The pressure waveform in the AMR grid is in good agreement with that in the regular Cartesian grid. Figure 9 shows the maximum pressure rise at r/L = 1 and the computational time according to the minimum grid spacing. The maximum pressure rise, caused by the front shock and compression waves, is

300

one of the most important parameters to evaluate sonic boom strength. Thus, an appropriate numerical scheme used in a full-field simulation is selected by investigating the relation between the maximum pressure rise and the computational time. As shown in Fig. 9 (a), the maximum pressure rise increases with decreasing the minimum grid spacing. For the MUSCL-Hancock scheme, the maximum pressure rise with the superbee limiter is higher than that for the minmod limiter. For the MPWENO scheme, the maximum pressure rise

305

in the 5th order version is higher than that at 3rd order. Furthermore, the maximum pressure rise of the MPWENO scheme is also higher than that of the MUSCL-Hancock scheme. As shown in Fig. 9(b), the computational time in the MUSCL-Hancock scheme does not change when using different limiter functions. In the MPWENO scheme, the computational time at 5th order is approximately 8% longer than that at 3rd 19

Figure 10: Computational time and total number of cells obtained by axisymmetric Euler analysis with the AMR grid in the near field around the NWM, where the minimum grid spacing is constant regardless of the maximum AMR level.

order. The computational time for the fine grid with grid spacing reduced by a factor 2 is 6-8 times longer 310

than that for the coarse grid, regardless of the numerical scheme. The computational time in the MPWENO scheme is 6-8 times longer than that of the MUSCL-Hancock scheme. Therefore, if the grid spacing in the MUSCL-Hancock scheme has double the spatial resolution of the MPWENO scheme, the shock wave from the MUSCL-Hancock scheme is a bit more highly resolved than that from the MPWENO scheme, for comparable computational time. In addition, it should be noted that the hydrostatic reconstruction

315

procedure for maintaining hydrostatic balance in high-order numerical schemes takes longer than that in lower-order numerical schemes.

This is because higher-order scheme generally use the larger number

of stencils at which hydrostatic reconstruction values and perturbations need to be calculated. For these reasons, the MUSCL-Hancock scheme with the super-bee limiter is most suitable for performing a full-field simulation of sonic boom in a reasonable computational time. All other simulations, in this paper, were 320

performed using the MUSCL-Hancock scheme with the superbee limiter. Figure 10 shows the computational time and the total number of cells used in AMR computations. The computational time with AMR grids, lmax = 3 was 47.4% of that with a uniform Cartesian grid; thus, the computational cost was significantly reduced by using AMR. Although the total number of cells decreased with increasing AMR level, the computational time increased when the AMR level was too high. Therefore, an appropriate maximum AMR

325

level needs to be chosen for performing efficient simulations. 4.2. Full-field simulation A full-field simulation for sonic boom propagation in the entire flow field, ranging from the near field around the NWM to the far field extending to the ground, was performed to reproduce the D-SEND#1 resutls. Accuracy of the numerical simulation was validated by comparison with previous computational 20

330

results [29], results of the WPM, and flight test data [40]. Comparison of simulation results in uniform and stratified atmospheres is helpful for investigating the effects of atmospheric stratification. Hence, the simulation was performed under the following two atmospheric conditions; uniform atmosphere at flight altitude and horizontally stratified atmosphere based on meteorological data. The governing equations in the uniform and stratified atmospheres were, respectively, the three-dimensional Euler equations and

335

those with the gravitational source term given by Eqs.(1)-(3). The three-dimensional Cartesian coordinates (x, y, z) were set as illustrated in Fig. 1. The computational domain ranged from h/L = 0 to 1078.4. The angle of computational domain was set to the Mach angle given by Eq.(22). When grid lines are aligned with shock waves, the grid resolution in the direction parallel to these waves can be relatively low. Thus, the anisotropic Cartesian AMR grid was set such that the maximum grid spacing in the x-direction was set

340

to ∆x = 1 m, and those in the y- and z-directions were set to ∆y = ∆z = 4 m; i.e., the aspect ratio of cells was 1:4:4. The maximum AMR level around the NWM was lmax = 7, based on the results of axisymmetric analysis in Sec. 4.1, and in the other regions was lmax = 2 or 3. To reduce the computational load as much as possible, the computational domain was split into multiple subdomains using extended regions as described in Sec. 2.5. Table 3 shows the split-domain length in the first and final subdomains with

345

the extended regions. Since the simulation was performed with subdomains that reached the ground, the number of subdomains was nmax = 360 for uniform atmosphere and nmax = 345 for stratified atmosphere. To further reduce the computational cost, the length of the fine-grid region in the y-direction was set to 28 m, which included the entire subdomain and a part of the extended region. The boundary condition at the surface of the body evaluated by the cut cell grid used a slip wall. The boundary at x = xmin was used

350

as a freestream condition, and those at x = xmax , y = ymax , and z = zmax were evaluated by 0th order extrapolation of the perturbation from a freestream condition. At y = ymin , a freestream condition was applied in the subdomain (0), and boundary conditions in the other (n) subdomains were set by 0th order extrapolation of the perturbation from a freestream condition obtained in the previous (n-1) subdomain. At z = zmin , a plane-symmetry boundary condition was applied. Shock-wave reflection at the ground was

355

ignored. The CFL number was set to 0.99, and the computation in subdomain (n) was conducted until the time variation in pressure at the boundary between subdomain (n) and subdomain (n+1) converged to less than 0.001∆pmax at the boundary. Table 4 shows the total number of cells in all subdomains and the computational time of the numerical simulation. Figure 11 shows the computational grid and pressure distribution in the near field around the NWM

360

under consideration of a stratified atmosphere. The computational results obtained in the subdomain (0) to the subdomain (5), excluding the extended regions, are depicted. The pressure distribution in the near field around the NWM is almost the same for both uniform and stratified atmospheres, although the atmospheric pressure in a stratified atmosphere increases as altitude decreases. As shown in the closeup of the near field around the NWM, the physical phenomena around the body are the same as those in the axisymmetric 21

Table 3: Split-domain length in the first (n = 0) and final (n = nmax ) subdomains with the extended regions used for full-field simulation using AMR with a maximum refinement level of lmax = 2 and 3 in the far field [m].

Direction

n

Uniform atmosphere

Stratified atmosphere

lmax = 2

lmax = 2

lmax = 3

0

lmax = 3

20

x nmax

24

26

0 y

24 (subdomain), 40 (extended region) nmax 0

32

z nmax

672

704

Table 4: Total number of cells in all subdomains used for full-field simulation using AMR with a maximum refinement level of lmax = 2 and 3 in the far field and computational time.

Total number of cells

Uniform atmosphere

Stratified atmosphere

lmax = 2

lmax = 3

lmax = 2

lmax = 3

121,275,872

574,690,672

102,104,008

520,857,080

10.7

84.4

22.1

126.1

in all subdomains Computational time, h

365

analysis shown in Fig. 5 (see Sec.4.1 for details). This is because atmospheric stratification has a negligible effect on sonic boom strength in the near field, although it has a significant effect on sonic boom propagation in the far field as described in the previous computational study [29]. Fine grids are generated only near shock waves and compression waves over the entire three-dimensional computational domain. Since the pressure distribution at the boundary between the subdomains is smoothly connected, this shows that

370

segmentation of the computational domain does not affect the computational accuracy. Figure 12 shows the pressure distributions at h = 3 km and the ground. In this figure, the aspect ratio is set to 4:1 to facilitate visualization of footprint shape. Since the shock wave propagates conically, its footprint shape is a hyperbola and the curvature decreases towards the ground. The results in Fig. 12 indicate that full-field simulation is a powerful tool for investigating the characteristics of three-dimensional sonic boom propagation, which

375

has presently been difficult to evaluate using conventional one-dimensional or two-dimensional analysis. Figure 13 shows near-field pressure waveforms at r/L = 1 and 10 on the symmetry plane beneath the body. The horizontal axis is the relative time calculated by the streamwise distance from the front shock wave and the freestream velocity. The relative time is set to zero when the pressure perturbation becomes one-tenth of the maximum pressure rise. The pressure waveform obtained in the axisymmetric simulation, 22

Figure 11: Computational grid (left) and pressure distribution (right) from the subdomain (0) to the subdomain (5) obtained by full-field simulation in the flow field around the NWM with a stratified atmosphere. The maximum refinement levels are 7 in the near field around the body and 3 in the far field.

380

described in section 4.1, is depicted in Fig. 13(a), and the three-dimensional simulation without segmentation of the computational domain is depicted in Fig. 13(b). For the simulation without segmentation, the length of the computational domain in the y-direction was 136 m, corresponding to the equivalent length covered by subdomain (0) to subdomain (3) and the extended region (3). The pressure waveforms at r/L = 1 and 10 in a stratified atmosphere agree well with those in a uniform atmosphere, because atmospheric stratification has

385

a negligible effect on the flow field near the NWM. The rapid pressure rise caused by the front and rear shock waves is well captured by the anisotropic AMR grid with lmax = 3. Consequently, the pressure waveform at r/L = 1 of the three-dimensional simulation is in good agreement with that of the axisymmetric simulation, as shown in Fig. 13(a); i.e., three-dimensional simulation with an anisotropic grid reaches a comparable level of accuracy as an axisymmetric simulation with a square grid. As shown in Fig. 13(b), the pressure

390

waveform at r/L = 10 is in agreement regardless of whether segmentation method of computational domain is applied; thus, an segmentation of computational domain did not affect the computational accuracy. The computational time in an equivalent region, with and without the segmentation methods, was 3.5 h and 12.3 h, respectively. Therefore, by using the segmentation method of the computational domain, the computational time was decreased by 72%, and this method allows for reduction of the computational cost 23

Figure 12: Distribution of pressure perturbations from atmospheric pressure at h = 3 km (upper) and ground level (lower) obtained by a full-field simulation over the entire flow field around the NWM, under consideration of a stratified atmosphere (4:1 aspect ratio). The AMR grid with a maximum refinement level of 3 in the far field is used. Shock-wave reflection at the ground is ignored.

395

of a full-field simulation over long distances. Figure 14 shows the far-field pressure waveform at h = 0.5 km on the symmetry plane beneath the body, obtained by a full-field simulation with AMR using lmax = 3.

In this figure, the previous computational

results obtained by the BFC grid [29], results of the WPM [22], and flight test data [40] are depicted. In the sonic boom prediction based on the WPM, the flight and atmospheric conditions were set to be the 400

same as those used in the simulation, and the near-field waveform at r/L = 10, obtained by AMR grids with lmax = 3, was extrapolated. According to the previous computational study that reproduced the DSEND#1 [29], the location of the extrapolated waveform had a negligible effect on the far-field waveform when r/L > 5. As shown in Fig. 14, the shock and compression waves emanating from the NWM coalesce into an N-shaped waveform (N-wave) in the uniform and stratified atmospheres, because of nonlinear effects

405

over fairly long propagation distances of the waves. The front shock wave is captured using only a few grid cells,

and pressure waveforms of the simulation

agree well with those of the WPM. Table 5 shows the

pressure rise of front shock wave and the Mark VII perceived level (PL) [47], calculated using Shepherd et al.’s method [48], at h = 0.5 km on the symmetry plane beneath the body. Since the PL has often been used for evaluating sonic boom strength [11-13], we adopt it as a measure to compare computational results with 410

the results of the previous computational study, the WPM, and flight test measurements. Note that the PL increases with increasing the pressure rise, whereas it decreases with increasing the rise time. The pressure rise of front shock wave, in the simulation assuming uniform properties of the atmosphere, is lower than that of the WPM, whereas the pressure rise in the simulation considering atmospheric stratification is higher 24

(a) r/L = 1

(b) r/L = 10 Figure 13: Near-field pressure waveform on the symmetry plane beneath the body obtained by full-field simulation for the flow field around the NWM in uniform (UA) and stratified (SA) atmospheres. The result of axisymmetric Euler analysis with a regular Cartesian grid of 1280×1280 cells is depicted in the top figures.

The results of three-dimensional analysis in the

segmented and unsegmented domains are depicted in the bottom figures.

than that of the WPM, as confirmed in the previous computational study [29]. The PL in the simulation 415

is lower than that in the WPM because shock wave thickness in the WPM is assumed to be zero. The computational results in a stratified atmosphere are in good agreement with previous results obtained by the BFC grid. Since the resolution of shock wave in the AMR grid is lower than that of the BFC grid, the PL in the AMR grid is lower than that in the BFC grid. However, the shock waves are well captured using only a few grid cells. Hence, the pressure rise of the front shock wave in the AMR grid was similar to that in

420

the BFC grid, and the shock wave was captured by the coarser grids in the AMR computation, relative to the BFC grid used in the previous simulation. The results above indicate that the full-field simulation using AMR grids reached comparable levels of accuracy as those using the BFC grid and atmospheric stratification 25

(a) uniform atmosphere

(b) stratified atmosphere Figure 14: Far-field pressure waveform at h = 0.5 km (r/L = 989.1) on the symmetry plane beneath the NWM obtained by full-field simulation using the AMR grid with a maximum refinement level of 3 in the far field. The previous computational result using the BFC grid [29], the result of the WPM, and flight test data [40] are depicted.

was achieved using well-balanced finite volume methods with AMR grids. The pressure waveform in the simulation is macroscopically similar to that of the flight test. This demonstrates that the real flight test 425

for sonic boom was successfully reproduced by the novel full-field simulation technique using the Cartesian AMR and cut cell grids in a stratified atmosphere. However, the rise time of sonic boom was observed in the flight test as shown in Fig. 14(b), and the pressure profile behind the front shock wave was fluctuated due to atmospheric turbulence as described in [7]. Consequently, the pressure rise in the simulation is higher than that in the flight test, whereas the PL in the simulation is lower than that in the flight test. Although

430

the rise time prediction has already been conducted using full-field simulation for the BFC grid [30], doing so with AMR, and the consideration of atmospheric turbulence, are left for future work. Figure 15 shows the variation in pressure rise of the front shock wave on the symmetry plane beneath 26

(b) stratified atmosphere

(a) uniform atmosphere

Figure 15: Variation in pressure rise of the front shock wave on the symmetry plane beneath the NWM obtained by full-field simulation using the AMR grid with a maximum refinement level of 3 in the far field. The previous computational results obtained using the BFC grid [29] and the results of the WPM are depicted. Table 5: Pressure rise of the front shock wave and perceived level at h = 0.5 km on the symmetry plane beneath the NWM obtained by full-field simulation using the AMR grid with a maximum refinement level of 3 in the far field.

Atmosphere

AMR grid

BFC grid [29]

WPM

Flight test [40]

Pressure rise [Pa] Uniform

18.0

-

18.7

-

Stratification

28.3

28.3

27.1

26.8

Perceived level [dB] Uniform

92.9

-

98.1

-

Stratification

97.7

100.7

101.6

99.2

the body, obtained by a full-field simulation with AMR using lmax = 3. The pressure rise in the near field decreases rapidly, whereas in the far field it decreases slowly. This is because the effects of geometrical 435

spreading of the shock wave propagation weaken as the propagation distance increases.

The pressure rise

in stratified atmosphere is higher than that in uniform atmosphere, because it depends on the atmospheric pressure which increases towards the ground. The attenuation characteristics of the pressure rise in the simulation are in good agreement with those of the WPM, regardless of the atmospheric condition. The pressure rise obtained using AMR also agrees well with that of the BFC grid as shown in Fig. 15 (b). These 440

results indicate that the full-field simulation using AMR grids reached similar levels of accuracy as the WPM, which is the representative prediction method for sonic boom, and the previous simulation using the BFC grid. Therefore, it can be seen that the sonic boom propagation through the real stratified atmosphere over 27

h/L = 1000 was successfully analyzed by full-field simulation using AMR and cut cell grids in a stratified atmosphere.

445

5. Conclusions We developed a novel full-field simulation method for modeling three-dimensional shock wave propagation through a real stratified atmosphere. This method included the following four numerical approaches: (i) hierarchical structured AMR, (ii) Cartesian cut cells, (iii) well-balanced finite volume methods in combination with a numerical correction method, and (iv) the segmentation method of the computational domain.

450

In particular, the 2nd order, well-balanced, MUSCL-Hancock scheme applied over the cut cell and AMR grids for a stratified atmosphere was newly formulated using a hydrostatic reconstruction within each cell and a perturbation from the freestream condition in a stratified atmosphere. The computational results of an oblique shock wave in a standard atmosphere were in good agreement with the exact solution. Thus, the well-balanced finite volume method used in Cartesian AMR and cut cell grids was effective for performing

455

stable simulations of the shock wave propagating through a stratified atmosphere. Considering the relation between the computational accuracy and the computational cost obtained in the axisymmetric analysis, a MUSCL-Hancock scheme with a superbee limiter was most suitable for predicting pressure signatures. The full-field simulation successfully reproduced the D-SEND #1 results, conducted by JAXA, in a reasonable computational time. The pressure waveform obtained in the full-field simulation was in good agreement

460

with those of the previous computational study, WPM, and flight test. The grid convergence study showed that the mesh size was fine enough to assess pressure signatures over a stratified atmosphere. The results presented in this work demonstrate that the shock wave propagation over h/L = 1000 can be accurately analyzed by a full-field simulation using AMR and cut cell grids in a stratified atmosphere. Full-field simulations using AMR and cut cell grids can be directly applied to analyze the flow field around

465

complex geometries such as three-dimensional aircraft configurations and ground topography. However, a supercomputing system may be required for performing the simulation in a reasonable computational time. Moreover, full-field simulations based on CFD analysis have the potential to consider accelerated flight, maneuvers, and atmospheric variability, such as atmospheric turbulence and wind. For these reasons, we believe that the presented methodology is a promising approach towards numerical flight tests for real-world

470

sonic booms.

Acknowledgments The first author was supported by the Japan Society for the Promotion of Science (JSPS) Overseas Research Fellowships. The authors would like to express their gratitude to the JAXA D-SEND project team, who provided the flight test data. 28

475

References [1] J. J. Alonso, M. R. Colonno, Multidisciplinary Optimization with Applications to Sonic-Boom Minimization, Annual Review of Fluid Mechanics, 44 (2012) 505-526. https://doi.org/10.1146/annurev-fluid-120710-101133 [2] M. A. Park, M. J. Aftosmis, R. L. Campbell, M. B. Carter, S. E. Cliff, L. S. Bangert, Summary of the 2008 NASA Fundamental Aeronautics Program Sonic Boom Prediction Workshop, AIAA J. Aircraft, 51(3) (2014), 987-1001.

480

https://doi.org/10.2514/1.C032589 [3] C. M. Darden, Sonic Boom Minimization with Nose-Bluntness Relaxation, NASA TP-1348, 1979. [4] S. K. Rallabhandi, E. J. Nielsen, B. Diskin, Sonic-Boom Mitigation through Aircraft Design and Adjoint Methodology, AIAA J. Aircraft, 51(2) (2014), 502-510. https://doi.org/10.2514/1.C032189 [5] D. J. Maglieri, P. J. Bobbitt, K. J. Plotkin, K. P. Shepherd, P. G. Coen, D. M. Richwine, Sonic Boom: Six Decades of

485

Research, NASA SP-2014-622 (2014). [6] J. Pawlowski, D. Graham, C. Boccadoro, P. Coen, D. Maglieri, Origins and Overview of the Shaped Sonic Boom Demonstration Program, AIAA Paper 2005-0005 (2005). https://doi.org/10.2514/6.2005-5 [7] M. Kanamori, T. Takahashi, Y. Makino, Y. Naka, H. Ishikawa, Comparison of Simulated Sonic Boom in Stratified Atmosphere with Flight Test Measurements, AIAA J., 56(7) (2018) 2743-2755. http://dx.doi.org/10.2514/1.J056155

490

[8] J. Page, K. Plotkin, C. Hobbs, V. Sparrow, J. Salamone, R. Cowart, K. Elmer, H. R. Welge, J. Ladd, D. Maglieri, A. Piacsek, Superboom Caustic Analysis and Measurement Program (SCAMP) Final Report, NASA CR-2015-218871 (2015). [9] L. J. Cliatt II, M. A. Hill, E. A. Haering Jr., Mach Cutoff Analysis and Results from NASA’s Farfield Investigation of No-boom Thresholds, AIAA Paper 2016-3011 (2016). https://doi.org/10.2514/6.2016-3011 [10] L. J. Cliatt II, E. A. Haering Jr., S. R. Arnac, M. A. Hill, Lateral Cutoff Analysis and Results from NASA’s Farfield

495

Investigation of No-boom Thresholds, NASA TM-2016-218850 (2016). [11] M. A. Park, J. M. Morgenstern, Summary and Statistical Analysis of the First AIAA Sonic Boom Prediction Workshop, AIAA J. Aircraft 53(2) (2016) 578-598, https://doi.org/10.2514/1.C033449 [12] M. A. Park, M. Nemec, Nearfield Summary and Statistical Analysis of the Second AIAA Sonic Boom Prediction Workshop, AIAA J. Aircraft, 56(3) (2019) 851-875. https://doi.org/10.2514/1.C034866

500

[13] S. K. Rallabhandi, A. Loubeau, Summary of Propagation Cases of the Second AIAA Sonic Boom Prediction Workshop, AIAA J. Aircraft, 56(3) (2019) 876-895. https://doi.org/10.2514/1.C034805 [14] J. A. Salamone III, V. W. Sparrow, K. J. Plotkin, Solution of the Lossy Nonlinear Tricomi Equation Applied to Sonic Boom Focusing, AIAA J, 51(7) (2013) 1745-1754. https://doi.org/10.2514/1.J052171 [15] H. Yamashita, S. Obayashi, Sonic Boom Variability Due to Homogeneous Atmospheric Turbulence, AIAA J. Aircraft,

505

46(6) (2009) 1886-1893. https://doi.org/10.2514/1.40215 [16] H. Yamashita, S. Obayashi, Global Sonic Boom Overpressure Variation from Seasonal Temperature, Pressure, and Density Gradients, AIAA J. Aircraft, 50(6) (2013) 1933-1938. https://doi.org/10.2514/1.C031339 [17] M. Kanamori, T. Takahashi, Y. Naka, Y. Makino, H. Takahashi, H. Ishikawa, Numerical Evaluation of Effect of Atmospheric Turbulence on Sonic Boom Observed in D-SEND#2 Flight Test, AIAA Paper 2017-0278 (2017).

510

https://doi.org/10.2514/6.2017-0278 [18] F. Alauzet, A. Loseille, High-order sonic boom modeling based on adaptive methods, J. Comput. Phys. 229(3) (2010) 561-593. https://doi.org/10.1016/j.jcp.2009.09.020 [19] K. J. Plotkin,

State of the art of sonic boom modeling,

J. Acoust. Soc. Am. 111(1) (2002) 530-536.

https://doi.org/10.1121/1.1379075 515

[20] O. A. Morris, D. S. Miller, Sonic-Boom Wind-Tunnel Testing Techniques at High Mach Numbers, AIAA J. Aircraft, 9(9) (1972) 664-667. https://doi.org/10.2514/3.59057

29

[21] S. E. Cliff, S. D. Thomas, Euler/Experiment Correlation of Sonic Boom Pressure Signatures, AIAA J. Aircraft, 30(5) (1993) 669-675. https://doi.org/10.2514/3.46396 [22] C. L. Thomas, Extrapolation of Sonic Boom Pressure Signatures by the Waveform Parameter Method, NASA TN D-6832 520

(1972). [23] A. R. Pilon,

Spectrally Accurate Prediction of Sonic Boom Signals,

AIAA J.,

45(9) (2007) 2149-2156.

https://doi.org/10.2514/1.28159 [24] S. K. Rallabhandi, Advanced Sonic Boom Prediction Using the Augmented Burgers Equation, J. Aircraft, 48(4) (2011) 1245-1253. https://doi.org/10.2514/1.C031248 525

[25] B. E. McDonald, High-angle formulation for the nonlinear progressive-wave equation model, Wave Motion, 31(2) (2000) 165-171. https://doi.org/10.1016/S0165-2125(99)00044-X [26] T. Auger, F. Coulouvrat, Numerical Simulation of Sonic Boom Focusing, AIAA J., 40(9) (2002) 1726-1734. https://doi.org/10.2514/2.1877 [27] A. Loseille,

530

R. L¨ ohner,

Anisotropic Adaptive Simulations in Aerodynamics,

AIAA Paper 2010-169 (2010).

https://doi.org/10.2514/6.2010-169 [28] R. Yamashita, K. Suzuki, Sonic Boom Analysis for Hypersonic Vehicle by Global Direct Simulation, Journal of the Japan Society for Aeronautical and Space Sciences, 62 (2014) 77-84. (in Japanese) https://doi.org/10.2322/jjsass.62.77 [29] R. Yamashita, K. Suzuki, Full-Field Sonic Boom Simulation in Stratified Atmosphere, AIAA J., 54(10) (2016) 3223-3231, https://doi.org/10.2514/1.J054581

535

[30] R. Yamashita, K. Suzuki, Rise Time Prediction of Sonic Boom by Full-Field Simulation, AIAA J. Aircraft, 55(3) (2018) 1305-1310, https://doi.org/10.2514/1.C034629 [31] R. Yamashita, K. Suzuki, Full-Field Simulation for Sonic Boom Cutoff Phenomena, Transactions of the Japan Society for Aeronautical and Space Sciences, 58(6) (2015) 327-336, https://doi.org/10.2322/tjsass.58.327 [32] R. Yamashita, K. Suzuki, Lateral cutoff analysis of sonic boom using full-field simulation, Aerospace Science and Tech-

540

nology, 88 (2019) 316-328. https://doi.org/10.1016/j.ast.2019.03.020 [33] M.J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53(3) (1984) 484-512. https://doi.org/10.1016/0021-9991(84)90073-1 [34] M. Wintzer, M. Nemec, M. Aftosmis, Adjoint-Based Adaptive Mesh Refinement for Sonic Boom Prediction, AIAA Paper 2008-6593 (2008). https://doi.org/10.2514/6.2008-6593

545

[35] G. R. Anderson, M. J. Aftosmis, M. Nemec, Cart3D Simulations for the Second AIAA Sonic Boom Prediction Workshop, AIAA J. Aircraft, 56(3) (2019) 896-911. https://doi.org/10.2514/1.C034842 [36] N. Botta, R. Klein, S. Langenberg, S. L¨ utzenkirchen, Well balanced finite volume methods for nearly hydrostatic flows, J. Comput. Phys. 196(2) (2004) 539-565. https://doi.org/10.1016/j.jcp.2003.11.008 [37] R. Klein, K. R. Bates, N. Nikiforakis, Well-balanced compressible cut-cell simulation of atmospheric flow, Phyl. Trans. R.

550

Soc. A 367(1907) (2009) 4559-4575. https://doi.org/10.1098/rsta.2009.0174 [38] R. K¨ appeli, S. Mishra, Well-balanced schemes for the Euler equations with gravitation, J. Comput. Phys. 259 (2014) 199-219. https://doi.org/10.1016/j.jcp.2013.11.028 [39] N. Gokhale, N. Nikiforakis, R. Klein, A dimensionally split Cartesian cut cell method for hyperbolic conservation laws, J. Comput. Phys. 364(1) (2018) 186-208. https://doi.org/10.1016/j.jcp.2018.03.005

555

[40] D-SEND Database, URL: http://d-send.jaxa.jp/d send e/index.html [cited 30 May 2012]. [41] Standard Atmosphere, International Organization for Standardization, ISO 2533:1975 (1975). [42] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, 2009. [43] B. Van Leer, On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe, SIAM J. Sci. Stat. Comput. 5(1) (1984) 1-20. https://doi.org/10.1137/0905001

30

560

[44] D. S. Balsara, Chi-Wang Shu, Monotonicity Preserving Weighted Essentially Non-oscillatory Schemes with Increasingly High Order of Accuracy, J. Comput. Phys. 160(2) (2000) 405-452. https://doi.org/10.1006/jcph.2000.6443 [45] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Rev. 43(2001) 89-112. https://doi.org/10.1137/S003614450036757X [46] H. M. Liepmann and A. Roshko, Elements of Gas dynamics, John Wiley & Sons, 1957.

565

[47] S. S. Stevens, Perceived Level of Noise by Mark VII and Decibels (E), J. Acoust. Soc. Am. 51(2) (1972) 575-601. https://doi.org/10.1121/1.1912880 [48] K. P. Shepherd and B. M. Sullivan, A Loudness Calculation Procedure Applied to Shaped Sonic Booms, NASA TP 3134 (1991).

31

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Conflict of Interest and Authorship Conformation Form Please check the following as appropriate: 

All authors have participated in (a) conception and design, or analysis and interpretation of the data; (b) drafting the article or revising it critically for important intellectual content; and (c) approval of the final version.



This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue.



The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript

•

The following authors have affiliations with organizations with direct or indirect financial interest in the subject matter discussed in the manuscript:

Author’s name

Affiliation