With the growing demand for innovative aircraft and aerial designs to combat climate change, there is an increasing focus on environmentally responsible commercial aviation, as reported by Coder and Maughmer (2014). The primary objectives are to enhance fuel efficiency and reduce aerodynamic drag, which can be achieved by maintaining laminar boundary layer flow. This approach offers a potential tenfold reduction in friction drag compared to turbulent boundary layers (Krishnan et al., 2017). Since friction drag can account for nearly 50% of the total drag experienced by aircraft during cruise, delaying the transition to turbulence in the boundary layer is crucial for developing fuel-efficient designs (Karpuk and Mosca, 2024). Natural laminar flow airfoils present a viable strategy for achieving higher fuel efficiency. Recent studies by Halila et al. (2020) have demonstrated that integrating natural laminar flow technology with flow separation control or control of shock waves in the transonic regime (Chakraborty et al., 2022) can reduce total drag by 15% or more for typical jetliners at cruise conditions. As a result, natural laminar flow designs, once conceived as strictly experimental and/or conceptual, have gained relevance in modern commercial aircraft, business jets (Fujino et al., 2003), and unmanned aerial vehicles. This resurgence of interest in natural laminar flow technology can be attributed to recent advancements in high-fidelity simulations that accurately predict laminar-turbulent transitions (Halila et al., 2019) and effectively capture unsteady shock structures (Sengupta et al., 2021) in transonic aerodynamics. Accurate forecasting of transition onset and shock locations significantly influences boundary layer development, flow separation, friction drag, and maximum lift coefficients (Campbell and Lynde, 2017), all of which are critical to the design and performance of aerodynamic bodies.

Shock boundary layer interactions take place when shock waves interact unfavorably with the boundary layer producing large undesirable local fluctuations in properties such as pressure, shear stress and the rate of heat transfer. This is particularly significant in the transonic flow regime, as these off-design events can substantially affect both aerodynamic and thermodynamic properties, altering the flow field in notable ways. Such modifications lead to changes in parameters like pressure distribution and boundary layer characteristics, resulting in increased unsteadiness and higher drag (Cole and Cook, 2012). A thorough investigation of this dynamic phenomenon requires solving the compressible Navier-Stokes equations to accurately capture critical parameters such as shock location and strength, unsteady aerodynamic forces, and potential strategies for reducing the effect of shock waves. While some canonical numerical (Larsson et al., 2013) and experimental (Barre et al., 1996) studies have explored the interaction between shock waves and turbulence, the flow over an airfoil at varying angles of attack presents a more complex scenario due to the presence of variable streamwise pressure gradients, even in laminar flows. Transonic flows are characterized by unsteady shock wave systems. To better understand this intricately time-dependent behavior, researchers have analyzed the effects of downstream periodic pressure perturbations on shock waves (Bruce and Babinsky, 2010), shedding light on the complex nature of shock-boundary layer interactions.

Previous studies by Touré and Schuelein (2018), Gross and Lee (2018), Quadros and Bernardini (2018) on transonic shock boundary layer interactions did not address efforts to modify or control the interactions between shock waves and the underlying boundary layer. Typically, these transient behaviors exhibit large-amplitude normal or near-normal shocks accompanied by low-frequency motion, as seen by (Dussauge and Piponniau, 2008). Such behavior is problematic due to the resulting unsteady pressure fluctuations on the airfoil, which can lead to increased aerodynamic loads (Giannelis et al., 2017). Another notable aspect of transonic shock boundary layer interactions is the presence of upstream-propagating Kutta waves interacting with the shock system, which includes both oblique and normal shocks (Lee, 2001). The low-frequency motion of these shock waves induces similar low-frequency pressure fluctuations on the airfoil surface, a phenomenon referred to as “transonic buffeting”. This buffeting can intensify airfoil vibrations and ultimately pose risks of structural failure (Giannelis et al., 2017).

