On Energy Redistribution for the Nonlinear Parabolized Stability Equations Method

Abstract

We identify and quantify a seemingly overlook mechanism for energy transfer between adjacent frequency disturbances in the Nonlinear Parabolized Stability Equations method. Physically, this energy transfer results from the finite-bandwidth nature of actual disturbance spectrums versus the common numerical assumption of a discrete spectrum representation. Both quiet wind tunnel and flight conditions are considered and it is found that, for Mack’s second-mode instability, the mechanism is most significant in the 0.1–1% disturbance amplitude range (based on normalized pressure) and is responsible for a 15–30% increase in predicted disturbance amplitude.

1. Introduction

Sustained and controlled hypersonic flight remains an open and active field of research with important implications to both private corporations and government agencies. To achieve vehicle control at flight conditions [1,2,3], accurate prediction of boundary layer laminar-turbulent transition is imperative; however, the fundamental, physical processes governing the transition are not fully understood [4,5]. A precursor to the highly nonlinear process of transition prediction is the accurate modeling and stability analyses of the mechanisms governing disturbance growth within high-speed boundary layers. The Nonlinear Parabolized Stability Equations (NPSE) technique is an often employed and powerful method for conducting such numerical analysis. Here, we point out a seemingly overlooked aspect of NPSE energetics that appears to be important (and potentially dominant) during the early stages of disturbance growth.

The modal growth scenario for boundary layer transition is as follows [6,7,8,9]. First, small disturbances introduced into a boundary layer will experience linear growth [10,11,12,13,14,15,16,17,18]. Next, the disturbances experience moderate nonlinearity as their amplitudes increase via harmonic generation and energy exchange with other disturbances [19]. Finally, highly nonlinear growth ensues via mean flow distortion and secondary instabilities [20] before transitioning to turbulence. The current goal of the hypersonics community is to compare numerical data with high-quality quiet wind-tunnel and flight experimental measurements, in order to elucidate the physical processes governing linear and nonlinear disturbance growth estimates. The implication of our findings extend to predictions directly applicable to hypersonic boundary layer stability and transition (BLST) flight tests such as BOLT [21,22,23,24,25] and HIFiRE [26,27,28]).

Numerical–experimental comparison usually starts with linear stability theory (LST). The LST remains an important tool to assess the quality of the numerical methodology applied (grid/shock resolutions, wall temperature and nose bluntness [29]). Comparison of LST with experimental data provides an important check on experimental methodologies (model alignment, flow quality, and sensor performance [30]). Only a decade ago such agreement was not to be taken for granted [31,32]. The next level of sophistication includes more advance techniques including linear parabolized stability equations (LPSE) or linearized Navier–Stokes equations solvers. Such techniques capture non-parallel flow effects and serve as a first estimate for amplitude-based comparison. Associated with such linear analysis, comes the concept of an N-factor: defined as $N=ln\left(A/A_0\right)$ where $A_0$ and A are the disturbance amplitudes at some initial and final location, respectively. By definition, the N-factor is a function of initial disturbance amplitude ($A_0$), and its use is usually restricted to linear calculations. Traditionally, the N-factor envelope is often considered to be an upper estimate on the amplitude a disturbance can reach, because as the amplitude of the disturbance grows nonlinear, saturation occurs. That is, it is commonly expected that nonlinear effects prevent disturbances from reaching their linear growth potential. For further details and exceptions, see [33,34,35]. Finally, more sophisticated tools account for nonlinear and multiple-mode interactions such as nonlinear parabolized stability equations (NPSE) or direct numerical simulations (DNS).

