Boundary layer stability and transition analysis using the nonlinear One-Way Navier-Stokes approach

Abstract

We extend the One-Way Navier Stokes (OWNS) approach to support nonlinear interactions between waves of different frequencies, which will enable nonlinear analysis of instability and transition. In OWNS, the linearized Navier-Stokes equations are parabolized and solved in the frequency domain as a spatial initial-value (marching) problem. OWNS yields a reduced computational cost compared to direct numerical simulation (DNS), while also conferring numerous advantages over the parabolized stability equations (PSE), despite its higher computational cost relative to PSE, that we seek to extend to nonlinear analysis. We validate the nonlinear OWNS (NOWNS) method by examining nonlinear evolution of two- and three-dimensional disturbances in a low-speed Blasius boundary layer compared to nonlinear PSE (NPSE) and DNS results from the literature. We demonstrate that NOWNS supports non-modal instabilities, is more robust to numerical noise, and converges for stronger nonlinearities, as compared to NPSE.

Nomenclature

$(x,y,z)$  =	streamwise, wall-normal, and spanwise coordinates
$\nu$  =	specific volume
$(u,v,w)$  =	velocity (streamwise, wall-normal, and spanwise)
$p$  =	pressure
$\mu$  =	dynamic viscosity
$a$  =	sound speed
$\mathrm{Ma}$  =	Mach number
$c_{p}$  =	specific heat capacity at constant pressure
$\gamma$  =	heat capacity ratio
$k$  =	thermal conductivity
$\delta_{0}^{*}$  =	inlet Blasius length scale, defined as $\delta_{0}^{*}=\sqrt{\frac{x_{0}^{*}\mu_{\infty}^{*}\nu_{\infty}^{*}}{U_{\infty}^{*}}}$
$Re$  =	Reynolds number in terms of $\delta_{0}^{*}$, defined as $Re=\frac{a_{\infty}^{*}\delta_{0}^{*}}{\nu_{\infty}^{*}\mu_{\infty}^{*}}$
$Pr$  =	Prandtl number, defined as $Pr=\frac{\delta_{0}^{*}\mu_{\infty}^{*}}{k^{*}}$
$\omega$  =	temporal frequency
$F$  =	dimensionless temporal frequency, defined as $F=\frac{\omega^{*}\nu_{\infty}^{*}\mu_{\infty}^{*}}{U_{\infty}^{*2}}$
$\beta$  =	spanwise wavenumber
$b$  =	dimensionless spanwise wavenumber, defined as $b=\frac{\beta^{*}\nu_{\infty}^{*}\mu_{\infty}^{*}}{U_{\infty}^{*}}$
$\alpha$  =	streamwise wave number (complex-valued, where $\alpha_{r}$ is the phase and $\alpha_{i}$ is the growth rate)
$\mathrm{Re}_{x}$  =	Reynolds number in terms of $x$, defined as $\mathrm{Re}_{x}=\frac{U_{\infty}^{*}x^{*}}{\nu_{\infty}^{*}\mu_{\infty}^{*}}$
$M$  =	denotes the number of temporal Fourier modes
$N$  =	denotes the number of spanwise Fourier modes
$(m,n)$  =	denotes Fourier mode with temporal frequency $m\omega$ and spanwise wavenumber $n\beta$
Subscripts
$\infty$  =	free-stream value
$mn$  =	a quantity associated with temporal frequency $m\omega$ and spanwise wavenumber $n\beta$
$+-0$  =	plus, minus, and zero characteristics
Superscripts
$*$  =	dimensional quantity
$\prime$  =	a disturbance to the equilibrium solution
$\ddagger$  =	a matrix or a vector containing the OWNS auxiliary values
Other symbols
$\bar{(\cdot)}$  =	a quantity associated with the equilibrium solution (e.g., Blasius base flow)
$\overline{(\cdot)}$  =	a complex conjugate value
$\hat{(\cdot)}$  =	a quantity in Fourier space

1 Introduction

Modern industry tools for transition prediction involve extrapolation from linear amplification through the $e^{N}$-method [1, 2] or the variable $N$-factor approach [3]. Fundamentally, transition can be studied using DNS [4, 5, 6, 7] or Large Eddy Simulation (LES) [8, 9], but these approaches are limited by their high computational cost. The NPSE entail a much lower computational cost [10], but convergence issues limit their use to the early stages of transition [11, 12]. Even in the linear case, however, PSE cannot accurately model non-modal and multi-modal instabilities [12], and its minimum step size for stability prevents capturing all length scales in the streamwise direction, which is particularly problematic when extending to the nonlinear regime.

OWNS overcomes these limitations, but at a greater computational cost. In particular, it properly removes all upstream-going waves so that the regularizations used by PSE are not required, thus allowing it to capture properly non-modal and multi-modal instabilities, and to take arbitrarily small steps to resolve all length scales [12]. The linear OWNS approach has previously been applied to two-dimensional (2D) and three-dimensional (3D) boundary-layer flows [13], and has proven useful hypersonic boundary layer flows [14, 15, 16] where classical stability analysis methods often fail.

The NOWNS methodology has been validated for Blasius boundary layer flows [17], and demonstrated to support non-modal instabilities, as well as stronger nonlinearities than NPSE [18]. In this paper, we refine the NOWNS methodology, validate it, and demonstrate its benefits with respect to the NPSE framework by showing that it is more robust to numerical noise and supports stronger nonlinearitie than NPSE, and that it supports non-modal instabilities.

2 Method

The OWNS procedure was first developed as a method for constructing well-posed one-way approximations of linear hyperbolic systems with slowly-varying coefficients in the direction of integration [19], allowing the equations to be solved as a spatial-initial-value problem, leading to a reduced computational cost. The framework has also been applied to the Navier-Stokes equations, linearized about a slowly-varying baseflow, to perform linear stability analysis of jets and boundary layer flows [19, 20, 14, 15, 16, 13]. Here we generalize the OWNS methodology to support nonlinear effects.

The non-dimensional, compressible Navier-Stokes equations for an ideal gas can be written as

$$\frac{D\nu}{Dt}-\nu(\nabla\cdot\bm{u}) =0,$$ (1a)

$$\frac{D\bm{u}}{Dt}+\nu\nabla p =\frac{1}{Re}\nu\nabla\cdot\tau,$$ (1b)

$$\frac{Dp}{Dt}+\gamma p(\nabla\cdot\bm{u}) =\frac{1}{Re}(\gamma-1)[(\nabla\bm{u}):\tau]+\frac{\gamma}{PrRe}(\nu\nabla^{2}p+2\nabla\nu\cdot\nabla p+p\nabla^{2}\nu),$$ (1c)

for the stress tensor

$$\tau=\mu\Big{(}\nabla\bm{u}+(\nabla\bm{u})^{T}\Big{)}-\Big{(}\frac{2}{3}\mu-\kappa\Big{)}(\nabla\cdot\bm{u})I,$$

where $\mu$ is the dynamic viscosity and $\kappa$ is the bulk viscosity. We take the bulk viscosity to be zero ($\kappa=0$) and we introduce the Reynolds number, the Prandtl number, and the Blasius length scale

$$Re=\frac{\delta_{0}^{*}a_{\infty}^{*}}{\nu_{\infty}^{*}\mu_{\infty}^{*}},\quad Pr=\frac{c_{p,\infty}^{*}\mu_{\infty}^{*}}{k_{\infty}^{*}},\quad\delta_{0}^{*}=\frac{\nu_{\infty}^{*}x_{0}^{*}}{U_{\infty}^{*}}.$$

Here $\nu_{\infty}^{*}$ denotes the free-stream specific volume, $\mu_{\infty}^{*}$ the free-stream dynamic viscosity, $k_{\infty}^{*}$ the free-stream thermal diffusivity, $c_{p,\infty}^{*}$ the free-stream specific heat capacity at constant pressure, $a_{\infty}^{*}$ the free-stream sound speed, while $x$, $y$, $z$, correspond to the streamwise, transverse and spanwise directions, respectively. At low Mach numbers, temperature fluctuations are minimal, so we take the fluid properties ($\mu$, $k$, and $c_{p}$) to be constant.

2.1 Linear OWNS

The NOWNS procedure follows closely the linear OWNS procedure, which we introduce first. We wish to study the evolution of a time-varying disturbance, $\bm{q}^{\prime}$, to a time-invariant equilibrium solution, $\bar{\bm{q}}$, so we decompose the flow into terms associated with $\bar{\bm{q}}$ and $\bm{q}^{\prime}$, as described in appendix A, yielding the stability equations (31), for the vector of primitive variables $\bm{q}=(\nu,u,v,w,p)$. To simplify the exposition, we temporarily neglect the streamwise diffusion terms by introducing the exogeneous forcing function $\bm{f}=\bm{f}_{\mathrm{forcing}}+\bm{f}_{\mathrm{viscous}}$ where

$$\bm{f}_{\mathrm{viscous}}=-B_{x}(\bar{\bm{q}})\frac{\partial\bm{q}^{\prime}}{\partial x}-B_{xx}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x^{2}}-B_{xy}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial y}-B_{xz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial z},$$

represents streamwise diffusion terms associated with the disturbance variable, while $\bm{f}_{\mathrm{forcing}}$ is an arbitrary forcing function. We then neglect the nonlinear term and re-write equation (31) to obtain the linear stability equation

$$A_{x}(\bar{\bm{q}})\frac{\partial\bm{q}^{\prime}}{\partial x}=\mathcal{L}(\bar{\bm{q}})\bm{q}^{\prime}+\bm{f}.$$ (2)

Previous work OWNS has often neglected $\bm{f}_{\mathrm{viscous}}$ because this simplifies the procedure without significantly impacting its accuracy [19, 13, 14], but we have found that this term impacts more significantly the nonlinear calculation, as discussed in appendix C, so we re-introduce it (approximately) in section 2.1.5.

In NOWNS, we consider a system of linear OWNS equations coupled together through the nonlinear term, which acts as an inhomogeneous forcing function taking the place of $\bm{f}$ in (2). The OWNS outflow (OWNS-O) approach [19] supports only homogeneous equations, so we must instead consider either the OWNS projection (OWNS-P) approach [20] or the OWNS recursive (OWNS-R) approach [21], which both support inhomogeneous equations. OWNS-R entails a reduced computational cost compared to OWNS-P, but we have found OWNS-P to be more robust, so we use it to develop the NOWNS approach.

2.1.1 Semi-discrete equations

Throughout this work, we assume that our disturbances are periodic in the spanwise direction with wavenumber $\beta$ and in time with frequency $\omega$ so that we can expand our disturbances as Fourier series

$$\bm{q}^{\prime}(x,y,z,t)=\sum_{m,n=-\infty}^{\infty}\hat{\bm{q}}_{mn}(x,y)e^{i(n\beta z-m\omega t)}.$$ (3)

