New version program summary
Program Title: Delft Harmonic Navier–Stokes Solver (DeHNSSo)
CPC Library link to program files:
Developer’s repository link: https://github.com/pcboldini/DeHNSSo
Licensing provisions: GPLv3
Programming language: MATLAB, C
Journal Reference of previous version: Westerbeek S., Hulshoff S., Schuttelaars H., Kotsonis M., DeHNSSo: The Delft Harmonic Navier–Stokes Solver for nonlinear stability problems with complex geometric features, Computer Physics Communications 302 (2024) 109250.
Does the new version supersede the previous version?: Yes
Reasons for the new version: Release of new features, extension to compressible regimes, and performance improvements.
Summary of revisions:
general curvilinear perturbation formulation and automatic body-fitted stability-grid generation;
nonuniform streamwise grids with fourth-order finite differences;
compressible primitive-variable formulation with thermal and far-field boundary conditions, freestream forcing, and parabolized outflow treatment;
linear-cost sparse assembly and reduced operator storage;
reusable per-mode factorizations with sparse LU, banded LU, and out-of-core block-Thomas back-ends;
Newton pseudo-arclength continuation in the inflow disturbance amplitude.
Nature of problem: DeHNSSo-v2 computes the linear and nonlinear spatial evolution of temporally and spanwise periodic disturbances in incompressible and compressible laminar boundary layers over general two-dimensional body-fitted geometries. It is designed for studies of receptivity, nonlinear modal interactions, transition, and interactions with localized surface features or flow-control devices.
Solution method: Perturbations to a steady base flow are expanded in temporal and spanwise Fourier modes and discretized using streamwise finite differences and wall-normal finite differences or Chebyshev collocation. The resulting sparse systems, which couple the whole two-dimensional domain, are solved by direct factorization, with nonlinear coupling treated by fixed-point iteration or Newton pseudo-arclength continuation.
E-mail address: pietro.c.boldini@gmail.com (P. C. Boldini)
Keywords: Boundary-layer transition, Nonlinear stability analysis, Compressible flows, Harmonic Navier–Stokes, Surface imperfections, Open source
Introduction
Reducing drag on transport aircraft and limiting the aerodynamic heating of high-speed vehicles are central objectives of aerodynamic design. Both depend strongly on whether the boundary layer over the vehicle surface remains laminar or becomes turbulent: a turbulent boundary layer produces substantially greater skin friction and wall heat flux than a laminar one. Predicting the transition location therefore underlies two distinct design problems. On subsonic transport aircraft, skin friction accounts for a large fraction of the total drag, so delaying transition is among the most effective aerodynamic means of reducing drag [1 , 2 ] . At high speed, transition governs the heat loads that size the thermal-protection system and is often driven by the first and second Mack modes [3 ] . Both regimes ultimately pose the same problem: accurately simulating how small disturbances are amplified within a laminar boundary layer and how that amplification is modified by wall geometry, thermal conditions, and Mach number, among other factors.
The subsonic swept wing illustrates the sensitivity of this process to surface geometry. At the low freestream-turbulence levels typical of cruise flight, transition is commonly driven by stationary crossflow instability [4 ] . Laminar-flow-control strategies targeting this instability can be compromised by the geometric imperfections carried by real wings, including steps and gaps at panel junctions and the bumps and waviness produced by manufacturing and structural loading [5 , 6 ] . Improvements in manufacturing have reduced distributed surface waviness, increasing the relative importance of localized steps and gaps, which modify the incoming instability and may displace the transition front upstream [7 , 8 , 9 ] . Conversely, a deliberately shaped smooth hump has recently been shown to delay swept-wing transition by approximately 14% of the chord when encountered by a low-amplitude crossflow disturbance, but to advance transition when the incoming amplitude is larger [10 , 11 ] . The effect of the same surface feature can thus change sign according to the amplitude of the incoming disturbance and hence the nonlinear state of the perturbation field. Predicting either response requires resolving the nonlinear interaction between the instability and the surface feature at the scale of the boundary layer. Closely related interactions underlie wall-based control devices, both active, such as suction and plasma actuators, and passive, such as Helmholtz resonators [12 ] and other surface metamaterials [13 ] . Through localized forcing, these devices can damp or reinforce selected instability components.
Several numerical frameworks are available for predicting disturbance evolution, as reviewed by Westerbeek et al. [14 ] . Local stability theory treats the flow at successive streamwise stations and does not retain the global streamwise coupling of the disturbance field. The parabolized stability equations (PSE) [15 , 16 ] , including their extensions to compressible boundary layers and supersonic oblique breakdown [17 , 18 , 19 ] , retain non-parallel development while neglecting the second-order streamwise derivatives of a slowly varying disturbance shape function. The regularization required to obtain a stable downstream march accurately retains a single slowly varying spatial component at each temporal frequency, but may distort other instability mechanisms, nonmodal responses, and acoustic components [20 ] . The slowly varying approximation itself may lose validity where localized geometry or forcing produces rapid streamwise variation, as occurs near steps and humps [21 , 22 ] , during the scattering of high-speed boundary-layer modes by wall roughness [23 ] , and near surface indentations, which have been treated with adaptive harmonic linearized formulations [24 , 25 ] .
The harmonic Navier–Stokes (HNS) framework [26 , 27 ] avoids the parabolizing approximation by retaining all streamwise derivatives and solving the perturbation equations as a global boundary-value problem. It can therefore represent upstream and downstream influence and the elliptic interaction between a disturbance and a localized surface feature. HNS exploits the periodic structure of the disturbance by expanding it into a finite set of temporal and spanwise Fourier harmonics. Its computational cost exceeds that of a parabolized marched method but can remain substantially below that of direct numerical simulation when the disturbance spectrum is represented accurately by a limited number of harmonics [14 ] .
Recent HNS developments include both linear and nonlinear formulations. Paredes et al. [28 ] solved the harmonic linearized Navier–Stokes (HLNS) equations about spanwise- or azimuthally varying base flows using a block-pentadiagonal direct factorization, while Zhao and Dong [29 ] coupled HLNS to a fitted bow shock to study hypersonic receptivity to freestream vortical, entropy, and acoustic disturbances. Both formulations remain linear, as do the open-source global tools BROADCAST [30 ] and LinStab2D [31 ] , which provide stability and resolvent analysis of compressible flows. Scholten et al. [32 ] extended HNS to nonlinear compressible boundary-layer disturbances, closing the coupled harmonic system implicitly [33 ] . A complementary route is provided by the one-way Navier–Stokes equations [34 ] . Recursive and carrier-wave-factored formulations reduce the cost of resolving multiple downstream-traveling components in wall-bounded and compressible flows [35 , 36 ] , and Sleeman et al. [37 ] extended the approach to nonlinear harmonic interactions, reporting greater robustness than nonlinear PSE under strong forcing. One-way methods exclude upstream influence by construction; global HNS remains appropriate when such influence or the elliptic interaction with localized geometry is significant.
The Delft Harmonic Navier–Stokes Solver (DeHNSSo) was introduced by Westerbeek et al. [14 ] , hereafter DeHNSSo-v1, as an open-source MATLAB implementation of nonlinear HNS for spanwise-invariant laminar base flows and temporally and spanwise periodic disturbances. It discretizes the wall-normal direction by Chebyshev collocation and the streamwise direction by high-order finite differences, assembles the resulting global sparse systems, and closes the nonlinear interactions by harmonic balancing. The original release was validated for Tollmien–Schlichting (TS) waves, swept-wing crossflow instability, including the interaction of the latter with a smooth hump and a forward-facing step, all within a strictly incompressible framework. To the authors’ knowledge, it remains the only openly available nonlinear HNS implementation for this class of boundary-layer problems. Open-source nonlinear alternatives based on parabolized formulations include glimPSE for incompressible flows [38 ] and Krypton for transonic flows [39 ] . Three problem classes fall outside the scope of the original release of DeHNSSo. The first concerns geometry. The swept wing considered above has a curved surface, as do airfoils, leading edges, and concave walls supporting Görtler-type instabilities [40 ] . DeHNSSo-v1 applied its user-tuned algebraic grids only to localized humps on flat plates, without the mixed second-order metrics, and provided neither a standalone preprocessing pipeline for external CFD data nor mesh smoothing. The second concerns the fluid model. The incompressible formulation excludes high-speed regimes, in which instabilities such as Mack modes can dominate, and the low-speed flows in which non-adiabatic wall conditions can be used to control boundary-layer stability [41 ] . Both regimes require the energy equation, which adds a fifth unknown variable, and temperature-dependent transport properties. The third concerns computational cost. DeHNSSo-v1 refactorizes each per-mode operator at every nonlinear iteration and solves the modes sequentially, and at larger resolutions the fill-in of the sparse factorization sets the peak memory. The fixed-point iteration that closes the nonlinear coupling also fails to converge at sufficiently high disturbance amplitudes. These limitations restrict practical use to modest grids and disturbance amplitudes. DeHNSSo-v2 addresses these three limitations. It generalizes the curvilinear formulation to include the complete metric terms required by the compressible operator and introduces a standalone grid-generation and preprocessing pipeline for base flows imported from external CFD solvers, with optional elliptic smoothing for curved and concave surfaces. It introduces a compressible primitive-variable formulation with an energy equation and temperature-dependent transport properties, following Scholten et al. [32 ] . It combines O ( N ξ ) O(N_{\xi}) O ( N ξ ) [14 ] , published PSE [15 , 19 , 42 ] , HNS [32 ] , and DNS [43 ] results, and the experimentally validated Helmholtz-resonator computations of Michelis et al. [12 ] . Section 7 reports computational performance, and Sec. 8 summarizes the principal findings. The coefficient matrices of the compressible operator are given in Appendix A, the elliptic smoother in Appendix B, and the parameters of all calculations in Appendix C.
Governing equations
This section formulates the harmonic Navier–Stokes operator discretized by DeHNSSo-v2. The formulation begins with the fully compressible Navier–Stokes equations in primitive variables, closed by the ideal-gas equation of state and temperature-dependent transport properties. The flow is then decomposed into a steady base flow and a perturbation, transformed to body-fitted coordinates, and expanded into temporal and spanwise Fourier harmonics. The resulting curvilinear operator and its boundary conditions are presented in Secs. 2.6–2.7.
Flow-conservation equations
The nondimensional, fully compressible Navier–Stokes equations in primitive variables are
∂ t ρ + ∇ ⋅ ( ρ u ) = 0 , (1) \partial_t \rho+ \nabla\cdot(\rho\mathbf{u}) = 0, \tag*{(1)} ∂ t ρ + ∇ ⋅ ( ρ u ) = 0 , (1) ρ ∂ t u + ρ ( u ⋅ ∇ ) u = − ∇ p + 1 Re ∇ ⋅ τ , (2) \rho\partial_t \mathbf{u} + \rho(\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \frac{1}{\operatorname{Re}}\nabla\cdot\boldsymbol{\tau}, \tag*{(2)} ρ ∂ t u + ρ ( u ⋅ ∇ ) u = − ∇ p + Re 1 ∇ ⋅ τ , (2) ρ ∂ t T + ρ ( u ⋅ ∇ ) T = − ( γ − 1 ) ρ T ∇ ⋅ u + γ Re Pr ∇ ⋅ ( κ ∇ T ) + γ ( γ − 1 ) Ma 2 Re Φ . (3) \rho\partial_t T + \rho(\mathbf{u}\cdot\nabla)T = -(\gamma-1)\rho T\nabla\cdot\mathbf{u} + \frac{\gamma}{\operatorname{Re}\operatorname{Pr}}\nabla\cdot(\kappa\nabla T) + \frac{\gamma(\gamma-1)\operatorname{Ma}^{2}}{\operatorname{Re}}\Phi. \tag*{(3)} ρ ∂ t T + ρ ( u ⋅ ∇ ) T = − ( γ − 1 ) ρT ∇ ⋅ u + Re Pr γ ∇ ⋅ ( κ ∇ T ) + Re γ ( γ − 1 ) Ma 2 Φ. (3) The physical coordinates ( x , y , z ) (x,y,z) ( x , y , z ) z z z u = ( u , v , w ) T \mathbf{u}=(u,v,w)^{\mathsf{T}} u = ( u , v , w ) T ρ \rho ρ p p p T T T q = ( u , v , w , ρ , T ) T \mathbf{q}=(u,v,w,\rho,T)^{\mathsf{T}} q = ( u , v , w , ρ , T ) T
The nondimensional viscous stress tensor and dissipation function are
τ = μ [ ∇ u + ( ∇ u ) T − 2 3 ( ∇ ⋅ u ) I ] , Φ = τ : ∇ u , (4) \boldsymbol{\tau}=\mu\left[\nabla\mathbf{u}+(\nabla\mathbf{u})^{\mathsf{T}}-\frac{2}{3}(\nabla\cdot\mathbf{u})\mathbf{I}\right],\qquad\Phi=\boldsymbol{\tau}:\nabla\mathbf{u}, \tag*{(4)} τ = μ [ ∇ u + ( ∇ u ) T − 3 2 ( ∇ ⋅ u ) I ] , Φ = τ : ∇ u , (4) where μ \mu μ κ \kappa κ I \mathbf{I} I λ = − 2 μ / 3 \lambda=-2\mu/3 λ = − 2 μ /3
The reference length, velocity, and time scales, and consequently the nondimensional frequency ω \omega ω β \beta β Re \operatorname{Re} Re [14 ] . The additional compressible variables are nondimensionalized according to
ρ = ρ ∗ ρ ∞ ∗ , T = T ∗ T ∞ ∗ , p = p ∗ ρ ∞ ∗ u ∞ ∗ 2 , (5) \rho=\frac{\rho^{*}}{\rho_{\infty}^{*}},\qquad T=\frac{T^{*}}{T_{\infty}^{*}},\qquad p=\frac{p^{*}}{\rho_{\infty}^{*}u_{\infty}^{*2}}, \tag*{(5)} ρ = ρ ∞ ∗ ρ ∗ , T = T ∞ ∗ T ∗ , p = ρ ∞ ∗ u ∞ ∗ 2 p ∗ , (5) μ = μ ∗ μ ∞ ∗ , κ = κ ∗ κ ∞ ∗ . \mu=\frac{\mu^{*}}{\mu_{\infty}^{*}},\qquad\kappa=\frac{\kappa^{*}}{\kappa_{\infty}^{*}}. μ = μ ∞ ∗ μ ∗ , κ = κ ∞ ∗ κ ∗ . A superscript ( ⋅ ) ∗ (\cdot)^{*} ( ⋅ ) ∗ ( ⋅ ) ∞ (\cdot)_{\infty} ( ⋅ ) ∞
Ma = u ∞ ∗ c ∞ ∗ , Pr = c p , ∞ ∗ μ ∞ ∗ κ ∞ ∗ , γ = c p , ∞ ∗ c v , ∞ ∗ . (6) \operatorname{Ma}=\frac{u_{\infty}^{*}}{c_{\infty}^{*}},\qquad\operatorname{Pr}=\frac{c_{p,\infty}^{*}\mu_{\infty}^{*}}{\kappa_{\infty}^{*}},\qquad\gamma=\frac{c_{p,\infty}^{*}}{c_{v,\infty}^{*}}. \tag*{(6)} Ma = c ∞ ∗ u ∞ ∗ , Pr = κ ∞ ∗ c p , ∞ ∗ μ ∞ ∗ , γ = c v , ∞ ∗ c p , ∞ ∗ . (6) The parameters Ma and Pr are the freestream Mach and Prandtl numbers, respectively, and γ \gamma γ c ∞ ∗ c_\infty^* c ∞ ∗ c p , ∞ ∗ c_{p,\infty}^* c p , ∞ ∗ c v , ∞ ∗ c_{v,\infty}^* c v , ∞ ∗
Equation of state and transport properties
DeHNSSo-v2 assumes a calorically perfect ideal gas. With the scaling in Eq. (5) , the thermal equation of state becomes
p = ρ T γ M a 2 . (7) p = \frac{\rho T}{\gamma\mathrm{Ma}^{2}}. \tag*{(7)} p = γ Ma 2 ρT . (7) The ratio of specific heats is constant, and the dimensional caloric equation of state is e ∗ = c v ∗ T ∗ e^* = c_v^*T^* e ∗ = c v ∗ T ∗
The nondimensional dynamic viscosity is evaluated using Sutherland’s law,
μ ( T ) = T 3 / 2 1 + S ∗ / T ∞ ∗ T + S ∗ / T ∞ ∗ . (8) \mu(T) = T^{3/2}\frac{1+S^*/T_\infty^*}{T+S^*/T_\infty^*}. \tag*{(8)} μ ( T ) = T 3/2 T + S ∗ / T ∞ ∗ 1 + S ∗ / T ∞ ∗ . (8) where S ∗ S^* S ∗ κ ( T ) = μ ( T ) \kappa(T)=\mu(T) κ ( T ) = μ ( T )
Perturbation decomposition
The flow is decomposed into a steady, spanwise-invariant base flow and a three-dimensional, time-dependent perturbation,
q ( x , y , z , t ) = q ‾ ( x , y ) + q ′ ( x , y , z , t ) . (9) \mathbf{q}(x,y,z,t)=\overline{\mathbf{q}}(x,y)+\mathbf{q}'(x,y,z,t). \tag*{(9)} q ( x , y , z , t ) = q ( x , y ) + q ′ ( x , y , z , t ) . (9) The pressure is similarly decomposed as p = p ‾ + p ′ p=\overline{p}+p' p = p + p ′ p ‾ = ρ ‾ T ‾ / ( γ M a 2 ) \overline{p}=\overline{\rho}\overline{T}/(\gamma\mathrm{Ma}^{2}) p = ρ T / ( γ Ma 2 ) ρ = ρ ‾ + ρ ′ \rho=\overline{\rho}+\rho' ρ = ρ + ρ ′ T = T ‾ + T ′ T=\overline{T}+T' T = T + T ′ (7) gives the exact pressure perturbation
p ′ = T ‾ ρ ′ + ρ ‾ T ′ + ρ ′ T ′ γ M a 2 . (10) p'=\frac{\overline{T}\rho'+\overline{\rho}T'+\rho'T'}{\gamma\mathrm{Ma}^{2}}. \tag*{(10)} p ′ = γ Ma 2 T ρ ′ + ρ T ′ + ρ ′ T ′ . (10) The first two terms are linear in the perturbations and enter the HNS operator, whereas the product ρ ′ T ′ \rho'T' ρ ′ T ′
The temperature dependence of the transport properties is linearized about the base-flow temperature,
μ ( T ‾ + T ′ ) ≈ μ ‾ + μ ‾ T T ′ , κ ( T ‾ + T ′ ) ≈ κ ‾ + κ ‾ T T ′ , (11) \mu(\overline{T}+T')\approx\overline{\mu}+\overline{\mu}_T T', \qquad\kappa(\overline{T}+T')\approx\overline{\kappa}+\overline{\kappa}_T T', \tag*{(11)} μ ( T + T ′ ) ≈ μ + μ T T ′ , κ ( T + T ′ ) ≈ κ + κ T T ′ , (11) where μ ‾ = μ ( T ‾ ) \overline{\mu}=\mu(\overline{T}) μ = μ ( T ) κ ‾ = κ ( T ‾ ) \overline{\kappa}=\kappa(\overline{T}) κ = κ ( T )
μ ‾ T = d μ d T ∣ T ‾ , κ ‾ T = d κ d T ∣ T ‾ . (12) \overline{\mu}_T=\left.\frac{\mathrm{d}\mu}{\mathrm{d}T}\right|_{\overline{T}}, \qquad\overline{\kappa}_T=\left.\frac{\mathrm{d}\kappa}{\mathrm{d}T}\right|_{\overline{T}}. \tag*{(12)} μ T = d T d μ T , κ T = d T d κ T . (12) The equality κ = μ \kappa=\mu κ = μ κ ‾ = μ ‾ \overline{\kappa}=\overline{\mu} κ = μ κ ‾ T = μ ‾ T \overline{\kappa}_T=\overline{\mu}_T κ T = μ T ∇ μ ‾ = μ ‾ T ∇ T ‾ \nabla\overline{\mu}=\overline{\mu}_T\nabla\overline{T} ∇ μ = μ T ∇ T ∇ μ ‾ T = μ ‾ T T ∇ T ‾ \nabla\overline{\mu}_T=\overline{\mu}_{TT}\nabla\overline{T} ∇ μ T = μ T T ∇ T μ ‾ T T = d 2 μ / d T 2 ∣ T ‾ \overline{\mu}_{TT}=\mathrm{d}^2\mu/\mathrm{d}T^2|_{\overline{T}} μ T T = d 2 μ / d T 2 ∣ T κ ‾ \overline{\kappa} κ
Substituting Eq. (9) into Eqs. (1) –(3) and subtracting the steady equations, which the base flow is assumed to satisfy, gives
L ( q ‾ ) q ′ = r ( 2 ) ( q ′ , q ′ ) + r ( 3 ) ( q ′ , q ′ , q ′ ) , (13) \mathcal{L}(\overline{\mathbf{q}})\mathbf{q}'=\mathbf{r}^{(2)}(\mathbf{q}',\mathbf{q}')+\mathbf{r}^{(3)}(\mathbf{q}',\mathbf{q}',\mathbf{q}'), \tag*{(13)} L ( q ) q ′ = r ( 2 ) ( q ′ , q ′ ) + r ( 3 ) ( q ′ , q ′ , q ′ ) , (13) where L \mathcal{L} L r ( 2 ) \mathbf{r}^{(2)} r ( 2 ) r ( 3 ) \mathbf{r}^{(3)} r ( 3 )
Coordinate transformation and first-order inverse metrics
The Cartesian derivatives appearing in both the linear operator and the nonlinear forcing are transformed to body-fitted coordinates. The transformation acts in the ( x , y ) (x,y) ( x , y ) z z z
Consider the forward mapping x = x ( ξ , η ) x=x(\xi,\eta) x = x ( ξ , η ) y = y ( ξ , η ) y=y(\xi,\eta) y = y ( ξ , η ) Figure 1 . Its first derivatives are denoted by a = x ξ a=x_\xi a = x ξ b = x η b=x_\eta b = x η c = y ξ c=y_\xi c = y ξ d = y η d=y_\eta d = y η J = a d − b c J=ad-bc J = a d − b c J J J
Figure 1. Physical-to-computational mapping used by DeHNSSo-v2. (a) Body-fitted grid in the physical ( x , y ) (x,y) ( x , y ) δ 99 \delta_{99} δ 99 ( ξ , η ) (\xi,\eta) ( ξ , η )
The first-order inverse metrics are
ξ x = d J , ξ y = − b J , η x = − c J , η y = a J . (14) \xi_x=\frac{d}{J}, \qquad\xi_y=-\frac{b}{J}, \qquad\eta_x=-\frac{c}{J}, \qquad\eta_y=\frac{a}{J}. \tag*{(14)} ξ x = J d , ξ y = − J b , η x = − J c , η y = J a . (14) These metrics are the coefficients of the chain rule that expresses the Cartesian derivatives of the perturbation in terms of computational ones,
∂ x = ξ x ∂ ξ + η x ∂ η , ∂ y = ξ y ∂ ξ + η y ∂ η . (15) \partial_x = \xi_x \partial_\xi+ \eta_x \partial_\eta,\qquad\partial_y = \xi_y \partial_\xi+ \eta_y \partial_\eta. \tag*{(15)} ∂ x = ξ x ∂ ξ + η x ∂ η , ∂ y = ξ y ∂ ξ + η y ∂ η . (15) Second-order inverse metrics
The second-order inverse metrics are obtained by differentiating the inverse-coordinate identities. In compact form, the physical Hessian of either inverse coordinate s ∈ { ξ , η } s \in\{\xi,\eta\} s ∈ { ξ , η }
H x , y ( s ) = − J f − T [ s x H ξ , η ( x ) + s y H ξ , η ( y ) ] J f − 1 , (16) \mathbf{H}_{x,y}(s)=-\mathbf{J}_f^{-\mathsf{T}}\left[s_x\mathbf{H}_{\xi,\eta}(x)+s_y\mathbf{H}_{\xi,\eta}(y)\right]\mathbf{J}_f^{-1}, \tag*{(16)} H x , y ( s ) = − J f − T [ s x H ξ , η ( x ) + s y H ξ , η ( y ) ] J f − 1 , (16) where H ξ , η \mathbf{H}_{\xi,\eta} H ξ , η J f \mathbf{J}_f J f a a a b b b c c c d d d J J J ξ x x \xi_{xx} ξ xx ξ x y \xi_{xy} ξ x y ξ y y \xi_{yy} ξ y y η x x \eta_{xx} η xx η x y \eta_{xy} η x y η y y \eta_{yy} η y y
For illustration purposes, the second x x x ξ \xi ξ
ξ x x = − d 3 x ξ ξ + b d 2 y ξ ξ − c 2 d x η η + b c 2 y η η + 2 c d 2 x ξ η − 2 b c d y ξ η J 3 . (17) \xi_{xx}=\frac{-d^3x_{\xi\xi}+bd^2y_{\xi\xi}-c^2d\,x_{\eta\eta}+bc^2y_{\eta\eta}+2cd^2x_{\xi\eta}-2bcd\,y_{\xi\eta}}{J^3}. \tag*{(17)} ξ xx = J 3 − d 3 x ξ ξ + b d 2 y ξ ξ − c 2 d x η η + b c 2 y η η + 2 c d 2 x ξ η − 2 b c d y ξ η . (17) The implementation evaluates the pure derivatives ξ x x \xi_{xx} ξ xx ξ y y \xi_{yy} ξ y y η x x \eta_{xx} η xx η y y \eta_{yy} η y y
ξ x y = ξ y ∂ ξ ξ x + η y ∂ η ξ x , η x y = ξ y ∂ ξ η x + η y ∂ η η x . (18) \xi_{xy}=\xi_y\partial_\xi\xi_x+\eta_y\partial_\eta\xi_x,\qquad\eta_{xy}=\xi_y\partial_\xi\eta_x+\eta_y\partial_\eta\eta_x. \tag*{(18)} ξ x y = ξ y ∂ ξ ξ x + η y ∂ η ξ x , η x y = ξ y ∂ ξ η x + η y ∂ η η x . (18) The second Cartesian derivatives transform according to
∂ x x 2 = ξ x 2 ∂ ξ ξ 2 + 2 ξ x η x ∂ ξ η 2 + η x 2 ∂ η η 2 + ξ x x ∂ ξ + η x x ∂ η , (19) \partial_{xx}^2=\xi_x^2\partial_{\xi\xi}^2+2\xi_x\eta_x\partial_{\xi\eta}^2+\eta_x^2\partial_{\eta\eta}^2+\xi_{xx}\partial_\xi+\eta_{xx}\partial_\eta, \tag*{(19)} ∂ xx 2 = ξ x 2 ∂ ξ ξ 2 + 2 ξ x η x ∂ ξ η 2 + η x 2 ∂ η η 2 + ξ xx ∂ ξ + η xx ∂ η , (19) ∂ y y 2 = ξ y 2 ∂ ξ ξ 2 + 2 ξ y η y ∂ ξ η 2 + η y 2 ∂ η η 2 + ξ y y ∂ ξ + η y y ∂ η , (20) \partial_{yy}^2=\xi_y^2\partial_{\xi\xi}^2+2\xi_y\eta_y\partial_{\xi\eta}^2+\eta_y^2\partial_{\eta\eta}^2+\xi_{yy}\partial_\xi+\eta_{yy}\partial_\eta, \tag*{(20)} ∂ y y 2 = ξ y 2 ∂ ξ ξ 2 + 2 ξ y η y ∂ ξ η 2 + η y 2 ∂ η η 2 + ξ y y ∂ ξ + η y y ∂ η , (20) ∂ x y 2 = ξ x ξ y ∂ ξ ξ 2 + ( ξ x η y + ξ y η x ) ∂ ξ η 2 + η x η y ∂ η η 2 + ξ x y ∂ ξ + η x y ∂ η . (21) \partial_{xy}^2=\xi_x\xi_y\partial_{\xi\xi}^2+(\xi_x\eta_y+\xi_y\eta_x)\partial_{\xi\eta}^2+\eta_x\eta_y\partial_{\eta\eta}^2+\xi_{xy}\partial_\xi+\eta_{xy}\partial_\eta. \tag*{(21)} ∂ x y 2 = ξ x ξ y ∂ ξ ξ 2 + ( ξ x η y + ξ y η x ) ∂ ξ η 2 + η x η y ∂ η η 2 + ξ x y ∂ ξ + η x y ∂ η . (21) Harmonic ansatz and curvilinear operator
The perturbation is expanded into a truncated ensemble of temporal and spanwise Fourier harmonics,
q ′ ( ξ , η , z , t ) = ∑ m = − M M ∑ n = − N N q ^ m , n ( ξ , η ) exp [ i ( β n z − ω m t ) ] , (22) \mathbf{q}'(\xi,\eta,z,t)=\sum_{m=-M}^{M}\sum_{n=-N}^{N}\hat{\mathbf{q}}_{m,n}(\xi,\eta)\exp\left[i\left(\beta_n z-\omega_m t\right)\right], \tag*{(22)} q ′ ( ξ , η , z , t ) = m = − M ∑ M n = − N ∑ N q ^ m , n ( ξ , η ) exp [ i ( β n z − ω m t ) ] , (22) where ω m = m ω \omega_m=m\omega ω m = mω β n = n β \beta_n=n\beta β n = n β
q ^ − m , − n = q ^ m , n † , (23) \hat{\mathbf{q}}_{-m,-n}=\hat{\mathbf{q}}_{m,n}^{\dagger}, \tag*{(23)} q ^ − m , − n = q ^ m , n † , (23) where † \dagger † ( − m , − n ) (-m,-n) ( − m , − n ) ( 2 M + 1 ) ( 2 N + 1 ) (2M+1)(2N+1) ( 2 M + 1 ) ( 2 N + 1 ) 1 2 [ ( 2 M + 1 ) ( 2 N + 1 ) + 1 ] \frac{1}{2}[(2M+1)(2N+1)+1] 2 1 [( 2 M + 1 ) ( 2 N + 1 ) + 1 ] [14 ] describe an additional reduction based on spanwise-reflection symmetry, but this requires w ˉ = 0 \bar{w}=0 w ˉ = 0
For each harmonic, the Fourier expansion replaces the temporal and spanwise derivatives by
∂ t → − i ω m , ∂ z → i β n , ∂ z z 2 → − β n 2 . (24) \partial_t\rightarrow-i\omega_m,\qquad\partial_z\rightarrow i\beta_n,\qquad\partial_{zz}^2\rightarrow-\beta_n^2. \tag*{(24)} ∂ t → − i ω m , ∂ z → i β n , ∂ z z 2 → − β n 2 . (24) Combining these substitutions with the coordinate transformations in Eqs. (15) and (19) –(21) gives, for each harmonic,
A m , n q ^ m , n + B m , n ∂ η q ^ m , n + C ∂ η η 2 q ^ m , n + D m , n ∂ ξ q ^ m , n + E ∂ ξ η 2 q ^ m , n + F ∂ ξ ξ 2 q ^ m , n = r m , n , (25) \mathbf{A}_{m,n}\hat{\mathbf{q}}_{m,n}+\mathbf{B}_{m,n}\partial_\eta\hat{\mathbf{q}}_{m,n}+\mathbf{C}\,\partial_{\eta\eta}^2\hat{\mathbf{q}}_{m,n}
+\mathbf{D}_{m,n}\partial_\xi\hat{\mathbf{q}}_{m,n}+\mathbf{E}\,\partial_{\xi\eta}^2\hat{\mathbf{q}}_{m,n}+\mathbf{F}\,\partial_{\xi\xi}^2\hat{\mathbf{q}}_{m,n}=\mathbf{r}_{m,n}, \tag*{(25)} A m , n q ^ m , n + B m , n ∂ η q ^ m , n + C ∂ η η 2 q ^ m , n + D m , n ∂ ξ q ^ m , n + E ∂ ξ η 2 q ^ m , n + F ∂ ξ ξ 2 q ^ m , n = r m , n , (25) where r m , n \mathbf{r}_{m,n} r m , n (28) . The inverse metrics are incorporated into the coefficient matrices, so Eq. (25) is discretized directly on the body-fitted ( ξ , η ) (\xi,\eta) ( ξ , η ) [44 ] . The base-flow gradients entering the matrices are precomputed in physical coordinates as described in § 3.2 .
The operator separates the wall-normal, streamwise, and mixed derivatives. The locality of the streamwise stencil makes the discretized operator block banded, which the direct solvers of § 4.2 exploit. The 5 × 5 5\times5 5 × 5 A m , n \mathbf{A}_{m,n} A m , n B m , n \mathbf{B}_{m,n} B m , n D m , n \mathbf{D}_{m,n} D m , n
Spatial discretization yields the per-mode system
M m , n q ^ m , n = r m , n , (26) \mathbf{M}_{m,n}\hat{\mathbf{q}}_{m,n}=\mathbf{r}_{m,n}, \tag*{(26)} M m , n q ^ m , n = r m , n , (26) where M m , n \mathbf{M}_{m,n} M m , n
M m , n = M ( 0 ) + ω m M ( ω ) + β n M ( β ) + β n 2 M ( β 2 ) + β n M ( β ξ ) . (27) \mathbf{M}_{m,n}=\mathbf{M}^{(0)}+\omega_m\mathbf{M}^{(\omega)}+\beta_n\mathbf{M}^{(\beta)}+\beta_n^2\mathbf{M}^{(\beta^2)}+\beta_n\mathbf{M}^{(\beta\xi)}. \tag*{(27)} M m , n = M ( 0 ) + ω m M ( ω ) + β n M ( β ) + β n 2 M ( β 2 ) + β n M ( β ξ ) . (27) The final term contains the compressible i β n ∂ ξ i\beta_n\partial_\xi i β n ∂ ξ ( 0 , 0 ) (0,0) ( 0 , 0 ) § 2.7 ). These blocks are constructed once and recombined for each retained harmonic, as described in § 4 .
Harmonic balancing resolves r ( 2 ) \mathbf{r}^{(2)} r ( 2 ) r ( 3 ) \mathbf{r}^{(3)} r ( 3 ) (13) into a forcing for each retained harmonic,
r m , n = r m , n ( 2 ) + r m , n ( 3 ) . (28) \mathbf{r}_{m,n}=\mathbf{r}_{m,n}^{(2)}+\mathbf{r}_{m,n}^{(3)}. \tag*{(28)} r m , n = r m , n ( 2 ) + r m , n ( 3 ) . (28) Both are obtained by convolution of retained Fourier modes, while contributions outside the selected truncation ensemble are discarded. The nonlinear forcing terms include the density–velocity coupling in continuity; the unsteady and convection terms weighted by the base-flow and perturbation densities in the momentum and energy equations; the product ρ ′ T ′ / ( γ M a 2 ) \rho'T'/(\gamma\mathrm{Ma}^2) ρ ′ T ′ / ( γ Ma 2 ) (11) .
As in DeHNSSo-v1 [14 ] , the nonlinear system is solved by default using a fixed-point (i.e., Picard) iteration. At iteration k k k r m , n ( k ) \mathbf{r}_{m,n}^{(k)} r m , n ( k ) { q ^ ( k ) } \{\hat{\mathbf{q}}^{(k)}\} { q ^ ( k ) } (25) while the individual harmonic systems are solved for { q ^ ( k + 1 ) } \{\hat{\mathbf{q}}^{(k+1)}\} { q ^ ( k + 1 ) }
σ m , n ( k ) = 1 N ξ max i A i ( k ) ∑ i = 1 i b ∣ A i ( k ) − A i ( k − 1 ) ∣ m , n , (29) \sigma_{m,n}^{(k)}=\frac{1}{N_{\xi}\max_i A_i^{(k)}}\sum_{i=1}^{i_b}\left|A_i^{(k)}-A_i^{(k-1)}\right|_{m,n}, \tag*{(29)} σ m , n ( k ) = N ξ max i A i ( k ) 1 i = 1 ∑ i b A i ( k ) − A i ( k − 1 ) m , n , (29) where A i A_i A i i i i σ m , n ( k ) \sigma_{m,n}^{(k)} σ m , n ( k ) [45 ] combines the last m A m_A m A = m A =m_A = m A q ^ m , n ( k + 1 ) = ( 1 − θ ) q ^ m , n ( k ) + θ M m , n − 1 r m , n ( k ) \hat{\mathbf{q}}_{m,n}^{(k+1)}=(1-\theta)\hat{\mathbf{q}}_{m,n}^{(k)}+\theta\mathbf{M}_{m,n}^{-1}\mathbf{r}_{m,n}^{(k)} q ^ m , n ( k + 1 ) = ( 1 − θ ) q ^ m , n ( k ) + θ M m , n − 1 r m , n ( k ) θ ∈ ( 0 , 1 ] \theta\in(0,1] θ ∈ ( 0 , 1 ] § 4.2 ). The coupled Newton alternative for cases beyond the convergence range of Picard iteration is described in § 4.3 .
In the limit M a → 0 \mathrm{Ma}\to0 Ma → 0 [14 ] . This asymptotic consistency is verified numerically in § 5.1 using a Cartesian benchmark at M a = 10 − 3 \mathrm{Ma}=10^{-3} Ma = 1 0 − 3
Boundary conditions
The treatment of no-slip and inhomogeneous wall conditions, freestream decay, inflow disturbances, outflow buffers, and embedded boundaries at sharp wall features follows Westerbeek et al. [14 ] . For the compressible state vector
q ^ m , n = ( u ^ m , n , v ^ m , n , w ^ m , n , ρ ^ m , n , T ^ m , n ) T , \hat{\mathbf{q}}_{m,n}=\left(\hat{u}_{m,n},\hat{v}_{m,n},\hat{w}_{m,n},\hat{\rho}_{m,n},\hat{T}_{m,n}\right)^{\mathrm{T}}, q ^ m , n = ( u ^ m , n , v ^ m , n , w ^ m , n , ρ ^ m , n , T ^ m , n ) T , DeHNSSo-v2 introduces thermal wall conditions, a density-based far-field condition, compressible and PSE inflow profiles, and a parabolized outflow treatment, as illustrated in Figure 1 (b).
Velocity and wall conditions. At the wall (η = 0 \eta=0 η = 0
u ^ m , n ( ξ , 0 ) = u ^ m , n w , (30) \hat{\mathbf{u}}_{m,n}(\xi,0)=\hat{\mathbf{u}}_{m,n}^{\mathrm{w}}, \tag*{(30)} u ^ m , n ( ξ , 0 ) = u ^ m , n w , (30) where u ^ m , n w \hat{\mathbf{u}}_{m,n}^{\mathrm{w}} u ^ m , n w § 3 ). The temperature perturbation satisfies either an isothermal or an adiabatic condition,
T ^ m , n ( ξ , 0 ) = 0 (isothermal) , \hat{T}_{m,n}(\xi,0)=0\qquad\text{(isothermal)}, T ^ m , n ( ξ , 0 ) = 0 (isothermal) , ∂ η T ^ m , n ( ξ , 0 ) = 0 (adiabatic) . (31) \partial_{\eta}\hat{T}_{m,n}(\xi,0)=0\qquad\text{(adiabatic)}. \tag*{(31)} ∂ η T ^ m , n ( ξ , 0 ) = 0 (adiabatic) . (31) The density perturbation ρ ^ m , n \hat{\rho}_{m,n} ρ ^ m , n [14 ] .
Density-based far-field condition. The homogeneous Dirichlet condition of the pressure-based formulation of Westerbeek et al. [14 ] may reflect acoustic disturbances at finite Mach number. At the upper boundary η = H \eta=H η = H [32 ] ,
ρ ^ m , n = u ^ m , n = w ^ m , n = T ^ m , n = 0 , \hat{\rho}_{m,n}=\hat{u}_{m,n}=\hat{w}_{m,n}=\hat{T}_{m,n}=0, ρ ^ m , n = u ^ m , n = w ^ m , n = T ^ m , n = 0 , v ^ m , n unconstrained , (32) \hat{v}_{m,n}\ \text{unconstrained}, \tag*{(32)} v ^ m , n unconstrained , (32) which admits a wall-normal mass-flux perturbation at the upper boundary. This generalizes to all ( m , n ) (m,n) ( m , n ) [14 ] to the MFD, while the incompressible operator retains the original condition. The height H H H [29 ] . Alternatively, a homogeneous Neumann condition may be imposed at η = H \eta=H η = H
Inflow condition. Dirichlet conditions are imposed on all perturbation components at the inflow. User-defined primary modes are prescribed there, while the MFD and remaining higher harmonics are set to zero and generated downstream by nonlinear interactions. As in DeHNSSo-v1, the inflow profiles may instead be supplied externally (LOAD) or left undisturbed (ZERO).
An inflow profile q ^ m , n ( η ) \hat{\mathbf{q}}_{m,n}(\eta) q ^ m , n ( η ) [14 ] . The base flow is assumed parallel at the inflow station, and the streamwise derivatives are replaced according to
∂ ξ → i α , ∂ ξ ξ 2 → − α 2 . (33) \partial_{\xi}\to\mathrm{i}\alpha,\qquad\partial_{\xi\xi}^{2}\to-\alpha^{2}. \tag*{(33)} ∂ ξ → i α , ∂ ξ ξ 2 → − α 2 . (33) For prescribed frequency and spanwise wavenumber, these substitutions yield the local eigenvalue problem for the complex streamwise wavenumber α \alpha α
DeHNSSo-v2 additionally provides a self-contained linear PSE inflow treatment, selected as PSE, which accounts for the nonparallel development upstream of the HNS domain; the compressible version follows Chang et al. [46 ] . For each primary mode, the Fourier coefficient of Eq. (22) is factorized as
q ^ m , n ( ξ , η ) = q ^ m , n ( ξ 0 , η ) exp [ i ∫ ξ 0 ξ α m , n ( ξ ′ ) d ξ ′ ] , (34) \hat{\mathbf{q}}_{m,n}(\xi,\eta)=\hat{\mathbf{q}}_{m,n}(\xi_{0},\eta)\exp\left[\mathrm{i}\int_{\xi_{0}}^{\xi}\alpha_{m,n}(\xi')\,\mathrm{d}\xi'\right], \tag*{(34)} q ^ m , n ( ξ , η ) = q ^ m , n ( ξ 0 , η ) exp [ i ∫ ξ 0 ξ α m , n ( ξ ′ ) d ξ ′ ] , (34) where α m , n ( ξ ) \alpha_{m,n}(\xi) α m , n ( ξ ) q ˇ m , n ( ξ , η ) \check{q}_{m,n}(\xi,\eta) q ˇ m , n ( ξ , η ) [46 ] , evaluated over the three velocity components. The calculation is initialized with a local eigenmode on a self-similar precursor and marched to the HNS inflow using a first-order backward-ξ \xi ξ α m , n \alpha_{m,n} α m , n Δ ξ > 1 / ∣ α r ∣ \Delta\xi> 1/|\alpha_r| Δ ξ > 1/∣ α r ∣ [47 ] .
Outflow parabolization. The amplitude-damping buffer of Westerbeek et al. [14 ] is retained over ξ ∈ [ ξ b , ξ L ] \xi\in[\xi_b,\xi_L] ξ ∈ [ ξ b , ξ L ] [48 ] ). As an alternative, Eq. (25) is gradually parabolized within the buffer [49 ] by ramping down the streamwise second-derivative contribution,
F ∂ ξ ξ 2 q ^ m , n ⟶ s ( ξ ) F ∂ ξ ξ 2 q ^ m , n , ξ ∈ [ ξ b , ξ L ] , (35) F \partial_{\xi\xi}^{2}\hat{q}_{m,n} \longrightarrow s(\xi)F \partial_{\xi\xi}^{2}\hat{q}_{m,n}, \qquad\xi\in[\xi_b,\xi_L], \tag*{(35)} F ∂ ξ ξ 2 q ^ m , n ⟶ s ( ξ ) F ∂ ξ ξ 2 q ^ m , n , ξ ∈ [ ξ b , ξ L ] , (35) where the ramp s ( ξ ) s(\xi) s ( ξ ) s ( ξ b ) = 1 s(\xi_b)=1 s ( ξ b ) = 1 s ( ξ L ) = 0 s(\xi_L)=0 s ( ξ L ) = 0
A m , n q ^ m , n + B m , n ∂ η q ^ m , n + C ∂ η η 2 q ^ m , n + D m , n ∂ ξ q ^ m , n + E ∂ ξ η 2 q ^ m , n = r m , n . (36) A_{m,n}\hat{q}_{m,n}+B_{m,n}\partial_{\eta}\hat{q}_{m,n}+C\partial_{\eta\eta}^{2}\hat{q}_{m,n}
+D_{m,n}\partial_{\xi}\hat{q}_{m,n}+E\partial_{\xi\eta}^{2}\hat{q}_{m,n}=r_{m,n}. \tag*{(36)} A m , n q ^ m , n + B m , n ∂ η q ^ m , n + C ∂ η η 2 q ^ m , n + D m , n ∂ ξ q ^ m , n + E ∂ ξ η 2 q ^ m , n = r m , n . (36) The buffer thus becomes first order in ξ \xi ξ
Stability-grid generation and numerical discretization
DeHNSSo-v2 introduces a stability-grid generator through the gridgen/ module. DeHNSSo-v1 already included user-tuned algebraic grid construction, whereas the new module provides a solver-independent pipeline for importing, resampling, and differentiating base flows defined on Cartesian, structured-curvilinear, or unstructured grids. It constructs the body-fitted mesh on which the perturbation operator is assembled, with the wall-normal points clustered by the mapping of Malik [50 ] , and evaluates the associated inverse metrics and base-flow derivatives with the selected discretization (§ 3.2 ). The input may be obtained from a numerical simulation or from an analytical solution sampled on its native grid. The processed grid, base flow, metric coefficients, derivatives, and forcing metadata are stored in a MATLAB .mat archive that can be read directly by the solver.
The grid-generation pipeline is summarized in Algorithm 1. ClassifyInput first identifies the input as Cartesian, structured-curvilinear, or unstructured, and the corresponding branch resamples it. For a Cartesian grid, ResampleCartesian interpolates the base flow onto the prescribed streamwise stations and the clustered wall-normal points, after DetectWallSolid has identified any solids resting on the wall, such as steps and blocks, from the points of zero velocity; these are passed to the embedded-boundary treatment of Westerbeek et al. [14 ] . For a structured-curvilinear input, ResampleCurvilinear takes the wall from the boundary row of the native mesh and resamples the geometry and base-flow fields along the columns of the native mesh and along the wall arc length. For unstructured input, ExtractWall identifies the wall from the low-velocity portion of the boundary of the point cloud, smooths it, and parameterizes it by arc length, and NormalMesh constructs an algebraic mesh along the local wall normals. The automatic wall detection assumes a no-slip wall. The algebraic mesh may optionally be smoothed by SmoothMesh (§ 3.1 ) before Interpolate maps the base flow onto it. Finally, BuildStabilityGrid evaluates the differentiation operators and inverse metrics on the resampled grid, together with the physical derivatives of the base flow (§ 3.2 ).
Algorithm 1 DeHNSSo-v2 stability-grid generation pipeline (gridgen/ module). Variable and function names follow the code; they are defined in the gridgen user guide distributed with the code.
procedure StabilityGridGen ( B F n a t , p a r a m s ) (\mathrm{BF}_{\mathrm{nat}},\mathrm{params}) ( BF nat , params )
c l a s s ← C l a s s i f y I n p u t ( B F n a t ) \mathrm{class} \leftarrow\mathrm{ClassifyInput}(\mathrm{BF}_{\mathrm{nat}}) class ← ClassifyInput ( BF nat )
f ^ ← ∅ \hat{f} \leftarrow\varnothing f ^ ← ∅
if params.forcing then
( x c , w , f ^ ) ← P r e p a r e F o r c i n g ( p a r a m s ) (x_c,w,\hat{f}) \leftarrow\mathrm{PrepareForcing}(\mathrm{params}) ( x c , w , f ^ ) ← PrepareForcing ( params )
params.mug ← x c \leftarrow x_c ← x c ← w \leftarrow w ← w
end if
if class = = = then
params.step ← D e t e c t W a l l S o l i d ( B F n a t ) \leftarrow\mathrm{DetectWallSolid}(\mathrm{BF}_{\mathrm{nat}}) ← DetectWallSolid ( BF nat )
B F ← R e s a m p l e C a r t e s i a n ( B F n a t , p a r a m s ) \mathrm{BF} \leftarrow\mathrm{ResampleCartesian}(\mathrm{BF}_{\mathrm{nat}},\mathrm{params}) BF ← ResampleCartesian ( BF nat , params )
else if class = = = then
B F ← R e s a m p l e C u r v i l i n e a r ( B F n a t , p a r a m s ) \mathrm{BF} \leftarrow\mathrm{ResampleCurvilinear}(\mathrm{BF}_{\mathrm{nat}},\mathrm{params}) BF ← ResampleCurvilinear ( BF nat , params )
else
( x w , y w , n ^ ) ← E x t r a c t W a l l ( B F n a t ) (x_w,y_w,\hat{\boldsymbol{n}}) \leftarrow\mathrm{ExtractWall}(\mathrm{BF}_{\mathrm{nat}}) ( x w , y w , n ^ ) ← ExtractWall ( BF nat )
( X , Y ) ← N o r m a l M e s h ( x w , y w , n ^ , p a r a m s ) (X,Y) \leftarrow\mathrm{NormalMesh}(x_w,y_w,\hat{\boldsymbol{n}},\mathrm{params}) ( X , Y ) ← NormalMesh ( x w , y w , n ^ , params )
if params.elliptic then
( X , Y ) ← S m o o t h M e s h ( X , Y , p a r a m s ) (X,Y) \leftarrow\mathrm{SmoothMesh}(X,Y,\mathrm{params}) ( X , Y ) ← SmoothMesh ( X , Y , params )
end if
B F ← I n t e r p o l a t e ( B F n a t , X , Y ) \mathrm{BF} \leftarrow\mathrm{Interpolate}(\mathrm{BF}_{\mathrm{nat}},X,Y) BF ← Interpolate ( BF nat , X , Y )
end if
S t a b G r i d ← B u i l d S t a b i l i t y G r i d ( B F , p a r a m s ) \mathrm{StabGrid} \leftarrow\mathrm{BuildStabilityGrid}(\mathrm{BF},\mathrm{params}) StabGrid ← BuildStabilityGrid ( BF , params )
StabGrid.forcing ← f ^ \leftarrow\hat{f} ← f ^
S a v e A r c h i v e ( S t a b G r i d ) \mathrm{SaveArchive}(\mathrm{StabGrid}) SaveArchive ( StabGrid )
end procedure
Wall forcing, namely the nonzero wall velocity u ^ m , n w \hat{u}_{m,n}^{w} u ^ m , n w (30) , is also prepared by the grid generator. A forcing patch is defined either by its center x c x_c x c w w w PrepareForcing uses the patch to refine the streamwise grid and returns a profile of unit peak, which is stored with the grid; the solver then scales this profile and applies it as the wall velocity. The same interface can also refine the grid around a wall-mounted device whose boundary condition is supplied separately, such as the resonator neck of Sec. 6.4.
Figure 2(a) shows a grid generated for the M3J swept wing of the TU Delft low-turbulence tunnel [10 ] , using the base flow of Moniripiri et al. [43 ] . The near-wall clustering resolves the boundary layer along the continuously curved surface. Panels (b) and (c) show grids over generic concave geometries, used to assess the robustness of the grid construction. Potential mesh singularities (i.e., folding) are mitigated using the optional elliptic smoothing described in § 3.1 .
Elliptic mesh smoothing
Over a strongly concave wall, neighboring wall normals may converge and intersect, causing the algebraic normal-ray mesh to fold; an example is shown in Figure 2 (c). A folded mesh contains cells for which the mapping Jacobian vanishes or changes sign, making the inverse metrics of § 2.4 singular or incorrectly oriented. The optional smoothing procedure starts from the algebraic mesh and, in its default Laplace mode, solves the elliptic inverse-grid equations [51 ] ,
Figure 2. Body-fitted stability grids generated by the gridgen/ module. (a) Grid for the clean M3J swept-wing geometry used in § 6.3 , with wall-normal distances enlarged four times, colored by the wall-tangential base-flow velocity normalized by its local edge value, u ˉ / u ˉ e \bar{u}/\bar{u}_e u ˉ / u ˉ e δ 99 \delta_{99} δ 99
g ξ ξ x η η − 2 g ξ η x ξ η + g η η x ξ ξ = 0 , g_{\xi\xi}x_{\eta\eta}-2g_{\xi\eta}x_{\xi\eta}+g_{\eta\eta}x_{\xi\xi}=0, g ξ ξ x η η − 2 g ξ η x ξ η + g η η x ξ ξ = 0 , g ξ ξ y η η − 2 g ξ η y ξ η + g η η y ξ ξ = 0 , (37) g_{\xi\xi}y_{\eta\eta}-2g_{\xi\eta}y_{\xi\eta}+g_{\eta\eta}y_{\xi\xi}=0, \tag*{(37)} g ξ ξ y η η − 2 g ξ η y ξ η + g η η y ξ ξ = 0 , (37) where
g ξ ξ = x ξ 2 + y ξ 2 , g ξ η = x ξ x η + y ξ y η , g η η = x η 2 + y η 2 . (38) g_{\xi\xi}=x_{\xi}^{2}+y_{\xi}^{2},\qquad g_{\xi\eta}=x_{\xi}x_{\eta}+y_{\xi}y_{\eta},\qquad g_{\eta\eta}=x_{\eta}^{2}+y_{\eta}^{2}. \tag*{(38)} g ξ ξ = x ξ 2 + y ξ 2 , g ξ η = x ξ x η + y ξ y η , g η η = x η 2 + y η 2 . (38) The equations are solved on a uniformly spaced auxiliary computational grid using pointwise successive over-relaxation, with the algebraic mesh as the initial guess; the discretization and the iterative scheme are given in Appendix B. To prevent the outer boundary itself from folding, it is translated along the mean outward wall normal rather than offset along the local normals. After convergence, the physical coordinates are redistributed onto the prescribed wall-normal points.
Discretization, inverse metrics, and base-flow derivatives
Two wall-normal schemes are available: Chebyshev collocation, as in DeHNSSo-v1, and finite differences of order 2, 4, or 6.
or 6. Both use the N η N_{\eta} N η [50 ] , which spans η ∈ [ 0 , H ] \eta\in[0,H] η ∈ [ 0 , H ] [52 ] , which accommodate nonuniform station spacing, both in the perturbation operator and for the metrics and base-flow derivatives. These options are summarized in Table 1 .
Table 1. Spatial discretization in DeHNSSo-v1 [14 ] and DeHNSSo-v2. FD denotes finite differences.
The fourth-order streamwise stencil couples each streamwise grid station to the two stations on either side of it, producing the block-pentadiagonal structure shown in Figure 3 . Chebyshev collocation gives global wall-normal coupling within each streamwise station, whereas the finite-difference alternative produces sparse blocks with local interior stencils, although its boundary rows remain dense. The wall-normal scheme therefore sets the block density and the factorization memory of the direct solvers of § 4.2 and of the Newton solver of § 4.3 .
Figure 3. Sparsity pattern of the mode-invariant block M ( 0 ) \mathbf{M}^{(0)} M ( 0 ) N ξ = 7 N_{\xi}=7 N ξ = 7 N η = 11 N_{\eta}=11 N η = 11 ξ \xi ξ η \eta η η \eta η ( u ^ , v ^ , w ^ , ρ ^ , T ^ ) (\hat{u},\hat{v},\hat{w},\hat{\rho},\hat{T}) ( u ^ , v ^ , w ^ , ρ ^ , T ^ ) x x x y y y z z z ( M x , M y , M z , C , E ) (M_x,M_y,M_z,C,E) ( M x , M y , M z , C , E )
With these operators, the grid generator evaluates the firstand second-order derivatives of the forward mapping and constructs the inverse metrics of § 2.4 and § 2.5 . The first-order metrics follow Eq. (14) . The pure second-order metrics are evaluated from the analytical relations represented by Eq. (16) , whereas the mixed metrics follow Eq. (18) . These coefficients enter the curvilinear chain rules that give the physical first and second derivatives of the base flow; the second derivatives and the mixed metrics are required only by the compressible operator. As in DeHNSSo-v1, compatible physical derivatives supplied with the input base flow may alternatively be retained.
Efficient linear algebra and nonlinear continuation
This section presents the linear-algebra and nonlinear-solver developments introduced in DeHNSSo-v2. They address four limitations of DeHNSSo-v1: the cost of assembling the global operators, the repeated factorization of the per-mode systems during Picard iteration, the memory required by the factors, and the loss of fixed-point convergence at large disturbance amplitudes.
In DeHNSSo-v1, the operator of each active mode is factorized again at every nonlinear iteration; for the nonlinear incompressible Blasius case of § 7.1.1 , this repeated factorize-and-solve operation accounts for 92% of the run time (see Figure 14 a). In addition, the operator assembly grows super-linearly with N ξ N_{\xi} N ξ O ( N ξ ) O(N_{\xi}) O ( N ξ ) § 4.1 ) and reusable per-mode factorizations (§ 4.2 ). A banded-LU back-end provides predictable factor storage and a faster factorization (§ 4.2.2 ), while an optionally out-of-core block-Thomas back-end permits calculations whose complete banded factor does not fit in memory (§ 4.2.3 ). Finally, the Picard iteration inherited from DeHNSSo-v1 is complemented by a Newton continuation method for large disturbance amplitudes (§ 4.3 ).
Figure 14. (a) Share of the wall-clock time per task in DeHNSSo-v1 (v1) [14 ] and DeHNSSo-v2 (v2) for the case of Table 4 , with the total wall-clock time of each calculation; in DeHNSSo-v1, factorization and solve form a single operation. (b) LHS-assembly wall-clock time against the streamwise resolution N ξ N_{\xi} N ξ N η = 40 N_{\eta}=40 N η = 40 M = 1 M=1 M = 1 ( m , 0 ) (m,0) ( m , 0 ) t ˉ f a c t \bar{t}_{\mathrm{fact}} t ˉ fact t ˉ c a c h e \bar{t}_{\mathrm{cache}} t ˉ cache
Sparse operator assembly
Each mode-invariant block in Eq. (27) has order N q N η N ξ N_qN_\eta N_\xi N q N η N ξ N q = 4 N_q = 4 N q = 4 N q = 5 N_q = 5 N q = 5 ( i , j , v ) (\mathbf{i},\mathbf{j},\mathbf{v}) ( i , j , v ) O ( N ξ 1.5 ) O(N_{\xi}^{1.5}) O ( N ξ 1.5 ) N ξ = 200 N_{\xi} = 200 N ξ = 200 3200 3200 3200 Figure 14 b). DeHNSSo-v2 instead stores the triplets of each streamwise station separately and concatenates them once, after all stations have been processed. This reduces the assembly cost to O ( N ξ ) O(N_{\xi}) O ( N ξ ) § 7.1.1 , and agreement with the DeHNSSo-v1 solution is shown in § 5.1 .
Direct solvers
For a fixed base flow and grid, the per-mode operator M m , n \mathbf{M}_{m,n} M m , n (27) depends only on ( ω m , β n ) (\omega_m,\beta_n) ( ω m , β n ) § 2.6 , in which only the forcing r m , n \mathbf{r}_{m,n} r m , n
Except for seq and none, the back-ends of Table 2 exploit this invariance by factorizing each operator once and reusing the factor at every iteration. The back-ends, selected through \texttt{Opt.lu_mode}, differ in the factorization method and its storage. Their agreement is verified in § 5.2 , and their time and memory requirements are compared in § 7.1.2 and § 7.1.3 . The default auto is suitable for most incompressible calculations. For the compressible operator, lapack_band is faster as long as its banded factors fit in memory; beyond that limit, block_thomas can store the factor on disk, which allows finer grids.
Table 2. Direct-solver back-ends selected through \texttt{Opt.lu_mode}. The first four settings use sparse-LU and differ in the number of resident per-mode factors.
Cached sparse-LU factorization
DeHNSSo-v1 passed each per-mode system to a new sparse direct solve at every Picard iteration, although M m , n \mathbf{M}_{m,n} M m , n M m , n \mathbf{M}_{m,n} M m , n Figure 14 c).
Sparse LU introduces fill-in: entries that are zero in M m , n \mathbf{M}_{m,n} M m , n [53 ] to the compressible operator, where it reduces the bandwidth of the block-banded structure, and the column approximate minimum degree ordering [54 ] , which targets fill-in directly, to the incompressible operator, for which it is leaner and faster. Fill-in remains particularly costly for the five-variable compressible system and under Chebyshev collocation, which couples all wall-normal points within each station block.
The sparse-LU settings of Table 2 provide different memory–time compromises. The full setting retains every per-mode factor and seq only the current one. The default auto keeps up to \texttt{Opt.lu_max_cache} factors, a limit that is at least the number of computed modes unless set lower, and, when this limit is reached, removes the oldest factor not needed by the active modes. The none setting keeps no factor and solves each system directly, as in DeHNSSo-v1. With full or auto, the per-mode solves can be distributed over parallel workers through \texttt{Opt.parfor}, which requires the MATLAB Parallel Computing Toolbox; the resulting speed-up is measured in § 7.1.3 .
The assembled blocks and cached factors may also be reused between solver calls through \texttt{Opt.rerun}, provided that all operator-defining quantities remain unchanged and only the right-hand side r m , n \mathbf{r}_{m,n} r m , n § 7.1.1 . Changes to the grid, base flow, frequency, or spanwise wavenumber require the affected operators to be rebuilt.
Banded LU factorization
The locality of the streamwise stencil confines every nonzero entry of the per-mode operator to a band about the diagonal, giving the block-pentadiagonal structure of Figure 3 . DeHNSSo-v2 stores only this band and factorizes it with the general-banded LU routines of LAPACK [55 ] , called through a C MEX interface. This avoids an additional MATLAB matrix representation and provides separate factorization and substitution calls, allowing the factor to be reused throughout the Picard iteration.
The band is set by the streamwise coupling and therefore has the same width for Chebyshev collocation and wall-normal finite differences, although it widens as N η N_{\eta} N η η \eta η Figure 3 ). For a matrix of dimension n n n b b b O ( n b ) O(nb) O ( nb ) O ( n b 2 ) O(nb^{2}) O ( n b 2 ) n n n b b b N q N_q N q O ( N q 2 ) O(N_q^{2}) O ( N q 2 ) O ( N q 3 ) O(N_q^{3}) O ( N q 3 ) § 7.1.2 ). As banded storage allocates memory for every position within the band, including positions that contain zeros, it does not benefit from a sparse η \eta η § 7.1.2 . The advantages of the banded path are predictable memory requirements and a substantially faster factorization, which compensate for the larger factor when solving the five-variable compressible system. The routines apply partial pivoting [56 ] , which is necessary because the HNS operator is not diagonally dominant and, in the incompressible formulation, the continuity rows carry no diagonal entry.
Partitioned block-Thomas factorization
The banded back-end of § 4.2.2 stores the complete factor in memory, which eventually limits the grid size. DeHNSSo-v2 therefore provides a block-Thomas back-end based on the dual-Schur [57 ] approach of Paredes et al. [28 ] , selected through Opt.lu_mode = ’block_thomas’. The streamwise domain is divided into n p n_p n p [28 ] , n p = 2 n_p=2 n p = 2 [58 ] ; storing their explicit inverses, as in Paredes et al. [28 ] , is available as an option (Opt.bt_apply). The factor is built once for each mode and reused across Picard iterations.
The factor can be kept in memory or streamed to disk with Opt.bt_ooc = ’on’. In this out-of-core mode, the dense factor blocks are written to disk as they are formed, and the only the interface system, the sparse operator blocks of the partitions, and the blocks of the current recursion remain in memory. Unlike banded storage, the resident memory therefore depends on the sparsity of the operator blocks and decreases with a finite-difference wall-normal discretization (Figure 3 ), as quantified in § 7.1.2 . This allows resolutions whose in-core factor would exceed the available memory.
The recursion between stations is unpivoted and therefore requires nonsingular diagonal station blocks. This holds for the compressible operator but not for the incompressible one, so the back-end is restricted to compressible calculations.
Newton continuation for high-amplitude nonlinear solutions
The Picard iteration inherited from DeHNSSo-v1 remains the default solver (see § 2.6 ), but at high amplitudes the lagged nonlinear forcing may no longer define a convergent fixed-point iteration. Newton’s method instead carries the coupling between harmonics in the Jacobian [33 , 32 ] , and pseudo-arclength continuation [59 , 60 , 61 ] provides successive initial guesses by advancing the disturbance field and the inflow amplitude A 0 A_0 A 0 Opt.solver = ’newton’.
Let Q ^ \hat{\mathbf{Q}} Q ^ q ^ m , n \hat{q}_{m,n} q ^ m , n (26) is written as
F ( Q ^ , Q ^ † , A 0 ) = M Q ^ − b ( A 0 ) − r ( Q ^ , Q ^ † ) = 0 , (39) \mathcal{F}(\widehat{\mathbf{Q}},\widehat{\mathbf{Q}}^\dagger,A_0)=\mathbf{M}\widehat{\mathbf{Q}}-\mathbf{b}(A_0)-\mathbf{r}(\widehat{\mathbf{Q}},\widehat{\mathbf{Q}}^\dagger)=0, \tag*{(39)} F ( Q , Q † , A 0 ) = M Q − b ( A 0 ) − r ( Q , Q † ) = 0 , (39) where M = diag ( M m , n ) \mathbf{M}=\operatorname{diag}(\mathbf{M}_{m,n}) M = diag ( M m , n ) r \mathbf{r} r b ( A 0 ) \mathbf{b}(A_0) b ( A 0 ) A 0 A_0 A 0 F \mathcal{F} F Q ^ † \widehat{\mathbf{Q}}^\dagger Q † Q ^ \widehat{\mathbf{Q}} Q Q ^ \widehat{\mathbf{Q}} Q F \mathcal{F} F x \mathbf{x} x f \mathbf{f} f J R = ∂ f / ∂ x \mathbf{J}_{R}=\partial\mathbf{f}/\partial\mathbf{x} J R = ∂ f / ∂ x r \mathbf{r} r Q ^ \widehat{\mathbf{Q}} Q Q ^ † \widehat{\mathbf{Q}}^\dagger Q †
Continuation treats A 0 A_0 A 0 ( x c , A 0 , c ) (\mathbf{x}_c,A_{0,c}) ( x c , A 0 , c ) ( x ˙ , A ˙ 0 ) (\dot{\mathbf{x}},\dot{A}_0) ( x ˙ , A ˙ 0 )
S ( x , A 0 ) = x ˙ T ( x − x c ) + A ˙ 0 ( A 0 − A 0 , c ) − Δ s = 0 , (40) \mathcal{S}(\mathbf{x},A_0)=\dot{\mathbf{x}}^{\mathsf{T}}(\mathbf{x}-\mathbf{x}_c)+\dot{A}_0(A_0-A_{0,c})-\Delta s=0, \tag*{(40)} S ( x , A 0 ) = x ˙ T ( x − x c ) + A ˙ 0 ( A 0 − A 0 , c ) − Δ s = 0 , (40) where Δ s \Delta s Δ s f \mathbf{f} f S \mathcal{S} S g = ∂ f / ∂ A 0 \mathbf{g}=\partial\mathbf{f}/\partial A_0 g = ∂ f / ∂ A 0 [60 ] avoids assembling and factorizing the bordered Jacobian by solving
J R z 1 = − f , J R z 2 = g , \mathbf{J}_{R}\mathbf{z}_1=-\mathbf{f},\qquad\mathbf{J}_{R}\mathbf{z}_2=\mathbf{g}, J R z 1 = − f , J R z 2 = g , Δ A 0 = − S − x ˙ T z 1 A ˙ 0 − x ˙ T z 2 , Δ x = z 1 − Δ A 0 z 2 . (41) \Delta A_0=\frac{-\mathcal{S}-\dot{\mathbf{x}}^{\mathsf{T}}\mathbf{z}_1}{\dot{A}_0-\dot{\mathbf{x}}^{\mathsf{T}}\mathbf{z}_2},\qquad\Delta\mathbf{x}=\mathbf{z}_1-\Delta A_0\mathbf{z}_2. \tag*{(41)} Δ A 0 = A ˙ 0 − x ˙ T z 2 − S − x ˙ T z 1 , Δ x = z 1 − Δ A 0 z 2 . (41) The iteration continues until ( ∥ f ∥ 2 2 + S 2 ) 1 / 2 < 10 − 6 (\lVert\mathbf{f}\rVert_2^2+\mathcal{S}^2)^{1/2}<10^{-6} (∥ f ∥ 2 2 + S 2 ) 1/2 < 1 0 − 6 z 2 \mathbf{z}_2 z 2 [ − z 2 T , 1 ] T [-\mathbf{z}_2^{\mathsf{T}},1]^{\mathsf{T}} [ − z 2 T , 1 ] T
By default, J R \mathbf{J}_{R} J R [62 , 63 ] , is carried across nearby continuation steps when the convergence history allows; it is refreshed after Opt.newton_reuse steps or when the corrector exceeds Opt.newton_it_thresh iterations, as also done within the corrector by Busquet et al. [61 ] . Most corrector iterations therefore require only a solve with the retained factor. The full Newton corrector and first-order predictor of [60 ] are available through Opt.newton_corrector and Opt.newton_predictor. The step Δ s \Delta s Δ s [64 ] . A final Newton solve at fixed amplitude returns the solution at the requested A 0 A_0 A 0 § 7.2 .
Verification
This section verifies four aspects of DeHNSSo-v2: recovery of the incompressible limit and grid convergence (§ 5.1 ), independence of the solution from the direct-solver back-end (§ 5.2 ), accuracy and coordinate invariance of the curvilinear formulation (§ 5.3 ), and high-amplitude Newton continuation (§ 5.4 ).
Sections 5.1, 5.2, and 5.4 use the canonical nonlinear TS-wave development case of Bertolotti et al. [15 ] and Westerbeek et al. [14 ] in a Blasius boundary layer. The local Reynolds number is R = R e x R=\sqrt{Re_x} R = R e x R 0 = 400 R_0=400 R 0 = 400 ω 1 = 0.0346 \omega_1=0.0346 ω 1 = 0.0346 F = 86 × 10 − 6 F=86\times10^{-6} F = 86 × 1 0 − 6 β = 0 \beta=0 β = 0 M a = 0.5 Ma=0.5 M a = 0.5 4.5 4.5 4.5
Incompressible-limit recovery and grid convergence
Figure 4 shows the amplitude evolution ∣ u ′ ∣ max , ( m , 0 ) ( R ) |u'|_{\max,(m,0)}(R) ∣ u ′ ∣ m a x , ( m , 0 ) ( R ) R = R e x R=\sqrt{Re_x} R = R e x ( 0 , 0 ) (0,0) ( 0 , 0 ) ( 1 , 0 ) (1,0) ( 1 , 0 ) ( 5 , 0 ) (5,0) ( 5 , 0 ) M = 5 M=5 M = 5 N = 0 N=0 N = 0 ( 1 , 0 ) (1,0) ( 1 , 0 ) ξ b = 90 % \xi_b=90\% ξ b = 90%
Figure 4. Incompressible-limit recovery: amplitude evolution ∣ u ′ ∣ max , ( m , 0 ) |u'|_{\max,(m,0)} ∣ u ′ ∣ m a x , ( m , 0 ) ( 0 , 0 ) (0,0) ( 0 , 0 ) ( 1 , 0 ) (1,0) ( 1 , 0 ) ( 5 , 0 ) (5,0) ( 5 , 0 ) ( R 0 = 400 (R_0=400 ( R 0 = 400 F = 86 × 10 − 6 F=86\times10^{-6} F = 86 × 1 0 − 6 M = 5 ) M=5) M = 5 ) M a = 10 − 3 Ma=10^{-3} M a = 1 0 − 3 1800 × 140 1800\times140 1800 × 140 2674 × 100 2674\times100 2674 × 100
The compressible and incompressible DeHNSSo-v2 solutions agree closely for all six modes, and both reproduce DeHNSSo-v1 computed on its own 2674 × 100 2674\times100 2674 × 100 [14 ] , previously verified against the NPSE reference case of Bertolotti et al. [15 ] . This agreement verifies recovery of the established incompressible solver in the low-Mach limit. At M a = 10 − 3 Ma=10^{-3} M a = 1 0 − 3 r m , n ( 3 ) \mathbf{r}_{m,n}^{(3)} r m , n ( 3 ) (28) is negligible and the quadratic forcing r m , n ( 2 ) \mathbf{r}_{m,n}^{(2)} r m , n ( 2 ) § 6.2 .
Grid convergence is measured using the mode-summed perturbation kinetic energy over the buffer-free domain,
E = ∑ ( m , n ) ∈ A ∫ ξ 0 ξ b ∫ 0 H 1 2 ( ∣ u ^ m , n ∣ 2 + ∣ v ^ m , n ∣ 2 + ∣ w ^ m , n ∣ 2 ) ∣ J ∣ d η d ξ , (42) E = \sum_{(m,n)\in\mathcal{A}} \int_{\xi_0}^{\xi_b} \int_{0}^{H} \frac{1}{2}\left(\lvert\hat{u}_{m,n}\rvert^2 + \lvert\hat{v}_{m,n}\rvert^2 + \lvert\hat{w}_{m,n}\rvert^2\right)\lvert J\rvert\,\mathrm{d}\eta\,\mathrm{d}\xi, \tag*{(42)} E = ( m , n ) ∈ A ∑ ∫ ξ 0 ξ b ∫ 0 H 2 1 ( ∣ u ^ m , n ∣ 2 + ∣ v ^ m , n ∣ 2 + ∣ w ^ m , n ∣ 2 ) ∣ J ∣ d η d ξ , (42) where A \mathcal{A} A ξ b \xi_b ξ b J J J § 2.4 . On the Cartesian grids used here, J J J R = 825 R = 825 R = 825
The resolution is assessed by varying N ξ = 500 N_{\xi} = 500 N ξ = 500 1000 1000 1000 N η = 140 N_{\eta} = 140 N η = 140 N η = 50 N_{\eta} = 50 N η = 50 100 100 100 N ξ = 1800 N_{\xi} = 1800 N ξ = 1800 1800 × 140 1800 \times140 1800 × 140 σ m , n = 10 − 5 \sigma_{m,n} = 10^{-5} σ m , n = 1 0 − 5 (29) . Figure 5 reports the relative difference ϵ = ∣ E − E r e f ∣ / ∣ E r e f ∣ \epsilon= \lvert E - E_{\mathrm{ref}}\rvert/\lvert E_{\mathrm{ref}}\rvert ϵ = ∣ E − E ref ∣ / ∣ E ref ∣ R = 825 R = 825 R = 825 N ξ = 1800 N_{\xi} = 1800 N ξ = 1800 N η = 70 N_{\eta} = 70 N η = 70 1800 × 140 1800 \times140 1800 × 140
Figure 5. Grid convergence for the Blasius TS calculation: relative difference in the domain-integrated energy (black) and local energy at R = 825 R = 825 R = 825 1800 × 140 1800 \times140 1800 × 140 N η = 140 N_{\eta} = 140 N η = 140 N ξ = 1800 N_{\xi} = 1800 N ξ = 1800
Direct-solver equivalence
A key requirement for robustness is that the converged solution must not depend on the direct-solver back-end. All settings of Tab. 2 are therefore applied to the Blasius TS case of § 5.1 on the reduced N ξ × N η = 800 × 40 N_{\xi} \times N_{\eta} = 800 \times40 N ξ × N η = 800 × 40 [14 ] , with only Opt.lu mode varied. The none setting, which factorizes each assembled operator at every solve as in DeHNSSo-v1, serves as the reference. The block_thomas setting is restricted to the compressible operator, so it is applied to the corresponding compressible case at M a = 10 − 3 \mathrm{Ma} = 10^{-3} Ma = 1 0 − 3
The back-ends are compared through the energy measures of § 5.1 . Every calculation is advanced for 80 Picard iterations, beyond the point where the tolerance σ m , n = 10 − 5 \sigma_{m,n} = 10^{-5} σ m , n = 1 0 − 5 (29) is reached, so that the comparison is not affected by independent stopping decisions. Table 3 lists the relative differences; the seq, auto, and full rows coincide because these settings use the same factorization and differ only in factor residency. For scale, continuing the none run from iteration 49, when the stopping criterion is first satisfied, to iteration 80 changes E E E 9.5 × 10 − 7 9.5 \times10^{-7} 9.5 × 1 0 − 7
Table 3. Direct-solver equivalence for the Blasius TS case of § 5.1 . Relative difference ϵ = ∣ E − E r e f ∣ / ∣ E r e f ∣ \epsilon= \lvert E - E_{\mathrm{ref}}\rvert/\lvert E_{\mathrm{ref}}\rvert ϵ = ∣ E − E ref ∣ / ∣ E ref ∣ (42) , evaluated over the buffer-free domain and locally at R = 825 R = 825 R = 825 M a = 10 − 3 \mathrm{Ma} = 10^{-3} Ma = 1 0 − 3
Curvilinear-operator accuracy and coordinate invariance
The curvilinear formulation is verified in two stages. First, the transformed derivative operators are applied to an analytical field to establish their formal order of accuracy. Second, coordinate invariance of the linear HNS solution is tested by comparing Cartesian and distorted-grid solutions for the M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 4.5 4.5 4.5 § 6.1 . The distortion activates all metric terms, including ξ y \xi_y ξ y η x y \eta_{xy} η x y
The grid is distorted with the sinusoidal mapping of Colella et al. [65 ] , in the normalization of Grimm-Strele et al. [66 ] :
x = ξ 0 + L [ s + ε sin ( 2 π s ) sin ( 2 π t ) ] , y = H [ t + ε sin ( 2 π s ) sin ( 2 π t ) ] , (43) \begin{aligned}
x &= \xi_0 + L\left[s + \varepsilon\sin(2\pi s)\sin(2\pi t)\right],\\
y &= H\left[t + \varepsilon\sin(2\pi s)\sin(2\pi t)\right], \tag*{(43)}
\end{aligned} x y = ξ 0 + L [ s + ε sin ( 2 π s ) sin ( 2 π t ) ] , = H [ t + ε sin ( 2 π s ) sin ( 2 π t ) ] , (43) with s = ( ξ − ξ 0 ) / L s = (\xi- \xi_0)/L s = ( ξ − ξ 0 ) / L t = η / H t = \eta/H t = η / H L L L H H H ε \varepsilon ε Figure 6 (a) shows the distorted M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 χ = ∣ r ξ ⋅ r η ∣ / ( ∣ r ξ ∣ ∣ r η ∣ ) \chi= \lvert\mathbf{r}_{\xi}\cdot\mathbf{r}_{\eta}\rvert/(\lvert\mathbf{r}_{\xi}\rvert\lvert\mathbf{r}_{\eta}\rvert) χ = ∣ r ξ ⋅ r η ∣ / (∣ r ξ ∣ ∣ r η ∣) χ = 0 \chi= 0 χ = 0 ε = 0.01 \varepsilon= 0.01 ε = 0.01 0.937 ≤ J ≤ 1.063 0.937 \le J \le1.063 0.937 ≤ J ≤ 1.063 χ max = 0.76 \chi_{\max} = 0.76 χ m a x = 0.76
Figure 6. Curvilinear-operator accuracy and coordinate invariance at ε = 0.01 \varepsilon= 0.01 ε = 0.01 χ \chi χ L 2 L_2 L 2 L ∞ L_\infty L ∞ L 2 L_2 L 2 M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 N ξ ≈ 10 N η N_\xi\approx10N_\eta N ξ ≈ 10 N η N η − 4 N_\eta^{-4} N η − 4
For the analytical test, a smooth field with nonzero first, second, and mixed derivatives is differentiated using the discrete metrics and matrices employed in the assembled operator. For each derivative D D D Figure 6 (b) reports the relative error ∥ D h f − D f ∥ 2 / ∥ D f ∥ 2 \lVert D_h f - D f \rVert_2/\lVert D f \rVert_2 ∥ D h f − D f ∥ 2 / ∥ D f ∥ 2 D h D_h D h (15) and (19) –(21) converge at fourth order over six grids, from 401 × 41 401 \times41 401 × 41 3040 × 306 3040 \times306 3040 × 306 ξ x y = η x y = 0 \xi_{xy} = \eta_{xy} = 0 ξ x y = η x y = 0 (18) .
The coordinate invariance of the complete HNS operator of Eq. (25) is then tested using the M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 § 6.1 . At each resolution, the problem is solved on the Cartesian grid ( ε = 0 ) (\varepsilon= 0) ( ε = 0 ) ( ε = 0.01 ) (\varepsilon= 0.01) ( ε = 0.01 ) Δ p w , r = ∥ p ^ w ( ε ) − p ^ w ( 0 ) ∥ L r / ∥ p ^ w ( 0 ) ∥ L r \Delta p_{w,r} = \lVert\hat{p}_{w}^{(\varepsilon)} - \hat{p}_{w}^{(0)} \rVert_{L_r}/\lVert\hat{p}_{w}^{(0)} \rVert_{L_r} Δ p w , r = ∥ p ^ w ( ε ) − p ^ w ( 0 ) ∥ L r / ∥ p ^ w ( 0 ) ∥ L r r ∈ { 2 , ∞ } r \in\{2,\infty\} r ∈ { 2 , ∞ } Figure 6 (c) plots both quantities against N η N_\eta N η L ∞ L_\infty L ∞ Δ p w , ∞ \Delta p_{w,\infty} Δ p w , ∞ L 2 L_2 L 2 Δ p w , 2 \Delta p_{w,2} Δ p w , 2 401 × 41 401 \times41 401 × 41 1201 × 121 1201 \times121 1201 × 121
The coordinate-invariance test is repeated for the M a = 4.5 \mathrm{Ma} = 4.5 Ma = 4.5 801 × 101 801 \times101 801 × 101 1201 × 145 1201 \times145 1201 × 145 L ∞ L_\infty L ∞ L 2 L_2 L 2 ξ x y \xi_{xy} ξ x y η x y \eta_{xy} η x y M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5
High-amplitude Newton continuation
The Newton continuation of § 4.3 is verified against the independent HNS implementation of Scholten et al. [32 ] using the nonlinear Blasius benchmark of § 5.1 . Starting from two converged lower-amplitude solutions, the second-order predictor and modified Newton corrector, here with the Jacobian refactorized at every step, advance the disturbance field and inflow amplitude A 0 A_0 A 0 A 0 A_0 A 0 u r m s , ( 1 , 0 ) = A 0 / 2 u_{\mathrm{rms},(1,0)} = A_0/\sqrt{2} u rms , ( 1 , 0 ) = A 0 / 2 u r m s , ( 1 , 0 ) = 4 × 10 − 3 u_{\mathrm{rms},(1,0)} = 4 \times10^{-3} u rms , ( 1 , 0 ) = 4 × 1 0 − 3 A 0 = 5.657 × 10 − 3 A_0 = 5.657 \times10^{-3} A 0 = 5.657 × 1 0 − 3 Figure 7 b).
Figure 7. High-amplitude Newton continuation for the Blasius TS case (M = 4 M = 4 M = 4 N = 0 N = 0 N = 0 N ξ × N η = 1050 × 60 N_{\xi} \times N_{\eta} = 1050 \times60 N ξ × N η = 1050 × 60 u r m s , ( m , 0 ) u_{\mathrm{rms},(m,0)} u rms , ( m , 0 ) u r m s , ( 1 , 0 ) = 4 × 10 − 3 u_{\mathrm{rms},(1,0)} = 4 \times10^{-3} u rms , ( 1 , 0 ) = 4 × 1 0 − 3 [32 ] (circles), and the linear fundamental (dashed line). Colors identify the MFD, the fundamental ( 1 , 0 ) (1, 0) ( 1 , 0 ) ( 2 , 0 ) (2, 0) ( 2 , 0 ) ( 4 , 0 ) (4, 0) ( 4 , 0 ) u r m s , ( 1 , 0 ) u_{\mathrm{rms},(1,0)} u rms , ( 1 , 0 )
Figure 7 (a) shows close agreement with the reference over the full range of the reference data for the MFD, the fundamental, and its three harmonics on a 1050 × 60 1050 \times60 1050 × 60 R ≈ 800 R \approx800 R ≈ 800 R ≈ 875 R \approx875 R ≈ 875 [32 ] fail to converge because of strong nonlinear coupling. DeHNSSo-v2 still converges in this regime, where the fundamental mode reaches u r m s , ( 1 , 0 ) ≈ 0.18 u_{\mathrm{rms},(1,0)} \approx0.18 u rms , ( 1 , 0 ) ≈ 0.18 R = 1100 R = 1100 R = 1100
Figure 7 (b) compares Picard iteration, supplemented by Anderson acceleration [45 ] , with Newton continuation along the same branch. Picard uses successive fixed-amplitude solves and converges up to an inflow amplitude u r m s , ( 1 , 0 ) ≈ 2.9 × 10 − 3 u_{\mathrm{rms},(1,0)} \approx2.9 \times10^{-3} u rms , ( 1 , 0 ) ≈ 2.9 × 1 0 − 3 u r m s , ( 1 , 0 ) = 4 × 10 − 3 u_{\mathrm{rms},(1,0)} = 4 \times10^{-3} u rms , ( 1 , 0 ) = 4 × 1 0 − 3 § 7 ).
Although Chebyshev collocation closely reproduces the modal evolution with N η = 40 N_\eta= 40 N η = 40 N η = 60 N_\eta= 60 N η = 60 Figure 3 ). Evaluated at zero disturbance on the 1050 × 40 1050 \times40 1050 × 40 1050 × 60 1050 \times60 1050 × 60
Validation
Four cases are used to validate the principal capabilities of DeHNSSo-v2 against independent reference data. The compressible formulation is assessed for linear and nonlinear disturbances in flat-plate boundary layers (Secs. 6.1 and 6.2), the curvilinear formulation for nonlinear stationary crossflow instability over a swept wing (Sec. 6.3), and the admittance wall condition for a TS wave interacting with a wall-embedded Helmholtz resonator (Sec. 6.4). Recovery of the incompressible limit is verified separately in Sec. 5.1.
Linear compressible flat-plate flow
The linear compressible operator, compressible local-stability inflow, and far-field treatment are validated against the PSE results of Zanus [42 ] and an in-house PSE solver validated by Barahona et al. [41 ] . Three adiabatic flat-plate cases span the subsonic, supersonic, and hypersonic regimes: M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 1.6 1.6 1.6 F = 40 × 10 − 6 F = 40 \times10^{-6} F = 40 × 1 0 − 6 M a = 4.5 \mathrm{Ma} = 4.5 Ma = 4.5 F = 153 × 10 − 6 F = 153 \times10^{-6} F = 153 × 1 0 − 6 ( M = 1 , N = 0 ) (M = 1, N = 0) ( M = 1 , N = 0 ) A 0 = 2 × 10 − 6 A_{0} = 2 \times10^{-6} A 0 = 2 × 1 0 − 6 800 × 81 800 \times81 800 × 81 M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 1.6 1.6 1.6 800 × 101 800 \times101 800 × 101 M a = 4.5 \mathrm{Ma} = 4.5 Ma = 4.5
Figure 8 (a) compares the spatial growth rate of the ( 1 , 0 ) (1, 0) ( 1 , 0 )
Figure 8. Linear compressible flat-plate validation. (a) Spatial growth rate − α i -\alpha_i − α i ( ρ u ) ′ (\rho u)' ( ρ u ) ′ M a = 0.5 \mathrm{Ma}=0.5 Ma = 0.5 1.6 1.6 1.6 u ′ u' u ′ M a = 4.5 \mathrm{Ma}=4.5 Ma = 4.5 R = R e x R=\sqrt{\mathrm{Re}_x} R = Re x M a = 0.5 \mathrm{Ma}=0.5 Ma = 0.5 1.6 1.6 1.6 F = 40 × 10 − 6 F=40\times10^{-6} F = 40 × 1 0 − 6 M a = 4.5 \mathrm{Ma}=4.5 Ma = 4.5 F = 153 × 10 − 6 F=153\times10^{-6} F = 153 × 1 0 − 6 [41 ] (filled diamonds), and the PSE data of Zanus [42 ] (circles). (b)–(d) Wall-normal profiles of the streamwise-velocity, temperature, and density perturbations for the second Mack mode at M a = 4.5 \mathrm{Ma}=4.5 Ma = 4.5 R ≈ 1486 R\approx1486 R ≈ 1486 ∣ u ^ ∣ max |\hat{u}|_{\max} ∣ u ^ ∣ m a x y / δ 99 y/\delta_{99} y / δ 99 [41 ] (black diamonds).
− α i = 1 2 ∣ q ′ ∣ max d ∣ q ′ ∣ max d R , (44) -\alpha_i = \frac{1}{2\lvert q'\rvert_{\max}} \frac{\mathrm{d}\lvert q'\rvert_{\max}}{\mathrm{d}R}, \tag*{(44)} − α i = 2 ∣ q ′ ∣ m a x 1 d R d ∣ q ′ ∣ m a x , (44) where ∣ q ′ ∣ max \lvert q'\rvert_{\max} ∣ q ′ ∣ m a x [42 ] , the streamwise mass flux ( ρ u ) ′ = ρ ‾ u ′ + u ‾ ρ ′ (\rho u)' = \overline{\rho}u' + \overline{u}\rho' ( ρ u ) ′ = ρ u ′ + u ρ ′ M a = 0.5 \mathrm{Ma} = 0.5 Ma = 0.5 1.6 1.6 1.6 u ′ u' u ′ M a = 4.5 \mathrm{Ma} = 4.5 Ma = 4.5 M a = 4.5 \mathrm{Ma} = 4.5 Ma = 4.5 Figure 8 (b)–(d), DeHNSSo-v2 also reproduces the wall-normal structure of the second Mack mode near its maximum growth rate at R ≈ 1486 R \approx1486 R ≈ 1486 [42 ] did not report these profiles, they are compared with the in-house PSE solver only. These comparisons validate both the streamwise amplification and the wall-normal structure of the compressible linear perturbations.
Nonlinear compressible oblique-mode breakdown
The compressible nonlinear forcing and its harmonic balancing are validated against the nonlinear PSE calculation of Chang and Malik [19 ] for oblique-mode breakdown over an adiabatic flat plate at M a = 1.6 \mathrm{Ma} = 1.6 Ma = 1.6 T ∞ ∗ = 300 K T_{\infty}^{*} = 300\ \mathrm{K} T ∞ ∗ = 300 K P r = 0.72 \mathrm{Pr} = 0.72 Pr = 0.72 ( 1 , ± 1 ) (1, \pm1) ( 1 , ± 1 ) R = R e x = 750 R = \sqrt{Re_x} = 750 R = R e x = 750 F = 2 × 10 − 5 F = 2 \times10^{-5} F = 2 × 1 0 − 5 β / R = ± 0.83 × 10 − 4 \beta/R = \pm0.83 \times10^{-4} β / R = ± 0.83 × 1 0 − 4 60 ∘ 60^\circ 6 0 ∘ A 0 = 8 × 10 − 4 A_{0} = 8 \times10^{-4} A 0 = 8 × 1 0 − 4 R = 1550 R = 1550 R = 1550 (35) begins at ξ b \xi_b ξ b R = 1500 R = 1500 R = 1500 1500 × 60 1500 \times60 1500 × 60 M = N = 4 M = N = 4 M = N = 4 θ = 0.4 \theta= 0.4 θ = 0.4 σ m , n < 7 × 10 − 4 \sigma_{m,n} < 7 \times10^{-4} σ m , n < 7 × 1 0 − 4
The pair of oblique modes is generated by the linear PSE inflow of Sec. 2.7, which develops both waves over a short upstream region and thereby avoids the transient that a locally parallel eigenfunction would introduce. The linear compressible operator was verified at these conditions beforehand: the N N N F = 2 × 10 − 5 F = 2 \times10^{-5} F = 2 × 1 0 − 5 β / R = 0.96 × 10 − 4 \beta/R = 0.96 \times10^{-4} β / R = 0.96 × 1 0 − 4 F = 4 × 10 − 5 F = 4 \times10^{-5} F = 4 × 1 0 − 5 [19 ] (not shown). All three modes are set to zero at the inflow, so the MFD and the higher harmonics arise from the nonlinear coupling alone.
Figure 9 shows the wall-normal maximum of the streamwise-velocity amplitude of every retained harmonic against R R R [19 ] . DeHNSSo-v2 reproduces the staggered harmonic cascade from the growing oblique pair through the MFD and the stationary streamwise streak ( 0 , 2 ) (0,2) ( 0 , 2 ) ( 0 , 2 ) (0,2) ( 0 , 2 ) ( 1 , 3 ) (1,3) ( 1 , 3 ) ( 3 , 1 ) (3,1) ( 3 , 1 ) ( 3 , 3 ) (3,3) ( 3 , 3 ) ( 2 , 4 ) (2,4) ( 2 , 4 ) R ≈ 1300 R\approx1300 R ≈ 1300 M = N = 4 M=N=4 M = N = 4 M = N = 9 M=N=9 M = N = 9
Figure 10 shows the instantaneous streamwise velocity at the constant scaled height y / l x = 2.3 y/l_x=2.3 y / l x = 2.3 l x l_x l x y = 2.3 R / R 0 y=2.3R/R_0 y = 2.3 R / R 0 [19 ] . As nonlinear interactions generate the stationary streamwise-streak mode ( 0 , 2 ) (0,2) ( 0 , 2 ) M = N = 4 M=N=4 M = N = 4
Figure 9. Nonlinear compressible validation for the oblique-mode breakdown of Chang and Malik [19 ] at M a = 1.6 \mathrm{Ma}=1.6 Ma = 1.6 F = 2 × 10 − 5 F=2\times10^{-5} F = 2 × 1 0 − 5 A 0 = 8 × 10 − 4 A_0=8\times10^{-4} A 0 = 8 × 1 0 − 4 R = R e x R=\sqrt{\mathrm{Re}_x} R = Re x M = N = 4 M=N=4 M = N = 4
Figure 10. Instantaneous streamwise velocity reconstructed from the DeHNSSo-v2 solution (with M = N = 4 M=N=4 M = N = 4 y / l x = 2.3 y/l_x=2.3 y / l x = 2.3 y = 2.3 R / R 0 y=2.3R/R_0 y = 2.3 R / R 0 x = R 2 / R 0 x=R^2/R_0 x = R 2 / R 0
x ≈ 3000 x\approx3000 x ≈ 3000 [19 ] .
Nonlinear crossflow instability over a curved swept wing
The curvilinear formulation of Secs. 2.4–2.5 and the grid-processing pipeline of Sec. 3 are validated through a simulation of nonlinear stationary crossflow instability over the M3J swept wing tested in the TU Delft low-turbulence tunnel [10 ] . DeHNSSo-v2 simulations are tested against the DNS of Moniripiri et al. [43 ] , itself designed to replicate the TU Delft experiment. The airfoil is a modified NACA 66018 at Λ = 45 ∘ \Lambda=45^\circ Λ = 4 5 ∘ 3 ∘ 3^\circ 3 ∘ R e c x = 2.38 × 10 6 \mathrm{Re}_{c_x}=2.38\times10^6 Re c x = 2.38 × 1 0 6 [43 ] reported good DNS agreement with the experimentally measured disturbance growth and transition to turbulence. Figure 11 presents the DeHNSSo-v2 stability grid and imported base flow for the hump configuration.
As shown in Figure 11 (a), the local airfoil surface is inclined by up to 20 ∘ 20^\circ 2 0 ∘ Figure 11 (b) is 1 1 1 1 : 50 1:50 1 : 50 x / c x = 0.15 x/c_x=0.15 x / c x = 0.15 0.155 ≤ x / c x ≤ 0.178 0.155\le x/c_x\le0.178 0.155 ≤ x / c x ≤ 0.178 Figure 11 (c)). The wall of the stability grid coincides with the DNS wall to within interpolation error. The wall-normal domain extends to η = 20 \eta= 20 η = 20 η = 60 \eta= 60 η = 60 δ ˉ r e f ∗ = 2.687 × 10 − 4 m \bar{\delta}^{*}_{\mathrm{ref}} = 2.687 \times10^{-4}\,\mathrm{m} δ ˉ ref ∗ = 2.687 × 1 0 − 4 m Q r e f = 23.496 m / s Q_{\mathrm{ref}} = 23.496\,\mathrm{m/s} Q ref = 23.496 m/s x / c x = 0.055 x/c_x = 0.055 x / c x = 0.055 R e 0 = Q r e f δ ˉ r e f ∗ / ν ∞ = 451 \mathrm{Re}_0 = Q_{\mathrm{ref}}\bar{\delta}^{*}_{\mathrm{ref}}/\nu_{\infty} = 451 Re 0 = Q ref δ ˉ ref ∗ / ν ∞ = 451
Figure 11. Stability domain, grid, and base flow for the hump configuration A1-H of the M3J swept wing of Moniripiri et al. [43 ] . (a) Body-fitted grid in the leading-edge-normal frame at true aspect ratio, with coordinates normalized by δ ˉ r e f ∗ \bar{\delta}^{*}_{\mathrm{ref}} δ ˉ ref ∗ η = 60 \eta= 60 η = 60 y ~ = ( y − y w , c ) / δ ˉ r e f ∗ \tilde{y} = (y-y_{w,c})/\bar{\delta}^{*}_{\mathrm{ref}} y ~ = ( y − y w , c ) / δ ˉ ref ∗ y w , c y_{w,c} y w , c x x x δ 99 \delta_{99} δ 99 − W c f / Q r e f -W_{cf}/Q_{\mathrm{ref}} − W c f / Q ref W c f = 0 W_{cf}=0 W c f = 0
In DeHNSSo-v2, only the stationary fundamental ( 0 , 1 ) (0,1) ( 0 , 1 ) β 1 = 0.225 \beta_1 = 0.225 β 1 = 0.225 λ z = 7.5 m m \lambda_z = 7.5\,\mathrm{mm} λ z = 7.5 mm A 0 = 2.62 × 10 − 3 A_0 = 2.62 \times10^{-3} A 0 = 2.62 × 1 0 − 3 10 − 1 10^{-1} 1 0 − 1 M = 0 M = 0 M = 0 N = 4 N = 4 N = 4 2428 × 60 2428 \times60 2428 × 60 Δ ξ ≈ 0.51 δ ˉ r e f ∗ \Delta\xi\approx0.51\bar{\delta}^{*}_{\mathrm{ref}} Δ ξ ≈ 0.51 δ ˉ ref ∗ ( 0 , 1 ) (0,1) ( 0 , 1 )
Figure 12(a,b) compares the chordwise development of the amplitude of every retained stationary mode, the MFD ( 0 , 0 ) (0,0) ( 0 , 0 ) ( 0 , 1 ) (0,1) ( 0 , 1 ) ( 0 , 4 ) (0,4) ( 0 , 4 ) [43 ] . On A1-C, all harmonics and the MFD follow the reference across the domain. On A1-H, the solver reproduces the fundamental mode stalling in the crossflow-reversal region (see Figure 12 b) while growing less than on A1-C beyond x / c x ≈ 0.2 x/c_x \approx0.2 x / c x ≈ 0.2 x / c x ≈ 0.185 x/c_x \approx0.185 x / c x ≈ 0.185 x / c x ≈ 0.35 x/c_x \approx0.35 x / c x ≈ 0.35 A ≈ 22 % A \approx22\% A ≈ 22% ( 0 , 3 ) (0,3) ( 0 , 3 ) ( 0 , 4 ) (0,4) ( 0 , 4 ) ( 0 , 1 ) (0,1) ( 0 , 1 )
Figure 12. Nonlinear stationary crossflow over the clean (A1-C) and hump (A1-H) configurations of the M3J swept wing of Moniripiri et al. [43 ] , retaining M = 0 M = 0 M = 0 N = 4 N = 4 N = 4 n n n u ~ ξ \tilde{u}_{\xi} u ~ ξ A = max η 2 ∣ u ~ ξ , ( 0 , n ) ∣ / Q r e f A = \max_{\eta} 2|\tilde{u}_{\xi,(0,n)}|/Q_{\mathrm{ref}} A = max η 2∣ u ~ ξ , ( 0 , n ) ∣/ Q ref x / c x = 0.36 x/c_x = 0.36 x / c x = 0.36 I P I_P I P [43 ] and normalized by Q r e f 3 Q_{\mathrm{ref}}^3 Q ref 3 I P = 0 I_P = 0 I P = 0 W c f = 0 W_{cf} = 0 W c f = 0
Figure 12(c–f) indirectly compares the shape of the perturbation rather than its amplitude: the integrand of the Reynolds–Orr production of the fundamental mode (see Moniripiri et al. [43 ] ), whose negative values indicate the reverse lift-up effect. On A1-C, kinetic energy is transferred from the base flow to the perturbation field. By contrast, a region of negative production appears near the wall downstream of the reversal region in A1-H, responsible for the reduced disturbance growth downstream of the hump. DeHNSSo-v2 accurately reproduces this near-wall region at the same location and with the same extent.
The A1-C and A1-H DeHNSSo-v2 calculations required 15.2 h and 11.5 h on system E (Appendix D), with four and two Jacobian factorizations, respectively. Each calculation used approximately 1.5 × 10 5 1.5 \times10^{5} 1.5 × 1 0 5 3 × 10 2 3 \times10^{2} 3 × 1 0 2 1.0 × 10 8 1.0 \times10^{8} 1.0 × 1 0 8 4 × 10 5 4 \times10^{5} 4 × 1 0 5 [43 ] .
Tollmien–Schlichting waves over a Helmholtz resonator
A passive wall-embedded Helmholtz resonator interacting with a linear TS wave is validated following the flow case investigated by Michelis et al. [12 ] , whose linearized Navier–Stokes computations were validated against wind-tunnel particle image velocimetry measurements. A linear two-dimensional TS wave develops in a Blasius boundary layer at M a = 0.06 \mathrm{Ma} = 0.06 Ma = 0.06 R e 0 = 371 \mathrm{Re}_0 = 371 Re 0 = 371 U ˉ r e f = 20.9 m / s \bar{U}_{\mathrm{ref}} = 20.9\,\mathrm{m/s} U ˉ ref = 20.9 m/s l ˉ r e f = 2.70 × 10 − 4 m \bar{l}_{\mathrm{ref}} = 2.70 \times10^{-4}\,\mathrm{m} l ˉ ref = 2.70 × 1 0 − 4 m x r = 0.55 m x_r = 0.55\,\mathrm{m} x r = 0.55 m f r = 295 H z f_r = 295\,\mathrm{Hz} f r = 295 Hz F = 64.4 × 10 − 6 F = 64.4 \times10^{-6} F = 64.4 × 1 0 − 6 Δ N \Delta N Δ N
Following Michelis et al. [12 ] , the resonator is represented by the complex wall admittance H B C H_{\mathrm{BC}} H BC v ^ m , n ∣ η = 0 = 0 \hat{v}_{m,n}|_{\eta=0} = 0 v ^ m , n ∣ η = 0 = 0 § 2.7 is therefore replaced by the wall row v ^ m , n − H B C p ^ m , n = 0 \hat{v}_{m,n} - H_{\mathrm{BC}}\hat{p}_{m,n} = 0 v ^ m , n − H BC p ^ m , n = 0 (25) , while the no-slip and isothermal conditions are retained. Unlike the prescribed wall forcing of Westerbeek et al. [14 ] , this condition specifies neither variable independently.
Figure 13 (a,Figure 13 b) compares the N N N N u ′ = ln ( A u ′ / A u ′ , I ) N_{u'} = \ln(A_{u'}/A_{u',I}) N u ′ = ln ( A u ′ / A u ′ , I ) N v ′ = ln ( A v ′ / A v ′ , I ) N_{v'} = \ln(A_{v'}/A_{v',I}) N v ′ = ln ( A v ′ / A v ′ , I ) A A A A I A_I A I Δ N ≈ 0.35 \Delta N \approx0.35 Δ N ≈ 0.35 Figure 13 (c) shows the wall-normal velocity perturbation of case R near the neck, where the resonator drives v ′ v' v ′ − 6.6 -6.6 − 6.6 δ r ∗ \delta_r^* δ r ∗ Figure 13 (d) compares the transmission coefficient ∣ T ∣ = exp [ Δ N ] |\mathcal{T}| = \exp[\Delta N] ∣ T ∣ = exp [ Δ N ] f T S / f r f_{\mathrm{TS}}/f_r f TS / f r
Figure 13. Tollmien–Schlichting wave interacting with a wall-embedded Helmholtz resonator at M a = 0.06 Ma = 0.06 M a = 0.06 R e 0 = 371 Re_0 = 371 R e 0 = 371 f T S = f r = 295 f_{\mathrm{TS}} = f_r = 295 f TS = f r = 295 F = 64.4 × 10 − 6 F = 64.4 \times10^{-6} F = 64.4 × 1 0 − 6 N N N − 1 -1 − 1 δ 99 \delta_{99} δ 99 ∣ T ∣ = e Δ N |T| = e^{\Delta N} ∣ T ∣ = e Δ N [12 ] ; C and R denote the baseline and resonator cases, respectively, and the two velocity components coincide in (d). In (a–c), the abscissa is the distance from the neck, x − x r x-x_r x − x r δ r ∗ \delta_r^* δ r ∗ x = x r x=x_r x = x r
Computational performance
The efficiency developments of § 4 are first evaluated using the nonlinear Tollmien–Schlichting case, verified in § 5 . On the same grid and hardware, DeHNSSo-v2 reaches convergence 18.9 times faster than DeHNSSo-v1. § 7.1.1 investigates the runtime distribution, LHS assembly, and factorization reuse. § 7.1.2 compares the time–memory trade-offs of the direct-solver back-ends, and § 7.1.3 examines scaling with resolution, hardware, and parallel workers. Finally, § 7.2 assesses the cost of Newton continuation. The hardware is described in Appendix D.
Linear-algebra performance
Runtime analysis
Table 4 compares the wall-clock time and memory of DeHNSSo-v1 and DeHNSSo-v2 for the nonlinear TS case, computed with the same setup, harmonic truncation, grid, and hardware. Both solver versions converge the nonlinear terms with the Picard iteration of § 2.6 , so that the comparison isolates the linear algebra. Figure 14 (a) shows how the runtime of each version is divided among its tasks. In DeHNSSo-v1, the repeated factorize-and-solve operations dominate, as noted in § 4 . In DeHNSSo-v2, factorization still takes the largest share, 52%, but it is performed only once per mode, when the mode enters the truncation, so that factorization takes about 21 s, against 691 s for the combined factorize-and-solve operations of DeHNSSo-v1. The two developments behind this reduction, linear-cost assembly and factorization reuse, are examined in turn.
Table 4. Wall-clock time and memory of DeHNSSo-v1 [14 ] and DeHNSSo-v2 for the nonlinear incompressible Blasius TS case (M = 5 M = 5 M = 5 N ξ = 800 N_{\xi} = 800 N ξ = 800 N η = 40 N_{\eta} = 40 N η = 40
The first development is the O ( N ξ ) O(N_{\xi}) O ( N ξ ) § 4.1 . Figure 14 (b) isolates the assembly over N ξ = 200 N_{\xi} = 200 N ξ = 200 3200 3200 3200 M = 1 M = 1 M = 1 O ( N ξ ) O(N_{\xi}) O ( N ξ ) O ( N ξ 1.5 ) O(N_{\xi}^{1.5}) O ( N ξ 1.5 ) N ξ = 3200 N_{\xi} = 3200 N ξ = 3200
The second development, the factorization reuse of § 4.2.1 , provides the larger gain. DeHNSSo-v1 refactorizes every active per-mode system at each nonlinear iteration, whereas DeHNSSo-v2 factorizes each mode once and reduces every later solve to a back-solve with the cached factor, which is 50 times faster (Table 4 ). Figure 14 (c) shows when these factorizations occur: in the first two iterations for modes ( 0 , 0 ) (0,0) ( 0 , 0 ) ( 2 , 0 ) (2,0) ( 2 , 0 ) ( 3 , 0 ) (3,0) ( 3 , 0 ) ( 5 , 0 ) (5,0) ( 5 , 0 ) Table 4 ). For two stationary-crossflow cases of Westerbeek et al. [14 ] , the speed-ups are 9.3 for nonlinear development and 2.0 for linear development over a hump, where the single solve precludes factorization reuse.
This gain requires the factors to remain in memory. As listed in Table 4 , the operator occupies only 112 MB, because its mode-invariant blocks are stored once and the MFD is retained as a 2 MB correction, whereas the six cached LU factors occupy 6.2 GB. Factorization reuse therefore saves time at the cost of memory, so that storing the factors, rather than computing them, becomes the main constraint.
The cache can also persist between solver calls through Opt.rerun (§ 4.2.1 ). When only the inflow amplitude A 0 A_0 A 0 Table 4 repeats in 15.3 s instead of 39.9 s, as it skips operator assembly and all six factorizations.
Direct-solver back-ends
Figure 15 compares the cached sparse-LU, banded-LU, and block-Thomas back-ends of § 4.2.1 –§ 4.2.3 by their factorization time and resident factor storage. The baseline case uses the 800 × 40 800 \times40 800 × 40 § 7.1.1 and one thread of system D, together with an 800 × 80 800 \times80 800 × 80 M a = 10 − 3 \mathrm{Ma}=10^{-3} Ma = 1 0 − 3 § 5.1 ) are included, so that the cost of the additional variable can be measured in the incompressible limit.
Figure 15. Resident factor storage versus factorization time for the direct-solver back-ends on one thread of system D (Appendix D). Marker shape denotes the back-end, color and label the operator and grid (I, incompressible; C, compressible), and fill the wall-normal discretization (filled, Chebyshev; open, fourth-order finite differences). For block Thomas, only the resident portion is shown; the remainder is streamed to disk.
On the baseline grid, banded LU exchanges memory for speed (see § 4.2.2 ). For the incompressible operator, banded LU factorizes 1.4 times faster than sparse LU, with factors 2.6 times larger. For the compressible operator, the factorization is 4.9 times faster and the factors are only 1.3 times larger. The estimates of § 4.2.2 predict that the fifth variable increases the banded storage and factorization time by ( 5 / 4 ) 2 = 1.56 (5/4)^2=1.56 ( 5/4 ) 2 = 1.56 ( 5 / 4 ) 3 = 1.95 (5/4)^3=1.95 ( 5/4 ) 3 = 1.95
Block Thomas, which is available only for the compressible operator and provides an out-of-core alternative (§ 4.2.3 ), keeps a resident factor of approximately one fifth of the sparse-LU factor and one seventh of the banded-LU factor on the baseline grid, with a factorization time between the two.
Fourth-order wall-normal finite differences (FD4), shown by the open markers, test the dependence on operator sparsity described in § 4.2.2 and § 4.2.3 . They reduce the operator nonzeros by factors of 5 and 9 on the 800 × 40 800 \times40 800 × 40 800 × 80 800 \times80 800 × 80 800 × 80 800 \times80 800 × 80
For the baseline compressible operator, one solve with a retained factor takes 0.18 s with sparse LU, 0.36 s with banded LU, and 1.00 s with out-of-core block Thomas, which does not alter the ordering of the back-ends for the iteration counts considered here. Banded LU is therefore generally fastest when sufficient memory is available, whereas block Thomas is preferable for memory-limited compressible calculations, particularly
with finite differences. These measurements support the choice of back-end described in § 4.2 .
Scaling with resolution and hardware
This subsection tests whether the conclusions of § 7.1.2 persist for the larger systems required by practical nonlinear-transition calculations. Unless stated otherwise, all tests use the compressible operator at M a = 10 − 3 \mathrm{Ma}=10^{-3} Ma = 1 0 − 3
Figure 16 examines scaling with both grid directions. In panels (a,c), increasing N ξ N_{\xi} N ξ N ξ N_{\xi} N ξ N η N_{\eta} N η N η 2 N_{\eta}^{2} N η 2 O ( N η 3 ) O(N_{\eta}^{3}) O ( N η 3 )
Figure 16. Resolution scaling for the compressible operator at Ma = 10 − 3 10^{-3} 1 0 − 3 N ξ N_{\xi} N ξ N η = 80 N_{\eta}=80 N η = 80 N η N_{\eta} N η N ξ = 800 N_{\xi}=800 N ξ = 800 N ξ N_{\xi} N ξ N η N_{\eta} N η 4 _4 4
In the FD4 _4 4 N ξ 0.5 N_{\xi}^{0.5} N ξ 0.5 N η 1.5 N_{\eta}^{1.5} N η 1.5
Hardware affects the back-ends according to their parallelism and memory demand. Sparse LU is effectively serial, whereas banded LU and block Thomas benefit from additional threads, so that at eight threads banded LU is 2.6 times faster on system E than on system D. Near the memory limit of the 32 GB system D, however, a solve with the 17.9 GB banded factor takes 7.5 s, against 1.1 s on system E, so that timings obtained near the memory limit are not transferable between machines.
Mode-level parallelism through Opt.parfor (Sec. 4.2.1) distributes the back-solves over parallel workers. For the incompressible case of Sec. 7.1.1 on system D, six workers accelerate the back-solve loop by 3.0 times but reduce the total runtime by only 1.22 times, because the loop accounts for about 29% of the calculation and parallelism saturates at the number of active modes. Factorization reuse, which is also available in serial, therefore provides the larger reduction in runtime.
Nonlinear-solver performance
The Picard iteration of Sec. 2.6 solves each harmonic mode separately, whereas Newton continuation includes their nonlinear coupling in the Jacobian of § 4.3 . This extends the convergence range, as verified in § 5.4 , at the cost of a larger factorization. For the Blasius TS continuation of Figure 7 (M = 4 M = 4 M = 4 1050 × 60 1050 \times60 1050 × 60 u r m s , ( 1 , 0 ) = 4 × 10 − 3 u_{\mathrm{rms},(1,0)} = 4 \times10^{-3} u rms , ( 1 , 0 ) = 4 × 1 0 − 3
Figure 17 shows how this cost grows with the harmonic truncation for both operators on the same 1050 × 60 1050 \times60 1050 × 60 ( 0 , 0 ) (0,0) ( 0 , 0 ) ( M , 0 ) (M,0) ( M , 0 ) M = 4 M = 4 M = 4 M a = 10 − 3 \mathrm{Ma} = 10^{-3} Ma = 1 0 − 3 M = 2 M = 2 M = 2 4 4 4 ( M + 1 ) 2 (M+1)^2 ( M + 1 ) 2 § 7.1.2 , with 3.2–3.3 times the factor storage and 29–42 times the factorization time. Since the in-core and out-of-core factorizations of the compressible M = 2 M = 2 M = 2
Figure 17. Newton Jacobian (a) factorization time and (b) factor size on disk versus harmonic truncation M M M M = 4 M=4 M = 4 ( M + 1 ) 2 (M+1)^2 ( M + 1 ) 2
The modified Newton corrector of § 4.3 avoids repeated factorizations by retaining the Jacobian factor across continuation steps. Over the first ten steps of the Blasius continuation, taken with identical arclength increments, reusing the initial factor throughout reduces the runtime from 299 min, with a factorization at every step, to 57 min, a factor of 5.2; the default reuse limit of eight steps gives 82 min. With reuse, the corrector count rises from two to five iterations per step, whereas it remains at two without reuse. The final inflow amplitudes u r m s , ( 1 , 0 ) u_{\mathrm{rms},(1,0)} u rms , ( 1 , 0 )
The swept-wing cases of § 6.3 show the cost of repeated solves on a larger grid. The N = 4 N = 4 N = 4 2428 × 60 2428 \times60 2428 × 60
Conclusions
We have presented DeHNSSo-v2, a new version of the open-source Delft Harmonic Navier–Stokes Solver [14 ] . The present release combines major updates to the governing formulation, grid preprocessing, nonlinear solution, and linear algebra. These additions extend the scope of DeHNSSo-v1 to general two-dimensional body-fitted geometries, compressible flows, and higher disturbance amplitudes, targeting boundary-layer stability and flow-control problems in which nonparallel development, nonlinear interactions, or localized geometry require global streamwise coupling.
A general curvilinear formulation and an automatic grid-generation and preprocessing pipeline replace the problem-specific grid construction of DeHNSSo-v1 and support convex and concave wall geometries. A compressible formulation in the primitive perturbations (u ′ , v ′ , w ′ , ρ ′ , T ′ u',v',w',\rho',T' u ′ , v ′ , w ′ , ρ ′ , T ′ O ( N ξ ) O(N_{\xi}) O ( N ξ ) N ξ 1.5 N_{\xi}^{1.5} N ξ 1.5
The curvilinear formulation preserves the fourth-order accuracy of the discretization on strongly distorted grids, and the compressible solver recovers the established incompressible solution in the low-Mach limit. Published comparisons span linear compressible instabilities up to M a = 4.5 \mathrm{Ma} = 4.5 Ma = 4.5 M a = 1.6 \mathrm{Ma} = 1.6 Ma = 1.6 § 4 , including their growth with the number of variables: banded LU is fastest when sufficient memory is available, whereas block Thomas with finite differences minimizes resident storage and supports out-of-core execution. Newton continuation extends the Picard limit by approximately 40%, reproduces the independent HNS solution of Scholten et al. [32 ] , including the approximately 45-fold amplification of the fundamental Tollmien–Schlichting wave, and provides the nonlinear crossflow solutions compared with direct numerical simulation. The coupled Newton Jacobian, particularly in compressible flows, remains the principal cost.
DeHNSSo-v2 provides a flexible software base for further development. Its present scope remains restricted to steady, spanwise-invariant base flows and, in the compressible formulation, to a calorically perfect ideal gas at constant Prandtl number. Future developments will target: (i) discrete adjoint and resolvent operators for sensitivity and receptivity analysis [67 ] , (ii) an HNS–NPSE–LPSE coupling for efficient treatment of nonparallel effects [16 ] , (iii) nonideal and high-temperature gas models [68 , 69 ] , and (iv) adapting scalable, memory-efficient iterative methods to three-dimensional HNS [70 ] .
The open-source version of DeHNSSo-v2 is available at the following link: https://github.com/pcboldini/DeHNSSo.
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
The European Research Council (grant no. ERC-2023-CoG-101125132, MetaWing) has supported this study. P. C. Boldini acknowledges Dr. M. Moniripiri (KTH Royal Institute of Technology) for providing the swept-wing DNS data. The authors thank the Flow Stability and Flow Control group at TU Delft for testing and debugging DeHNSSo-v2.
Compressible LHS coefficient matrices
The 5 × 5 5 \times5 5 × 5 (25) are listed below. Rows correspond to the streamwise, wall-normal, and spanwise momentum equations, the continuity equation, and the energy equation; columns to the state components (u ′ u^\prime u ′ v ′ v^\prime v ′ w ′ w^\prime w ′ ρ ′ \rho^\prime ρ ′ T ′ T^\prime T ′ (15) and (19) –(21) are absorbed into the coefficients through the building blocks defined below. Derivatives of the base flow are precomputed physical quantities and are therefore written with physical derivatives. The convention follows Scholten et al. [32 ] . Only nonzero entries are listed.
Notation. The building blocks that collect the inverse-metric factors are defined as follows. The curvilinear advection velocities are
U ξ ≡ u ˉ ξ x + v ˉ ξ y , U η ≡ u ˉ η x + v ˉ η y , (45) U_{\xi} \equiv\bar{u}\xi_x+\bar{v}\xi_y,\qquad U_{\eta}\equiv\bar{u}\eta_x+\bar{v}\eta_y, \tag*{(45)} U ξ ≡ u ˉ ξ x + v ˉ ξ y , U η ≡ u ˉ η x + v ˉ η y , (45) and the metric weights that distribute a viscous Laplacian a ∂ x x 2 + b ∂ y y 2 a\partial_{xx}^{2}+b\partial_{yy}^{2} a ∂ xx 2 + b ∂ y y 2
Λ a , b ξ ξ = a ξ x 2 + b ξ y 2 , Λ a , b η η = a η x 2 + b η y 2 , \Lambda_{a,b}^{\xi\xi}=a\xi_x^2+b\xi_y^2,\qquad
\Lambda_{a,b}^{\eta\eta}=a\eta_x^2+b\eta_y^2, Λ a , b ξ ξ = a ξ x 2 + b ξ y 2 , Λ a , b η η = a η x 2 + b η y 2 , Λ a , b ξ η = 2 ( a ξ x η x + b ξ y η y ) , Λ a , b ξ = a ξ x x + b ξ y y , \Lambda_{a,b}^{\xi\eta}=2(a\xi_x\eta_x+b\xi_y\eta_y),\qquad
\Lambda_{a,b}^{\xi}=a\xi_{xx}+b\xi_{yy}, Λ a , b ξ η = 2 ( a ξ x η x + b ξ y η y ) , Λ a , b ξ = a ξ xx + b ξ y y , Λ a , b η = a η x x + b η y y , (46) \Lambda_{a,b}^{\eta}=a\eta_{xx}+b\eta_{yy}, \tag*{(46)} Λ a , b η = a η xx + b η y y , (46) such that a ∂ x x 2 + b ∂ y y 2 = Λ a , b ξ ξ ∂ ξ ξ 2 + Λ a , b η η ∂ η η 2 + Λ a , b ξ η ∂ ξ η 2 + Λ a , b ξ ∂ ξ + Λ a , b η ∂ η a\partial_{xx}^{2}+b\partial_{yy}^{2}=\Lambda_{a,b}^{\xi\xi}\partial_{\xi\xi}^{2}+\Lambda_{a,b}^{\eta\eta}\partial_{\eta\eta}^{2}+\Lambda_{a,b}^{\xi\eta}\partial_{\xi\eta}^{2}+\Lambda_{a,b}^{\xi}\partial_{\xi}+\Lambda_{a,b}^{\eta}\partial_{\eta} a ∂ xx 2 + b ∂ y y 2 = Λ a , b ξ ξ ∂ ξ ξ 2 + Λ a , b η η ∂ η η 2 + Λ a , b ξ η ∂ ξ η 2 + Λ a , b ξ ∂ ξ + Λ a , b η ∂ η ∂ x = ξ x ∂ ξ + η x ∂ η \partial_x=\xi_x\partial_\xi+\eta_x\partial_\eta ∂ x = ξ x ∂ ξ + η x ∂ η ∂ y = ξ y ∂ ξ + η y ∂ η \partial_y=\xi_y\partial_\xi+\eta_y\partial_\eta ∂ y = ξ y ∂ ξ + η y ∂ η μ ˉ \bar{\mu} μ ˉ
S ˉ x x = 4 3 ∂ x u ˉ − 2 3 ∂ y v ˉ , S ˉ y y = 4 3 ∂ y v ˉ − 2 3 ∂ x u ˉ , \bar{S}_{xx}=\frac{4}{3}\partial_x\bar{u}-\frac{2}{3}\partial_y\bar{v},\qquad
\bar{S}_{yy}=\frac{4}{3}\partial_y\bar{v}-\frac{2}{3}\partial_x\bar{u}, S ˉ xx = 3 4 ∂ x u ˉ − 3 2 ∂ y v ˉ , S ˉ y y = 3 4 ∂ y v ˉ − 3 2 ∂ x u ˉ , S ˉ x y = ∂ y u ˉ + ∂ x v ˉ , S ˉ z z = − 2 3 ( ∂ x u ˉ + ∂ y v ˉ ) , \bar{S}_{xy}=\partial_y\bar{u}+\partial_x\bar{v},\qquad
\bar{S}_{zz}=-\frac{2}{3}(\partial_x\bar{u}+\partial_y\bar{v}), S ˉ x y = ∂ y u ˉ + ∂ x v ˉ , S ˉ z z = − 3 2 ( ∂ x u ˉ + ∂ y v ˉ ) , S ˉ x z = ∂ x w ˉ , S ˉ y z = ∂ y w ˉ . (47) \bar{S}_{xz}=\partial_x\bar{w},\qquad
\bar{S}_{yz}=\partial_y\bar{w}. \tag*{(47)} S ˉ x z = ∂ x w ˉ , S ˉ y z = ∂ y w ˉ . (47) Base-flow viscous Laplacians acting on u ˉ \bar{u} u ˉ v ˉ \bar{v} v ˉ w ˉ \bar{w} w ˉ T ˉ \bar{T} T ˉ
L ˉ u = 4 3 ∂ x x 2 u ˉ + ∂ y y 2 u ˉ + 1 3 ∂ x y 2 v ˉ , \bar{L}_{u}=\frac{4}{3}\partial_{xx}^{2}\bar{u}+\partial_{yy}^{2}\bar{u}+\frac{1}{3}\partial_{xy}^{2}\bar{v}, L ˉ u = 3 4 ∂ xx 2 u ˉ + ∂ y y 2 u ˉ + 3 1 ∂ x y 2 v ˉ , L ˉ v = ∂ x x 2 v ˉ + 4 3 ∂ y y 2 v ˉ + 1 3 ∂ x y 2 u ˉ , \bar{L}_{v}=\partial_{xx}^{2}\bar{v}+\frac{4}{3}\partial_{yy}^{2}\bar{v}+\frac{1}{3}\partial_{xy}^{2}\bar{u}, L ˉ v = ∂ xx 2 v ˉ + 3 4 ∂ y y 2 v ˉ + 3 1 ∂ x y 2 u ˉ , L ˉ w = ∂ x x 2 w ˉ + ∂ y y 2 w ˉ , \bar{L}_{w}=\partial_{xx}^{2}\bar{w}+\partial_{yy}^{2}\bar{w}, L ˉ w = ∂ xx 2 w ˉ + ∂ y y 2 w ˉ , L ˉ T = − γ μ ˉ T Pr Re ( ∂ x x 2 T ˉ + ∂ y y 2 T ˉ ) − γ μ ˉ T T Pr Re ∣ ∇ T ˉ ∣ 2 . (48) \bar{L}_{T}=-\frac{\gamma\bar{\mu}_{T}}{\operatorname{Pr}\operatorname{Re}}\left(\partial_{xx}^{2}\bar{T}+\partial_{yy}^{2}\bar{T}\right)-\frac{\gamma\bar{\mu}_{TT}}{\operatorname{Pr}\operatorname{Re}}\lvert\nabla\bar{T}\rvert^{2}. \tag*{(48)} L ˉ T = − Pr Re γ μ ˉ T ( ∂ xx 2 T ˉ + ∂ y y 2 T ˉ ) − Pr Re γ μ ˉ T T ∣ ∇ T ˉ ∣ 2 . (48) Base-flow normalized dissipation Φ ˉ 0 ≡ Φ ˉ / μ ˉ \bar{\Phi}_{0}\equiv\bar{\Phi}/\bar{\mu} Φ ˉ 0 ≡ Φ ˉ / μ ˉ
Φ ˉ 0 = S ˉ x x ∂ x u ˉ + S ˉ y y ∂ y v ˉ + S ˉ x y 2 + S ˉ x z ∂ x w ˉ + S ˉ y z ∂ y w ˉ . (49) \bar{\Phi}_{0}=\bar{S}_{xx}\partial_x\bar{u}+\bar{S}_{yy}\partial_y\bar{v}+\bar{S}_{xy}^{2}+\bar{S}_{xz}\partial_x\bar{w}+\bar{S}_{yz}\partial_y\bar{w}. \tag*{(49)} Φ ˉ 0 = S ˉ xx ∂ x u ˉ + S ˉ y y ∂ y v ˉ + S ˉ x y 2 + S ˉ x z ∂ x w ˉ + S ˉ y z ∂ y w ˉ . (49) The base-flow advective operator and divergence:
u ˉ ⋅ ∇ q ˉ i ≡ u ˉ ∂ x q ˉ i + v ˉ ∂ y q ˉ i , D ˉ u ≡ ∂ x u ˉ + ∂ y v ˉ . (50) \bar{\mathbf{u}}\cdot\nabla\bar{q}_i\equiv\bar{u}\partial_x\bar{q}_i+\bar{v}\partial_y\bar{q}_i,\qquad
\bar{D}_{u}\equiv\partial_x\bar{u}+\partial_y\bar{v}. \tag*{(50)} u ˉ ⋅ ∇ q ˉ i ≡ u ˉ ∂ x q ˉ i + v ˉ ∂ y q ˉ i , D ˉ u ≡ ∂ x u ˉ + ∂ y v ˉ . (50) Dissipation prefactor:
c Φ ≡ 2 γ ( γ − 1 ) Ma 2 Re . (51) c_{\Phi}\equiv\frac{2\gamma(\gamma-1)\operatorname{Ma}^{2}}{\operatorname{Re}}. \tag*{(51)} c Φ ≡ Re 2 γ ( γ − 1 ) Ma 2 . (51) The shorthand κ ˉ = μ ˉ \bar{\kappa}=\bar{\mu} κ ˉ = μ ˉ κ ˉ T = μ ˉ T \bar{\kappa}_{T}=\bar{\mu}_{T} κ ˉ T = μ ˉ T κ ˉ T T = μ ˉ T T \bar{\kappa}_{TT}=\bar{\mu}_{TT} κ ˉ T T = μ ˉ T T § 2.2 .
Viscous-diffusion matrices F \mathbf{F} F C \mathbf{C} C E \mathbf{E} E (coefficients of ∂ ξ ξ 2 q ^ \partial_{\xi\xi}^{2}\hat{\mathbf{q}} ∂ ξ ξ 2 q ^ ∂ η η 2 q ^ \partial_{\eta\eta}^{2}\hat{\mathbf{q}} ∂ η η 2 q ^ ∂ ξ η 2 q ^ \partial_{\xi\eta}^{2}\hat{\mathbf{q}} ∂ ξ η 2 q ^ Writing M \mathbf{M} M F , C , E \mathbf{F},\mathbf{C},\mathbf{E} F , C , E Λ ∘ \Lambda^{\circ} Λ ∘ Λ ξ ξ \Lambda^{\xi\xi} Λ ξ ξ F \mathbf{F} F Λ η η \Lambda^{\eta\eta} Λ η η C \mathbf{C} C Λ ξ η \Lambda^{\xi\eta} Λ ξ η E \mathbf{E} E
M 1 , 1 = − μ ˉ Re Λ 4 / 3 , 1 ∘ , M 2 , 2 = − μ ˉ Re Λ 1 , 4 / 3 ∘ , M_{1,1}=-\frac{\bar{\mu}}{\operatorname{Re}}\Lambda^{\circ}_{4/3,1},\qquad
M_{2,2}=-\frac{\bar{\mu}}{\operatorname{Re}}\Lambda^{\circ}_{1,4/3}, M 1 , 1 = − Re μ ˉ Λ 4/3 , 1 ∘ , M 2 , 2 = − Re μ ˉ Λ 1 , 4/3 ∘ , M 3 , 3 = − μ ˉ Re Λ 1 , 1 ∘ , M 5 , 5 = − γ μ ˉ Pr Re Λ 1 , 1 ∘ . (52) M_{3,3}=-\frac{\bar{\mu}}{\operatorname{Re}}\Lambda^{\circ}_{1,1},\qquad
M_{5,5}=-\frac{\gamma\bar{\mu}}{\operatorname{Pr}\operatorname{Re}}\Lambda^{\circ}_{1,1}. \tag*{(52)} M 3 , 3 = − Re μ ˉ Λ 1 , 1 ∘ , M 5 , 5 = − Pr Re γ μ ˉ Λ 1 , 1 ∘ . (52) together with the Stokes cross-couplings
F 1 , 2 = F 2 , 1 = − μ ˉ 3 Re ξ x ξ y , C 1 , 2 = C 2 , 1 = − μ ˉ 3 Re η x η y , F_{1,2}=F_{2,1}=-\frac{\bar{\mu}}{3\operatorname{Re}}\xi_x\xi_y,\qquad
C_{1,2}=C_{2,1}=-\frac{\bar{\mu}}{3\operatorname{Re}}\eta_x\eta_y, F 1 , 2 = F 2 , 1 = − 3 Re μ ˉ ξ x ξ y , C 1 , 2 = C 2 , 1 = − 3 Re μ ˉ η x η y , E 1 , 2 = E 2 , 1 = − μ ˉ 3 Re ( ξ x η y + ξ y η x ) . (53) E_{1,2}=E_{2,1}=-\frac{\bar{\mu}}{3\operatorname{Re}}(\xi_x\eta_y+\xi_y\eta_x). \tag*{(53)} E 1 , 2 = E 2 , 1 = − 3 Re μ ˉ ( ξ x η y + ξ y η x ) . (53) Matrices D \mathbf{D} D and B \mathbf{B} B coefficients of ∂ ξ q ^ \partial_{\xi}\hat{\mathbf{q}} ∂ ξ q ^ and ∂ η q ^ \partial_{\eta}\hat{\mathbf{q}} ∂ η q ^ (advection, pressure gradient, viscous-shear couplings). Each splits as D = D ( 0 ) + i β D ( β ) \mathbf{D}=\mathbf{D}^{(0)}+i\beta\mathbf{D}^{(\beta)} D = D ( 0 ) + i β D ( β ) B = B ( 0 ) + i β B ( β ) \mathbf{B}=\mathbf{B}^{(0)}+i\beta\mathbf{B}^{(\beta)} B = B ( 0 ) + i β B ( β ) β \beta β D \mathbf{D} D
D 1 , 1 ( 0 ) = ρ ˉ U ξ − μ ˉ Re Λ 4 / 3 , 1 ξ − μ ˉ T Re ( 4 3 ξ x ∂ x T ˉ + ξ y ∂ y T ˉ ) , D_{1,1}^{(0)}=\bar{\rho}U_{\xi}-\frac{\bar{\mu}}{\operatorname{Re}}\Lambda_{4/3,1}^{\xi}-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\frac{4}{3}\xi_x\partial_x\bar{T}+\xi_y\partial_y\bar{T}\right), D 1 , 1 ( 0 ) = ρ ˉ U ξ − Re μ ˉ Λ 4/3 , 1 ξ − Re μ ˉ T ( 3 4 ξ x ∂ x T ˉ + ξ y ∂ y T ˉ ) , D 1 , 2 ( 0 ) = − μ ˉ 3 Re ξ y − μ ˉ T Re ( ξ x ∂ y T ˉ − 2 3 ξ y ∂ x T ˉ ) , D_{1,2}^{(0)}=-\frac{\bar{\mu}}{3\operatorname{Re}}\xi_y-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\xi_x\partial_y\bar{T}-\frac{2}{3}\xi_y\partial_x\bar{T}\right), D 1 , 2 ( 0 ) = − 3 Re μ ˉ ξ y − Re μ ˉ T ( ξ x ∂ y T ˉ − 3 2 ξ y ∂ x T ˉ ) , D 1 , 4 ( 0 ) = T ˉ γ Ma 2 ξ x , D_{1,4}^{(0)}=\frac{\bar{T}}{\gamma\operatorname{Ma}^{2}}\xi_x, D 1 , 4 ( 0 ) = γ Ma 2 T ˉ ξ x , D 1 , 5 ( 0 ) = μ ˉ γ Ma 2 ξ x − μ ˉ T Re ( S ˉ x x ξ x + S ˉ x y ξ y ) , D_{1,5}^{(0)}=\frac{\bar{\mu}}{\gamma\operatorname{Ma}^{2}}\xi_x-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\bar{S}_{xx}\xi_x+\bar{S}_{xy}\xi_y\right), D 1 , 5 ( 0 ) = γ Ma 2 μ ˉ ξ x − Re μ ˉ T ( S ˉ xx ξ x + S ˉ x y ξ y ) , D 2 , 1 ( 0 ) = − μ ˉ 3 Re ξ y + μ ˉ T Re ( 2 3 ξ x ∂ y T ˉ − ξ y ∂ x T ˉ ) , D_{2,1}^{(0)}=-\frac{\bar{\mu}}{3\operatorname{Re}}\xi_y+\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\frac{2}{3}\xi_x\partial_y\bar{T}-\xi_y\partial_x\bar{T}\right), D 2 , 1 ( 0 ) = − 3 Re μ ˉ ξ y + Re μ ˉ T ( 3 2 ξ x ∂ y T ˉ − ξ y ∂ x T ˉ ) , D 2 , 2 ( 0 ) = ρ ˉ U ξ − μ ˉ Re Λ 1 , 4 / 3 ξ − μ ˉ T Re ( ξ x ∂ x T ˉ + 4 3 ξ y ∂ y T ˉ ) , D_{2,2}^{(0)}=\bar{\rho}U_{\xi}-\frac{\bar{\mu}}{\operatorname{Re}}\Lambda_{1,4/3}^{\xi}-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\xi_x\partial_x\bar{T}+\frac{4}{3}\xi_y\partial_y\bar{T}\right), D 2 , 2 ( 0 ) = ρ ˉ U ξ − Re μ ˉ Λ 1 , 4/3 ξ − Re μ ˉ T ( ξ x ∂ x T ˉ + 3 4 ξ y ∂ y T ˉ ) , D 2 , 4 ( 0 ) = T ˉ γ Ma 2 ξ y , D_{2,4}^{(0)}=\frac{\bar{T}}{\gamma\operatorname{Ma}^{2}}\xi_y, D 2 , 4 ( 0 ) = γ Ma 2 T ˉ ξ y , D 2 , 5 ( 0 ) = μ ˉ γ Ma 2 ξ y − μ ˉ T Re ( S ˉ x y ξ x + S ˉ y y ξ y ) , D_{2,5}^{(0)}=\frac{\bar{\mu}}{\gamma\operatorname{Ma}^{2}}\xi_y-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\bar{S}_{xy}\xi_x+\bar{S}_{yy}\xi_y\right), D 2 , 5 ( 0 ) = γ Ma 2 μ ˉ ξ y − Re μ ˉ T ( S ˉ x y ξ x + S ˉ y y ξ y ) , D 3 , 3 ( 0 ) = ρ ˉ U ξ − μ ˉ Re Λ 1 , 1 ξ − μ ˉ T Re ( ξ x ∂ x T ˉ + ξ y ∂ y T ˉ ) , D_{3,3}^{(0)}=\bar{\rho}U_{\xi}-\frac{\bar{\mu}}{\operatorname{Re}}\Lambda_{1,1}^{\xi}-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\xi_x\partial_x\bar{T}+\xi_y\partial_y\bar{T}\right), D 3 , 3 ( 0 ) = ρ ˉ U ξ − Re μ ˉ Λ 1 , 1 ξ − Re μ ˉ T ( ξ x ∂ x T ˉ + ξ y ∂ y T ˉ ) , D 3 , 5 ( 0 ) = − μ ˉ T Re ( ξ x ∂ x w ˉ + ξ y ∂ y w ˉ ) , D_{3,5}^{(0)}=-\frac{\bar{\mu}_{T}}{\operatorname{Re}}\left(\xi_x\partial_x\bar{w}+\xi_y\partial_y\bar{w}\right), D 3 , 5 ( 0 ) = − Re μ ˉ T ( ξ x ∂ x w ˉ + ξ y ∂ y w ˉ ) , D 4 , 1 ( 0 ) = ρ ˉ ξ x , D 4 , 2 ( 0 ) = ρ ˉ ξ y , D 4 , 4 ( 0 ) = U ξ , D_{4,1}^{(0)}=\bar{\rho}\xi_x,\qquad
D_{4,2}^{(0)}=\bar{\rho}\xi_y,\qquad
D_{4,4}^{(0)}=U_{\xi}, D 4 , 1 ( 0 ) = ρ ˉ ξ x , D 4 , 2 ( 0 ) = ρ ˉ ξ y , D 4 , 4 ( 0 ) = U ξ , D 5 , 1 ( 0 ) = ( γ − 1 ) ρ ˉ T ˉ ξ x − c Φ μ ˉ ( S ˉ x x ξ x + S ˉ x y ξ y ) , D_{5,1}^{(0)}=(\gamma-1)\bar{\rho}\bar{T}\xi_x-c_{\Phi}\bar{\mu}\left(\bar{S}_{xx}\xi_x+\bar{S}_{xy}\xi_y\right), D 5 , 1 ( 0 ) = ( γ − 1 ) ρ ˉ T ˉ ξ x − c Φ μ ˉ ( S ˉ xx ξ x + S ˉ x y ξ y ) , D 5 , 2 ( 0 ) = ( γ − 1 ) ρ ˉ T ˉ ξ y − c Φ μ ˉ ( S ˉ x y ξ x + S ˉ y y ξ y ) , D_{5,2}^{(0)}=(\gamma-1)\bar{\rho}\bar{T}\xi_y-c_{\Phi}\bar{\mu}\left(\bar{S}_{xy}\xi_x+\bar{S}_{yy}\xi_y\right), D 5 , 2 ( 0 ) = ( γ − 1 ) ρ ˉ T ˉ ξ y − c Φ μ ˉ ( S ˉ x y ξ x + S ˉ y y ξ y ) , D 5 , 3 ( 0 ) = − c Φ μ ˉ ( ξ x ∂ x w ˉ + ξ y ∂ y w ˉ ) , D_{5,3}^{(0)}=-c_{\Phi}\bar{\mu}\left(\xi_x\partial_x\bar{w}+\xi_y\partial_y\bar{w}\right), D 5 , 3 ( 0 ) = − c Φ μ ˉ ( ξ x ∂ x w ˉ + ξ y ∂ y w ˉ ) , D 5 , 5 ( 0 ) = ρ ˉ U ξ − γ μ ˉ Pr Re Λ 1 , 1 ξ − 2 γ μ ˉ T Pr Re ( ξ x ∂ x T ˉ + ξ y ∂ y T ˉ ) . (54) D_{5,5}^{(0)}=\bar{\rho}U_{\xi}-\frac{\gamma\bar{\mu}}{\operatorname{Pr}\operatorname{Re}}\Lambda_{1,1}^{\xi}-\frac{2\gamma\bar{\mu}_{T}}{\operatorname{Pr}\operatorname{Re}}\left(\xi_x\partial_x\bar{T}+\xi_y\partial_y\bar{T}\right). \tag*{(54)} D 5 , 5 ( 0 ) = ρ ˉ U ξ − Pr Re γ μ ˉ Λ 1 , 1 ξ − Pr Re 2 γ μ ˉ T ( ξ x ∂ x T ˉ + ξ y ∂ y T ˉ ) . (54) The i β i\beta i β
D 1 , 3 ( β ) = D 3 , 1 ( β ) = − μ ˉ 3 Re ξ x , D 2 , 3 ( β ) = D 3 , 2 ( β ) = − μ ˉ 3 Re ξ y . (55) D_{1,3}^{(\beta)}=D_{3,1}^{(\beta)}=-\frac{\bar{\mu}}{3\operatorname{Re}}\xi_x,\qquad
D_{2,3}^{(\beta)}=D_{3,2}^{(\beta)}=-\frac{\bar{\mu}}{3\operatorname{Re}}\xi_y. \tag*{(55)} D 1 , 3 ( β ) = D 3 , 1 ( β ) = − 3 Re μ ˉ ξ x , D 2 , 3 ( β ) = D 3 , 2 ( β ) = − 3 Re μ ˉ ξ y . (55) The matrix B \mathbf{B} B B ( 0 ) \mathbf{B}^{(0)} B ( 0 ) B ( β ) \mathbf{B}^{(\beta)} B ( β ) D \mathbf{D} D ξ \xi ξ η \eta η U ξ → U η U_{\xi}\to U_{\eta} U ξ → U η ξ x → η x \xi_x\to\eta_x ξ x → η x ξ y → η y \xi_y\to\eta_y ξ y → η y ξ x y → η x y \xi_{xy}\to\eta_{xy} ξ x y → η x y Λ a , b ξ → Λ a , b η \Lambda_{a,b}^{\xi}\to\Lambda_{a,b}^{\eta} Λ a , b ξ → Λ a , b η η \eta η ξ \xi ξ
Matrix A \mathbf{A} A coefficient of q ^ \hat{\mathbf{q}} q ^ (algebraic block). The algebraic block lumps the no-derivative contributions of the linearized perturbation equations together with the ( − i ω ) (-i\omega) ( − iω ) i β i\beta i β β 2 \beta^2 β 2
A = A ( 0 ) + ( − i ω ) A ( ω ) + i β A ( β ) + β 2 A ( β β ) . (56) \mathbf{A}=\mathbf{A}^{(0)}+(-i\omega)\mathbf{A}^{(\omega)}+i\beta\mathbf{A}^{(\beta)}+\beta^2\mathbf{A}^{(\beta\beta)}. \tag*{(56)} A = A ( 0 ) + ( − iω ) A ( ω ) + i β A ( β ) + β 2 A ( β β ) . (56) The time-derivative mass matrix is
A ( ω ) = diag ( ρ ˉ , ρ ˉ , ρ ˉ , 1 , ρ ˉ ) . (57) \mathbf{A}^{(\omega)}=\operatorname{diag}(\bar{\rho},\bar{\rho},\bar{\rho},1,\bar{\rho}). \tag*{(57)} A ( ω ) = diag ( ρ ˉ , ρ ˉ , ρ ˉ , 1 , ρ ˉ ) . (57) The β 2 \beta^2 β 2
A ( β β ) = diag ( μ ˉ R e , μ ˉ R e , 4 μ ˉ 3 R e , 0 , μ ˉ P r R e ) . (58) \mathbf{A}^{(\beta\beta)}=\operatorname{diag}\left(\frac{\bar{\mu}}{\mathrm{Re}},\frac{\bar{\mu}}{\mathrm{Re}},\frac{4\bar{\mu}}{3\mathrm{Re}},0,\frac{\bar{\mu}}{\mathrm{Pr}\,\mathrm{Re}}\right). \tag*{(58)} A ( β β ) = diag ( Re μ ˉ , Re μ ˉ , 3 Re 4 μ ˉ , 0 , Pr Re μ ˉ ) . (58) The i β i\beta i β A ( β ) \mathbf{A}^{(\beta)} A ( β )
A 1 , 1 ( β ) = ρ ˉ w ˉ , A 1 , 3 ( β ) = 2 μ ˉ T ∂ x T ˉ 3 R e , A 1 , 5 ( β ) = − μ ˉ T ∂ x w ˉ R e , A 2 , 2 ( β ) = ρ ˉ w ˉ , A 2 , 3 ( β ) = 2 μ ˉ T ∂ y T ˉ 3 R e , A 2 , 5 ( β ) = − μ ˉ T ∂ y w ˉ R e , A 3 , 1 ( β ) = − μ ˉ T ∂ x T ˉ R e , A 3 , 2 ( β ) = − μ ˉ T ∂ y T ˉ R e , A 3 , 3 ( β ) = ρ ˉ w ˉ , A 3 , 4 ( β ) = T ˉ γ M a 2 , A 3 , 5 ( β ) = ρ ˉ γ M a 2 + 2 μ ˉ T D ˉ u 3 R e , A 4 , 3 ( β ) = ρ ˉ , A 4 , 5 ( β ) = w ˉ , A 5 , 1 ( β ) = − c Φ μ ˉ S ˉ x z , A 5 , 2 ( β ) = − c Φ μ ˉ S ˉ y z , A 5 , 3 ( β ) = ( γ − 1 ) ρ ˉ T ˉ − c Φ μ ˉ S ˉ z z , A 5 , 5 ( β ) = ρ ˉ w ˉ . (59) \begin{aligned}
A_{1,1}^{(\beta)}&=\bar{\rho}\bar{w},&
A_{1,3}^{(\beta)}&=\frac{2\bar{\mu}_{T}\,\partial_x\bar{T}}{3\mathrm{Re}},\\
A_{1,5}^{(\beta)}&=-\frac{\bar{\mu}_{T}\,\partial_x\bar{w}}{\mathrm{Re}},&
A_{2,2}^{(\beta)}&=\bar{\rho}\bar{w},\\
A_{2,3}^{(\beta)}&=\frac{2\bar{\mu}_{T}\,\partial_y\bar{T}}{3\mathrm{Re}},&
A_{2,5}^{(\beta)}&=-\frac{\bar{\mu}_{T}\,\partial_y\bar{w}}{\mathrm{Re}},\\
A_{3,1}^{(\beta)}&=-\frac{\bar{\mu}_{T}\,\partial_x\bar{T}}{\mathrm{Re}},&
A_{3,2}^{(\beta)}&=-\frac{\bar{\mu}_{T}\,\partial_y\bar{T}}{\mathrm{Re}},\\
A_{3,3}^{(\beta)}&=\bar{\rho}\bar{w},&
A_{3,4}^{(\beta)}&=\frac{\bar{T}}{\gamma\mathrm{Ma}^{2}},\\
A_{3,5}^{(\beta)}&=\frac{\bar{\rho}}{\gamma\mathrm{Ma}^{2}}+\frac{2\bar{\mu}_{T}\bar{D}_{u}}{3\mathrm{Re}},&
A_{4,3}^{(\beta)}&=\bar{\rho},\\
A_{4,5}^{(\beta)}&=\bar{w},&
A_{5,1}^{(\beta)}&=-c_{\Phi}\bar{\mu}\,\bar{S}_{xz},\\
A_{5,2}^{(\beta)}&=-c_{\Phi}\bar{\mu}\,\bar{S}_{yz},&
A_{5,3}^{(\beta)}&=(\gamma-1)\bar{\rho}\bar{T}-c_{\Phi}\bar{\mu}\,\bar{S}_{zz},\\
A_{5,5}^{(\beta)}&=\bar{\rho}\bar{w}. \tag*{(59)}
\end{aligned} A 1 , 1 ( β ) A 1 , 5 ( β ) A 2 , 3 ( β ) A 3 , 1 ( β ) A 3 , 3 ( β ) A 3 , 5 ( β ) A 4 , 5 ( β ) A 5 , 2 ( β ) A 5 , 5 ( β ) = ρ ˉ w ˉ , = − Re μ ˉ T ∂ x w ˉ , = 3 Re 2 μ ˉ T ∂ y T ˉ , = − Re μ ˉ T ∂ x T ˉ , = ρ ˉ w ˉ , = γ Ma 2 ρ ˉ + 3 Re 2 μ ˉ T D ˉ u , = w ˉ , = − c Φ μ ˉ S ˉ y z , = ρ ˉ w ˉ . A 1 , 3 ( β ) A 2 , 2 ( β ) A 2 , 5 ( β ) A 3 , 2 ( β ) A 3 , 4 ( β ) A 4 , 3 ( β ) A 5 , 1 ( β ) A 5 , 3 ( β ) = 3 Re 2 μ ˉ T ∂ x T ˉ , = ρ ˉ w ˉ , = − Re μ ˉ T ∂ y w ˉ , = − Re μ ˉ T ∂ y T ˉ , = γ Ma 2 T ˉ , = ρ ˉ , = − c Φ μ ˉ S ˉ x z , = ( γ − 1 ) ρ ˉ T ˉ − c Φ μ ˉ S ˉ z z , (59) The β \beta β A ( 0 ) \mathbf{A}^{(0)} A ( 0 )
A 1 , 1 ( 0 ) = ρ ˉ ∂ x u ˉ , A 1 , 2 ( 0 ) = ρ ˉ ∂ y u ˉ , A 1 , 4 ( 0 ) = u ˉ ⋅ ∇ u ˉ + ∂ x T ˉ γ M a 2 , A 1 , 5 ( 0 ) = ∂ x ρ ˉ γ M a 2 − μ ˉ T L ˉ u R e − μ ˉ T T R e ( S ˉ x x ∂ x T ˉ + S ˉ x y ∂ y T ˉ ) , A 2 , 1 ( 0 ) = ρ ˉ ∂ x v ˉ , A 2 , 2 ( 0 ) = ρ ˉ ∂ y v ˉ , A 2 , 4 ( 0 ) = u ˉ ⋅ ∇ v ˉ + ∂ y T ˉ γ M a 2 , A 2 , 5 ( 0 ) = ∂ y ρ ˉ γ M a 2 − μ ˉ T L ˉ v R e − μ ˉ T T R e ( S ˉ y y ∂ y T ˉ + S ˉ x y ∂ x T ˉ ) , A 3 , 1 ( 0 ) = ρ ˉ ∂ x w ˉ , A 3 , 2 ( 0 ) = ρ ˉ ∂ y w ˉ , A 3 , 4 ( 0 ) = u ˉ ⋅ ∇ w ˉ , A 3 , 5 ( 0 ) = − μ ˉ T L ˉ w R e − μ ˉ T T R e ( ∂ y w ˉ ∂ y T ˉ + ∂ x w ˉ ∂ x T ˉ ) , A 4 , 1 ( 0 ) = ∂ x ρ ˉ , A 4 , 2 ( 0 ) = ∂ y ρ ˉ , A 4 , 4 ( 0 ) = D ˉ u , A 5 , 1 ( 0 ) = ρ ˉ ∂ x T ˉ , A 5 , 2 ( 0 ) = ρ ˉ ∂ y T ˉ , A 5 , 4 ( 0 ) = u ˉ ⋅ ∇ T ˉ + ( γ − 1 ) T ˉ D ˉ u , A 5 , 5 ( 0 ) = ( γ − 1 ) ρ ˉ D ˉ u + L ˉ T − γ ( γ − 1 ) M a 2 R e μ ˉ T Φ ˉ 0 . (60) \begin{aligned}
A_{1,1}^{(0)}&=\bar{\rho}\,\partial_x\bar{u},&
A_{1,2}^{(0)}&=\bar{\rho}\,\partial_y\bar{u},\\
A_{1,4}^{(0)}&=\bar{\mathbf{u}}\cdot\nabla\bar{u}+\frac{\partial_x\bar{T}}{\gamma\mathrm{Ma}^{2}},\\
A_{1,5}^{(0)}&=\frac{\partial_x\bar{\rho}}{\gamma\mathrm{Ma}^{2}}-\frac{\bar{\mu}_{T}\bar{L}_{u}}{\mathrm{Re}}
-\frac{\bar{\mu}_{TT}}{\mathrm{Re}}\left(\bar{S}_{xx}\partial_x\bar{T}+\bar{S}_{xy}\partial_y\bar{T}\right),\\
A_{2,1}^{(0)}&=\bar{\rho}\,\partial_x\bar{v},&
A_{2,2}^{(0)}&=\bar{\rho}\,\partial_y\bar{v},\\
A_{2,4}^{(0)}&=\bar{\mathbf{u}}\cdot\nabla\bar{v}+\frac{\partial_y\bar{T}}{\gamma\mathrm{Ma}^{2}},\\
A_{2,5}^{(0)}&=\frac{\partial_y\bar{\rho}}{\gamma\mathrm{Ma}^{2}}-\frac{\bar{\mu}_{T}\bar{L}_{v}}{\mathrm{Re}}
-\frac{\bar{\mu}_{TT}}{\mathrm{Re}}\left(\bar{S}_{yy}\partial_y\bar{T}+\bar{S}_{xy}\partial_x\bar{T}\right),\\
A_{3,1}^{(0)}&=\bar{\rho}\,\partial_x\bar{w},&
A_{3,2}^{(0)}&=\bar{\rho}\,\partial_y\bar{w},\\
A_{3,4}^{(0)}&=\bar{\mathbf{u}}\cdot\nabla\bar{w},\\
A_{3,5}^{(0)}&=-\frac{\bar{\mu}_{T}\bar{L}_{w}}{\mathrm{Re}}
-\frac{\bar{\mu}_{TT}}{\mathrm{Re}}\left(\partial_y\bar{w}\,\partial_y\bar{T}+\partial_x\bar{w}\,\partial_x\bar{T}\right),\\
A_{4,1}^{(0)}&=\partial_x\bar{\rho},&
A_{4,2}^{(0)}&=\partial_y\bar{\rho},&
A_{4,4}^{(0)}&=\bar{D}_{u},\\
A_{5,1}^{(0)}&=\bar{\rho}\,\partial_x\bar{T},&
A_{5,2}^{(0)}&=\bar{\rho}\,\partial_y\bar{T},\\
A_{5,4}^{(0)}&=\bar{\mathbf{u}}\cdot\nabla\bar{T}+(\gamma-1)\bar{T}\,\bar{D}_{u},\\
A_{5,5}^{(0)}&=(\gamma-1)\bar{\rho}\bar{D}_{u}+\bar{L}_{T}
-\frac{\gamma(\gamma-1)\mathrm{Ma}^{2}}{\mathrm{Re}}\bar{\mu}_{T}\bar{\Phi}_{0}. \tag*{(60)}
\end{aligned} A 1 , 1 ( 0 ) A 1 , 4 ( 0 ) A 1 , 5 ( 0 ) A 2 , 1 ( 0 ) A 2 , 4 ( 0 ) A 2 , 5 ( 0 ) A 3 , 1 ( 0 ) A 3 , 4 ( 0 ) A 3 , 5 ( 0 ) A 4 , 1 ( 0 ) A 5 , 1 ( 0 ) A 5 , 4 ( 0 ) A 5 , 5 ( 0 ) = ρ ˉ ∂ x u ˉ , = u ˉ ⋅ ∇ u ˉ + γ Ma 2 ∂ x T ˉ , = γ Ma 2 ∂ x ρ ˉ − Re μ ˉ T L ˉ u − Re μ ˉ T T ( S ˉ xx ∂ x T ˉ + S ˉ x y ∂ y T ˉ ) , = ρ ˉ ∂ x v ˉ , = u ˉ ⋅ ∇ v ˉ + γ Ma 2 ∂ y T ˉ , = γ Ma 2 ∂ y ρ ˉ − Re μ ˉ T L ˉ v − Re μ ˉ T T ( S ˉ y y ∂ y T ˉ + S ˉ x y ∂ x T ˉ ) , = ρ ˉ ∂ x w ˉ , = u ˉ ⋅ ∇ w ˉ , = − Re μ ˉ T L ˉ w − Re μ ˉ T T ( ∂ y w ˉ ∂ y T ˉ + ∂ x w ˉ ∂ x T ˉ ) , = ∂ x ρ ˉ , = ρ ˉ ∂ x T ˉ , = u ˉ ⋅ ∇ T ˉ + ( γ − 1 ) T ˉ D ˉ u , = ( γ − 1 ) ρ ˉ D ˉ u + L ˉ T − Re γ ( γ − 1 ) Ma 2 μ ˉ T Φ ˉ 0 . A 1 , 2 ( 0 ) A 2 , 2 ( 0 ) A 3 , 2 ( 0 ) A 4 , 2 ( 0 ) A 5 , 2 ( 0 ) = ρ ˉ ∂ y u ˉ , = ρ ˉ ∂ y v ˉ , = ρ ˉ ∂ y w ˉ , = ∂ y ρ ˉ , = ρ ˉ ∂ y T ˉ , A 4 , 4 ( 0 ) = D ˉ u , (60) Elliptic mesh smoother
The Laplace mesh equations (37) are discretized by second-order central differences on the uniform auxiliary grid (Δ ξ = Δ η = 1 \Delta\xi=\Delta\eta=1 Δ ξ = Δ η = 1 i i i ξ \xi ξ j j j η \eta η X i , j X_{i,j} X i , j Y i , j Y_{i,j} Y i , j (38) are evaluated by central differences on the current mesh, so the discrete system is nonlinear. Solving the discrete equation for the self-term gives the Gauss–Seidel update
X i , j G S = 1 2 ( g ξ ξ + g η η ) [ g ξ ξ ( X i , j + 1 + X i , j − 1 ) + g η η ( X i + 1 , j + X i − 1 , j ) − g ξ η 2 ( X i + 1 , j + 1 − X i − 1 , j + 1 − X i + 1 , j − 1 + X i − 1 , j − 1 ) ] . (61) X_{i,j}^{\mathrm{GS}}=\frac{1}{2(g_{\xi\xi}+g_{\eta\eta})}\left[g_{\xi\xi}(X_{i,j+1}+X_{i,j-1})+g_{\eta\eta}(X_{i+1,j}+X_{i-1,j})-\frac{g_{\xi\eta}}{2}(X_{i+1,j+1}-X_{i-1,j+1}-X_{i+1,j-1}+X_{i-1,j-1})\right]. \tag*{(61)} X i , j GS = 2 ( g ξ ξ + g η η ) 1 [ g ξ ξ ( X i , j + 1 + X i , j − 1 ) + g η η ( X i + 1 , j + X i − 1 , j ) − 2 g ξ η ( X i + 1 , j + 1 − X i − 1 , j + 1 − X i + 1 , j − 1 + X i − 1 , j − 1 ) ] . (61) and analogously for Y Y Y X i , j ← ( 1 − ω S O R ) X i , j + ω S O R X i , j G S X_{i,j}\leftarrow(1-\omega_{\mathrm{SOR}})X_{i,j}+\omega_{\mathrm{SOR}}X_{i,j}^{\mathrm{GS}} X i , j ← ( 1 − ω SOR ) X i , j + ω SOR X i , j GS ω S O R = 1.7 \omega_{\mathrm{SOR}}=1.7 ω SOR = 1.7 η \eta η § 3.1 are held fixed. The sweeps are repeated until max ∣ Δ X ∣ + max ∣ Δ Y ∣ \max|\Delta X|+\max|\Delta Y| max ∣Δ X ∣ + max ∣Δ Y ∣ ξ \xi ξ
Summary of calculation parameters
Tables C.5 and C.6 collect the physical and numerical parameters of every calculation presented in this paper.
Computational hardware
Wall-clock measurements use two systems, continuing the nomenclature of Westerbeek et al. [14 ] . System D is a laptop with an Apple M1 Pro processor (6 performance and 2 efficiency cores, up to 3.2 GHz) and 32 GB of unified memory. System E is a workstation with an Intel Core i7-14700K processor (8 performance and 12 efficiency cores, up to 5.6 GHz) and 188 GB of DDR5 memory available. The thread count used for each measurement is stated in § 7 . Factor sizes are identical between runs of the same calculation, whereas factorization times vary slightly between repeated runs on the same system, owing to operating-system scheduling and memory pressure.
Table C.5. Physical parameters of the calculations. R = R e x R = \sqrt{Re_x} R = R e x γ \gamma γ ( M , N ) (M,N) ( M , N ) A 0 A_0 A 0 2 ∣ u ^ ∣ 2|\hat{u}| 2∣ u ^ ∣ F = ω 1 / R e 0 F = \omega_1/Re_0 F = ω 1 / R e 0
Table C.6. Numerical settings of the calculations. Entries in italic are configurable code keywords. Solver: the Picard iteration of Sec. 2.6, with Anderson acceleration (AA, m A = 5 m_A = 5 m A = 5 θ \theta θ η \eta η FD 4 = fourth-order finite difference, Ch = Chebyshev); IC: inflow condition (CLST , ILST = compressible, incompressible local stability solution; PSE = embedded parabolized stability solution; LOAD = externally supplied profile, here the reference DNS); BC: top boundary (n n n H _ D R H\_DR H _ D R isoth. = isothermal); buffer: outflow treatment (para = parabolization, on = amplitude damping) with onset ξ b \xi_b ξ b