The nonlinear techniques described above require an estimate of the initial disturbance amplitude and structure [36]. Such small initial amplitude disturbances are quite challenging to measure experimentally, and may even be close to the noise floor of the some computation. To overcome this challenge, the aerospace community has borrowed from the dynamical systems and other fluid dynamics communities [37,38,39,40,41] to rapidly develop adjoint-based and input–output/resolvent based tools to solve this receptivity problem (i.e., how freestream disturbances are mapped into boundary layer disturbances). Since it is extremely difficult to obtain experimental data for such small amplitude disturbances, the approach is to correlate small disturbance measurements with linear growth calculations under the inherent assumption that the dynamics at small amplitudes are sufficiently linear (i.e., small in amplitude). Based on the second-mode experimental work of Marineau [42,43,44], Chynoweth and Schneider [45,46], and others [28,47,48], a reasonable lower bound on observable disturbance amplitudes is around 0.5–1% based on wall pressure measurements normalized by the freestream pressure. Naturally, details concerning tunnel facilities, measurement techniques, and data processing are important. Time harmonic disturbances in this amplitude range have harmonics on the order of 0.0025–0.01% amplitude, which are less than 1% of the primary disturbance. Hence, in this amplitude range one considers linear analysis to apply at least for the initial growth phase of these observable disturbances. However, Marineau et al. (2017 [44], Figures 7–13 ) showed a case where LPSE analysis underpredicts disturbance growth by approximately 10–20%.

While Marineau’s results are just a single example of this phenomena, the meticulous experimental and computation treatment of their study caught our attention and motivated this study. The discrepancy between computed and observed disturbance growth, particularly in the early stages of growth, being the primary motivation for the current study. Further, this is exacerbated by the traditional view described above, that suggests LPSE analysis should overpredict disturbance growth. To address this discrepancy, we pose the question: How “receptive” are downstream disturbances to upstream disturbances? That is, the traditional interpretation of boundary layer stability experimental and numerical data follows the very linear hypothesis that disturbances are spectrally decomposed and that neighboring disturbances in frequency space do not significantly interact. While in the limit of vanishing amplitude, this would appear to hold, in this manuscript, we show that there is a seemingly overlooked mechanism by which energy from high-frequency upstream disturbances transfer to neighboring frequencies, particularly lower-frequency downstream disturbances. Then, we estimate the potential size of this effect for typical hypersonic wind tunnel experimental and flight conditions.

2. Basic States

Our stability analysis begins with calculations of high-fidelity, laminar basic states. Based on previous work [49,50], we choose to focus our investigation on: (1) A 1 m long, 7° half opening angle, 1mm circular nose bluntness, flared cone (with a flare radius of 4 m) at flow conditions consistent with those expected in AFOSR-Notre Dame Large Mach 6 Quiet Tunnel [51]. (2) A 1 m long, 7° half opening angle, 1mm circular nose bluntness, straight cone at flow conditions consistent with those expected in AFOSR-Notre Dame Large Mach 6 Quiet Tunnel (Figure 1). (3) A 1 m long, 7° half opening angle, 1 mm circular nose bluntness, straight cone at flow conditions consistent with the HIFiRE flight experiments [27]. Flow conditions are provided in

Table 1. Test conditions.

	Mach	Re/m	${\mathit{\rho }}_{\mathit{\infty }}$ [kg/m${}^3$]	T${}_{\mathit{\infty }}$ [K]	u${}_{\mathit{\infty }}$ [m/s]	T${}_{\mathbf{wall}}$ [K]
Tunnel	6	11.0 × 10${}^6$	0.0432	53.0488	875.9795	300
Flight	5.3	13.42 × 10${}^6$	0.1190	201.4	1509.2075	393.4

US3D inflow free-stream conditions for Re/m = 11.0 × 10${}^6$, corresponding to AFOSR-Notre Dame Large Mach 6 Quiet Tunnel conditions and flight conditions from HIFiRE experiment 21.5 s data [27].

and calculated using standard isentropic relationships.

Basic state calculations are generated using US3D, a CFD software package specifically designed for high-speed flows [52,53,54]. US3D uses the finite volume method to calculate fluxes across cells, where each grid has 1000 streamwise points and 690 wall normal points. The grid and numerical convergence of our computational methodology has been established in the literature [30,31], but is not presented here. Wall normal basic state boundary layer profiles are extracted using in-house scripts to prepare the data for stability analysis. The 4 m flared cone is chosen as an idealized example to isolate the effect of energy transfer among disturbances that is anticipated to be experimentally reproducible. Note, most of the investigation presented here is confined to the initial 0.5 m of the cones considered, which is consistent with other quiet tunnel experimental entries.

3. Stability Analysis

3.1. Jokher