In the linear case, all Fourier modes evolve independently, and we consider a single disturbance of the form $\bm{q}^{\prime}(x,y,z,t)=\hat{\bm{q}}(x,y)e^{i(n\beta z-m\omega t)}$, while we discretize in the wall-normal direction using a 4th-order central finite differences which we represent using $D\approx\partial/\partial y$. Our semi-discrete linear operator is then

$$L(\bar{\bm{q}})=i\omega I-[A_{y}(\bar{\bm{q}})+B_{y}(\bar{\bm{q}})]D-i\beta[A_{z}(\bar{\bm{q}})+B_{z}(\bar{\bm{q}})]-C(\bar{\bm{q}})-B_{yy}(\bar{\bm{q}})D^{2}+\beta^{2}B_{zz}(\bar{\bm{q}})-i\beta B_{yz}(\bar{\bm{q}})=0,$$ (4)

and we obtain

$$A_{x}(\bar{\bm{q}})\frac{\partial\bm{q}^{\prime}}{\partial x}=L(\bar{\bm{q}})\bm{q}^{\prime}+\bm{f},$$ (5)

a system of ODEs in $x$ comprising $N_{v}$ variables, where $N_{v}=5N_{y}$ ($N_{v}=4N_{y}$) in 3D (2D), for the $N_{y}$ grid points in $y$.

2.1.2 Parabolization using the OWNS projection procedure

The above ODEs in $x$ contain components that propagate both upstream and downstream, and cannot be stably integrated without further intervention. Therefore, we remove upstream effects using a projection operator that we apply in the characteristic variables, $\bm{\phi}=T\bm{q}^{\prime}$, where $T$ are the eigenvectors of $A_{x}$, while $\tilde{A}_{x}=TA_{x}T^{-1}$ are the eigenvalues. Here we have dropped the argument $\bar{\bm{q}}$ for brevity. We transform our equations to characteristic variables as

$$\tilde{A}_{x}\frac{\partial\hat{\bm{\phi}}}{\partial x}=\tilde{L}\hat{\bm{\phi}}+\hat{\bm{f}}_{\phi}.$$ (6)

with $\tilde{L}=TLT^{-1}-TA_{x}\frac{\partial T^{-1}}{\partial x}$ and $\hat{\bm{f}}_{\phi}=T\hat{\bm{f}}$, and re-organize the diagonal matrix $\tilde{A}_{x}$

$$\tilde{A}_{x}=\begin{bmatrix}\tilde{A}_{++}&0&0\\ 0&\tilde{A}_{--}&0\\ 0&0&\tilde{A}_{00}\end{bmatrix},$$ (7)

for the $N_{+}$ positive eigenvalues $\tilde{A}_{++}>0$, the $N_{-}$ negative eigenvalues $\tilde{A}_{--}<0$, and the $N_{0}$ zero eigenvalues $\tilde{A}_{00}=0$, where $N_{0}+N_{-}+N_{+}=N_{v}$. We further define

$$\tilde{A}_{\pm\pm}=\begin{bmatrix}\tilde{A}_{++}&0\\ 0&\tilde{A}_{--}\end{bmatrix},\quad\tilde{L}_{\pm\pm}=\begin{bmatrix}\tilde{L}_{++}&\tilde{L}_{+-}\\ \tilde{L}_{-+}&\tilde{L}_{--}\end{bmatrix},\quad\tilde{L}_{\pm 0}=\begin{bmatrix}\tilde{L}_{+0}\\ \tilde{L}_{-0}\end{bmatrix},\quad\tilde{L}_{0\pm}=\begin{bmatrix}\tilde{L}_{0+}&\tilde{L}_{0-}\end{bmatrix},\quad\hat{\bm{\phi}}=\begin{bmatrix}\hat{\bm{\phi}}_{\pm}\\ \hat{\bm{\phi}}_{0}\end{bmatrix},\quad\hat{\bm{f}}_{\phi}=\begin{bmatrix}\hat{\bm{f}}_{\phi,\pm}\\ \hat{\bm{f}}_{\phi,0}\end{bmatrix},$$

so that our equations become

$$\tilde{A}_{\pm\pm}\frac{\partial\hat{\bm{\phi}}_{\pm}}{\partial x} =\tilde{L}_{\pm\pm}\hat{\bm{\phi}}_{\pm}+\tilde{L}_{\pm 0}\hat{\bm{\phi}}_{0}+\hat{\bm{f}}_{\phi,\pm},$$ (8a)

$$0 =\tilde{L}_{0\pm}\hat{\bm{\phi}}_{\pm}+\tilde{L}_{00}\hat{\bm{\phi}}_{0}+\hat{\bm{f}}_{\phi,0},$$ (8b)

which is a differential algebraic equation (DAE) of index 1. We can use the algebraic constraint (8b) to obtain $\hat{\bm{\phi}}_{0}=-\tilde{L}_{00}^{-1}[\tilde{L}_{0\pm}\hat{\bm{\phi}}_{\pm}+\hat{\bm{f}}_{\phi,0}]$ so that we can re-write our DAE as an ODE

$$\frac{\partial\hat{\bm{\phi}}_{\pm}}{\partial x}=M\hat{\bm{\phi}}_{\pm}+\hat{\bm{g}},$$ (9)

for $M=\tilde{A}_{\pm\pm}^{-1}[\tilde{L}_{\pm\pm}-\tilde{L}_{\pm 0}\tilde{L}_{00}^{-1}\tilde{L}_{0\pm}]$, and $\hat{\bm{g}}=\tilde{A}_{\pm\pm}^{-1}[\hat{\bm{f}}_{\phi,\pm}-\tilde{L}_{\pm 0}L_{00}^{-1}\hat{\bm{f}}_{\phi,0}]$.

The upstream- and downstream-going modes of (9) can be determined based on the eigenvalues of $M$, according to Brigg’s criterion [22], which can then be used to introduce well-posed one-way equations, according to the criterion of Kreiss [23], as described in Towne and Colonius (2015) [19] and Towne et al. (2022) [20]. We expand the solution as a linear combination of the eigenvectors $\hat{\bm{\phi}}_{\pm}=V\hat{\bm{\psi}}=\sum_{k=1}^{N}\bm{v}^{(k)}\hat{\psi}_{k},$ where $V$ are the eigenvectors $M=VDV^{-1}$. From Brigg’s criterion, we know that $M$ has $N_{+}$ downstream- and $N_{-}$ upstream-going modes, so that we can further partition $V$ into its downstream- ($V_{+}$) and upstream-going ($V_{-}$) components yielding the downstream- and upstream-going solutions $\hat{\bm{\phi}}_{\pm}^{\prime}=V_{+}\hat{\bm{\psi}}_{+}$ and $\hat{\bm{\phi}}_{\pm}^{\prime\prime}=V_{-}\hat{\bm{\psi}}_{-}$, respectively. We then use the eigensystem of $M$ to define the projection operator

$$P=V\begin{bmatrix}I_{++}&0\\ 0&0\end{bmatrix}V^{-1},$$

which retains the downstream-going modes associated with $\bm{\psi}_{+}$ while removing the upstream-going modes associated with $\bm{\psi}_{-}$. We apply $P$ to our ODE (9) and use linearity to assert that

$$\frac{\partial\hat{\bm{\phi}}^{\prime}_{\pm}}{\partial x}+\frac{\partial\hat{\bm{\phi}}_{\pm}^{\prime\prime}}{\partial x}=P[M\hat{\bm{\phi}}_{\pm}+\hat{\bm{g}}]+[I-P][M\hat{\bm{\phi}}_{\pm}+\hat{\bm{g}}],$$

yielding the equation for the downstream-going solution

$$\frac{\partial\hat{\bm{\phi}}_{\pm}^{\prime}}{\partial x}=P[M\hat{\bm{\phi}}_{\pm}^{\prime}+\hat{\bm{g}}],$$ (10)

where we used that $P$ and $M$ commute, as shown by Towne et al. [20].

2.1.3 Approximate projection operator

Explicitly using Brigg’s criterion to identify upstream- and downstream-going modes in a rigorous way would be computationally expensive and prone to numerical error, so we instead apply the projection operator approximately using a recursive filter [19, 20].

We define the residuals $\hat{\bm{r}}_{\pm}(\bm{\phi})=\tilde{A}_{x,\pm\pm}^{-1}[\tilde{L}_{\pm\pm}\hat{\bm{\phi}}_{\pm}+\tilde{L}_{\pm 0}\hat{\bm{\phi}}_{0}+\hat{\bm{f}}_{\phi,\pm}],$ and $\hat{\bm{r}}_{0}(\bm{\phi})=\tilde{L}_{0\pm}\hat{\bm{\phi}}_{\pm}+\tilde{L}_{00}\hat{\bm{\phi}}_{0}+\hat{\bm{f}}_{\phi,0}$ based on (8), and following the approach of Towne et al. [20], we can apply $P$ approximately to the residual using the recursive filter

$$\hat{\bm{r}}_{+}^{(-N_{b})} =0$$ (11a)

$$(\tilde{L}-i\beta_{-}^{(j)}\tilde{A})\hat{\bm{r}}^{(-j)}-(\tilde{L}-i\beta_{+}^{(j)}\tilde{A})\hat{\bm{r}}^{(-j-1)} =0,\quad j=1,\dots,N_{b}-1$$ (11b)

$$(\tilde{L}-i\beta_{-}^{(0)}\tilde{A})\hat{\bm{r}}^{(0)}-(\tilde{L}-i\beta_{+}^{(0)}\tilde{A})\hat{\bm{r}}^{(-1)} =(\tilde{L}-i\beta_{-}^{(0)}\tilde{A})\hat{\bm{r}}$$ (11c)

$$\hat{\bm{r}}_{0}^{(0)} =\hat{\bm{r}}_{0},$$ (11d)

$$(\tilde{L}-i\beta_{+}^{(j)}\tilde{A})\hat{\bm{r}}^{(j)}-(\tilde{L}-i\beta_{-}^{(j)}\tilde{A})\hat{\bm{r}}^{(j+1)} =0,\quad j=0,\dots,N_{b}-1$$ (11e)

$$\hat{\bm{r}}_{-}^{(N_{b})} =0,$$ (11f)

where $\hat{\bm{r}}_{0}=0$ when the algebraic constraint (8b) is satisfied. Here, $\{\beta_{\pm}^{(j)}\}_{j=0}^{N_{b}-1}$ are termed the recursion parameters, while $\{\hat{\bm{r}}^{(j)}\}_{j=-N_{b}}^{N_{b}}$ are termed the auxiliary variables, where $N_{\mathrm{aux}}=(2N_{v}+1)N_{b}$. We introduce the vector auxiliary variables $\hat{\bm{r}}_{\mathrm{aux}}\in\mathcal{C}^{N_{\mathrm{aux}}}$, and the approximate projection operators $P_{1}\in\mathbb{C}^{N_{\mathrm{aux}}\times N_{v}}$, $P_{2}\in\mathbb{C}^{N_{\mathrm{aux}}\times N_{\mathrm{aux}}}$, $P_{3}\in\mathbb{R}^{N_{\pm}\times N_{\mathrm{aux}}}$, where $P_{1}\hat{\bm{r}}$ and $P_{2}\hat{\bm{r}}_{\mathrm{aux}}$ give the right- and left-hand sides of (11), respectively, while $P_{3}$ extracts $\hat{\bm{r}}_{\pm}^{(0)}$ from $\hat{\bm{r}}_{\mathrm{aux}}$ as $\hat{\bm{r}}_{\pm}^{(0)}=P_{3}\hat{\bm{r}}_{\mathrm{aux}}$. The action of the approximate projection operator on the DAE (8) can expressed compactly as

