Numerical Simulation of Disturbance Evolution in the Supersonic Boundary Layer over an Expansion Corner

The linear and nonlinear stages of disturbance development in the supersonic boundary layer over a 10° expansion corner is investigated numerically within the framework of Navier—Stokes equations for Mach number 3. The effect of sudden flow expansion on the disturbance evolution is analyzed. The flow stabilization effect observable in the aerodynamic experiment is also discussed.


NOMENCLATURE
x, y, z = longitudinal, vertical, and lateral Cartesian coordinates composing the right-hand triple u, , w = Cartesian coordinates of the velocity vector p = pressure l = coordinate along a surface: l = x at x < 0 and l = x/cos(ε) at x ≥ 0 d w = distance to the wall = pressure disturbance on the wall M = Mach number Re ∞L = Reynolds number based on the scale length L β = lateral wavenumber δ = boundary layer thickness according to the u(δ) = 0.99u e criterion ω = cyclic frequency ε = expansion corner in degrees Subscripts e = at the outer edge of the boundary layer ∞ = in the freestream Abbreviations FP = flat plate, ε = 0°E C = expansion corner, ε = 10°W P = wave packet (linear regime) TS = turbulent spot (essentially nonlinear regime) Flow past the elements of supersonic flight vehicles is associated with the formation of flow acceleration and deceleration zones. Within these zones the boundary layer can interact with shocks, be separated, and form zones of enhanced heat transfer on reattachment to the surface. The boundary layer turbulization can considerably increase this effect. In practice, flow acceleration zones with a favorable (negative) pressure gradient are encountered not less than deceleration zones. Despite this fact, the fundamental problem of the influence of a favorable pressure gradient accelerating the flow on the laminar boundary v ' w p layer stability and the return of turbulent flow to the laminar state (relaminarization, see, for example [1]) has been studied considerably less.
The problem of relaminarization has been studied from the last mid-century. An investigation of subsonic turbulent flows in the presence of a large negative pressure gradient showed the possibility of complete boundary layer relaminarization [2][3][4][5], which is connected with the curved nature of streamlines and favorable longitudinal and normal pressure gradients which lead to a rapid diminishing of the turbulent fluctuation scale in the flow acceleration region [6]. The longitudinal gradient of the static pressure and the freestream Reynolds number were noted as the basic flow parameters influencing the generation and development of the relaminarization process.
In the supersonic flow regime relaminarization leads to a considerable decrease in viscous friction forces and heat transfer from the warm gas to the surface in a flow. It is known that in the outer region of the boundary layer the compressibility effects are predominant over the other effects (see, for example, [7,8]). The compressibility effects include a weak decay and enlargement of large-scale eddy structures in an expansion flow at Mach number 3 [7]. As noted in [7], the Reynolds shear stresses are considerably suppressed with the result that large eddies grow weak downstream. It was also noted that small-scale structures are considerably suppressed immediately behind a rarefaction fan. These inferences were confirmed in experiments at Mach number 4.9 [9,10], where it was found that the intermittence region of the boundary layer shortens in an acceleration flow being displaced toward the boundary layer edge. This fact makes difficult the mixing of the gas from the external flow in the boundary layer.
In [11] the possibility of partial boundary layer relaminarization was demonstrated for Mach number 4 (in its wall region consisting up to 40% of its total thickness). It was noted that an increase in the Reynolds number leads mainly to the stretching of the relaminarized flow region, while the greater flow acceleration (greater negative pressure gradient) in the interaction region results in a progressive decrease of turbulent fluctuations. It was noted that the relaminarization criteria obtained at subsonic velocities can also be applied at large supersonic velocities.
The experiments [12] performed in a shock tunnel at Mach numbers from 5 to 8 confirm indirectly the stabilizing effect of the favorable pressure gradient on the "ogive-cone-cone-cylinder" body of revolution and also indicate the relaminarization of a turbulent wedge behind an isolated surface roughness.
The experimental studies, such as [10], provide a good foundation for developing and testing the computational models of flows of this type, which can be applied together with the Reynolds equations and large eddy simulation. However, owing to computational complexity, there are only few publications concerning direct numerical simulation of laminar flow stability and turbulent flow relaminarization at supersonic velocities. Of these few studies we can note [13] performed at the Mach number M = 2.9 and [14] at M = 2.7. In these two studies turbulent flow around an expansion corner was considered. The two-layer structure of the accelerated flow near the corner was found out; it confirms experimental observations. The flow in the upper layer is characterized by strong suppression of turbulent fluctuations, which are slowly recovered further downstream. In the lower layer the fluctuations are suppressed only in a small vicinity of the flow turn region and are rapidly recovered further downstream.
The mechanisms of disturbance growth and flow relaminarization control the laminar-turbulent transition process. It is known that in the case of a low level of external disturbances, which is inherent in flight conditions, the transition in a supersonic boundary layer over a smooth surface includes successive formation, growth, and coalescence of individual turbulent spots. A spot is formed at the nonlinear stage of wave packet development, whose evolution is determined by the undisturbed flow parameters, the disturbance background, and the temperature of the surface in a flow. In accordance with the linear stability theory, it is the oblique first-mode waves that are predominant over a hot (thermally insulated) surface, while the plane second-mode waves predominate over a cooled surface at large Mach numbers. Although a turbulent spot is an essentially nonlinear object, its development can be influenced by the linear stability mechanism. The calculations [15] establish the relation between the characteristics of turbulent spots and linear wave packets. There are certain experimental and numerical investigations that confirm this assertion (their detailed review is given in [15]).
It should be noted that there are only few calculated data on the development of wave packets and turbulent spots in supersonic boundary layers over expansion corners. In this study, the evolution of a wave packet (linear regime) and a turbulent spot (nonlinear regime) in the supersonic boundary layer over an expansion corner is investigated at Mach number 3. The numerical approaches used are described in Section 1 and an analysis of the numerical results is made in Section 2, where the mean flowfield, the disturbance development in the physical plane, and the spectral characteristics of the disturbances are presented. The conclusions are given in the Summary.