The preceding discussion underscores the necessity for a comprehensive assessment of the aerodynamic performance of a natural laminar flow airfoil under both design and off-design conditions. This assessment will serve as a benchmark for future validation efforts and establish a foundation for strategies aimed at controlling flow separation and shocks to minimize friction drag. Additionally, a thorough investigation into boundary layer characteristics and unsteady separation will contribute to the design of optimized natural laminar flow airfoils. In this study, we will simulate a range of free-stream Mach numbers, Reynolds numbers, and operational conditions for the SHM1 airfoil, which is integral to the design of the Honda business jet (Fujino et al., 2003). To accurately capture pressure waves, shock structures, and boundary layer interactions, we will employ dispersion relation-preserving compact schemes (Sagaut et al., 2023), which are effective in resolving both temporal and spatial scales within the flow. We ensure the integrity of our numerical approach by implementing an error-free non-overlapping parallelization strategy, which maintains the same level of accuracy as sequential computing (Sundaram et al., 2023) through global spectral analysis. This strategy has previously demonstrated its capability in effectively capturing shock boundary layer interactions in compressible transonic flow, as in Chakraborty et al. (2022), Sengupta et al. (2022), with results validated against flight test data from Fujino et al. (2003) as well as benchmark wind tunnel results of Harris (1981).

The structure of the paper is as follows: The next section outlines the problem formulation for simulating flow around the SHM1 airfoil, detailing the governing equations and test cases. Section 3 describes the numerical methods employed and the validation efforts with experimental data. In Section 4, we present the results and discussions, including instantaneous Schlieren visualizations under various design and off-design conditions, vorticity spectra, and coefficients of pressure and skin friction. Finally, we evaluate the aerodynamic performance by comparing lift and drag coefficients across the different operating conditions. The paper concludes with a summary and final remarks in Section 5.

Figure 1 depicts the schematic of the computational domain, which employs an O-grid topology. This grid was generated using a hyperbolic technique in Pointwise, featuring 1,251 points in the azimuthal ( ξ ) direction and 401 points in the wall-normal ( η ) direction. The corresponding Δ y + values are in the range of 0.09–0.25 for the test cases reported here. The zoomed view of the grid including wall-normal spacing has been reported in a previous work Sengupta et al. (2021), which used the same grid. Similarly, the mesh independence and validation efforts for the adopted mesh on the SHM1 airfoil has been reported (Sundaram et al., 2023), prompting us to use the same grid here. A cut is introduced along the ξ direction, as shown in Figure 1, to simplify the domain. This cut allows for natural periodicity of all variables along its boundaries. The outer boundary of the computational domain extends to 16 c , ensuring a nearly uniform grid in the wake region in both directions. A shock-free mesh is crucial for accurately capturing the dynamics of traveling waves in compressible flow around airfoils. The grid points are concentrated near the wall in the η direction, as well as near the leading edge, trailing edge, and the region ahead of the separation point on the upper surface in the ξ direction. The numerical simulations solve the two-dimensional compressible Navier-Stokes equations, formulated by established notations, as detailed in (Sengupta et al., 2022; Sengupta and Shandilya, 2024). ∂ Q ̂ ∂ t + ∂ E ̂ ∂ x + ∂ F ̂ ∂ y = ∂ E v ̂ ∂ x + ∂ F v ̂ ∂ y (1) where the conservative variables are given as, Q ̂ = [ ρ ρ u ρ v e t ] T . The convective flux vectors are similarly given as, E ̂ = [ ρ u ρ u 2 + p ρ u v ( ρ e t + p ) u ] T F ̂ = [ ρ v ρ u v ρ v 2 + p ( ρ e t + p ) v ] T

and the viscous flux vectors are given as, E v ̂ = [ 0 τ x x τ x y ( u τ x x + v τ x y − q x ) ] T F v ̂ = [ 0 τ y x τ y y ( u τ y x + v τ y y − q y ) ] T

In the given expressions, ρ , e t , and p denote dimensionless values of density, total specific energy, and pressure, respectively. These physical variables are normalized with respect to the free-stream density ( ρ ∞ ) , free-stream velocity ( U ∞ = M ∞ γ R T ∞ ) , free-stream temperature ( T ∞ ) , free-stream dynamic viscosity ( μ ∞ ) , length scale ( c ) , and the time scale, ( c / U ∞ ) . Here, γ = 1.4 represents the specific heat capacity ratio or the adiabatic index for air. The dimensionless parameters, namely, the Prandtl number ( P r ) , free-stream Reynolds number ( R e ∞ ) , and free-stream Mach number, ( M ∞ ) , are defined as follows: P r = μ C p k ̄ ; R e ∞ = ρ ∞ U ∞ c μ ∞ ; M ∞ = U ∞ a ∞ where a ∞ is the free-stream speed of sound. The heat conduction terms involved are given by q x = − μ P r R e ∞ γ − 1 M ∞ 2 ∂ T ∂ x q y = − μ P r R e ∞ γ − 1 M ∞ 2 ∂ T ∂ y