$$\frac{\partial\hat{\bm{\phi}}_{\pm}^{\prime}}{\partial x} =P_{3}\hat{\bm{r}}_{\mathrm{aux}},$$ (12a)

$$P_{2}\hat{\bm{r}}_{\mathrm{aux}} =P_{1}\begin{bmatrix}\tilde{A}_{\pm\pm}^{-1}[\tilde{L}_{\pm\pm}\hat{\bm{\phi}}^{\prime}_{\pm}+\tilde{L}_{\pm 0}\hat{\bm{\phi}}^{\prime}_{0}+\hat{\bm{f}}_{\phi,\pm}]\\ \tilde{L}_{0\pm}\hat{\bm{\phi}}^{\prime}_{\pm}+\tilde{L}_{00}\hat{\bm{\phi}}^{\prime}_{0}+\hat{\bm{f}}_{\phi,0}\end{bmatrix},$$ (12b)

$$0 =\tilde{L}_{0\pm}\hat{\bm{\phi}}^{\prime}_{\pm}+\tilde{L}_{00}\hat{\bm{\phi}}^{\prime}_{0}+\hat{\bm{f}}_{\phi,0}.$$ (12c)

2.1.4 Fully-discrete equations

We define $A^{\ddagger},L^{\ddagger}\in\mathbb{R}^{(N_{v}+N_{\mathrm{aux}})\times(N_{v}+N_{\mathrm{aux}})}$ and $\hat{\bm{\phi}}^{\ddagger},\hat{\bm{f}}_{\phi}^{\ddagger}\in\mathbb{C}^{N_{v}+N_{\mathrm{aux}}}$ such that

$$A^{\ddagger}=\begin{bmatrix}I_{\pm\pm}&0&0\\ 0&0&0\\ 0&0&0\end{bmatrix},\quad L^{\ddagger}=\begin{bmatrix}0&0&\Delta xP_{3}\\ P_{1}\begin{bmatrix}\tilde{A}_{\pm\pm}^{-1}L_{\pm\pm}\\ L_{0\pm}\end{bmatrix}&P_{1}\begin{bmatrix}\tilde{A}_{\pm\pm}^{-1}L_{\pm 0}\\ L_{00}\end{bmatrix}&-P_{2}\\ L_{0\pm}&L_{00}&0\end{bmatrix},\quad\hat{\bm{\phi}}^{\ddagger}=\begin{bmatrix}\bm{\phi}_{\pm}^{\prime}\\ \bm{\phi}_{0}\\ \bm{r}_{\mathrm{aux}}\end{bmatrix},\quad\hat{\bm{f}}_{\phi}^{\ddagger}=\begin{bmatrix}0\\ P_{1}\begin{bmatrix}\tilde{A}_{\pm\pm}^{-1}\hat{\bm{f}}_{\phi,\pm}\\ \hat{\bm{f}}_{\phi,0}\end{bmatrix}\\ \hat{\bm{f}}_{\phi,0}\\ \end{bmatrix},$$

and apply an $s$-order BDF scheme to obtain the linear OWNS equation

$$[c^{(0)}A^{\ddagger}-L^{\ddagger(k+1)}]\hat{\bm{\phi}}^{\ddagger(k+1)}=-\sum_{l=1}^{s-1}c^{(l)}A^{\ddagger}\hat{\bm{\phi}}^{\ddagger(k+1-l)}+\hat{\bm{f}}_{\phi}^{\ddagger(k+1)},$$ (13)

as has been done in previous work on linear OWNS [13, 14].

2.1.5 Streamwise diffusion terms

Following discretization in the wall-normal direction and trasnformation to characteristic variables, our streamwise diffusion terms become

$$\hat{\bm{f}}_{\phi,\mathrm{viscous}}=-\tilde{B}_{xx}\frac{\partial^{2}\hat{\bm{\phi}}}{\partial x^{2}}-\tilde{B}_{x}\frac{\partial\hat{\bm{\phi}}}{\partial x}-\tilde{B}\hat{\bm{\phi}}.$$

for $\tilde{B}_{xx}=TB_{xx}T^{-1}$, $\tilde{B}_{x}=2TB_{xx}\frac{\partial T^{-1}}{\partial x}+T[B_{x}+B_{xy}D+i\beta B_{xz}]T^{-1}$, and $\tilde{B}=TB_{xx}\frac{\partial^{2}T^{-1}}{\partial x^{2}}+T[B_{x}+B_{xy}D+i\beta B_{xz}]\frac{\partial T^{-1}}{\partial x}$. We discretize the second-derivative using a second-order backward difference, while we discretize the first-derivative using the BDF scheme to obtain

$$\hat{\bm{f}}^{(k+1)}_{\phi,\mathrm{viscous}}=-\tilde{B}_{xx}\frac{\hat{\bm{\phi}}^{(k+1)}-2\hat{\bm{\phi}}^{(k)}+\hat{\bm{\phi}}^{(k-1)}}{(\Delta x)^{2}}-\Big{(}\sum_{l=0}^{s-1}c^{(l)}\tilde{B}_{x}\hat{\bm{\phi}}^{(k+1-l)}\Big{)}-\tilde{B}\hat{\bm{\phi}}^{(k+1)},$$

and we add this term back into the fully-discrete OWNS equations (13).

2.2 Nonlinear OWNS

Whereas infinitesimal disturbances (in the linear case) evolve independently from each other so that each Fourier mode can be considered separately, finite amplitude disturbances (in the nonlinear case) are coupled through the nonlinear term. Since it is not feasible to consider an infinite number of Fourier modes, we truncate the Fourier (3) series

$$\bm{q}^{\prime}(x,y,z,t)=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\hat{\bm{q}}_{mn}e^{i(n\beta z-m\omega t)},$$ (14)

resulting in $(2M+1)\times(2N+1)$ Fourier modes. However, we require that $\bm{q}^{\prime}$ be real-valued, which yields the constraint $\hat{q}_{mn}=\overline{\hat{q}_{-mn}}$, so that we need only track $(M+1)\times(2N+1)$ Fourier modes. We can additionally introduce a spanwise symmetry condition, as described in section 2.2.5, which further reduces the number of modes to $(M+1)\times(N+1)$.

We use our assumption of periodicity to obtain the Fourier series

$$\tilde{L}(\bar{\bm{q}})\bm{\phi}=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\hat{L}_{mn}\hat{\bm{\phi}}_{mn}e^{i(n\beta z-m\omega t)},\quad\tilde{\bm{F}}(\bm{\phi})=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\hat{\bm{F}}_{mn}(\bm{\phi})e^{i(n\beta z-m\omega t)},$$

where the transformation of $\bm{F}(\bm{\phi})$ to characteristic variables, $\tilde{\bm{F}}(\bm{\phi})$, mimics the transformation of $L(\bar{\bm{q}})$. The Fourier modes are mutually orthogonal, yielding the following system of equations

$$\tilde{A}\frac{\partial\hat{\bm{\phi}}_{mn}}{\partial x}=\hat{L}_{mn}\hat{\bm{\phi}}_{mn}+\hat{\bm{F}}_{mn}(\bm{\phi})+\hat{\bm{f}}_{\phi,mn},\quad\forall m\in\mathbb{Z}_{M},\quad\forall n\in\mathbb{Z}_{N},$$ (15)

for $\mathbb{Z}_{M}\equiv\{x\in\mathbb{Z}|-M\leq x\leq M\}$ and $\mathbb{Z}_{N}\equiv\{x\in\mathbb{Z}|-N\leq x\leq N\}$, where $\mathbb{Z}$ is the set of all integers. We follow a procedure that mimics the linear OWNS approach to obtain

$$\frac{\partial\bm{\phi}_{\pm,mn}^{\prime}}{\partial x}=\hat{P}_{mn}[\hat{M}_{mn}\hat{\bm{\phi}}^{\prime}_{\pm,mn}+\hat{\bm{g}}_{mn}(P\bm{\phi}_{\pm}^{\prime})],\quad\forall m\in\mathbb{Z}_{M},\quad\forall n\in\mathbb{Z}_{N}.$$ (16)

In the linear case, $P$ and $M$ commute so that two one-way parabolic equations can be solved to recover the full elliptic solution [20]. However, this property does not apply in the nonlinear case because the $P$ does not commute with the nonlinear term ($P\bm{g}(\bm{\phi}_{\pm})\neq P\bm{g}(P\bm{\phi}_{\pm})$ in general), so that (16) removes the upstream effect, $\bm{\phi}^{\prime\prime}$, from the nonlinear term, and when we sum the upstream- and downstream-going equations together, we do not recover the elliptic equation. This is a reasonable choice for convective instabilities, where the disturbances travel primarily in one direction, and we verify a posteriori that we match closely DNS results in the literature.

2.2.1 Fully-discrete equations

We define $L_{mn}^{\ddagger}\in\mathbb{R}^{(N_{v}+N_{\mathrm{aux}})\times(N_{v}+N_{\mathrm{aux}})}$ and $\hat{\bm{\phi}}_{mn}^{\ddagger},\hat{\bm{f}}_{mn,\phi}^{\ddagger}\in\mathbb{C}^{N_{v}+N_{\mathrm{aux}}}$, as in the linear case, for all $m\in\mathbb{Z}_{M}$ and $n\in\mathbb{Z}_{N}$, and further introduce the nonlinear term $\hat{\bm{F}}_{mn}^{\ddagger}\in\mathbb{C}^{N_{v}+N_{\mathrm{aux}}}$ such that it mimics the definition of the forcing function $\hat{\bm{f}}_{mn,\phi}^{\ddagger}$. Then we discretize the first-derivatives in both the linear and nonlinear terms using the BDF scheme yielding the fully-discrete nonlinear system of equations

$$\sum_{l=0}^{s-1}c^{(l)}A^{\ddagger}\hat{\bm{\phi}}_{mn}^{\ddagger(k+1-l)}=\hat{L}_{mn}^{\ddagger(k+1)}\hat{\bm{\phi}}_{mn}^{\ddagger(k+1)}+\hat{\bm{F}}_{mn}^{\ddagger(k+1)}+\hat{\bm{f}}_{\phi,mn}^{\ddagger(k+1)},\quad\forall m\in\mathbb{Z}_{M},\quad\forall n\in\mathbb{Z}_{N}.$$ (17)