METHODOLOGY OF THE CALCULATIONS
The numerical simulation of the flowfields was performed using the in-house software package HSFlow [16], in which the Navier-Stokes equations are discretized by means of the second-order finitevolume TVD-scheme. The reconstruction of the convective flux terms on the cell sides is performed using the WENO-3 scheme. The mean stationary flowfield is obtained using a second-order method of steadystate attainment with application of structured multiblock computation grids. The calculations were carried out for a perfect gas (air) with a constant adiabatic exponent γ = 1.4 and the Prandtl number Pr = 0.72. The viscosity coefficient is calculated from the Sutherland formula with the characteristic temperature T = 110.4 K: μ = (1 + T μ )/(T + T μ )T 1.5 . The flow velocity, temperature, density, and pressure are nondimensionalized as follows: (u, , w) = (u*, *, w*)/u , T = T*/T , ρ = ρ*/ρ , and p = p*/(ρ u ); all the dimensional parameters are marked by the asterisk as a superscript. The no-slip condition is imposed on the wall: (u, *, w) = (0, 0, 0); the wall temperature is constant and is equal to the flow recovery temperature calculated basing on the given stagnation temperature T = 290 K for the freestream Mach number M ∞ = 3: T = T = T (1 + Pr 0.5 0.5(γ -1)M ) ≈ 261.8 K; T = T (1 + 0.5(γ -1)M ) -1 ≈ 103.57 K; T w /T e ≈ 2.53. The simplified boundary condition of extrapolation is preassigned for the pressure: (∂p/∂n) w = 0. The dimensional quantities can be assessed using the scale length L* = 0.1 m and the time scale τ* = L*/u ≈ 1.64 × 10 -4 . The dimensionless frequency and wave vector components are ω = ω*L*/u and (α, β) = (α, β)L*, respectively. We will consider flow over a flat plate (ε = 0°, abbreviated FP) and over an expansion corner (ε = 10°, abbreviated EC). We consider the development of a wave packet (linear regime, abbreviated WP) and a turbulent spot (nonlinear regime, abbreviated TS) over both geometries, four cases, altogether. The expansion corner is located at x = l = 0, where l is the coordinate along the surface (x = l at x < 0 and x = lcos(ε) at x ≥ 0).
The calculation procedure consists of five steps and ensures the same initial disturbance fields within the WP and TS groups. First, the mean flow at the plate at -7.5 ≤ x ≤ 0.2 is calculated using the method of steady-state attainment. The freestream conditions are imposed on the entry boundary (left and upper sides), while on the exit (backward) boundary all the variables of the problem are linearly extrapolated from within the computation domain.
Secondly, a subdomain is separated out under the bow shock wave generated by viscous-inviscid interaction; its left entry boundary is located at l = -7.2, 0 ≤ y ≤ ∼0.12. The flow is fixed on the new entry boundaries basing on the data obtained at step 1. Then the mean flow is refined again using the steadystate attainment method, up to the moment, when the calculated flow parameters change within 10 -11 for the unit of the computational time.
Thirdly, the mean flowfield is doubled in the third direction for 0 ≤ z ≤ 0.7; in this direction, the boundary conditions of symmetry are applied. Disturbances are introduced in the boundary layer using a generator of the "injection-suction" type through a square orifice on the wall (-6.9923 = x s ≤ x ≤ x e = -6.6623, -0.036 = z s ≤ z ≤ z e = 0.036), which is modeled by means of a time-dependent boundary condition for the gas flow rate in the direction normal to the surface: Here, ω 0 = 10.021 is the baseline generator frequency and A is its amplitude. To form the linear wave packet we take A = 5 × 10 -7 and to generate the turbulent spot we take A = 5 × 10 -4 . The generator excites disturbances near the lower branch of the neutral curve for the fundamental harmonic ω 0 . The generator works during the time period 0 < t < π/ω 0 , which provides a near-uniform initial disturbance spectrum for ω < ω 0 , perceived by the boundary layer. The computation grid in the subdomain has an excessive resolution N x × N y × N z = 2023 × 262 × 330 corresponding to approximately 90 points on the longitudinal wavelength of the fundamental harmonic ω 0 . The temporal resolution amounts to about 125 points on the fundamental harmonic period 2π/ω 0 . On this grid the disturbances are calculated up to the moment t = 7.5, when their leading front approaches the section l = 0.
Fourthly, in accordance with the procedure described above, the mean flow is calculated in a domain extended in both x and z: 0 ≤ z ≤ 1.5; the length along the wall l ≤ 6; the left boundary of the subdomain begins at l = -3.5, 0 ≤ y ≤ 0.55; the grid dimensions are N x × N y × N z = 1317 × 262 × 470. The extended grid resolves the baseline disturbance wavelength on 45 points. The transverse boundary layer resolution, normal to the surface, does not change, as compared with the original grid, and amounts to about 100- cos π cos π (sin(ω ) sin (2ω )).
120 grid lines. At l > 6 the extended computation domain is closed by a buffer zone with the cells strongly enlarged in both x and y intended for suppressing disturbances arriving to the exit boundary and capable to cause the numerical method instability in the case of an intense turbulent spot.
Fifthly, the disturbance field is determined by means of subtracting the mean field from the disturbed field obtained at steps 1 to 3. It is added to the mean flow over the flat plate or the expansion corner in the extended domain obtained at step 4 with application of first-order interpolation to transfer the disturbances between noncoinciding computation grids. Thus, the initial disturbance fields in the extended computation domain are the same for all the cases considered in this study. It should be noted that the procedure of the disturbance transfer from one grid to another does not introduce any spurious disturbances into the solution. The calculations are continued until the disturbances leave the computation domain, that is, until the greatest correction to the dependent variables of the problem on a time step amounts the values of 10 -7 (t max = 20) for the wave packets and 10 -4 (t max = 26 for the expansion corner and t max = 30.5 for the flat plate) for turbulent spots throughout the entire flowfield. Thereupon, the pressure disturbance field p (t, l, z) on the surface is analyzed.
The extended computation domain and the steady flowfield in it are presented in Fig. 1 for the case of the expansion corner in the full computation domain (the detailed discussion is given in Section 2).
The experience of the previous calculations of the authors has shown that the grid resolution of 90 (correspondingly, 47) points over a wavelength leads to a smaller than 0.01% (correspondingly, 0.27%) reduction in the disturbance amplitude on one length of the monochromatic acoustic wave with respect to the natural level of its viscous decay. A preliminary analysis within the framework of the linear stability theory has shown that this numerical decay is small compared with the physical disturbance growth or decay in boundary layers. The total error of the calculations of the disturbance field for the case of the wave packet can be assessed as less than 10% at the end of the computation domain. In the case of the turbulent spot the generated disturbances have even smaller scales. They are subjected to a greater numerical decay; it is supposed that this fact cannot influence the conclusions made in this paper.