The JoKHeR Parabolized Stability Equations (PSE) package by Kuehl et al. [19,33] was developed in collaboration with Dr. Helen Reed at Texas A&M as part of the efforts of the National Center for Hypersonic Laminar-Turbulent Transition Research. The code employs a quasi-3D, compressible, ideal gas, primitive variable formulation, that is, it marches disturbances along a predefined path with the assumption of uniformity in the perpendicular direction. The package consists of Linear Stability Theory (LST), Linear Parabolized Stability Equations (LPSE), and Nonlinear Parabolized Stability Equations (NPSE) codes. These codes have been extensively validated against experimental [19,30,31,55,56] and numerical [29,31,56,57,58,59] datasets. A unique feature of JoKHeR is that it employs a nonlinear wave packet formulation for NPSE implementation which allows for the modeling of finite bandwidth disturbances [19,60]. Thus, it accounts for spectral broadening and low-frequency content generation, which is important for the accurate prediction of nonlinear energy exchanges [33,60,61,62]. Note, in addition to the reference just cited, this manuscript represents another example of such nonlinear energy exchange.

3.1.1. Linear Stability Theory

LST considers a linear, steady, and parallel basic flow state determined from separate CFD simulations, and solves the disturbance equations via substitution of Equation (1) into the Navier–Stokes equations. The disturbance is assumed of the form indicated by Equation (2), and is substituted into the disturbance equations. This leads to the generalized eigenvalue problem with $\alpha$, $\beta$ and $\omega$ being the streamwise wave number, spanwise wave number, and the frequency, respectively. The resulting eigenvalues are used to determine instability, and the corresponding eigenvectors represent the shape of the disturbance in the wall normal direction.

$$\varphi (x,y,z,t)=\underset{basicstate}{\underbrace{\overline{\varphi }\left(y\right)}}+\underset{disturbance}{\underbrace{{\varphi }^{\prime }(x,y,z,t)}}$$

(1)

$${\varphi }^{\prime }=\widehat{\varphi }\left(y\right)e^{i(\alpha x+\beta z-\omega t)}$$

(2)

3.1.2. Parabolized Stability Equations

Originally identified by Herbert [63] and Bertolotti [64], during a critical review of Gaster [65] early nonparallel work, the parabolized stability equations have been developed as an efficient and powerful tool for studying the stability of advection-dominated laminar flows. Excellent introductions to the PSE method and summary of its early development were provided by Herbert [63]. During the early stages of both linear and nonlinear development of this technique, much was established related to basic marching procedures, curvature, normalization conditions and numerical stability of the method itself [64,66,67,68,69]. In a relatively short time, the field rapidly expanded [55,59,70,71,72,73,74,75,76,77,78] to include complex geometries, compressible flow, and finite-rate thermodynamics.

PSE is similar to the Fourier/Laplace transform, where an initial-value problem is considered. However, the slowly varying basic state assumption is made in the streamwise direction where a slow variable $\overline{x}=\frac{x}{Re}$ is introduced. Ultimately, disturbances are assumed of the form

$$F\left[{\varphi }^{\prime }\right]=\underset{\mathrm{shape}}{\underbrace{\tilde{\varphi }(\overline{x},y)}}\underset{\mathrm{wave}}{\underbrace{\mathsf{\Phi }(x,t)}}$$

where the wave part satisfies

$$\begin{array}{ccc}\hfill \frac{\partial \mathsf{\Phi }}{\partial x}& =& i\alpha \left(\overline{x}\right)\mathsf{\Phi }\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \frac{\partial \mathsf{\Phi }}{\partial t}& =& -i\omega \mathsf{\Phi },\hfill \end{array}$$

(4)

and $\mathit{Re}=\frac{U_e{\delta }_r}{{\nu }_e}$ is a Reynolds number based on characteristic values of edge velocity ($U_e$), edge kinematic viscosity (${\nu }_e$), and reference boundary-layer length scale (${\delta }_r$). Thus, PSE considers disturbances of the form

$${\varphi }^{\prime }=\underset{-\infty }{\overset{\infty }{\int }}\underset{\mathrm{shape}}{\underbrace{\tilde{\varphi }(\overline{x},y,\omega )}}\underset{\mathrm{wave}}{\underbrace{A(\overline{x},\omega )e^{-i\omega t}}}d\omega$$