2.2.2 Pseudo-spectral method

We employ a pseudo-spectral method whereby we solve our equations in Fourier space, while we evaluate the nonlinear terms in physical space. Given the Fourier coefficients $\hat{\bm{\phi}}_{mn}$ for $m=0,\dots,2M$ and $n=0,\dots,N$, we can use the inverse discrete Fourier transform (iDFT) to compute the solution in physical space as

$$\bm{\phi}_{mn}=\frac{1}{2M+1}\frac{1}{2N+1}\sum_{k=0}^{2M}\sum_{l=0}^{2N}\hat{\bm{\phi}}_{kl}e^{i2\pi mk/(2M+1)}e^{i2\pi nl/(2N+1)}.$$ (18)

Then, we can evaluate the nonlinear terms in physical space, $\tilde{\bm{F}}_{mn}=\tilde{\bm{F}}(\bm{\phi}_{mn})$, and employ the discrete Fourier transform (DFT) to compute the Fourier components of the nonlinear terms as

$$\hat{\bm{F}}_{kl}=\sum_{m=0}^{2M}\sum_{n=0}^{2N}\tilde{\bm{F}}_{mn}e^{-i2\pi mk/(2M+1)}e^{-i2\pi nl/(2N+1)}.$$ (19)

In practice, the DFT and iDFT are performed using fast Fourier transform (FFT) libraries. At each step, we require the Euclidean norm of the residual for (17) be converged in both an absolute and relative sense to a tolerance of $10^{-10}$.

2.2.3 Nonlinear solution procedure and computational cost

We can take the derivative of the NOWNS equations with respect to the disturbance variable to obtain the Newton iteration (in 2D with $M=1$)

$$\begin{bmatrix}(c^{(0)}A^{\ddagger}-\hat{L}_{0}^{\ddagger})+\hat{J}^{\ddagger}_{0}&\overline{\hat{J}^{\ddagger}_{1}}&\hat{J}^{\ddagger}_{1}\\ \hat{J}^{\ddagger}_{1}&(c^{(0)}A^{\ddagger}-\hat{L}^{\ddagger}_{1})+\hat{J}^{\ddagger}_{0}&\overline{\hat{J}^{\ddagger}_{1}}\\ \overline{\hat{J}^{\ddagger}_{1}}&\hat{J}^{\ddagger}_{1}&(c^{(0)}A^{\ddagger}-\hat{L}_{-1}^{\ddagger})+\hat{J}^{\ddagger}_{0}\end{bmatrix}\begin{bmatrix}\Delta\hat{\bm{\phi}}_{0}^{\ddagger}\\ \Delta\hat{\bm{\phi}}_{1}^{\ddagger}\\ \overline{\Delta\hat{\bm{\phi}}_{1}}^{\ddagger}\end{bmatrix}=\begin{bmatrix}\hat{r}_{0}^{\ddagger}\\ \hat{r}_{1}^{\ddagger}\\ \overline{\hat{r}_{1}}^{\ddagger}\end{bmatrix},$$ (36)

as discussed in appendix B, where $J_{p-m}^{\ddagger}=-\partial\hat{\bm{F}}^{\ddagger}_{m}/\partial\hat{\bm{\phi}}_{p}^{\ddagger}$. If we neglect $J_{m}^{\ddagger}$ for $m\neq 0$, then we obtain

$$\begin{bmatrix}(c^{(0)}A^{\ddagger}-\hat{L}_{0}^{\ddagger})+\hat{J}_{0}^{\ddagger}&0&0\\ 0&(c^{(0)}A^{\ddagger}-\hat{L}_{1}^{\ddagger})+\hat{J}_{0}^{\ddagger}&0\\ 0&0&(c^{(0)}A^{\ddagger}-\hat{L}_{-1}^{\ddagger})+\hat{J}_{0}^{\ddagger}\end{bmatrix}\begin{bmatrix}\Delta\hat{\bm{\phi}}_{0}^{\ddagger}\\ \Delta\hat{\bm{\phi}}_{1}^{\ddagger}\\ \overline{\Delta\hat{\bm{\phi}}_{1}}^{\ddagger}\end{bmatrix}=\begin{bmatrix}\hat{r}_{0}^{\ddagger}\\ \hat{r}_{1}^{\ddagger}\\ \overline{\hat{r}_{1}}^{\ddagger}\end{bmatrix}.$$ (20)

An extension to 3D for Newton’s method (36) is straightforward, while system (21) can be written in 3D as

$$[c^{(0)}A^{\ddagger}-\hat{L}_{mn}^{\ddagger}+\hat{J}_{0}^{\ddagger}]\hat{\bm{\phi}}_{mn}^{\ddagger(k+1)}=-\sum_{l=1}^{s-1}c^{(l)}A^{\ddagger}\hat{\bm{\phi}}_{mn}^{\ddagger(k+1-l)}+\hat{\bm{F}}_{mn}^{\ddagger(k+1)}+\hat{\bm{f}}_{mn}^{\ddagger(k+1)},$$ (21)

for $m=-M,\dots,M$ and $n=-N,\dots,N$. If we neglect $J^{\ddagger}_{m}$ entirely, then we obtain

$$[c^{(0)}A^{\ddagger}-\hat{L}_{mn}^{\ddagger}]\hat{\bm{\phi}}_{mn}^{\ddagger(k+1)}=-\sum_{l=1}^{s-1}c^{(l)}A^{\ddagger}\hat{\bm{\phi}}_{mn}^{\ddagger(k+1-l)}+\hat{\bm{F}}_{mn}^{\ddagger(k+1)}+\hat{\bm{f}}_{mn}^{\ddagger(k+1)},$$ (22)

for $m=-M,\dots,M$ and $n=-N,\dots,N$. To solve (21) or (22), we mimic the NPSE solution procedure [24, 25, 11]: we take the lower-upper (LU) decomposition of $[c^{(0)}A^{\ddagger}-\hat{L}_{mn}^{\ddagger}]$ to solve for $\hat{\bm{\phi}}_{m}^{\ddagger(k+1)}$, then we update the nonlinear term, and repeat until the residual is converged. We typically prefer solving (21) over (22) because it reduces the number of iterations to convergence, without increasing the computational cost. Alternatively, if we accept a larger computational cost, then we can further reduce the number of iterations by solving (36).

The linear OWNS system (13) comprises $N_{v}+N_{\mathrm{aux}}=N_{v}+(2N_{v}+1)N_{b}$ equations, where $N_{v}$ scales with $N_{y}$, so that the computational cost to solve this system using a direct multifrontal solver (LU decomposition) scales as $\mathcal{O}(N_{y}^{a}N_{b}^{a})$, where $a$ is a problem dependent coefficients that depend on the sparsity pattern. Theoretically, $1<a\leq 3$ and we typically observe $a\approx 1.5$ for 2D problems and $a\approx 2$ for 3D ones, see Towne et al. (2022) [20] for further details. The cost to integrate over $N_{x}$ stations is then $\mathcal{O}(N_{x}N_{y}^{a}N_{b}^{a})$, while a global method entails a cost of $\mathcal{O}(N_{x}^{a}N_{y}^{a})$, so that OWNS is more efficient for $N_{b}\ll N_{x}$. The nonlinear OWNS system (17) comprises $(M+1)\times(2N+1)\times(N_{v}+N_{\mathrm{aux}})$ equations, and can be solved using Newton’s method (36), which entails a cost of $\mathcal{O}(M^{a}N^{a}N_{y}^{a}N_{b}^{a})$. However, this cost can be reduced by instead solving (21) or (22), since we can perform the LU decomposition of $[c^{(0)}A^{\ddagger}-\hat{L}_{m}^{\ddagger}]$ separately for each Fourier mode, yielding a cost of $\mathcal{O}(MNN_{y}^{a}N_{b}^{a})$. Therefore, the cost to integrate the NOWNS equations is $\mathcal{O}(N_{x}MNN_{y}^{a}N_{b}^{3})$ using (21) or (22), while it increases to $\mathcal{O}(N_{x}M^{a}N^{a}N_{y}^{a}N_{b}^{3})$ for Newton’s method (36), as compared to the cost $\mathcal{O}(N_{x}^{a}M^{a}N^{a}N_{y}^{a})$ for nonlinear global methods.

Although (21) and (22) entail a lower computational cost, these methods fail for strong nonlinearities and we must instead employ Newton’s method (36). In practice we implement a hybrid approach whereby we first solve (21) to harness it’s reduced computational cost, and then if more than $\mathcal{O}(100)$ iterations have elapsed, we switch to Newton’s method to harness it’s better convergence properties. We compare the performance of these three methods in appendix E.

To reduce the computational cost of Newton’s method, we re-use the LU factors from the first iteration as a pre-conditioner for the Generalized Minimal Residual Method (GMRES), and we note that a similar procedure was performed for linear OWNS in Araya et al. [26]. We also tested the block-Jacobi relaxation method, which entails a cost of $\mathcal{O}(MNN_{y}^{a}N_{b}^{a})$, and although we found this approach converged quickly in the early stages of the march, it failed as the nonlinearity grew stronger and we did not pursue relaxation methods further.

2.2.4 Special treatment of the zero-frequency modes

We have three options for obtaining a stable march for the zero-frequency modes, as depicted in table 2. In the first approach, which mimics how the zero-frequency modes are handled by NPSE [25, 11], we neglect the streamwise presssure gradient, $\partial p_{0n}/\partial x$, and the streamwise diffusion terms, $\partial^{2}q_{0n}/\partial x^{2}$, associated with the zero-frequency modes. However, in general we would prefer to avoid neglecting terms, and to instead parabolize these equations using the OWNS approach. In the second approach, we include the streamwise diffusion terms (but exclude the streamwise pressure gradient) associated with the zero-frequency modes, and parabolize the equations for all modes using OWNS. This approach yields a stable spatial march that agrees well with DNS, but offers no advantages over the first approach: for all of the cases examined in this paper, neglecting the streamwise diffusion terms associated with the zero-frequency modes does not change substantially the results of the NOWNS calculation. In the third approach, we include both the streamwise diffusion terms and the streamwise pressure gradient associated with the zero-frequency modes. Although the OWNS approach yields a stable spatial march, we have found that this approach produces inferior comparisons to DNS solutions from the literature, as discussed in appendix D.

In summary, including the streamwise pressure gradient for the zero-frequency modes worsens agreement with DNS, while including the streamwise diffusion terms for the zero-frequency modes has nearly no impact on the solution. Since neglecting these terms is more computationally efficient while still providing excellent agreement with DNS, we recommend that these terms be neglected.

	Include $\partial p_{0n}/\partial x$?	Include $\partial^{2}q_{0n}/\partial x^{2}$?	Parabolized using OWNS?	Agreement with DNS?
1	No	No	No	Yes
2	No	Yes	Yes	Yes
3	Yes	Yes	Yes	No