Mean Flow
The calculated mean flowfield over the expansion corner is presented in Fig. 1. The inviscid approximation of this class of flows is known as the Prandtl-Meyer flow [17]. The flow parameters change rapidly in a small vicinity of the corner l = 0; viscosity smoothens the distributions of the gasdynamic parameters in this region. The pressure on the surface varies on the scale of the boundary layer thickness calculated directly ahead of the bend. The temperature and velocity profiles are recovered on a greater distance. The boundary layer thickness increases rapidly on several scales of δ ahead of the corner.
In Fig. 2 the calculated profiles of the absolute value of the velocity, the Mach number M, and the static temperature T are plotted ahead of and behind the corner; here, d w is the distance to the wall. Ahead of the corner the profiles are in good agreement with those for the case of the flat plate and begin experience distortions only in the immediate proximity of the corner, l > -0.1, that is, the upstream effect of the ' w  With the passage through the corner the flow velocity increases only slightly; the new absolute value of the velocity at the boundary layer edge V e ≈ 1.056 is about 6% greater than the corresponding value upstream of the corner. The flow acceleration, regarded as an increase in the Mach number, is realized at the expense of flow cooling (Fig. 2b). As can be seen in Fig. 2, the flow parameters at the outer edge of the boundary layer become steady even at l < 0.3, which is of the order of 10δ. In this case, the expansion fan is located fairly high above the surface and the boundary layer is formed beneath it. With increase in the Mach number the fan becomes more gently sloping, with the result that the interaction region behind the corner lengthens.
It should be noted that the flow parameters at the boundary layer edge behind the corner are in good agreement with the values predicted by Prandtl-Meyer theory: M e ≈ 3.58 and T e ≈ 0.79. For this reason, this theory can be used for rapidly estimating the parameters of the flow over the expansion corner. For example, we will assess the ratio of the boundary layer thicknesses over the corner and the flat plate under the assumption of a constant gas flow rate through the boundary layer (2.1) Substituting M 1 = M ∞ = 3 and M 2 = M e2 ≈ 3.58 we obtain δ 2 /δ 1 ≈ 1.7. In the direct calculations of Navier-Stokes equations the corresponding value is 1.5.
It should be noted that the rapid boundary layer broadening must be accompanied by a rapid variation in its stability characteristics. In particular, assuming that the phase velocities of the disturbances for the FP and EC cases are similar in value, we obtain that the unstable frequency range must be scaled together with the boundary layer thickness: f 2 /f 1 ≈ λ 1 /λ 2 ≈ δ 1 /δ 2 , where f is the frequency of the predominant disturbance harmonic. Because of this, it might be expected that the disturbances growing ahead of the corner will transform into disturbances decaying downstream of the corner. This supposition is discussed below.