(5)

where $A(\overline{x},\omega )=e^{i\int \alpha (\overline{x},\omega )dx}$. The dependence of the shape function ($\tilde{\varphi }$) and amplitude function (A) on $\omega$ are made explicit. The shape and amplitude functions are essentially the Fourier transform of the disturbance. Upon expansion of the streamwise derivatives

$$\begin{array}{ccc}\hfill \frac{\partial {\varphi }^{\prime }}{\partial x}& =& \underset{-\infty }{\overset{\infty }{\int }}\left(\frac{1}{\mathit{Re}}\frac{\partial \tilde{\varphi }}{\partial \overline{x}}+i\alpha \tilde{\varphi }\right)A{e}^{-i\omega t}d\omega \hfill \\ \hfill \frac{{\partial }^2{\varphi }^{\prime }}{\partial x^2}& =& \underset{-\infty }{\overset{\infty }{\int }}\left(\frac{1}{{\mathit{Re}}^2}\frac{{\partial }^2\tilde{\varphi }}{\partial {\overline{x}}^2}+\frac{2i\alpha }{\mathit{Re}}\frac{\partial \tilde{\varphi }}{\partial \overline{x}}+\frac{i\tilde{\varphi }}{\mathit{Re}}\frac{\partial \alpha }{\partial \overline{x}}-{\alpha }^2\tilde{\varphi }\right)A{e}^{-i\omega t}d\omega ,\hfill \end{array}$$

it is found that the second spatial derivative $\frac{{\partial }^2\tilde{\varphi }}{\partial {\overline{x}}^2}$ is of the highest order, and a perturbation expansion is consistently truncated, thus the term is neglected. This leaves the disturbance equations nearly parabolized, and an efficient marching solution is sought. JoKHeR implements a wave packet formulation that improves representation of energy transfer among modes in a nonlinear calculation, compared to the traditional discrete mode formulation [33]. Ultimately, in the Quasi-3D formulation, the disturbance is discretely represented as ${\varphi }^{\prime }={\sum }_k\tilde{\varphi }{(\overline{x},y)}_{k}A{\left(\overline{x}\right)}_{k}W{\left(\omega \right)}_k{e}^{-i{\omega }_{k}t},$ and a frequency content for each mode is assumed of the form $W_0=\frac{1}{{\sigma }_0\sqrt{2\pi }}{e}^{-\frac{{\left(\omega -{\omega }_0\right)}^2}{2{\sigma }_0^2}}$. Whereas, the standard, discrete mode NPSE formulation uses $W_0=\delta \left(\omega -{\omega }_0\right)$. The bandwidth of harmonics obeys ${\sigma }_i=\sqrt{i+1}{\sigma }_0$, and harmonic balancing is used to calculate nonlinear interactions [33]. The finite bandwidth representation of spectral energy is crucial for modeling certain nonlinear phenomena, such as spectral broadening and low frequency content generation, which is observed in experiments [31,32].

4. Linear Results

Our stability analysis begin with a LST investigation, which is summarized in Figure 2. Much experimental hypersonic second-mode dominated transition work focuses on flared cones (i.e., Purdue 3 m flared cone) as a way to isolate a particular frequency instability via controlling the boundary layer height. This methodology achieves maximum growth in the limited streamwise length available in most quiet hypersonic tunnels. Consistent with the hypothesis of this manuscript, such results have led to the neglect of energy transfer among neighboring frequencies, while still achieving good experimental–numerical agreement. Notice that in both straight cone cases, there is an appreciable span of unstable frequencies, i.e., ≈400–150 kHz at tunnel conditions and ≈1300–500 kHz at flight conditions. For our 4m flared case, we have isolated the shifting frequency region (≈525–325 kHz) to the front 0.35 m of the cone, leaving the latter part of the cone with a nearly constant boundary layer height.