Table 2: Three approaches to parabolizing the zero-frequency modes

2.2.5 Spanwise symmetry

If the disturbances are symmetric, then we can enforce a symmetry condition to reduce the number of equations from $(M+1)\times(2N+1)$ to $(M+1)\times(N+1)$. All variables have even-symmetry ($\hat{\nu}_{m,-n}^{\prime}=\hat{\nu}_{m,n}^{\prime}$, $\hat{u}_{m,-n}^{\prime}=\hat{u}_{m,n}^{\prime}$, $\hat{v}_{m,-n}^{\prime}=\hat{v}_{m,n}^{\prime}$, $\hat{p}_{m,-n}^{\prime}=\hat{p}_{m,n}^{\prime}$), with the exception of the $w$-velocity which has odd-symmetry ($\hat{w}_{m,-n}^{\prime}=-\hat{w}_{m,n}^{\prime}$).

2.2.6 Boundary conditions

At the wall, we impose no-slip isothermal boundary conditions ($u^{\prime}=v^{\prime}=w^{\prime}=T^{\prime}=0$) and solve for the specific volume, $\nu^{\prime}$, using the (nonlinear) continuity equation. At the far-field boundary, we impose 1D (in $y$) inviscid Thompson characteristic boundary conditions to minimize spurious numerical reflections [27], which we implement using the linearized boundary-layer flow equations.

Some previous work on NPSE has used similar characteristic far-field boundary conditions [11]. Chang et al. (1993) [25] used the far-field boundary condition $\hat{\bm{q}}_{mn}^{\prime}(y_{\max})=0$ for $m,n\neq 0$. As the boundary layer must be allowed to grow in the wall-normal direction (due to nonlinear interactions), they used $\partial\hat{v}_{00}^{\prime}/\partial y=0$ at $y_{\max}$ for the mean-flow distortion (MFD). The characteristic far-field boundary conditions are advantageous because they allow us to use the same boundary conditions for all Fouriers modes, instead of handling the MFD as a separate case.

2.2.7 Effects of the mean-flow distortion

In NOWNS, the disturbances interact to excite the MFD, so that the corrected mean flow, $\bar{\bm{q}}+\hat{\bm{q}}_{00}$, differs from the baseflow, $\bar{\bm{q}}$. We have experimented with linearizing about both the baseflow and the corrected mean flow and found that it does not have a large impact on the NOWNS calculation, as discussed in appendix F. Linearizing about the corrected mean flow increases the computational cost of NOWNS because the projection operators change between iterations, since the MFD changes, so the LU factorization must be updated. On the other hand, linearizing about the baseflow allows us to use the same LU factorization at each iteration because the baseflow is not affected by changes in the MFD. Since it is more computationally efficient to linearize only about the baseflow, we choose this approach moving forward.

2.2.8 Recursion parameters

The choice of recursion parameters is described in Sleeman et al. (2024) [18], and matches the recursion parameters used by Rigas et al. (2017) [13], which are based on the recursion parameters originally developed by Towne and Colonius (2015) [19]. We briefly discuss our strategy for selecting the recursion parameters when $\omega=0$ in appendix G.

3 Validation

We validate NOWNS by applying it to 2D and 3D Blasius boundary layer flows for which there are existing DNS and NPSE results in the literature. We choose a Mach number of $\mathrm{Ma}_{\infty}=0.1$ to study flows near the incompressible limit. In what follows, we use the following dimensionless quantities:

$$\mathrm{Re}_{x}=\frac{U_{\infty}^{*}x^{*}}{\nu_{\infty}^{*}\mu_{\infty}^{*}},\quad y=\frac{y^{*}}{\delta_{0}^{*}},\quad F=\frac{\omega^{*}\nu_{\infty}^{*}\mu_{\infty}^{*}}{U_{\infty}^{*2}},\quad b=\frac{\beta^{*}\nu_{\infty}^{*}\mu_{\infty}^{*}}{U_{\infty}^{*}},$$

where $\mathrm{Re}_{x}$ is the streamwise coordinate, $F$ is the temporal frequency, and $b$ is the spanwise wavenumber. We refer to modes according to their temporal frequency and their spanwise wave number as $(m,n)$, where $m$ refers to the frequency $\omega_{m}=m\omega$ and $n$ refers to the spanwise wave number $\beta_{n}=n\beta$. To be consistent with previous literature, we measure the amplitude of disturbances as

$$u^{\prime(m,n)}_{\max}(x)=c_{m,n}\ \max_{y}|u_{m,n}^{\prime}(x,y)|,\qquad c_{m,n}=\begin{cases}1&\mbox{$m=n=0$},\\ \sqrt{2}&\mbox{$m=0,n\neq 0$; $n=0,m\neq 0$},\\ 2&\mbox{otherwise}.\end{cases}$$ (23)

3.1 2D evolution of a Tollmien-Schlichting wave

We consider the test case developed by Bertolotti et al. [10] which has been widely used in the literature as a validation case for NPSE [28, 29, 30]. This case examines the evolution of a Tollmien-Schlichting (TS) wave excited at the inlet at a frequency $F=86\times 10^{-6}$ and amplitude $u^{\prime(1)}_{\max}(x_{0})=0.25\%$. All other Fourier components initially have zero-amplitude and are generated through nonlinear interactions with the TS wave. The grid extends over the domain $\mathrm{Re}_{x}\in[1.6\times 10^{5},10^{6}]$ and $y\in[0,75]$ with $4000$ stations evenly spaced in $x$ and 150 grid points in $y$, with the majority of the grid points clustered towards the wall, while the Fourier series is truncated at $M=5$ temporal modes.

Figure 1 compares NOWNS to DNS and NPSE for $u^{\prime(m,n)}_{\max}(x)$, while the $u$- and $v$-velocity profiles for the MFD and TS waves are shown in figure 2. Excellent agreement is obtained; the discrepancy for the MFD of $v$ can be attributed to the Dirichlet boundary conditions used in the DNS.

Figure 1: Amplitude of $u^{\prime}$ v.s. streamwise coordinate, $\mathrm{Re}_{x}$, for 2D evolution of TS wave.

(a) $u$-velocity for MFD

(b) $u$-velocity for TS wave

(c) $v$-velocity for MFD

(d) $v$-velocity for TS wave

Figure 2: $u$-and $v$-velocity profiles at streamwise coordinate $\mathrm{Re}_{x}=7.80\times 10^{5}$ for 2D evolution of TS wave.

3.2 3D oblique-wave breakdown

Next we consider the oblique-wave breakdown case studied by Joslin et al. (1993) [29], where transition is initiated by two oblique waves with opposite wave angle. They defined two cases–with small and large initial amplitude. In this section, for validation purposes, we consider the small amplitude case for which NPSE was previously successful. Oblique-wave breakdown has also been studied using both experiment and spatial DNS by Berlin et al. (1999) [31], while it was studied using NPSE for compressible flows by Chang and Malik (1994) [32]. We further note that whereas fundamental and subharmonic transition can be studied using Herbert’s secondary stability theory [33], no such theory exists for oblique-wave breakdown, so that either experiment or numerical simulation is necessary to study this transition scenario [29, 31].

The oblique waves have amplitude $u_{\max}^{\prime(1,1)}(x_{0})=\sqrt{2}\times 10^{-3}$ at the inlet at a frequency $F=86\times 10^{-6}$ and spanwise wavenumber $b=2/9\times 10^{-3}$, while the grid extends over the domain $\mathrm{Re}_{x}\in[2.74\times 10^{5},6.08\times 10^{5}]$ and $y\in[0,75]$ with $2000$ stations evenly spaced in $x$ and 150 grid points in $y$, while the Fourier series is truncated at $M=3$ and $N=4$. Figure 3 compares NOWNS to DNS and NPSE for $u_{\max}^{\prime(m,n)}(x)$, while figure 4 compares the $u$-velocity profiles at $\mathrm{Re}_{x}=4.69\times 10^{5}$; we see that we have excellent agreement for between the DNS and NOWNS results for all for modes.

Figure 3: Amplitude of $u^{\prime}$ v.s. streamwise coordinate, $\mathrm{Re}_{x}$, for the small-amplitude oblique-wave breakdown case at frequency $F=86\times 10^{-6}$, spanwise wavenumber $b=2/9\times 10^{-3}$, with initial amplitude of $u_{\max}^{\prime(1,1)}(x_{0})=\sqrt{2}\times 10^{-3}$.

(a) MFD

(b) Oblique mode

(c) Vortex mode

(d) First harmonic in $z$

Figure 4: $u$-velocity profiles at streamwise coordinate $\mathrm{Re}_{x}=4.69\times 10^{5}$ for small-amplitude oblique-wave breakdown.

4 Advantages of the NOWNS procedure

In this section, we apply NOWNS in three scenarios where NPSE is known to break down.

4.1 High amplitude oblique-wave breakdown

It is well-known that NPSE can fail for sufficiently strong nonlinearities [11, 12], as demonstrated in Joslin et al. (1993) for the oblique-wave breakdown case [29]. We modify the setup of section 3.2 by increasing the amplitude of the oblique wave at the inlet by a factor of 10 to $u_{\max}^{\prime(1,1)}(x_{0})=\sqrt{(}2)\times 10^{-2}$, while we change the streamwise domain so that it comprises 2000 grid points over the domain $\mathrm{Re}_{x}\in[2.73529\times 10^{5},4.9\times 10^{5}]$, and the Fourier series are truncated at $M=N=7$.

We plot the amplitudes of the $u$-velocity as a function of streamwise station in figure 5. First we note that we have excellent agreement between the amplitudes predicted by NOWNS, NPSE, and DNS for the early stages of transition. However, NPSE begins to fail towards the end of the domain, while NOWNS is able to march all the way to the end of the DNS calculation. The DNS calculation is under-resolved near the end of the domain [29], which may explain the discrepancy between the DNS and NOWNS calculations. We further note that although the oblique-wave, $(1,1)$, is initially the dominant instability, it is rapidly overtaken by the vortex mode, $(0,2)$, leading to streaks, which we observe in the contour plots of the $u$-velocity ($u=\bar{u}+u^{\prime}$), shown figure 6. The maximum amplitude of the oblique-wave occurs at $y=2.55$ (figure 6(a)), while the maximum amplitude of the vortex mode occurs at $y=1.31$ (figure 6(b)).

Figure 5: Amplitude of $u^{\prime}$ v.s. streamwise coordinate, $\mathrm{Re}_{x}$, for high amplitude oblique-wave breakdown.

(a) $y=2.55$

(b) $y=1.31$

Figure 6: Contours of the instantaneous $u$-velocity for high amplitude oblique-wave breakdown.

4.2 Low amplitude oblique-wave breakdown with random noise