Disturbance Evolution in the Physical Plane
We will consider the disturbance evolution with reference to the example of a pressure disturbance p on the wall for several successive moments of time (Fig. 3).
In accordance with the linear stability theory for the flat-plate boundary layer, for the flow parameters under consideration there exists only one unstable mode, namely, the first mode according to Mack's terminology [18,19]. The first-mode disturbances are oblique waves propagating with the front inclination  As can be seen in Figs. 3a-3d, the wave packet FP enhances downstream monotonically, increasing in proportion in the longitudinal and lateral directions. The wave packet EC develops in the same fashion up to the bend of the surface at l = 0 and decays monotonically and rapidly further downstream. It should be noted that, as the EC packet passes over the corner, any new disturbances do not appear in the p (x, z) field in the boundary layer. The quantitative comparison of the amplitudes of the wave packets FP and EC is given in Fig. 4a, where the distributions of the greatest-in-z quantity p (x, z) are presented in each grid section x = const. The packet FP amplitude increases downstream almost exponentially, whereas the EC packet decays exponentially at l > 0.
The turbulent spot on the flat plate evolves in a different way (Figs. 3e-3h): it enlarges monotonically downstream but its amplitude remains on an approximately same level (Fig. 4b). On the spot periphery, particularly, in its forward region, there is a region of lowered pressure, while the spot core is the region of elevated pressure. In the vicinity of the spot, where the nonlinear interaction is almost absent, oblique wave packets are formed; their geometric parameters are characteristic of the first-mode wave packet (Figs. 3a-3d). As distinct from the case of the wave packet, the turbulent spot EC is attenuated when moving over the expansion corner but this attenuation is only local in nature and the spot continues to grоw, similarly to the case of FP at l > 1. The presence of the expansion corner delays the spot development. This observation is confirmed in Fig. 4b: the level of the maximum disturbance amplitude in the spot decreases sharply behind the corner and is slowly recovered further downstream approaching the corresponding disturbance level over the flat plate.