The aforementioned second-mode frequency isolation or “tuning” [60,79] is further clarified from the LPSE investigation shown in Figure 3. The LPSE N-factor plots are usually interpreted via the envelope method, in which one hypothesizes a transitional N-factor, for a given set of conditions, then determines which disturbance first reaches that N-factor. For example, in Figure 3 (lower panel), a transitional N-factor of 10 would correspond to a 750 kHz disturbance transitioning at about 0.5 m. It is also evident from the upper panel, that once the boundary layer height stabilizes, the transition is dominated by a much narrower frequency band. Conversely, the straight cones show a strong dependence of transition on disturbance frequency. Thus, there exists a fundamental difference between the commonly used flared cones in wind-tunnel experiments, and the canonical straight cone geometries typical in flight. If energy transfers among neighboring frequencies, as speculated by Batista [50], and is consistent with the calculations of Khan [49]), one anticipates the effects to be more significant for typical flight geometries.

5. Nonlinear Energy Transfer Mechanism

The finite bandwidth disturbance representation is necessary to properly represent spectral energy redistribution as nonlinearities become important [30,45]. Graphically, this situation is illustrated in Figure 4. The top panel depicts quadratic nonlinearity, i.e., the convolution of the spectral weighting functions of a 300 kHz second-mode (bandwidth of ${\sigma }_0$) with itself. This nonlinear interaction results in a 600 kHz first harmonic (middle top panel) of bandwidth $\sqrt{2}{\sigma }_0$. Depicted in the middle lower panel, this first harmonic feeds back onto the primary second-mode via another quadratic interaction with the complex conjugate of the primary second-mode. Notice that this nonlinear feedback term is forcing a larger bandwidth than that of the primary mode. Not all of the forcing is applied directly to the primary mode, but instead ≈40% of the forcing is projected onto “side lobes”. The side lobe frequencies are represented discretely in the NPSE code, and are depicted as vertical bars at 240 kHz and 360 kHz in the lower panel. Thus, in addition to the linear growth a particular frequency disturbance experiences, there is a secondary forcing applied to disturbances who have neighbors in frequency space of sufficient amplitude. Later we show that a “sufficient” amplitude is around 0.5–1.0% of normalized pressure. In contrast to the traditional "envelope" method, an interesting consequence of such energy redistribution provides a pathway along which energy shifts frequencies as a disturbance propagates downstream.

Re-consider the stability diagram for the 4 m flared cone at tunnel conditions (Figure 2 upper) and the corresponding LPSE results (Figure 3 upper). Eventually, the boundary layer is dominated by a ≈340 kHz frequency disturbance. Now, assume that the transitional N-factor in this case is estimated to the exceptionally large value of 30. Traditionally, one would initialize this 340 kHz disturbance at its N1 neutral point (at ≈0.15 m) with some small initial amplitude and evolve the disturbance downstream. However, it is often neglected that this initial amplitude is itself a function of the upstream instabilities.

Meaning, through the nonlinear spectral broadening mechanism described above, upstream disturbance evolution influences the initial amplitude condition for disturbances farther downstream. As an example, considering Figure 2 (upper), a 460 kHz disturbance begins to grow around 0.075 m from the nose tip. Despite the small amplitude achieved due to the limited extent of its streamwise growth, some energy will be transferred to its side lobe frequencies, particularly the frequencies slightly below 460 kHz. When the 460-$\mathsf{\Delta }\omega$ kHz disturbance begins linear growth, at its neutral point, it already has experienced slight forcing due to the presence of the 460 kHz mode. This energy transfer will repeat itself from the 460-$\mathsf{\Delta }\omega$ kHz disturbance to the 460-2$\mathsf{\Delta }\omega$ kHz disturbance, and the cycle continues to cascade energy from the higher frequency upstream disturbances to the lower frequency downstream disturbances (or vice versa depending on boundary layer height evolution). Each step builds upon the previous, increasing the initial amplitude until transition. Thus, this scenario is valid in the weak to medium nonlinear regime (i.e., prior to the onset of high nonlinearity) and suggests the envelope method should serve as a lower bound to disturbance growth, however the effect of nonlinear saturation may reduce the disturbance amplitude below the envelope prediction level before transition.

6. NPSE Results