Here we demonstrate for the low amplitude oblique wave breakdown case of section 3.2 that NOWNS is robust to random noise applied to the inlet boundary condition. Given the eigenfunction from the locally parallel linear stability theory, $\bm{q}_{\mathrm{LST}}$, we add random noise, $\bm{q}_{\mathrm{noise}}$, to obtain the inlet condition $\bm{q}_{\mathrm{LST}}+\varepsilon\tilde{\bm{q}}_{\mathrm{noise}}$. We choose complex random noise such that

$$\tilde{\bm{q}}_{\mathrm{noise}}=\tilde{\bm{q}}_{\mathrm{noise},r}+i\tilde{\bm{q}}_{\mathrm{noise},i},\quad\tilde{\bm{q}}_{\mathrm{noise},r},\tilde{\bm{q}}_{\mathrm{noise},i}\sim\mathcal{U}_{[0,1]},$$

where $\mathcal{U}_{[a,b]}$ represents the uniform distribution over the interval $[a,b]$. We then normalize the noise to obtain $q_{\mathrm{noise}}$, such that the maximum amplitude of the $u$-velocity noise is equal to the free-stream $u$-velocity, $U_{\infty}$. We recall that the amplitude of $\bm{q}_{\mathrm{LST}}$ is $\sqrt{2}\times 10^{-3}$, and we choose $\varepsilon=\sqrt{2}\times 10^{-5}$. In figure 7, we plot the profile of the $u$- and $v$-velocities, as well as the thermodynamic variables. We see that the random noise has a relatively small effect on the $u$-velocity, and a slightly more pronounced effect on the $v$-velocity, while it has a larger impact on the thermodynamic variables. Like the $u$-velocity profile, the random noise has a relatively small effect on the $w$-velocity profile, so we omit this plot.

The NOWNS march succeeds even if we introduce large disturbances to the velocity fields, but the $u^{\prime}$ amplitudes differ substantially from the noise-less case due to the large perturbations. Therefore, we instead introduce relatively large random disturbances to the thermodynamic variables, and relatively small ones to the velocity field. In figure 8, we see that the amplitudes predicted by NOWNS for the noisy inlet condition agree closely with those without noise. We also see that although NPSE initially is able to accurately predict the evolution of the $u$-velocity amplitudes, it eventually becomes inaccurate (especially for $u_{(0,0)}^{\prime}$ and $u_{(1,1)}^{\prime}$) before failing. We also plot the contours of the real part of the $u$- and $v$-velocities of the oblique wave, with and without noise in figure 9. We see that despite the noisy inlet condition, NOWNS evolves the $u$-velocity of the oblique wave such that it matches closely the case without noise. On the other hand, the $v$-velocity is more affected by the numerical noise, yet we still obtain good qualitative agreement. The other modes (e.g., the vortex mode) are evolved accurately by the NOWNS calculation with noisy inlet condition, and the contour plots with and without the numerical noise are indistinguishable from each other, and so are not plotted here.

Figure 7: Oblique-wave inlet boundary condition with and without random noise.

Figure 8: Amplitude of $u^{\prime}$ v.s. streamwise coordinate, $\mathrm{Re}_{x}$, for low amplitude oblique-wave breakdown with and without random noise.

Figure 9: Contour plots of the $u$- and $v$ velocities of the oblique wave, with and without random noise.

4.3 Blowing/suction strip

Blowing/suction strips are frequently used to study laminar-turbulent boundary layer transition in low-speed boundary layer flows [4, 5, 6, 9, 34]. Here, we introduce disturbances by specifying a non-zero wall-normal velocity such that $v(y=0)=f(x,z,t)$, for some function $f(x,z,t)$ that is periodic in $t$ and $z$. NPSE does not support blowing/suction strips because they introduce non-modal disturbances, demonstrating an advantage of NOWNS. However, we note that Herbert’s second stability theory yields an inlet boundary condition so that NPSE (and NOWNS) can be used to study fundamental (K-type) and subharmonic (H-type) transition [33, 35].

Rist and Fasel (1995) used DNS with a blowing/suction strip to study K-type transition [5], while similar studies were performed by Sayadi et al. (2013) using DNS [6], and by Rigas et al. (2021) using a harmonic balance method (HBM) [34]. The blowing/suction strip is given by

$$f(x,z,t)=5\times 10^{-3}\sin(\omega t)v_{a}(x)+1.3\times 10^{-4}\cos(\beta z)v_{s}(x),$$ (24)

where

$$v_{a}(\mathrm{Re}_{x})=\left\{\begin{array}[]{cc}0,&\mathrm{Re}_{x}\leq\mathrm{Re}_{x}(x_{1})\\ 15.1875\xi^{5}-35.4375\xi^{4}+20.25\xi^{3},&\mathrm{Re}_{x}(x_{1})<\mathrm{Re}_{x}\leq\mathrm{Re}_{x}(x_{m})\\ -v_{a}\big{(}2\mathrm{Re}_{x}(x_{m})-\mathrm{Re}_{x}\big{)},&\mathrm{Re}_{x}(x_{m})<\mathrm{Re}_{x}\leq\mathrm{Re}_{x}(x_{2})\\ 0,&\mathrm{Re}_{x}(x_{2})<\mathrm{Re}_{x}\end{array}\right.$$ (25a)

$$v_{s}(\mathrm{Re}_{x})=\left\{\begin{array}[]{cc}0,&R\leq\mathrm{Re}_{x}(x_{1})\\ -3\xi^{4}+4\xi^{3},&\mathrm{Re}_{x}(x_{1})<\mathrm{Re}_{x}\leq\mathrm{Re}_{x}(x_{m})\\ v_{s}\big{(}2\mathrm{Re}_{x}(x_{m})-\mathrm{Re}_{x}\big{)},&\mathrm{Re}_{x}(x_{m})<\mathrm{Re}_{x}\leq\mathrm{Re}_{x}(x_{2})\\ 0,&\mathrm{Re}_{x}(x_{2})<\mathrm{Re}_{x}\end{array}\right.$$ (25b)

for $\mathrm{Re}_{x}(x_{1})=1.3438\times 10^{5}$, $\mathrm{Re}_{x}(x_{2})=1.5532\times 10^{5}$, $x_{m}=(x_{1}+x_{2})/2$, and $\xi=(\mathrm{Re}_{x}-\mathrm{Re}_{x}(x_{1}))/(\mathrm{Re}_{x}(x_{m})-\mathrm{Re}_{x}(x_{1}))$. We choose $F=110\times 10^{-6}$ and $b=0.423\times 10^{-3}$ with $M=N=4$, while the grid extends over $\mathrm{Re}_{x}\in[1.33956\times 10^{5},2.72\times 10^{5}]$ and $y\in[0,60]$, with 1300 stations in $x$ and 100 grid points in $y$.

Figure 10 shows excellent agreement between the $u^{\prime}$ amplitudes of NOWNS and the DNS of Rist and Fasel (1995) [5]. We note that in the early stages of the march there is disagreement between the DNS and NOWNS calculations because the blowing/suction strip causes upstream effects that NOWNS neglects by construction. However, these disturbances are convective in nature, and the amplitudes predicted by NOWNS rapidly converge to those predicted by DNS as the march progresses downstream.

Figure 10: Amplitude of $u^{\prime}$ v.s. streamwise coordinate, $\mathrm{Re}_{x}$, for K-type transition

5 Conclusion

We have extended the OWNS approach, a fast marching procedure previously developed for solving linear flow disturbance equations, to support nonlinear interactions. We have demonstrated that it is effective for 2D and 3D disturbances in a wall-bounded flow by comparing against the spatial DNS and NPSE studies of Joslin et al. (1993) [29]. Moreover, we have demonstrated that NOWNS is successful in cases where NPSE fails. In particular, we have shown that NOWNS can handle stronger nonlinearities and continues marching past the point where NPSE fails, and that NOWNS is robust to numerical noise. We further demonstrated that NOWNS supports non-modal instabilities in the form of a blowing/suction strip, which would not be supported by NPSE. Future work will extend the NOWNS procedure to support high-speed flows and will concentrate on studying optimal transition mechanisms, seeking to apply NOWNS in cases where spatial DNS and other global methods are too computationally expensive to be feasible.

Appendix A Navier-Stokes Equations

Given the non-dimensional Navier-Stokes equations 1, we define the vector $\bm{q}=(\nu,u,v,w,p)$ and write

$$\begin{split}\frac{\partial\bm{q}}{\partial t}&+[A_{x}(\bm{q})+B_{x}(\bm{q})]\frac{\partial\bm{q}}{\partial x}+[A_{y}(\bm{q})+B_{y}(\bm{q})]\frac{\partial\bm{q}}{\partial y}+[A_{z}(\bm{q})+B_{z}(\bm{q})]\frac{\partial\bm{q}}{\partial z}\\ &+B_{xx}(\bm{q})\frac{\partial^{2}\bm{q}}{\partial x^{2}}+B_{yy}(\bm{q})\frac{\partial^{2}\bm{q}}{\partial y^{2}}+B_{zz}(\bm{q})\frac{\partial^{2}\bm{q}}{\partial z^{2}}+B_{xy}(\bm{q})\frac{\partial^{2}\bm{q}}{\partial x\partial y}+B_{xz}(\bm{q})\frac{\partial^{2}\bm{q}}{\partial x\partial z}+B_{yz}(\bm{q})\frac{\partial^{2}\bm{q}}{\partial y\partial z}=0.\end{split}$$ (26)

where $B$ denotes viscous terms while $A$ denotes inviscid terms. Next we decompose the flow into a time-invariant equilibrium solution, $\bar{\bm{q}}$, and a time-varying disturbance variable, $\bm{q}^{\prime}$, such that $\bm{q}=\bar{\bm{q}}+\bm{q}^{\prime}$, which we use to obtain

$$\begin{split}\frac{\partial\bm{q}^{\prime}}{\partial t}&+[A_{x}(\bar{\bm{q}})+B_{x}(\bar{\bm{q}})]\frac{\partial\bm{q}^{\prime}}{\partial x}+[A_{y}(\bar{\bm{q}})+B_{y}(\bar{\bm{q}})]\frac{\partial\bm{q}^{\prime}}{\partial y}+[A_{z}(\bar{\bm{q}})+B_{z}(\bar{\bm{q}})]\frac{\partial\bm{q}^{\prime}}{\partial z}+C(\bar{\bm{q}})\bm{q}^{\prime}\\ &+B_{xx}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x^{2}}+B_{yy}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial y^{2}}+B_{zz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial z^{2}}+B_{xy}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial y}+B_{xz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial z}+B_{yz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial y\partial z}=\bm{F}(\bm{q}^{\prime}),\end{split}$$ (27)

where we have defined $C(\bar{\bm{q}})$ such that