Frequency-Amplitude Analysis of Disturbances
We will consider the disturbance evolution in the spectral plane. For this purpose, we will take the field of the two-dimensional fast Fourier transform p (ω, β) for the quantity p (t, x, z) with respect to the variables t and z; here, x is a parameter. In the case of the wave packet the spectra contain two symmetric maxima determining the angle of inclination of the wave front (Fig. 5a). Further downstream the frequencies and the wave numbers of these maxima decrease slowly, which is in agreement with the results of the linear stability theory and will not be discussed here in detail; as for the maximum values, they increase downstream. Ahead of the corner (l < 0) the wave packet spectra are identical for the FP and EC cases. Behind the corner the disturbance spectrum in the EC packet decays monotonically and uniformly throughout the entire spectral range, except from the near vicinity of zero: β ∼ 2, ω ∼ 10. Here, a weak growth on the level of the background noise is observable. Supposedly, this is a new wave packet which is generated from background disturbances in the restructured boundary layer. Owing to restricted dimensions of the computation domain, a more detailed analysis on the basis of the calculations performed seems purposeless.
The corresponding spectra for the turbulent spot are qualitatively different. Due to a large amplitude of the disturbance generator, the initial spectrum of disturbances generated in the boundary layer is not simple, although two maxima of the first-mode disturbances can again be separated out; they correspond to the maxima in Fig. 6. The nonlinear interaction leads to the appearance of harmonics with multiple frequencies and wave numbers; the spectrum is rapidly filled and breaks into small parts; the (ω, β) = (0, 0) harmonic appears and enhances, which indicates an increasing variation of the mean flow within the spot. Again, the spectra turn out to be identical for l < 0, while differences appear directly at l ≥ 0. The main difference is that the EC spectrum amplitude is rapidly reduced throughout the entire frequency-wave range and the spectrum becomes less filled but further downstream it is slowly filled and recovered with respect to the amplitude. Nevertheless, the saturation region is not recovered and remains smaller than in the FP case (Figs. 6e and 6f).
The spectral behavior of the wave packet and the turbulent spot described above is clearly illustrated by the distributions of the maximum Fourier amplitude on the surface presented in Fig. 7. Clearly that the growth of the strongest harmonic of the packet is near-exponential and is the same for the FP and EC cases at l < 0; the decay of this harmonic behind the expansion corner is also near-exponential. The growth of the strongest harmonic in the turbulent spot, (ω, β) ≈ (0, 0), proceeds most actively at the spot formation stage, l < -2, and then decelerates considerably. Directly behind the expansion corner this harmonic attenuates jumpwise at 0 ≤ l ≤ 0.5 but its growth is renewed, the growth rate (slope of the curve) being similar in value with that for the flat plate.