To illustrate the integrated effect of this energy transfer mechanism, we first consider the 4 m flared cone case at Notre Dame tunnel conditions (Figure 5). In this case, the most unstable frequency shifts from around 525 kHz down to around 350 kHz within the first 0.3 m. NPSE calculations are presented which span the frequency range from 460–370 kHz at 15 kHz intervals. The traditional discrete frequency methodology (blue), which neglects finite bandwidth effects is compared to the wavepacket methodology (black), which accounts for finite bandwidth effects. The energy shifting mechanism described above is readily observed and it is found that the discrete methodology underpredicts disturbance amplitude by approximately 15% within the first $\approx 0.28$ m of the cone. Note, this effect becomes significant in the disturbance amplitude range from about 0.5–1% amplitude based on normalized pressure. Initial disturbance amplitudes of $6.5\times {10}^{-5}$ are chosen to illustrate the effect.

The 4 m flared cone geometry was chosen to isolate this seemingly overlooked energy transfer mechanism, and may provide an experimentally viable test case. However, actual flight vehicles are more appropriately represented by straight cone geometries. As such, 1m long straight cones at Notre Dame quiet tunnel conditions and typical flight conditions are considered to determine the extent of the “energy frequency shifting” effects on disturbance amplitude prediction. Figure 6 illustrates the comparison between discrete and finite bandwidth methodologies for Notre Dame quite tunnel conditions, while Figure 7 considers flight conditions, specifically those experienced during a HIFiRE flight [27]. In both cases, we again see the energy shifting mechanism is active.

For the Notre Dame tunnel conditions, a frequency range between 200–300 kHz (with a $\mathsf{\Delta }\omega =20$) kHz was considered and the disturbances were initialized with an amplitude of $9\times {10}^{-5}$ at their respective neutral points. As the disturbances grow, the energy transfer becomes noticeable once the disturbances reach an amplitude of roughly 0.5% and becomes significant as the disturbances reach 1% amplitude. Ultimately, the finite bandwidth formulation results in a 35% larger disturbance amplitude than the discrete formulation at 0.65 m from the nose tip. Despite the nonlinear nature of these calculations, and for convenience, the lower panel of Figure 6 shows an N-factor comparison between the discrete and finite bandwidth calculations. The energy frequency shifting has resulted in a N-factor difference of approximately 0.3 at 0.65 m. Or, in other words, the energy frequency shifting mechanism is responsible for a 35% increase in disturbance amplitude prediction compared with the traditional discrete formulation.

For the HIFiRE flight conditions, a frequency range between 725–850 kHz, with a $\mathsf{\Delta }\omega =25$ kHz is considered. The disturbances are initiated with an amplitude of $2.5\times {10}^{-5}$ at their respective neutral points. Again, as the disturbances grow, the energy transfer becomes noticeable once the disturbances reach an amplitude of roughly 0.5%, and becomes significant in this case as the disturbances reach 1–2% amplitude. Ultimately, the finite bandwidth formulation results in a 22% larger disturbance amplitude than the discrete formulation near the middle of the cone (0.49 m from the nose tip). The energy frequency shifting has resulted in a N-factor difference of approximately 0.2 by the mid point of the cone, meaning the energy frequency shifting mechanism is responsible for a 22% increase in disturbance amplitude.

7. Conclusions

In this manuscript, we have begun to address the question: How “receptive” are downstream disturbances to upstream disturbances? This question was motivated by experimental observations which indicated a numerical under-prediction of second-mode wave amplitudes in a hypersonic boundary layer. The wave packet NPSE formulations was compared to the traditional discrete methodology, to verify a seemingly overlooked mechanism by which energy is transferred between neighboring disturbances in frequency space. Such a mechanism appears to have been unaccounted for, due to an emphasis on the linear interpretation of experimental and numerical datasets. This mechanism appears most active (and possibly dominant) in the 0.1–2% disturbance amplitude range (based on normalized pressure measurements of second-mode instability waves), and is important for accurate receptivity calculations. The results suggest that this amplitude range is also near the limit of experimental observations, so knowledge of this mechanism is important when interpreting numerical-experimental data comparisons. It was shown that the “energy-frequency-shifting” mechanism is geometry and flow condition dependent, but in general is expected to account for a 20–30% underprediction in disturbance amplitude, over a 0.5 m streamwise distance, when using the traditional (discrete) PSE formulation.