$$\begin{split}C(\bar{\bm{q}})\bm{q}^{\prime}&=[A_{x}(\bm{q}^{\prime})+B_{x}(\bm{q}^{\prime})]\frac{\partial\bar{\bm{q}}}{\partial x}+[A_{y}(\bm{q}^{\prime})+B_{y}(\bm{q}^{\prime})]\frac{\partial\bar{\bm{q}}}{\partial y}+[A_{z}(\bm{q}^{\prime})+B_{z}(\bm{q}^{\prime})]\frac{\partial\bar{\bm{q}}}{\partial z}\\ &+B_{xx}(\bm{q}^{\prime})\frac{\partial^{2}\bar{\bm{q}}}{\partial x^{2}}+B_{yy}(\bm{q}^{\prime})\frac{\partial^{2}\bar{\bm{q}}}{\partial y^{2}}+B_{zz}(\bm{q}^{\prime})\frac{\partial^{2}\bar{\bm{q}}}{\partial z^{2}}+B_{xy}(\bm{q}^{\prime})\frac{\partial^{2}\bar{\bm{q}}}{\partial x\partial y}+B_{xz}(\bm{q}^{\prime})\frac{\partial^{2}\bar{\bm{q}}}{\partial x\partial z}+B_{yz}(\bm{q}^{\prime})\frac{\partial^{2}\bar{\bm{q}}}{\partial y\partial z},\end{split}$$ (28)

and the nonlinear term

$$\begin{split}\bm{F}(\bm{q}^{\prime})&=-[A_{x}(\bm{q}^{\prime})+B_{x}(\bm{q}^{\prime})]\frac{\partial\bm{q}^{\prime}}{\partial x}-[A_{y}(\bm{q}^{\prime})+B_{y}(\bm{q}^{\prime})]\frac{\partial\bm{q}^{\prime}}{\partial y}-[A_{z}(\bm{q}^{\prime})+B_{z}(\bm{q}^{\prime})]\frac{\partial\bm{q}^{\prime}}{\partial z}\\ &-B_{xx}(\bm{q}^{\prime})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x^{2}}-B_{yy}(\bm{q}^{\prime})\frac{\partial^{2}\bm{q}^{\prime}}{\partial y^{2}}-B_{yy}(\bm{q}^{\prime})\frac{\partial^{2}\bm{q}^{\prime}}{\partial z^{2}}-B_{xy}(\bm{q}^{\prime})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial y}-B_{xz}(\bm{q}^{\prime})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial z}-B_{yz}(\bm{q}^{\prime})\frac{\partial^{2}\bm{q}^{\prime}}{\partial y\partial z}.\end{split}$$ (29)

We further define the linear operator

$$\begin{split}\mathcal{L}(\bar{\bm{q}})=&-\frac{\partial}{\partial t}-[A_{y}(\bar{\bm{q}})+B_{y}(\bar{\bm{q}})]\frac{\partial\bm{q}^{\prime}}{\partial y}-[A_{z}(\bar{\bm{q}})+B_{z}(\bar{\bm{q}})]\frac{\partial\bm{q}^{\prime}}{\partial z}-C(\bar{\bm{q}})\bm{q}^{\prime}\\ &-B_{yy}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial y^{2}}-B_{zz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial z^{2}}-B_{yz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial y\partial z}=0,\end{split}$$ (30)

yielding

$$A_{x}(\bar{\bm{q}})\frac{\partial\bm{q}^{\prime}}{\partial x}=\mathcal{L}(\bar{\bm{q}})\bm{q}^{\prime}+\bm{F}(\bm{q}^{\prime})-B_{x}(\bar{\bm{q}})\frac{\partial\bm{q}^{\prime}}{\partial x}-B_{xx}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x^{2}}-B_{xy}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial y}-B_{xz}(\bar{\bm{q}})\frac{\partial^{2}\bm{q}^{\prime}}{\partial x\partial z}.$$ (31)

The operators for the first derivatives (without viscous terms), $A_{x}$, $A_{y}$, and $A_{z}$, are given by

$$A_{x}(\bm{q})=\left[\begin{array}[]{ccccc}u&-\nu&0&0&0\\ 0&u&0&0&\nu\\ 0&0&u&0&0\\ 0&0&0&u&0\\ 0&\gamma p&0&0&u\end{array}\right],\quad A_{y}(\bm{q})=\left[\begin{array}[]{ccccc}v&0&-\nu&0&0\\ 0&v&0&0&0\\ 0&0&v&0&\nu\\ 0&0&0&v&0\\ 0&0&\gamma p&0&v\end{array}\right],\quad A_{z}(\bm{q})=\left[\begin{array}[]{ccccc}w&0&0&-\nu&0\\ 0&w&0&0&0\\ 0&0&w&0&0\\ 0&0&0&w&\nu\\ 0&0&0&\gamma p&w\end{array}\right],$$

while the operators for the first derivatives (with viscous terms) $B_{x}$, $B_{y}$, and $B_{z}$, are given by

$$B_{x}(\bm{q})=-\frac{2\gamma}{PrRe}\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ \frac{\partial p}{\partial x}&0&0&0&\frac{\partial\nu}{\partial x}\end{array}\right],\quad B_{y}(\bm{q})=-\frac{2\gamma}{PrRe}\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ \frac{\partial p}{\partial y}&0&0&0&\frac{\partial\nu}{\partial y}\end{array}\right],$$

and

$$B_{z}(\bm{q})=-\frac{2\gamma}{PrRe}\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ \frac{\partial p}{\partial z}&0&0&0&\frac{\partial\nu}{\partial z}\end{array}\right],$$

For the second derivatives $B_{xx}$ and $B_{yy}$ we have

$$B_{xx}(\bm{q})=\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&-\frac{4}{3}\frac{\nu}{Re}&0&0&0\\ 0&0&-\frac{\nu}{Re}&0&0\\ 0&0&0&-\frac{\nu}{Re}&0\\ -\frac{\gamma p}{RePr}&0&0&0&-\frac{\gamma\nu}{RePr}\end{array}\right],\quad B_{yy}(\bm{q})=\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&-\frac{\nu}{Re}&0&0&0\\ 0&0&-\frac{4}{3}\frac{\nu}{Re}&0&0\\ 0&0&0&-\frac{\nu}{Re}&0\\ -\frac{\gamma p}{RePr}&0&0&0&-\frac{\gamma\nu}{RePr}\end{array}\right],$$

while for $B_{zz}$ and $B_{xy}$ we have

$$B_{zz}(\bm{q})=\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&-\frac{\nu}{Re}&0&0&0\\ 0&0&-\frac{\nu}{Re}&0&0\\ 0&0&0&-\frac{4}{3}\frac{\nu}{Re}&0\\ -\frac{\gamma p}{RePr}&0&0&0&-\frac{\gamma\nu}{RePr}\end{array}\right],\quad B_{xy}(\bm{q})=\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&0&-\frac{\nu}{3Re}&0&0\\ 0&-\frac{\nu}{3Re}&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\end{array}\right],$$

and for $B_{xz}$ and $B_{yz}$ we have

$$B_{xz}(\bm{q})=\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&0&0&-\frac{\nu}{3Re}&0\\ 0&0&0&0&0\\ 0&-\frac{\nu}{3Re}&0&0&0\\ 0&0&0&0&0\end{array}\right],\quad B_{yz}(\bm{q})=\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&-\frac{\nu}{3Re}&0\\ 0&0&-\frac{\nu}{3Re}&0&0\\ 0&0&0&0&0\end{array}\right].$$

Finally, the operator $C$ is given by

$$C(\bm{q})=\left[\begin{array}[]{ccccc}-\nabla\cdot\bm{u}&\frac{\partial\nu}{\partial x}&\frac{\partial\nu}{\partial y}&\frac{\partial\nu}{\partial z}&0\\ \frac{\partial p}{\partial x}-\frac{1}{Re}\nabla^{2}u-\frac{1}{3Re}[\partial_{xx}u+\partial_{xy}v+\partial_{xz}w]&\frac{\partial u}{\partial x}&\frac{\partial u}{\partial y}&\frac{\partial u}{\partial z}&0\\ \frac{\partial p}{\partial y}-\frac{1}{Re}\nabla^{2}v-\frac{1}{3Re}[\partial_{xy}u+\partial_{yy}v+\partial_{yz}w]&\frac{\partial v}{\partial x}&\frac{\partial v}{\partial y}&\frac{\partial v}{\partial z}&0\\ \frac{\partial p}{\partial z}-\frac{1}{Re}\nabla^{2}w-\frac{1}{3Re}[\partial_{xz}u+\partial_{yz}v+\partial_{zz}w]&\frac{\partial w}{\partial x}&\frac{\partial w}{\partial y}&\frac{\partial w}{\partial z}&0\\ -\frac{\gamma}{RePr}\nabla^{2}p&\frac{\partial p}{\partial x}&\frac{\partial p}{\partial y}&\frac{\partial p}{\partial z}&\gamma\nabla\cdot\bm{u}-\frac{\gamma}{RePr}\nabla^{2}\nu\end{array}\right],$$

while the nonlinear term is a vector comprising the following components

$$F_{1}(\bm{q}) =u\frac{\partial\nu}{\partial x}+v\frac{\partial\nu}{\partial y}+w\frac{\partial\nu}{\partial z}-\nu\frac{\partial u}{\partial x}-\nu\frac{\partial v}{\partial y}-\nu\frac{\partial w}{\partial z},$$

$$F_{2}(\bm{q}) =\nu\frac{\partial p}{\partial x}-u\frac{\partial u}{\partial x}-v\frac{\partial u}{\partial y}-w\frac{\partial u}{\partial z}+\frac{1}{Re}\nu\nabla^{2}u+\frac{1}{3Re}\Big{(}\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}v}{\partial x\partial y}+\frac{\partial^{2}w}{\partial x\partial z}\Big{)},$$

$$F_{3}(\bm{q}) =\nu\frac{\partial p}{\partial y}-u\frac{\partial v}{\partial x}-v\frac{\partial v}{\partial y}-w\frac{\partial v}{\partial z}+\frac{1}{Re}\nu\nabla^{2}u+\frac{1}{3Re}\Big{(}\frac{\partial^{2}u}{\partial x\partial y}+\frac{\partial^{2}v}{\partial y^{2}}+\frac{\partial^{2}w}{\partial y\partial z}\Big{)},$$

$$F_{4}(\bm{q}) =\nu\frac{\partial p}{\partial z}-u\frac{\partial w}{\partial x}-v\frac{\partial w}{\partial y}-w\frac{\partial w}{\partial z}+\frac{1}{Re}\nu\nabla^{2}u+\frac{1}{3Re}\Big{(}\frac{\partial^{2}u}{\partial x\partial z}+\frac{\partial^{2}v}{\partial y\partial z}+\frac{\partial^{2}w}{\partial z^{2}}\Big{)},$$

$$F_{5}(\bm{q}) =\frac{\gamma}{RePr}[\nu\nabla^{2}p+2\nabla p\cdot\nabla\nu+p\nabla^{2}\nu]-\bm{u}\cdot\nabla p-\gamma p\nabla\cdot{\bm{u}}.$$

Appendix B Jacobian of the NOWNS equations