Delay of the Turbulent Spot
As noted above, the presence of the expansion corner delays the development of the turbulent spot. We will illustrate this thesis. In Fig. 8a the turbulent spot is visualized on the plate, when it is completely  located at l < 0 (that is, coincides with the EC case) and after the passage through the line l < 0 for both the plate and the expansion corner. In both cases (Figs. 8b and 8c) the spot shape is near-triangular. However, in the FP case (Fig. 8b) the spot is greater than in the EC case (Fig. 8c), although the spot structures seen in profile differ only slightly. It might be expected that due to viscous friction the contribution of the greater spot will be greater.
The contribution of the spot can be calculated, as follows: where S is the area of the surface in a flow and Δc f, x is the excessive friction coefficient, compared to the case of undisturbed laminar flow. The friction coefficient is determined in the standard way where ∂V/∂n is the derivative of the absolute value of the velocity vector with respect to the normal to the surface. The location of the application of the excessive force ΔF v,x can be determined from simple geometric considerations, similarly to the location of the body center-of-mass (2.4) In Fig. 9 the evolution of the turbulent spot contribution to the viscous drag force is presented. On passage through the corner (EC case) the quantity l c < 0. This is due to the fact that a part of the turbulent spot that has penetrated into the l > 0 region loses rapidly its intensity, compared with the other part of the spot. When the EC spot has penetrated completely into the l > 0 region, its contribution into the viscous friction force begins to increase again. In this case, the growth rates are similar in value for the FP and EC cases, which confirms the supposition about the delaying effect of the expansion corner on the turbulent     Fig. 9; it is approximately unity.

SUMMARY
Within the framework of Navier-Stokes equations the development of wave packets and turbulent spots in a supersonic boundary layer (Mach number 3) over a 10° expansion corner is investigated numerically. The occurrence of the expansion corner leads to flow stabilization. It is shown that the wave packets of the first unstable mode decay monotonically behind the expansion corner. This is due to the sudden flow restructuring, the boundary layer broadening, and, as a consequence, a variation in the boundary layer stability characteristics. The unstable region is scaled with respect to the frequency together with the boundary layer thickness, whose variation may be estimated on the basis of Prandtl-Meyer theory.
It is also shown that the localized regions of turbulent flow, or turbulent spots, are suppressed on passage through the expansion corner only locally. The occurrence of the expansion corner delays the spot development on a scale of the order of 40 thicknesses of the undisturbed boundary layer ahead of the corner; the disturbance amplitude decreases and the frequency-wave spectrum becomes less filled. Downstream of this region the turbulent spot grows similarly to the spot the gradient-free flow over the flat plate. For this reason the turbulent flow relaminarization effect observable in the experimental heat-transfer patterns can actually indicate only local attenuation of turbulence rather than its complete suppression.
An investigation of the second unstable mode of the boundary layer, which appears at large Mach numbers, and the behavior of viscous friction and heat transfer to the expansion corner surface is the theme of further studies of the authors.

DECLARATION OF CONFLICTING INTERESTS
The Authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.