The components of the symmetric Newtonian viscous stress tensors, τ x x , τ x y , τ y x , τ y y , are defined as τ x x = 1 R e ∞ 4 3 μ + κ μ ∞ ∂ u ∂ x + − 2 3 μ + κ μ ∞ ∂ v ∂ y τ y y = 1 R e ∞ 4 3 μ + κ μ ∞ ∂ v ∂ y + − 2 3 μ + κ μ ∞ ∂ u ∂ x τ x y = τ y x = μ R e ∞ ∂ u ∂ y + ∂ v ∂ x

The equations from the Cartesian space ( x , y ) are transformed to body-fitted computational grid ( ξ , η ) using the following relations: ξ = ξ ( x , y ) and η = η ( x , y ) . The transformed plane equations in strong conservation form are given as, ∂ Q ∂ t + ∂ E ∂ ξ + ∂ F ∂ η = ∂ E v ∂ ξ + ∂ F v ∂ η (2) with the state variables and flux vectors, given as Q = Q ̂ / J E = ξ x E ̂ + ξ y F ̂ / J F = η x E ̂ + η y F ̂ / J E v = ξ x E v ̂ + ξ y F v ̂ / J F v = η x E v ̂ + η y F v ̂ / J

Here J is the Jacobian of the grid transformation given by J = 1 x ξ y η − x η y ξ

The grid metrics, ξ x , ξ y , η x , and η y , are computed during the creation of the O-grid using a hyperbolic grid generation technique in Pointwise. These are expressed as follows: ξ x = J y η ; ξ y = − J x η ; η x = − J y ξ ; η y = J x ξ . This grid transformation ensures that the solid airfoil boundary aligns with one of the grid lines ( η = 0 ) . Using these transformations, the heat conduction terms in the transformed plane are given by, q x = − μ P r R e ∞ γ − 1 M ∞ 2 ξ x T ξ + η x T η q y = − μ P r R e ∞ γ − 1 M ∞ 2 ξ y T ξ + η y T η

The viscous stress components in the transformed plane are given by τ x x = 1 R e ∞ 4 3 μ + κ μ ∞ ξ x u ξ + η x u η − − 2 3 μ + κ μ ∞ ξ y v ξ + η y v η τ y y = 1 R e ∞ 4 3 μ + κ μ ∞ ξ y v ξ + η y v η − − 2 3 μ + κ μ ∞ ξ x v ξ + η x u η τ x y = τ y x = μ R e ∞ ξ y u ξ + η y u η + ξ x v ξ + η x v η

On the airfoil surface, a no-slip wall boundary condition is applied to the velocity components, specified as u = v = 0 . An adiabatic wall condition is enforced on the heat conduction terms, ensuring no heat transfer occurs across the airfoil surface: ( q x ) wall = ( q y ) wall = 0 . In the far field, the flow is approximated as one-dimensional in the η -direction, where a non-reflective characteristic-like boundary condition (Pulliam, 1986) derived from the Euler equations is implemented. These conditions are based on the signs of the eigenvalues of the linearized one-dimensional Euler equations, which depend on whether the flow is subsonic or supersonic at the far field. At the domain’s inflow, characteristic variables are set to the free-stream values corresponding to the specified Mach number. At the outflow, variables are extrapolated from the interior when the local Mach number is supersonic. To prevent spurious acoustic wave reflections that could compromise the physical integrity of the domain and to ensure accurate implementation of the far-field boundary conditions, the outer boundary is positioned at approximately 16 c .

The current simulations utilize highly accurate dispersion relation preserving compact schemes for the spatial discretization and time integration of the governing equations (Sagaut et al., 2023). Convective flux derivatives in Equations 1, 2 are computed using an optimized upwind compact scheme, O U C S 3 , with explicit boundary closures (Sagaut et al., 2023). To maintain the isotropic nature of viscous flux derivatives, a second-order central difference method (CD2) is employed for discretization. A novel parallelization strategy (Sundaram et al., 2022) is implemented without overlap points at the sub-domain boundaries, minimizing errors associated with parallelization by computing derivatives using the interior compact scheme with global spectral analysis. Time integration is conducted using a fourth-order, four-stage Runge-Kutta method (RK4) with a non-dimensional time step of 2.5 × 1 0 − 6 for all the reported simulations, to ensure all spatial and temporal scales in the flow are captured. This numerical framework has been validated in a prior work (Sundaram et al., 2023) by comparing the pressure distribution on the SHM1 airfoil with experimental results from (Fujino et al., 2003) for the following non-dimensional parameters: R e ∞ = 13.6 × 1 0 6 , M s = 0.62 and α = 0.27 ° . Additionally, the same numerical methods were used Sengupta et al. (2013) to validate against benchmark wind tunnel results reported by Harris (1981) for the NACA0012 airfoil. The large wind tunnel employed (Harris, 1981) mitigated issues associated with testing in transonic wind tunnels (Garbaruk et al., 2003), such as wall interference and three-dimensional effects.