Here we derive the Jacobian of the fully-discrete system of NOWNS equations (17) for 2D disturbances. However, this analysis readily extends to 3D disturbances as well. We first introduce the NOWNS residual

$$\hat{\bm{r}}^{\ddagger(k+1)}_{m}=\sum_{l=0}^{s-1}c^{(l)}A^{\ddagger}\hat{\bm{\phi}}_{m}^{\ddagger(k+1-l)}-\hat{L}_{m}^{\ddagger}\hat{\bm{\phi}}_{m}^{\ddagger(k+1)}-\hat{\bm{F}}_{m}^{\ddagger(k+1)}-\hat{\bm{f}}_{m}^{\ddagger(k+1)}=0,\quad m=0,\dots,2M,$$ (32)

and take it’s derivative with respect to $\hat{\bm{\phi}}_{p}^{\ddagger(k+1)}$ to obtain

$$\frac{\partial\hat{\bm{r}}^{\ddagger(k+1)}_{m}}{\partial\hat{\bm{\phi}}_{p}^{\ddagger(k+1)}}=(c^{(0)}A^{\ddagger}-\hat{L}_{m}^{\ddagger})\delta_{mp}-\frac{\partial\hat{\bm{F}}_{m}^{\ddagger(k+1)}}{\partial\hat{\bm{\phi}}_{p}^{\ddagger(k+1)}},\quad m,p=0,\dots,2M,$$ (33)

Using the definition DFT, it can be shown (for 2D disturbances) that

$$\frac{\partial\hat{\bm{F}}^{\ddagger(k+1)}_{m}}{\partial\hat{\bm{\phi}}^{\ddagger(k+1)}_{p}}=\frac{1}{2M+1}\sum_{l=0}^{2M}\frac{\partial\tilde{\bm{F}}^{\ddagger(k+1)}}{\partial\bm{\phi}^{\ddagger(k+1)}}\Big{|}_{\bm{\phi}=\bm{\phi}_{l}}e^{i2\pi l\frac{p-m}{2M+1}},$$ (34)

which allows us to write

$$\hat{J}_{p-m}^{\ddagger}=-\frac{\partial\hat{\bm{F}}_{m}^{\ddagger(k+1)}}{\partial\hat{\bm{\phi}}_{p}^{\ddagger(k+1)}}=-\begin{bmatrix}0&0&0\\ P_{1,m}\begin{bmatrix}\tilde{A}_{\pm\pm}^{-1}\hat{J}_{\pm\pm,p-m}\\ \hat{J}_{0\pm,p-m}\end{bmatrix}&P_{1,m}\begin{bmatrix}\tilde{A}_{\pm\pm}^{-1}\hat{J}_{\pm 0,p-m}\\ \hat{J}_{00,p-m}\end{bmatrix}&0\\ \hat{J}_{0\pm,p-m}&\hat{J}_{00,p-m}&0\end{bmatrix},\quad m,p=0,\dots,2M.$$ (35)

We drop the $(k+1)$ superscript for simplicity, and for $M=1$ we can write

$$\begin{bmatrix}(c^{(0)}A^{\ddagger}-\hat{L}_{0}^{\ddagger})+\hat{J}^{\ddagger}_{0}&\overline{\hat{J}^{\ddagger}_{1}}&\hat{J}^{\ddagger}_{1}\\ \hat{J}^{\ddagger}_{1}&(c^{(0)}A^{\ddagger}-\hat{L}^{\ddagger}_{1})+\hat{J}^{\ddagger}_{0}&\overline{\hat{J}^{\ddagger}_{1}}\\ \overline{\hat{J}^{\ddagger}_{1}}&\hat{J}^{\ddagger}_{1}&(c^{(0)}A^{\ddagger}-\hat{L}_{-1}^{\ddagger})+\hat{J}^{\ddagger}_{0}\end{bmatrix}\begin{bmatrix}\Delta\hat{\bm{\phi}}_{0}^{\ddagger}\\ \Delta\hat{\bm{\phi}}_{1}^{\ddagger}\\ \overline{\Delta\hat{\bm{\phi}}_{1}}^{\ddagger}\end{bmatrix}=\begin{bmatrix}\hat{r}_{0}^{\ddagger}\\ \hat{r}_{1}^{\ddagger}\\ \overline{\hat{r}_{1}}^{\ddagger}\end{bmatrix}.$$ (36)

where we have used $\hat{J}^{\ddagger}_{m}=\hat{J}^{\ddagger}_{m+2M+1}$ and $\hat{J}^{\ddagger}_{m}=\overline{\hat{J}^{\ddagger}_{-m}}$ for $m=0,\dots,M$.

Appendix C Streamwise diffusion terms

Here we demonstrate for the 2D validation case discussed in section 3.1 that although streamwise diffusion effects have a minimal impact on the linear calculation (figure 11(a)), the impact is more pronounced in the nonlinear case (figure 11(b)). We see that prior to the second neutral stability point, the calculations with and without the streamwise diffusion terms yield similar amplitudes, but that downstream of this point, the amplitudes of the calculation including the streamwise viscous terms has much lower amplitudes (particularly for the higher harmonics in the nonlinear case).

(a) Linear calculation

(b) Nonlinear calculation

Figure 11: $u$-velocity amplitudes with and without streamwise diffusion terms for the 2D evolution of the TS wave.

Appendix D Inclusion of the streamwise pressure gradient

Unlike the NPSE approach, the NOWNS approach has no minimum step size requirement for the march to be numerically stable, but it is still necessary to neglect the streamwise pressure gradient for the zero-frequency modes. If the streamwise pressure gradient for these modes is retained, the NOWNS march remains stable but becomes inaccurate, as demonstrated in figure 12. We plot the $u$-velocity amplitudes with and without $\partial_{x}p_{0}$ in figure 12(a), where we see that we have reasonable agreement for $m\neq 0$, but disagreement for $m=0$. In figure 12(b) we plot the profile of the $v_{0}$, and we notice that the profile predicted by NOWNS when $\partial_{x}p_{0}$ is included is substantially different from the profiles predicted by NPSE and DNS.

If we neglect $\partial_{x}p_{0}$, but include the streamwise diffusion terms $\partial_{xx}q_{0}$, then we must project the MFD and we find that we have good agreement between the DNS and NOWNS calculations. Therefore, we can conclude that the recursion parameters we are using for the zero-frequency modes are valid. However, as discussed above, the calculation remains stable but becomes inaccurate when we include $\partial_{x}p_{0}$. Although our march remains stable, there including $\partial_{x}p_{0}$ leads to inaccuracies in the march for unknown reasons.

In figure 13(a), we plot the $u$-velocity amplitudes computed by NOWNS with and without the pressure gradient for the zero-frequency modes. Mode (1,0) is tracked reasonably accurately, but the other modes are not. In particularly, the (1,2), (1,1), and (0,1) modes and the MFD have higher amplitudes than they should. In figure 13(b), we see that we have better qualitative agreement in the early stages of the march between the NOWNNS and DNS calculations when the streamwise pressure gradient terms are include for the zero-frequency modes. However, we have worse quantitative agreement in the later stages of the march, which in turn causes the larger amplitudes of the (1,2), (1,1) and (0,1) modes observed in figure 13(a).

(a) $u$-velocity amplitudes

(b) $v$-velocity profile

Figure 12: 2D validation case with and without streamwise pressure gradient for zero-frequency modes.

(a) $u$-velocity amplitudes

(b) $u$-velocity amplitudes for MFD

Figure 13: K-type transition case with and without streamwise pressure gradient for zero-frequency modes.

Appendix E Comparison of nonlinear solution procedures

We have discussed three procedures for solving the nonlinear system of equations: (i) Newton’s method, (ii) a quasi-Newton method that includes part of the nonlinear Jacobian, and (iii) a quasi-Newton method that excludes the Jacobian of the nonlinear term. We plot the iterations to converge as a function of streamwise station in figure 14. For the 2D validation case, we compare the iterations to convergence for the quasi-Newton method and Newton’s method in figure 14(a), which shows that Newton’s method converges in fewer iterations. In figure 14(b), we make the same comparison for linear and nonlinear quasi-Newton methods for K-type transition. We see that farther downstream, where the nonlinearity is stronger, the nonlinear quasi-Newton method converges in fewer iterations than the linear version. However, for sufficiently strong nonlinearities, it is necessary to use Newton’s method.

(a) Comparison of Newton’s method and the quasi-Newton method for the 2D validation case

(b) Comparison of the linear and nonlinear quasi-Newton methods for K-type transition

Figure 14: Comparison of solution procedures for the nonlinear system of equations

Appendix F Linearization about the baseflow vs. the corrected mean flow

Whereas in linear stability analysis, the mean flow is determined only by the baseflow, $\bar{\bm{q}}$, in the nonlinear case the disturbances interact to excite the MFD, $\bm{q}_{00}^{\prime}$, which yields the corrected mean flow, $\bar{\bm{q}}+\bm{q}_{00}^{\prime}$. In linear OWNS, we linearize the projection operators about the baseflow, while in nonlinear OWNS, we can choose to instead linearize about the corrected mean flow. For the K-type transition case of section 4.3, we perform the NOWNS calculation again, but we instead linearize about the corrected mean flow. In figure 15, we compare the amplitudes computed by the two approaches to NOWNS and we see that in the early stages of the march, when the amplitudes are small, the two calculations are nearly identical, but that they begin to differ slightly as the disturbance amplitudes increase. The differences are small, but are most notable for the MFD, $\bm{q}_{00}^{\prime}$, and the vortex mode, $q_{01}^{\prime}$, after $\mathrm{Re}_{x}=2.5\times 10^{5}$, which corresponds to an MFD amplitude of roughly 3% of $U_{\infty}$. The two NOWNS marches continue to yield similar results until the end of the calculation at $\mathrm{Re}_{x}=2.74\times 10^{5}$, where the MFD amplitude is 11% of $U_{\infty}$.

Linearizing about the baseflow is more computationally efficient because the projection operators do not change between iterations. Moreover, the choice to linearize about the corrected mean flow instead of the baseflow does not appear to substantially affect the stability calculations, for MFD amplitudes less than roughly 10% of $U_{\infty}$. Therefore, we recommend that NOWNS calculations be linearized about the baseflow rather than the corrected mean flow.

Figure 15: Amplitude of $u^{\prime}$ v.s. streamwise coordinate, $\mathrm{Re}_{x}$, for fundamental breakdown. Compare linearizing NOWNS about $\bar{\bm{q}}$ vs. $\bar{\bm{q}}+\bm{q}^{\prime}_{00}$.

Appendix G Recursion parameter sets

The recursion parameter sets for non-zero frequency modes match those used in [18]. They depend on $\omega$ through the streamwise wave number $k=\omega/\bar{c}$, and when $k=0$, the recursion parameters associated with the vortical and prpagating acoustic waves go to zero, leaving only the evanscent acoustic waves.