The study includes six test cases, detailed in Table 1, each initiated with free-stream conditions relevant to various design (Cases 1–3) and off-design (Cases 4–6) scenarios for the SHM1 airfoil (Fujino et al., 2003). The simulated free-stream Mach numbers, Reynolds numbers and angles of attack have been chosen carefully to replicate the flow conditions from Fujino et al. (2003). The free-stream temperature and density are set as: T ∞ = 288.15   K and ρ ∞ = 1.2256   k g / m 3 . Case-1 represents conditions for testing the SHM1 in a low-speed wind tunnel. Case-2 corresponds to a climb condition with M ∞ = 0.31 , where low profile drag is desirable. Case-3 reflects a cruise condition for the SHM1 at M ∞ = 0.62 . Cases 4 through 6 operate in the transonic regime, where shock boundary layer interactions occur. In Case-4, a small normal shock is observed in the Schlieren visualizations (Sengupta et al., 2021). Case-5 simulates the drag divergence Mach number, resulting in a strong normal shock on the suction surface (Fujino et al., 2003). Conversely, Case-6 features transonic shock boundary layer interactions that induce separation on the suction surface, leading to a λ -shock and wedge-shaped shock structure (Sengupta and Shandilya, 2024). The R e ∞ and angles of attack, α are selected to replicate conditions from experimental low-speed and transonic wind tunnel tests, as well as flight test data (Fujino et al., 2003). Since Case-5 and 6 have been computed for the same R e ∞ and α , but for different free-stream Mach numbers, the free-stream temperature of air considered in these two simulations is different. Here, Cases 3–6 represent the validation cases which are compared with the corresponding experimental visualizations and load distributions reported by Fujino et al. (2003).

Figure 2 compares the numerically computed time-averaged coefficient of pressure, C p distributions for a design condition ( M ∞ = 0.62) and an off-design condition ( M ∞ = 0.72 ) with experimental flight test data from Fujino et al. (2003). The results show good agreement between the simulated and experimental C p on both the suction and pressure surfaces. The C p distribution reveals a plateau, characteristic of natural laminar airfoils (Somers, 1992). The SHM1 is designed to maintain a favorable pressure gradient on the suction surface up to x / c = 0.42 , followed by a concave pressure recovery, which is seen in Figure 2A for M ∞ = 0.62 . In contrast, the pressure surface exhibits a favorable pressure gradient up to x / c = 0.63 to enhance drag reduction (Fujino et al., 2003). A steeper pressure recovery is evident on the pressure surface under the design condition, as shown in Figure 2A. The simulation also accurately captures the shock location, as depicted in Figure 2B, occurring at x / c ≈ 0.45 on the suction surface (Sengupta et al., 2021) for the off-design transonic condition ( M ∞ = 0.72 ) . The steep pressure rise induced by the shock is frequently observed upstream of the point where the shock interacts with the suction surface, as the shock transmits a “pressure signal” in the upstream direction, within the subsonic inner part of the boundary layer (Délery et al., 1986).

In this section, we compare the numerical Schlieren for the different operating conditions outlined in Table 1, emphasizing the varying flow physics and associated features. We also discuss the dominant time scales for these various design and off-design conditions, along with the aerodynamic performance of the SHM1 airfoil. These insights serve as supplementary benchmark datasets to the experimental data provided by Fujino et al. (2003).

The present study investigates the aerodynamic performance of the SHM1 airfoil (depicted in Figure 1) across various design and off-design conditions, with a focus on understanding flow physics and associated dominant time-scales. For the off-design cases, in particular, the effects of transonic shock boundary layer interactions as a function of increasing free-stream Mach number is also explored. Utilizing dispersion relation preserving highly accurate compact schemes for discretization of the governing compressible Navier-Stokes equations, six implicit large eddy simulations are performed for a range of free-stream Mach numbers. The numerics are validated against the experiments of (Fujino et al., 2003), revealing a good agreement of the computed pressure distribution with that in the flight test in Figure 2. The instantaneous flow features, spectra, pressure distributions, skin friction coefficients, and overall aerodynamic efficiency provide critical insights into the operation of the SHM1 for the following design scenarios: (i) low speed operation, (ii) climb conditions, (iii) cruise conditions; and the following off-design operation points: (i) first appearance of the near-normal transonic shock, (ii) drag-divergence condition, and (iii) shock-induced separation.

The various flow features are traced via density gradients in the numerical Schlieren plots in Figures 3–5. For low-speed operation, vortex shedding from the trailing edge follows the characteristic shedding frequency of Kelvin-Helmholtz eddies (Sengupta et al., 2024). Upon increasing the free-stream Mach number to that in the climb condition, in addition to the trailing edge vortices, upstream propagating pressure pulses are noted in the inviscid part of the flow, termed as Kutta waves (Lee, 2001). These undergo nonlinear interaction with each other and the underlying boundary layer at the cruise condition, leading to a perfect symmetric arrangement on the pressure and suction surfaces of the SHM1 airfoil. For a free-stream Mach number of 0.72, a near-normal shock (Zierep, 2003) is obtained on the suction surface which locally alters the supersonic flow to a subsonic one. At the drag-divergence Mach number, a stronger normal shock with a larger spread across the suction surface is observed. A weaker shock is obtained on the pressure surface also. For the shock-induced separation, λ shock waves are obtained due to merging of the normal shock waves with the oblique ones on the suction surface. Strong, nonlinear shock waves are obtained on the pressure surface, which are expected to contribute detrimentally to the skin friction drag at this operational regime.

The dominant time-scales in the flow are identified by examining the spectrum of vorticity for the six computed cases in Figures 6, 7. While the low-speed and climb conditions reveal distinct peaks in the frequency plane at sub-harmonics of the Kelvin-Helmholtz shedding frequency, the transonic shock boundary layer interactions evoke multi-periodicity, which follows the typical chaotic behavior displayed by other fluid dynamical systems exhibiting transition to turbulence. The associated spectra are chaotic with no discernible peak along the frequency plane. The time-averaged streamwise pressure coefficient computed in Figure 8 demonstrate characteristic plateaus along the suction surface indicative of natural laminar airfoils. Pressure spikes were identified at shock locations for off-design conditions along both suction and pressure surfaces (for free-stream Mach numbers of 0.73 and 0.78). A concave pressure recovery is noted on the pressure surface beyond x / c = 0.63 , matching the design requirements of the SHM1 (Fujino et al., 2003). The steepest spikes in the pressure distribution are noted for the shock-induced separation due to the presence of stronger shocks spread over a wider section of the suction surface. A delayed separation, with stronger shocks are noted with an increase in the free-stream Mach number. The analysis of skin-friction coefficients in Figure 9 revealed flow separation behavior, with shorter separation bubbles at lower Mach numbers and larger, delayed separation regions as free-stream Mach numbers increased. These findings underscore the complex interactions between Kutta waves, the boundary layer and the shock structures.

The aerodynamic performance of the SHM1 airfoil was quantified in Table 2 by computing the integrated coefficients of lift and drag. The ratio of lift and drag coefficient is a metric for the aerodynamic efficiency, which has also been compared for the six test cases. It has been highlighted that Case-6 (shock-induced separation) has the least drag but also suffers from a significant loss of lift. Conversely, Case-5 (drag-divergence) exhibited the highest lift, although this was offset by substantial drag, resulting in low aerodynamic efficiency. The cruise condition emerged as the most aerodynamically efficient, achieving the best ratio of C l to C d .

The present work suggests that the operating conditions of the Honda jet determine the range of load and pressure distribution on the suction and pressure surfaces, which in turn affect the aerodynamic performance. The free stream Mach number, Reynolds number, and angle of attack, are important design parameters which demarcate the range of operation into design and off-design regimes. These are furthermore subdivided into subsonic and transonic flow regimes, with the best aerodynamic efficiency observed at transonic conditions. For the Honda jet, specifically designed as a natural laminar flow airfoil, efficient transonic aerodynamic performance can be achieved with higher C l / C d and minimal flow separation.

In conclusion, this work enhances our understanding of the SHM1 airfoil’s performance across different operating conditions. The results provide valuable benchmarks for future optimization efforts, particularly in managing shock-induced separations and enhancing overall aerodynamic efficiency in transonic flight regimes. The insights gained from this study are essential for improving the design and operational strategies of lightweight business jets and other aircraft utilizing similar airfoil configurations.