On Numerical Simulations of Turbulent Flows Over a Bluff-Body with Aerodynamic Flow Control Based on Trapped Vortex Cells: Effects of Compressibility

1. Introduction

Aerodynamic flow control (AFC) using Trapped Vortex Cells (TVC), strategically distributed across various surfaces of an obstacle, represents a promising alternative technology in the aerospace industry. By harnessing naturally occurring vortical structures and stabilizing them within engineered cavities, TVC systems enable enhanced lift generation, delayed flow separation, and improved overall aerodynamic efficiency – without relying on traditional moving parts or energy-intensive actuators. This semi-passive control approach offers significant potential for applications in aircraft design, unmanned aerial vehicles, and high-performance aerospace platforms, particularly where compactness, reliability, and energy savings are critical. Ongoing research continues to explore the scalability of TVC systems, their integration with advanced materials, and their performance under varying flow regimes, including compressible and transonic conditions.

A historical perspective on the AFC technology based on trapped vortex cell (TVC) is provided in the foundational works of Sedda et al. [1] and Lysenko et al. [2]. The concept of vortex trapping was first introduced by Ringleb [3], and later advanced by Kasper [4], who patented a glider design incorporating a trapped vortex structure to enhance lift. This innovation culminated in the first documented flight experiment using TVC technology – the tailless Kasper Bekas BKB-1 glider [5]. In the early 1990s, Savitsky et al. [6] expanded on these principles with the development and testing of the EKIP aircraft, a blended wing-body configuration equipped with multiple TVCs. Between 1995 and 2005, both experimental and numerical studies on thick airfoils (with chord thicknesses ranging from $40\%$ to $50\%$) confirmed the potential for achieving nearly unseparated flow, leading to enhanced aerodynamic performance and minimal energy losses from flow control systems [7].

Further systematic investigations were conducted on thick airfoils such as MQ1 CIRA and Göttingen (with $37\%$ chord thickness), integrated with TVCs. These efforts were part of the EU’s 6th Framework Programme (VortexCell2050) and research initiatives funded by the Russian Foundation for Basic Research (2005–2009). These studies deepened the understanding of vortex dynamics and flow control mechanisms, revealing critical Mach number thresholds where significant reorganization of flow structures within TVCs occurs – typically accompanied by increased drag and reduced lift. Comprehensive analyses of these phenomena are available in the works of Donelli et al. [8], Sedda et al. [1], Isaev [9], and Lysenko [2,10].

A novel bluff-body configuration featuring aerodynamic flow control (AFC) based on trapped vortex cells (TVCs) was recently introduced, supported by pioneering large-eddy simulations (LES) [2]. The design comprises a half-cylinder main body with diameter and span length D, complemented by two semi-spherical attachments. Despite a degree of abstraction in the AFC implementation, the system effectively maintained an almost fully unseparated flow over the obstacle. This configuration represents a preliminary concept, where the engine-driven air suction mechanism was emulated through boundary conditions imposing fixed negative mass flow rates. Aerodynamic performance, quantified by the lift-to-drag ratio (K), showed a marked improvement: the AFC-equipped bluff body achieved $K=1.5$, in contrast to the baseline case without AFC, which exhibited negative lift and a ratio of $K=-0.9$. Energy losses attributed to the AFC system were estimated at approximately $10\%$ of the total drag force, indicating a favorable trade-off between control effort and aerodynamic gain. Figure 1 illustrates conceptual flow visualization using time-averaged streamlines over bluff bodies with and without AFC, respectively.

The next logical advancement in AFC technology was its integration into a fixed-wing configuration. For this purpose, a Clark Y airfoil was combined with a semicircular cylinder, serving as the fuselage, and outfitted with a system of three trapped vortex cells positioned on its upper surface. This concept is schematically illustrated in Figure 2,a. A numerical study was conducted using LES at baseline flow conditions: Mach number $M=0.1$, Reynolds number $Re=130,000$ and zero angle of attack. Figure 2 presents a schematic overview of the concept and its flow visualizations using the Q-criterion and line integral convolution (LIC) based on the time-averaged velocity field in the central vertical cross-section.

Despite demonstrating strong aerodynamic performance, with a lift-to-drag ratio of approximately $K\approx 7$, the simulations revealed several limitations of the AFC system. Specifically, additional losses were observed within the internal channels, along with the formation of secondary vortex structures inside the vortex cells (Figure 2,c). Another notable design feature identified during modeling was pronounced flow circulation near the lateral surfaces of the vortex cells, as shown in Figure 2,b. This effect is attributed to circulation interference caused by geometric discontinuities between the Clark Y airfoil and the semicircular fuselage. To mitigate these interactions in future designs, the addition of vertical ridges is recommended. This would suppress cross-flow along the aerodynamic surfaces and help maintain independent pressure distributions between the wing and fuselage.

Despite significant improvements in aerodynamic performance, the previously developed flow control system had not yet achieved fully unseparated flow. Building on this foundation, further advancements were made in the work by Lysenko [10], aiming to eliminate flow separation entirely. The bluff body was redesigned as a thick airfoil with a relative thickness of $50\%$ (equivalent to a semi-circular cylinder) and integrated with two trapped vortex cells. To streamline development and reduce interference from secondary vortex structures, the semi-spherical blocks used in earlier designs were removed. This simplification helped avoid complex interactions between unwanted eddies and the primary flow. Numerical simulations were conducted at a free-stream Reynolds number of $Re=130,000$, based on diameter, and a Mach number of $0.05$, with a zero angle of attack.

From a numerical modeling standpoint, both Reynolds-averaged Navier-Stokes (RANS) and large-eddy-simulation (LES) approaches were employed by Lysenko [10]. Each method confirmed the presence of fully unseparated flow over the airfoil equipped with AFC, with aerodynamic coefficient discrepancies limited to approximately $5\%$. Ultimately, the refined design achieved nearly complete unseparated flow, generating a lift force equivalent to roughly $80\%$ of the theoretical maximum value ($2\pi$), and delivering four times the lift-to-drag ratio of the original configuration [2]. Additionally, the influence of viscosity (Reynolds number effects) was examined. Results indicated that the aerodynamic characteristics of the improved concept were largely insensitive to viscosity variations. Across a broad Reynolds number range ($75,000\le Re\le 1,000,000$), a modest increase in lift and aerodynamic efficiency was observed, with the drag coefficient remaining nearly constant and the performance trending linearly toward an asymptotic limit.

This study builds upon previous research [2,10] by numerically investigating the influence of compressibility by varying Mach numbers on the aerodynamic performance of the proposed concept. A series of numerical simulations are conducted using the classical RANS approach for Mach numbers $M=0.05$, $0.1$, $0.2$ and $0.3$, with a fixed Reynolds number of $Re=130,000$. Additionally, the system’s behavior is evaluated across a range of angles of attack, $\alpha =\pm {15}^0$.

The specific goals and objectives of the study are as follows:

• Perform flow simulations around the concept using RANS methodology at Mach numbers $M=0.05$, $0.1$, $0.2$ and $0.3$, $Re=130,000$ and zero angle of attack;
• Conduct cross-platform validation for the limiting case at $M=0.3$, when the flow over airfoil’s upper surface becomes locally transonic, by comparing different numerical solvers for the compressible Navier-Stokes equations – specifically pressure-based, density-based, and hybrid approaches;
• Assess the impact of turbulence modeling by comparing Spalart-Allmaras (SA) and Realizable k-$ϵ$ (RKE) models across the selected Mach number range; analyze the limiting case assuming the fluid inviscid.
• Highlight the significance of accurately specifying thermodynamic properties of air, such as viscosity, heat capacity, and thermal conductivity, and their influence on integral aerodynamic characteristics;
• Analyze the aerodynamic performance of the system at various angles of attack within the range of $\pm {15}^0$;
• Ultimately, the study revealed that the aerodynamic performance of the system remains largely unaffected by Mach numbers up to $M=0.3$, consistently maintaining the fully unseparated flow. Under cruise conditions at $M=0.2$, this behavior is associated with a high lift coefficient of $C_l\approx 6$, which corresponds to approximately $95\%$ of the theoretical maximum value of $2\pi$. The resulting lift-to-drag ratio is around 24.

The article is presented in five parts. The first two sections are devoted to the problem statement and aspects of numerical and mathematical modeling. Then the main results are presented. The fifth part provides a brief discussion and critical review. At the end, the main conclusions are drawn.

2. Problem Statement and Theoretical Considerations

In this section a brief problem statement and key design features of the aerodynamic flow control system are presented. A brief mathematical background of the Reynolds-averaged Navier-Stokes approach accompanied by the Spalart-Allmaras and k-$ϵ$ realizable turbulence models are discussed as well.

The turbulent flow is characterized by the following Reynolds and Mach numbers: $Re={\rho }_{\infty }{U}_{\infty }D/{\mu }_{\infty }$ and $M=U_{\infty }/a_{\infty }$, respectively. Here $\rho$ is the density, U is the far-field velocity, a is the speed of sound, $\mu$ represents the dynamic molecular viscosity, D is the diameter of the thick airfoil (semi-circular cylinder). The subscript ∞ indicates the flow parameters at the inlet boundary. The drag and lift coefficients are defined as $C_d={\displaystyle \frac{2F_D}{{\rho }_{\infty }{U}_{\infty }^{2}D}}$ and $C_l={\displaystyle \frac{2F_L}{{\rho }_{\infty }{U}_{\infty }^{2}D}}$, where $F_D$ and $F_L$ are forces acting in the stream-wise and vertical directions. The pressure coefficient is specified as $C_p=2{\displaystyle \frac{(p-p_{\infty })}{\left({\rho }_{\infty }{U}_{\infty }^2\right)}}$, where p being the static pressure. All numerical simulations were performed at the fixed Reynolds number, $Re=130,000$, and the following Mach numbers, $M=0.05$, $0.1$, $0.2$ and $0.3$.

2.1. Overview of the Aerodynamic Flow Control Concept Based on the Trapped Vortex Cells

The primary design dimensions of the thick airfoil, both with and without an aerodynamic flow control system, are depicted in Figure 3. A bluff-body (BB hereafter) is composed of a main block shaped as a half-cylinder with a diameter D and a span length of D. The diameter of the half-cylinder, D, serves as the linear reference scale. The overall width, height, and span lengths of the bluff body are $L_x=D$, $L_z=0.5\times D$, and $L_y=1\times D$, respectively. BB with aerodynamic flow control incorporates a system of the two trapped vortex cells, with diameters of $0.2\times D$ and $0.12\times D$ respectively, positioned on its suction (upper) side. Each TVC has a span length equal to D. Two vortex cells are connected by an arch-shaped channel. The outlet section of the channel, with a surface area of $0.05\times D\times D$, is used for air suction (Figure 4,a).

The present concept has been redesigned compared to initial study [2], resulting in nearly unseparated flow past the airfoil. It is important to highlight that this system was developed after several dozen optimization iterations. This concept also abstracts the aerodynamic flow control system for air suction by imposing a fixed pressure drop assuming a passive method is employed. As initially reported by Lysenko [10], a pressure drop of $-2,000$ Pa was sufficient to sustain fully unseparated flow over a thick airfoil at a low Mach number of $M=0.05$. However, as the free-stream velocity increases toward Mach $0.3$, a greater static pressure drop, up to $-10,000$ Pa, is required to maintain similar flow behavior. In the current model, the consumed air is withdrawn inside BB, which is not a realistic scenario for real-life applications. However, in a broader context, aircraft dynamics require a feedback control loop in which air suction forcing is dynamically adjusted to achieve optimal system performance. In this context, the present flow control method can be regarded as active. Preliminary studies indicate that increasing the forcing up to $50\%$ allows the flow over the obstacle to remain unseparated while preserving similar large-scale vortex structures in the TVCs. However, this is accompanied by penalties in the form of increased drag force. While further detailed studies are necessary to model realistic air suction using an aircraft engine system, the authors believe that this design serves as a satisfactory proof-of-concept and holds significant potential for further development.

The physics and topology of the flow over BB without trapped vortex cells are qualitatively similar to the flow over a semi-circular cylinder at $Re=50,000$, as discussed in detail previously [11]. At a Reynolds number of $Re=130,000$, the separation of the laminar boundary layers occurs in the sub-critical flow regime. This condition results in complex, nonlinear interactions between near and far wakes characterized by two types of flow instabilities (Figure 4, b). The Kelvin-Helmholtz instability of the separated shear layer, along with vortex shedding (Bénard/von Kármán instability), dominates the wake. Figure 4, c shows the schematic flow over the BB with integrated vortex cells. Unlike the initial study [2], this design completely rearranges the flow over the airfoil, making it predominantly unseparated.

To avoid any confusion, all numerical simulations presented in this work were conducted in two-dimensional ($2D$) space.

2.2. Basic Governing Equations

In this subsection the governing equations for the turbulent compressible flows for the Reynolds-averaged simulations are presented. The Cartesian tensor notation, where t and $x_i$ are independent variables representing time and spatial coordinates of a Cartesian coordinate system $\mathbf{x}$, is used. The three components of the velocity vector $\mathbf{u}$ are denoted $u_i$ ($i=1,2,3$). The summation convection over repeated indices applies unless otherwise is noted.

The Favre-averaged (i.e. mass-density weighted) equations for mass, momentum, and energy for the turbulent compressible flows are:

$${\displaystyle \frac{\partial \overline{\rho }}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }{\tilde{\mathbf{u}}}_j\right)=0,$$

(1)

$${\displaystyle \frac{\partial }{\partial t}}\left(\overline{\rho }{\tilde{\mathbf{u}}}_i\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }{\tilde{\mathbf{u}}}_i{\tilde{\mathbf{u}}}_j\right)=-{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{{\tau }_{ij}}-\overline{\rho }\widetilde{{\mathbf{u}}_i^{\prime \prime }{\mathbf{u}}_j^{\prime \prime }}\right),$$

(2)

$${\displaystyle \frac{\partial }{\partial t}}\left(\overline{\rho }\tilde{h}\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }\tilde{h}{\tilde{\mathbf{u}}}_j\right)=-{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }{\alpha }_t{\displaystyle \frac{\tilde{h}}{\partial x_j}}-\overline{\rho }\widetilde{{\mathbf{u}}_j^{\prime \prime }{h}^{\prime \prime }}\right).$$

(3)

Here, the overbar denotes Reynolds averaging, while the tilde denotes Favre averaging: $\rho$ denotes the density, $\mathbf{u}$ is the velocity, p is the pressure, $h=\int C_{p}dT$ is the enthalpy, T is the temperature, $C_p$ is the heat capacity and ${\alpha }_t$ is the thermal diffusivity.

Here, $C_p$ is calculated as a function of temperature from a set of coefficients taken from NIST-JANAF thermochemical tables [12] or technical report NASA/TP-2002-211556 [13]. The thermal diffusivity is modeled as

$$\overline{\rho }{\alpha }_t=\mu {\displaystyle \frac{C_v}{C_{pBar}}}\left(1.32+1.77{\displaystyle \frac{R}{C_v}}\right),$$

where $\mu$ is the molecular viscosity (calculated according to the Sutherland’s law), R is the gas constant, $C_v$ represents heat capacity at constant volume. $C_{pBar}$ is the specific heat averaged over the temperature interval from the reference temperature to the actual mean temperature, $C_{pBar}={\displaystyle \frac{\Delta h}{\Delta T}}$, where $\Delta T=\tilde{T}-T_{ref}$ and $\Delta h=\tilde{h}\left(T\right)-\tilde{h}\left(T_{ref}\right)$. It worth noting when equation of state is resented by the ideal gas law, the alternative method based on the kinetic theory can be applied to compute air thermodynamic properties (i.e., viscosity, thermal conductivity and specific heat).

The turbulence flux $\widetilde{{\mathbf{u}}_j^{\prime \prime }{h}^{\prime \prime }}$ is derived according to the gradient hypothesis

$$-\overline{\rho }\widetilde{{\mathbf{u}}_j^{\prime \prime }{h}^{\prime \prime }}\approx {\displaystyle \frac{{\mu }_t}{{Pr}_t}}{\displaystyle \frac{\tilde{h}}{\partial x_j}},$$

where ${\mu }_t$ is the turbulence viscosity and ${Pr}_t$ is a turbulence Prandtl number (here ${Pr}_t=0.8$).

The stress tensor for a Newtonian fluid is expressed as

$$\overline{{\tau }_{ij}}=\mu \left({\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\mu {\delta }_{ij}{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_k}{\partial x_k}}.$$

The Reynolds stresses are modeled according to

$$\overline{\rho }\widetilde{{\mathbf{u}}_i^{\prime \prime }{\mathbf{u}}_j^{\prime \prime }}=-{\mu }_t\left({\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_j}{\partial x_i}}-{\displaystyle \frac{2}{3}}{\delta }_{ij}{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_k}{\partial x_k}}\right)+{\displaystyle \frac{2}{3}}\overline{\rho }\tilde{k}{\delta }_{ij}.$$

In order to close Eqs. 3-3 the model for the turbulence viscosity shall be provided. In the present study a one equation model by Spalart-Allmaras [14] and the realizable k-$ϵ$ turbulence model of Shih [15] are employed.

2.3. The Spalart-Allmaras Turbulence Model

The Spalart-Allmaras turbulence model [14] (SA here and after) solves a modeled transport equation for the kinematic eddy (turbulent) viscosity ($\tilde{\nu }$) as

$${\displaystyle \frac{\partial }{\partial t}}\left(\overline{\rho }\tilde{\nu }\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }\tilde{\nu }{\tilde{\mathbf{u}}}_j\right)={\displaystyle \frac{1}{{\sigma }_{\tilde{\nu }}}}\left[{\displaystyle \frac{\partial }{\partial x_j}}\left(\left(\mu +\overline{\rho }\tilde{\nu }\right){\displaystyle \frac{\partial \tilde{\nu }}{\partial x_j}}\right)+C_{b2}\overline{\rho }{\left({\displaystyle \frac{\partial \tilde{\nu }}{\partial x_j}}\right)}^2\right]+G_{\nu }-Y_{\nu },$$

(4)

where $G_{\nu }$ and $Y_{\nu }$ are production and destruction of the turbulence viscosity that occurs in the near-wall region due to wall blocking and viscous damping, respectively. ${\sigma }_{\tilde{\nu }}=2/3$ and $C_{b2}=0.622$ are the constants.

The turbulent viscosity, ${\mu }_t$, is computed as follows:

$${\mu }_t=\overline{\rho }\tilde{\nu }{f}_{\nu 1},f_{\nu 1}={\displaystyle \frac{{\chi }^3}{{\chi }^3+C_{\nu 1}^3}},$$

where $\chi \equiv \tilde{\nu }/\nu$, $C_{\nu 1}=7.1$ – constant and $\nu$ is the kinematic molecular viscosity.

The production term, $G_{\nu }$, is modeled as

$$G_{\nu }=C_{b1}\overline{\rho }\tilde{S}\tilde{\nu },\tilde{S}\equiv \widehat{S}+{\displaystyle \frac{\tilde{\nu }}{{\kappa }^2{d}^2}}{f}_{\nu 2},f_{\nu 2}=1-{\displaystyle \frac{\chi }{1+\chi f_{\nu 1}}},$$

where d is the distance from the wall, $\kappa =0.41$ and $C_{b1}=0.1355$ are the constants and $\widehat{S}$ is a scalar to measure of the deformation tensor. In the present study, $\widehat{S}$ is computed according to the modification proposed by Dacles-Mariani et al. [16] to combine both vorticity and the strain tensors:

$$\widehat{S}\equiv \left|{\Omega }_{ij}\right|+2min\left(0,\sqrt{2S_{ij}{S}_{ij}}-\left|{\Omega }_{ij}\right|,\right)$$

where

$$\left|{\Omega }_{ij}\right|=\sqrt{2{\Omega }_{ij}{\Omega }_{ij}},S_{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_j}{\partial x_i}}\right),{\Omega }_{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_i}{\partial x_j}}-{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_j}{\partial x_i}}\right).$$

The destruction term, $Y_{\nu }$ is modeled as

$$Y_{\nu }=C_{w1}\overline{\rho }{f}_w{\left({\displaystyle \frac{\tilde{\nu }}{d}}\right)}^2,f_w=g{\left[{\displaystyle \frac{1+C_{w3}^6}{g^6+C_{w3}^6}}\right]}^{1/6},g=r+C_{w2}\left(r^6-r\right),r\equiv {\displaystyle \frac{\tilde{\nu }}{\widehat{S}{\kappa }^2{d}^2}}.$$

Here, $C_{w1}={\displaystyle \frac{C_{b1}}{{\kappa }^2}}+{\displaystyle \frac{1+C_{b2}}{{\sigma }_{\tilde{\nu }}}}$, $C_{w2}=0.3$, $C_{w3}=2.0$ are model constants.

2.4. The Realizable k-$ϵ$ Turbulence Model

The realizable k-$ϵ$ turbulence model of Shih [15] (RKE here and after) is based on the turbulence kinetic energy ($\tilde{k}$) and its dissipation rate ($\widetilde{ϵ}$), where the turbulence viscosity is defined as

$${\mu }_t=C_{\mu }\overline{\rho }{\tilde{k}}^2/\widetilde{ϵ}.$$

(5)

The modeled transport equations are:

$${\displaystyle \frac{\partial }{\partial t}}\left(\overline{\rho }\tilde{k}\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }\tilde{k}{\tilde{\mathbf{u}}}_j\right)={\displaystyle \frac{\partial }{\partial x_j}}\left(\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right){\displaystyle \frac{\tilde{k}}{\partial x_j}}\right)+G-\overline{\rho }\widetilde{ϵ}-Y_M,$$

(6)

$${\displaystyle \frac{\partial }{\partial t}}\left(\overline{\rho }\widetilde{ϵ}\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\overline{\rho }\widetilde{ϵ}{\tilde{\mathbf{u}}}_j\right)={\displaystyle \frac{\partial }{\partial x_j}}\left(\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_{ϵ}}}\right){\displaystyle \frac{\widetilde{ϵ}}{\partial x_j}}\right)+\overline{\rho }{C}_{1}S\widetilde{ϵ}-\overline{\rho }{C}_2{\displaystyle \frac{{\widetilde{ϵ}}^2}{\tilde{k}+\sqrt{\nu \widetilde{ϵ}}}},$$

(7)

where

$$C_1=max\left[0.43,{\displaystyle \frac{\eta }{\eta +5}}\right],\eta =S{\displaystyle \frac{\tilde{k}}{\widetilde{ϵ}}},S=\sqrt{2S_{ij}{S}_{ij}}.$$

Here, G represents production of the turbulence kinetic energy due to the mean velocity gradients, $Y_M$ is contribution of the fluctuating dilatation in compressible turbulence, $C_2=1.9$ is a constant and ${\sigma }_k=1.0$ and ${\sigma }_{ϵ}=1.2$ are the turbulence Prandtl numbers for $\tilde{k}$ and $\widetilde{ϵ}$, respectively.

The production of the turbulence kinetic energy is modeled as

$$G=-\overline{\rho }\widetilde{{\mathbf{u}}_i^{\prime \prime }{\mathbf{u}}_j^{\prime \prime }}{\displaystyle \frac{\partial {\tilde{\mathbf{u}}}_j}{\partial x_i}}={\mu }_t{S}^2.$$

The coefficient $C_{\mu }$ required to compute the turbulence viscosity in Eq. 5 is no longer constant and defined as the following equation:

$$C_{\mu }={\displaystyle \frac{1}{A_0+A_s\frac{\tilde{k}{U}^{*}}{\widetilde{ϵ}}}},U^{*}\equiv \sqrt{S_{ij}{S}_{ij}+{\tilde{\Omega }}_{ij}{\tilde{\Omega }}_{ij}}.$$

Under some assumptions the term ${\tilde{\Omega }}_{ij}$ can be computed as the mean rate-of-rotation tensor viewed in a moving reference frame with the angular velocity ${\omega }_k$, ${\tilde{\Omega }}_{ij}\approx {\overline{\Omega }}_{ij}$. The model constants are $A_0=4.04$ and $A_s=\sqrt{6}cos\varphi$, where

$$\varphi ={\displaystyle \frac{1}{3}}{cos}^{-1}\left(\sqrt{6}W\right),W={\displaystyle \frac{S_{ij}{S}_{jk}{S}_{ki}}{{\tilde{S}}^3}},\tilde{S}=\sqrt{S_{ij}{S}_{ij}}.$$

Effects of compressibility for high-Mach-number flows in the k-$ϵ$ model is accounted by including the dilatation dissipation term $Y_M$, which is modeled according to proposal by Sarkar [17]:

$$Y_M=2\overline{\rho }\widetilde{ϵ}{M}_t^2,M_t=\sqrt{{\displaystyle \frac{\tilde{k}}{a^2}}},$$

where $M_t$ represents the turbulent Mach number and a is the speed of sound.

2.5. Basic Governing Equations in Conservative Form

Sometimes, it is useful to provide the governing equations in the conservation form:

$${\displaystyle \frac{\partial \mathbf{Q}}{\partial t}}+{\displaystyle \frac{\partial {\mathbf{F}}_j}{\partial x_j}}+{\displaystyle \frac{\partial {\mathbf{V}}_j}{\partial x_j}}=0,$$

(8)

where $\mathbf{Q}$ is the vector of conservative variables, ${\mathbf{F}}_j$ and ${\mathbf{V}}_j$ are the convective and viscous fluxes vectors, respectively:

$$\mathbf{Q}=\left[\begin{array}{c}\rho \\ \rho {\mathbf{u}}_i\\ \rho E\end{array}\right],{\mathbf{F}}_j=\left[\begin{array}{c}\rho {\mathbf{u}}_j\\ \rho {\mathbf{u}}_i{\mathbf{u}}_j+{\delta }_{ij}p\\ \rho {\mathbf{u}}_{i}H\end{array}\right],{\mathbf{V}}_j=\left[\begin{array}{c}0\\ -{\tau }_{ij}\\ -{\tau }_{ik}{\mathbf{u}}_k+{\mathbf{q}}_i\end{array}\right].$$

Here, the stare equations are defined as follow: $E=e+\rho {\displaystyle \frac{1}{2}}{\mathbf{u}}_j^2$ represents the total energy per unit mass, $e=RT/(\gamma -1)$ is the internal energy per unit mass, $H=E+p/\rho$ is the total enthalpy, R is the gas constant, $\gamma$ is the specific heat ratio, ${\tau }_{ij}$ is the viscous stress tensor and ${\mathbf{q}}_i$ is the heat flux vector.

3. Brief Aspects of Numerical Simulations

This section provides a concise overview of the numerical method, along with descriptions of the computational domain, grids, boundary and initial conditions.

3.1. Overview of the Numerical Methodology

Several numerical approaches to solve the Favre-averaged Navier-Stokes (RANS) equations are used in the present study. Two of them are so-called pressure- and density-based solvers implemented in the CFD platform Ansys Fluent 2021R1 [18] (hereinafter, AF). And an alternative so-called hybrid solver implemented in the open-source OpenFOAM [19] (hereinafter, OF) framework. The term ’solver’ here refers to a generalized numerical method used to solve the fluid dynamics equations.

• Pressure-based approach belongs to a general class of methods called the projection method [20,21], where the continuity of the velocity filed is achieved by solving pressure-correction equation. The pressure-based solver is implemented as a finite volume method (FVM) utilizing a second-order upwind scheme (SOU) for all convective terms and an implicit Euler method (bounded BDF-2) for time integration [22]. The time integration step is chosen to ensure that the local Courant number is less than ten, $CFL<10$. The so-called pressure-based coupled algorithm is used, which solves a coupled systems of equations comprising the momentum and pressure-based continuity equations. The system of linear algebraic equations is solved using the algebraic multigrid method (AMG) accompanied by the additive correction strategy [23] and the classical iterative Gauss-Seidel procedure.
• Density-based approach is a Godunov-type FVM solver. In his classical paper of 1959, Godunov [24] presented a conservative extension of the first–order upwind scheme to non-linear systems of hyperbolic conservation laws. The key ingredient of the scheme is the solution of the Riemann problem [25], which is approximated by the advection upstream splitting method (AUSM) [26] in the present study. Other numerical aspects like time integration, discretization schemes and linear solvers were set up in the same spirit as for the pressure-based approach.
• Hybrid approach is alternative method to integrate the mass, momentum and energy conservative equations by coupling the PIMPLE (combination of PISO/SIMPLE) algorithm [27] with a corresponding Kurganov-Noelle-Petrova scheme (KNP) [28] for non-oscillating discretization of the convective and diffusive terms. The hybrid method based on the coupled PIMPLE algorithm and KNP scheme meets the requirements on monotonicity and allows to simulate turbulent flows in a wide range of a Mach number, $0\le M\le 3$. The hybrid approach was implemented as a so-called `pimpleCentralFoam’ solver based on FVM supported by the OpenFOAM platform [19] with specific details provided by Krpaposhin et al. [29]. A stabilized preconditioned (bi)-conjugate gradients method (PBiCGStab) accompanied by the diagonal incomplete-LU preconditioner was used to solve the linear systems of the algebraic equations with the absolute tolerance of ${10}^{-12}$. A limited advection scheme, so called `vanLeer’ [30], which belong to a family of TVD schemes, was utilized to approximate convective and diffusive fluxes computed by KNP. Integration in time was handled by the first-order implicit Euler method with a restriction of $CFL<0.2$.

3.2. Computational Grids

The planar, unstructured, hexahedral-dominated baseline grid (designated as A0) was constructed as follows. The computational domain was defined as a rectangular block with dimensions $25D\times 25D$, initially seeded with $60\times 60$ nodes along the x and y directions. To enhance resolution, the region spanning $y/D=-4$ to $y/D=4$ was adapted using a $2\times 2\times 2$ refinement coefficient. Within this, the central block from $y/D=-2.5$ to $y/D=2.5$ was further refined to a uniform cell size of $0.025/D$. At the core of the domain, a cylindrical region with a radius of $0.13/D$, centered at the origin of the Cartesian coordinate system, was meshed with a fine uniform cell size of $0.00125/D$. A viscous sub-grid was applied around the obstacle, featuring an average wall-normal spacing of $y^{+}\simeq 1$, an expansion factor of $1.1$, and a total of five layers. For the second grid (referred to as A1), the cylindrical region of radius $0.13/D$ was re-meshed with an even finer uniform cell size of $0.000833/D$. The finest grid (referred to as A2) was then generated by adapting the cylindrical region of the A1 grid using an additional $2\times 2$ refinement. Details of the A0 and A3 grids are illustrated in Figure 5.

3.3. A Near-Wall Modeling Treatment

A wall treatment combining a two-layer model with enhanced wall functions is used for the realizable k-$ϵ$ turbulence model [18]. This method maintains the accuracy of the standard two-layer approach for fine near-wall meshes while minimizing accuracy loss for wall-function meshes. The two-layer approach assumes that the computational domain is divided into viscosity-affected and fully turbulent regions. The demarcation is determined by a wall-distance-based turbulent Reynolds number, ${Re}_y$ with the wall-normal distance, y, calculated at the cell centers. In the fully turbulent region (${Re}_y$>200), the standard RKE model is utilized. In the viscosity-affected near-wall region, the one-equation model by Wolfshtein [31] is applied, where the computation of turbulent viscosity is modified by introducing a length scale parameter proposed by Chen and Patel [32]. Finally, the blending between the two regions is computed using the method proposed by Jongen [33].

The Spalart-Allmaras turbulence model implemented in AF has been extended with a $y^{+}$ - insensitive wall treatment, which allows to apply model independent of the near-wall $y^{+}$ resolution and blends automatically from a viscous sublayer formulation to a logarithmic formulation [18]. At the same time on the intermediate grids, $1<y^{+}<30$, its continues to handle integrity and provides consistent wall shear stresses. For the hybrid solver implemented in OF, the SA turbulence model was used accompanied by the `nutUBlendedWallFunction’ boundary condition, which provides a wall constraint on the turbulent viscosity, i.e. $\tilde{\nu }$, based on velocity. In this case a binomial-type wall-function blending method was applied to approximate $\tilde{\nu }$ between the viscous and inertial sublayers [34].

3.4. Boundary and Initial Conditions

Boundary conditions were defined to ensure that the free-stream parameters matched the target Reynolds and Mach numbers, set at $Re=130,000$ and $M=0.05-0.3$, respectively. Table 1 outlines the free-stream values used to determine the Mach number at the inlet boundary, adjusted according to bluff body altitude to reflect atmospheric variations with elevation above sea level (h). Notably, the Reynolds number remained constant across all Mach number cases. Compressibility effects were modeled using the ideal gas law. Thermodynamic properties like dynamic viscosity, thermal conductivity, and heat capacity were computed either via kinetic gas theory or temperature-dependent polynomial functions. Free-stream turbulence was characterized by a turbulence intensity of $5\%$ and a hydraulic diameter equal to D. At the inlet boundary, fixed values were assigned for free-stream velocity and total temperature, while a fixed static pressure was imposed at the outlet. The body surface was treated as adiabatic with no-slip wall conditions. Air suction from the vortex cell system was simulated by applying a static pressure drop ranging from $-2,000$ to $-11,000$ Pa, depending on the Mach number. The initial conditions corresponded to a sudden stop of the body in the fluid flow, meaning the inlet boundary conditions were initially extended to the entire computational domain.

3.5. Critical Remark on RANS Methodology Validation and Verification

It is worth noting that the development of the AFC technology in previous studies [2,10] primarily focused on large-eddy simulation (LES), with the classical RANS approach serving as a supplementary tool. A parallel objective was the rigorous evaluation of the numerical platform’s capability to predict turbulent separated flows, encompassing both internal and external aerodynamic configurations. To support this effort, the subgrid-scale incompressibility assumption was frequently employed, restricting free-stream Mach numbers to $M=0.3$. The Reynolds number was similarly constrained to $Re=130,000$, under the assumption that the boundary layer along bluff body walls remains laminar. Key reference cases for AFC development included flow studies over semicircular [11] and circular [35] cylinders, conducted at $Re=50,000$ and $Re=130,000$, respectively. These cases provided foundational insights into AFC performance under the above flow conditions.

As an initial step in further investigating the Mach-number dependence of the AFC technology, this paper focuses on the application of the RANS methodology. It is important to note that over the past two decades, extensive research has been performed to validate and verify the RANS approach in the context of turbulent separated flows. Several fundamental studies addressing critical components of the evolving AFC framework are particularly noteworthy:

• Isaev et al. [36,37] examined the lid-driven cavity flow using an incompressible formulation. The analogy to trapped vortex cells (TVCs) is evident, as the primary vortex core can be considered inviscid with constant vorticity, closely linked to boundary layer development along channel walls. Additionally, assuming a steady-state flow regime is reasonable for moderate to high Reynolds numbers in this configuration.
• Circular and semicircular cylinders serving as baseline geometries for AFC development were investigated by Isaev et al. [38,39] and Lysenko et al. [40,41], providing key insights into flow control mechanisms over curved surfaces.
• Mach 3 wind channel with a step, modeled under the inviscid assumption, was studied by Isaev and Lysenko [42,43], contributing to the understanding of high-speed flow behavior in AFC-relevant geometries.

One of the remaining challenges not previously addressed is the formation of local transonic flow ($M\approx 1$) and the associated interaction of local shock waves with laminar or turbulent boundary layers. To investigate these phenomena, this paper extends the verification of numerical methods to the following sub-problem:

• Shock-boundary layer interaction (SBLI): Verification was performed using the well-documented experiment by Garnier [44], conducted in a supersonic wind tunnel at a nominal Mach number of $M=2.3$. An oblique shock is generated by a deflector mounted on the upper wall of the tunnel at a flow deflection angle of $8^0$, which reflects off the lower wall. At this deflection angle, boundary layer separation is observed.

To enhance the readability of the article, numerical results for the SBLI problem is presented in Appendix A. Several critical observations should be noted. This test case was modeled in a two-dimensional space, despite the inherently three-dimensional nature of the corresponding physical experiment. The flow poses significant challenges for the RANS methodologies, which rely on the assumption of equilibrium turbulence and lack mechanisms to capture low-frequency instabilities. As a result, RANS models are unable to resolve phenomena such as the Kelvin-Helmholtz instability that arises in shear layers during boundary layer separation from bluff body surfaces. Predicting the behavior of large recirculation zones (both attached and detached) remains a known limitation of RANS [45]. To assess numerical robustness, three algorithms for solving the compressible Navier-Stokes equations were tested within RANS: pressure-based, density-based, and hybrid solvers. These were implemented using two state-of-the-art computational platforms: a commercial solver, Ansys Fluent, and an open-source alternative, OpenFOAM. Additionally, two turbulence models, Spalart-Allmaras and realizable k-$ϵ$, were evaluated without the use of wall functions. Overall, this study demonstrated satisfactory agreement between numerical simulations and experimental data. Moreover, the convergence of results across different solvers and turbulence models suggests consistent performance within the bounds of their respective assumptions.

4. Results

The results section is organized into three distinct parts. The first part presents a critical commentary on the role of thermodynamic parameter modeling in accurately predicting integral aerodynamic characteristics. The second part discusses the primary findings obtained at Mach numbers of $0.05$, $0.1$, $0.2$, and $0.3$ under zero angle of attack conditions. The final part focuses on numerical simulation results for varying angles of attack: $\alpha =\pm 5^0$, $\pm {10}^0$ and $\pm {15}^0$.

4.1. Critical Remark on modeling Aspects of Thermodynamical Properties

It is important to note that the results reported in the previous study [10] were obtained using a numerical setup in which thermodynamic properties such as dynamic viscosity, specific heat, and thermal conductivity were treated as constants. While this assumption may be acceptable for low-Mach-number simulations where the flow is considered `incompressible’ ($M<0.1$), it becomes less valid under conditions involving significant local acceleration. In the present configuration, featuring a thick airfoil based on a semicircular cylinder equipped with AFC, the flow experiences substantial acceleration on the upper surface due to suction effects. For an incoming flow at $M<0.05$, the local Mach number on the upper surface reaches approximately $0.15$, and within the TVC channel, it increases to around $0.2$. Under these conditions, assuming constant thermodynamic properties may lead to overly conservative estimates.

Figure 6 illustrates isocontours of Mach number, density, dynamic viscosity, and thermal conductivity for flow over the thick airfoil with integrated TVCs at $M_{\infty }=0.05$. The resulting aerodynamic coefficients differ significantly, by approximately $25\%$, from those obtained using constant-property assumptions: $C_d=0.57$ and $C_l=5.67$ versus $C_d=0.82$ and $C_l=4.86$ (run AFC-RANS-A1 by Lysenko [10]), respectively. Therefore, all subsequent results presented in this study incorporate variable thermodynamic properties of air, modeled either through kinetic gas theory or temperature-dependent polynomial functions.

4.2. Effects of the Mach Number

The main results are summarized in Table 2 and Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. Figure 7 illustrates the dependence of the lift coefficient ($C_l$) and the lift-to-drag ratio ($K=C_l/C_d$)) on Mach number. The iso-contours of the pressure coefficient and its Mach number dependence, shown in Figure 8, reveal a pronounced rarefaction region above the upper surface of the airfoil and within the vortex cells. This rarefaction physically explains the observed increase in lift. As the Mach number increases from $0.05$ to $0.2$, both $C_l$ and K show a clear upward trend. Specifically, the lift-to-drag ratio rises from 10 at $M=0.05$ to a peak value of 24 at $M=0.2$. Notably, the results obtained using two different turbulence models, SA and RKE, exhibit strong agreement, both qualitatively and quantitatively. Figure 9, Figure 10 and Figure 11 present LIC visualizations of the flow over a thick airfoil equipped with AFC across Mach numbers $0.05$, $0.1$, $0.2$, and $0.3$. In all cases, the flow remains fully unseparated, with no evidence of vortex structure breakdown within the AFC cells. Figure 11 displays isocontours of vorticity magnitude, revealing a pronounced increase in vorticity intensity with rising Mach number.

The developed thick airfoil design, featuring two trapped vortex cells, demonstrates effective scalability up to a Mach number of $0.2$ and, as shown in previous work [10], exhibits Reynolds number invariance. At $M=0.2$, the configuration achieves a maximum lift-to-drag ratio of 24. A key characteristic of the design is the pronounced flow acceleration over the upper surface, driven by the air suction system integrated through the vortex cells. Despite the acceleration, the local flow remains subsonic, with Mach numbers ranging between $0.5$ and $0.6$ (see Figure 10). RANS simulations confirm the absence of secondary vortices or separation zones, both within the suction system and inside the vortex cells, indicating stable and efficient flow control performance.

At a Mach number of $0.3$ (here referred to as the critical value) both the lift coefficient and the lift-to-drag ratio exhibit a slight decline. This reduction is attributed to pronounced local flow acceleration over the airfoil’s upper surface, where transonic conditions emerge with local Mach numbers reaching approximately $1.0-1.1$, as shown in Figure 10. Despite these conditions, the flow over the airfoil equipped with AFC remains fully unseparated yet (Figure 9). The trapped vortex cells continue to facilitate air circulation and effectively drive flow through the AFC system. Notably, the performance of the air suction mechanism approaches its operational limit, with the outlet flow nearing sonic conditions ($M\approx 1$). Identifying this critical Mach number represents a key milestone in the ongoing development of the AFC-integrated airfoil design, offering valuable insight into its aerodynamic limits and scalability.

4.3. Effects of Angle of Attack

Further development of the concept includes investigating its aerodynamic characteristics across a range of angles of attack. Numerical simulations were conducted at the base Reynolds number of $Re=130,000$ and Mach number $M=0.2$. The angle of attack was varied between $-{15}^0$ and $+{15}^0$ in increments of 5 degrees, with the bluff-body rotating about an axis passing through the geometric center of the coordinate system (as shown in Figure 3). The key results, summarized in Table 2 (cases R19-R24) and Figure 12 and Figure 13, are compared to the baseline case R06 where the angle of attack is zero. It is evident that as the angle of attack increases from $0^0$ to ${15}^0$, the lift coefficient $C_l$ rises, reaching a maximum value of $7.57$ at $\alpha ={15}^0$. Throughout this range, the flow remains essentially unseparated, except for a small attached vortex near the leading edge on the upper surface. Conversely, when the angle of attack decreases from $0^0$ to $-{15}^0$, the lift coefficient drops to a minimum of $4.82$. Despite this reduction, the flow over the upper surface remains unseparated, and the loss in lift is primarily due to the growth of an attached separation zone on the lower surface. Notably, even at $\alpha =-{15}^0$, the system maintains operational stability without any abrupt decline in the lift-to-drag ratio, when the lift coefficient decreases by approximately $20\%$, indicating robust aerodynamic performance under adverse angles.

5. Discussion

This section presents several key observations, focusing on grid convergence behavior, cross-platform validation of results obtained using different numerical algorithms and a brief analysis of the influence of turbulence models. Finally, the findings of this study are compared with the recent work by Isaev et al. [46], who reported RANS results using the k-$\omega$ SST turbulence model for an alternative configuration involving a thick airfoil equipped with two vortex cells.

5.1. Impact of Numerical Methods and Turbulence Models

To support cross-platform validation of the results, this study examines the influence of various numerical algorithms used to solve the Navier-Stokes equations. The following approaches were employed: a pressure-based solver, classified under generalized projection methods; a Godunov-type (density-based) solver; and a hybrid method combining the PIMPLE algorithm with the KLP scheme. A detailed overview of these solvers is provided in Sec. Section 3.1. The analysis focuses on the limiting case of flow over a bluff body equipped with the AFC system at $M=0.3$, where local transonic conditions develop on the airfoil’s upper surface. The impact of turbulence modeling was also assessed by comparing results obtained using the Spalart-Allmaras (SA) and realizable k-$ϵ$ (RKE) models, alongside a limiting case assuming inviscid flow.

Figure 14, Figure 15, along with Table 2, summarize the key findings. All numerical solvers produced a consistent flow pattern over the AFC-equipped airfoil, characterized by the fully unseparated flow, high lift generation and strong aerodynamic performance. Minor discrepancies were noted between SA and RKE turbulence models, particularly in the vortex structure within the second trapped vortex cell. These differences are reflected quantitatively in the lift coefficient: $C_l\approx 0.56$ for RKE versus $C_l\approx 6.2$ for SA.

A comparison of local flow parameters within the primary trapped vortex cell (TVC) reveals key insights. Figure 15 illustrates the distributions of velocity and vorticity magnitudes along two orthogonal directions intersecting the vortex cell center. Due to the internal air circulation, the velocity magnitude at the cell center approaches zero, with the velocity deficit gradually recovering near the cell walls. All curves in Figure 15a,b exhibit satisfactory collapse. In the inviscid case, velocity recovery occurs more rapidly, attributed to the absence of viscous friction. The vorticity magnitude remains nearly constant at the cell center and begins to vary as the flow approaches the channel boundaries. Convergence is satisfactory across all cases, except for the inviscid scenario, where a more intense vorticity core is observed.

5.2. Remark on the Grid Independence Study

To assess grid sensitivity, the most demanding case at $M=0.3$ was selected for analysis. RANS simulations were performed using both Spalart-Allmaras and Realizable k-$ϵ$ turbulence models across three grids: A0, A1, and A2. Figure 16 shows the probability density functions of $y^{+}$ for each grid, clearly illustrating a shift in the mean $y^{+}$ value corresponding to grid refinement. The primary aerodynamic coefficients computed for $M=0.3$ exhibited strong convergence, with relative differences below $3\%$ across the grids (see Table 2). Additionally, the vorticity magnitude $\Omega$ was sampled at the center of the largest trapped vortex cell (TVC) and summarized in Table 3. These findings confirm that grid independence was satisfactory.

5.3. Comparison of Present Results with an Alternative Configuration

Recently, Isaev et al. [46] presented computational results for an alternative thick airfoil configuration featuring two trapped vortex cells, but employing a fundamentally different channel system for air suction. A schematic comparison of the two configurations is provided in Figure 7c-d. Their simulations were conducted using the classical RANS approach with the modified k-$\omega$ shear stress transport (SST hereafter) turbulence model [47], implemented via the VP2/3 multi-block CFD framework. The calculations were performed at a Reynolds number of ${10}^5$ and across a Mach number range from $0.01$ to $0.4$.

Results by Isaev et al. [46] have been incorporated into Table 2. Figure 7 presents a comparative analysis of the lift coefficient and lift-to-drag ratio distributions obtained in the present study and those reported by Isaev et al. [46]. A qualitative agreement is observed for both parameters. Notably, the lift coefficient differs by approximately $15\%$ between the two AFC configurations, despite the use of distinct computational methodologies. The lift coefficient also exhibits a consistent trend with increasing Mach number up to the value of $M=0.2$. For lower Mach numbers ($M\le 0.05$), the lift-to-drag ratios are nearly identical across both configurations. However, beginning at $M=0.1$, discrepancies in the lift-to-drag ratio become more pronounced, primarily due to variations in the computed drag coefficients.

Overall, the study demonstrates that incorporating an air suction system is essential for the effective operation of a two-cell flow control (AFC) configuration. In both examined cases (Figure 7,c and d), the feasibility of maintaining fully unseparated flow over a thick airfoil while achieving high lift coefficients and enhanced aerodynamic efficiency is clearly established. These findings are supported by numerical results obtained from three independent CFD platforms (Ansys Fluent, OpenFoam and VP2/3), employing diverse algorithms for solving the Navier-Stokes equations and a range of turbulence modeling approaches, including SA, RKE, SST, and LES.

6. Conclusions

This paper explores a novel aerodynamic flow control (AFC) concept featuring two trapped vortex cells integrated into a thick airfoil modeled as a semi-circular cylinder. Flow regimes were investigated at Mach numbers of $0.05$, $0.1$, $0.2$, and $0.3$, with the Reynolds number fixed at $Re=130,000$.

A previous study [10] focused on the influence of viscosity (specifically the Reynolds number) at a near-zero Mach number ($M=0.05$). It demonstrated that across a broad Reynolds number range (${75}^4\le Re\le {10}^6$), the AFC system maintains flow stability, characterized by fully unseparated flow and a lift coefficient reaching approximately $80\%$ of the theoretical maximum ($2\pi$). Numerical results obtained using both RANS and LES approaches showed good agreement, with discrepancies within $5\%$, and did not indicate the formation of secondary vortex structures inside the trapped vortex cells.

The present results obtained for Mach numbers of $0.05$, $0.1$, $0.2$, and $0.3$ demonstrate that the developed AFC concept remains robust against compressibility effects. Aerodynamic efficiency improves with increasing Mach number, rising from $K=10$ at $M=0.05$ to $K=24$ at $M=0.2$. The lift coefficient also increases, reaching approximately $90\%$ of the theoretical maximum ($2\pi$) at $M=0.2$. Across all cases, the flow over the AFC-equipped bluff body remains fully unseparated. From a practical standpoint, $M=0.2$ appears optimal for selecting the cruising flight mode of the proposed system, as the flow remains subsonic while maximizing aerodynamic performance. In contrast, $M=0.3$ represents a critical threshold for the two vortex cell design, where localized transonic flow emerges on the airfoil’s upper surface. Notably, even under these conditions, the flow remains unseparated. The limitations observed at $M=0.3$ are primarily attributed to the air suction system, specifically, the restricted area of the suction slot. Additionally, the high volume of air being extracted leads to significant flow acceleration over the upper surface of the semi-circular airfoil, resulting in localized transonic conditions.

The sensitivity of the AFC concept to angle of attack was investigated over a range of $\pm {15}^0$ in $5^0$ increments. Across all tested conditions, the system remained stable, consistently delivering high lift and favorable lift-to-drag ratios. Peak aerodynamic performance was observed at a positive angle of attack of $5^0$, yielding a lift coefficient of $C_l=6.64$ and a lift-to-drag ratio of $K=35$. At negative angles of attack, a substantial attached recirculation zone developed along the flat lower surface of the airfoil, reducing the lift-to-drag ratio to approximately 20. This behavior is attributed to the airfoil’s geometric configuration – a semicircular cylinder with a flat underside. Further optimization of the lower surface geometry holds promise for mitigating or eliminating the formation of this recirculation zone.

As part of the numerical modeling effort, a comprehensive cross-platform verification of results was performed using two state-of-the-art CFD platforms: the commercial Ansys Fluent and the open-source framework OpenFOAM. These tools are widely adopted across both industry and academia. Within the Reynolds-averaged Navier-Stokes (RANS) framework, three solver strategies were employed: pressure-based, density-based, and hybrid approaches. For the case of $M=0.3$, all solvers demonstrated satisfactory convergence in both integral and local flow characteristics, reinforcing the robustness of the modeling approach.

The next milestones in the development of this AFC concept include:

• Further validation of the large-eddy simulations for flow regimes at Mach numbers of $0.2$ and $0.3$;
• Preliminary conceptual design of an aircraft featuring aerodynamic boundary layer control, integrating the fuselage with the thick airfoil (bluff-body) and AFC system.

Appendix A.1. Brief Description of the Wind Tunnel and Experimental Conditions

The experimental facility was consist of the the rectangular shaped nozzle (height 120 mm, width 170 mm) designed for a nominal Mach number of $2.3$. A deflector as a flat plane with a sharp leading edge is held by the wind tunnel upper wall to generate an oblique shock wave which reflects on the lower wall. With the deviation angle of $8^0$ the interaction between the shock and the lower wall boundary layer is strong and the boundary layer fully separates. Figure A1 displays a description of the nozzle in two-dimensional space and associated distribution of the Mach number. Figure A2 shows a zoom near deflector and prime features of the flow topology: incident and reflected shock waves, incoming boundary layer, separation and relaxation zones with expansion fan. The boundary layer is fully turbulent and develops on a lower wall with nearly adiabatic constant wall temperature. Just upstream of the interaction ($x=240$ mm), it has thickness of 10 mm. The stagnation pressure and temperature being respectively equal to $50K$ Pa and 300 K, the Reynolds number based on the momentum thickness is about $4,600$.

Figure A1. Description of the wind tunnel with installed deflector and isocontours of the Mach number

Figure A1. Description of the wind tunnel with installed deflector and isocontours of the Mach number

Figure A2. Sketch of the turbulent shock wave/boundary layer interaction with description of main flow features

Figure A2. Sketch of the turbulent shock wave/boundary layer interaction with description of main flow features

Both Particle Image Velocimetry (PIV hereafter) and Laser Doppler Velocimetry (LDV hereafter) were used as independent measurement methods. LDV was used mainly for comparisons with PIV, and by cross- check, to validate PIV measurements. The flow was seeded with particles of incense smoke. Test through a shock wave showed that their size is about 1 micron, and provides sufficient resolution. A PIV/LDA cross- check in the incoming boundary layer showed a good agreement of mean velocity, even at distances as small as $y_{+}=40$. This comparison suggests that the longitudinal velocity is measured with a typical accuracy of $1\%$. The level of background turbulence in the outer flow is essentially due to aerodynamic noise radiated by the boundary layers. Its level is less than $1\%$ for velocity turbulence intensity. For these experiments we can consider two-dimensionality of the flow, since it was found that a the span-wise variations of the velocity are less than $\pm 5$ m/s, i.e. less than $1\%$ of the bulk velocity and therefore at the limit of the accuracy of measurements [48]. Further details and experimental data are available online via the ERCOFTAC knowledge base Wiki [49].

Appendix A.2. Brief Description of the Computational Grids, Boundary Condition and Run Matrix for Numerical Exercises

The following grid strategy was applied. The computational domain was the entire wind tunnel including the nozzle and the deflector. The initial unstructured grid (A0) was constructed from homogeneous, predominantly quad elements with a characteristic linear size of $0.007$ mm and viscous subgrid layers attached to the walls (with the following parameters: the height of the first cell is $1^{-5}$ mm, the number of layers is 20, and the growth factor of $1.15$). The second grid (A1) was constructed in the same way, when the characteristic linear cell size was adapted to $0.0035$ mm, and the parameters of the viscous subgrid were: the height of the first cell is $5^{-6}$ mm, the number of layers is 20, and the growth factor is $1.15$). It should be noted that such parameters provided already for the initial grid the $y^{+}$ distribution within about 1. In addition, for these two grids a fairly good grid dependence of the solution was achieved, when the variance for local velocity distributions was only a few percent.

The initial and boundary conditions were set as follows. At the wind tunnel inlet, the total pressure and temperature were set, and the turbulence intensity was low (about $0.1\%$). The outlet conditions are not critical for this problem, since supersonic flow is observed at the outlet. The initial conditions simulate the situation when the shutter from the chamber with parameters $P_0$ and $T_0$ is suddenly opened. The run matrix is given in Table A1. We used all three numerical approaches to solve the viscous, compressible Navier-Stokes equations given in Sec. Section 3.1 and two turbulence models, SA and RKE, respectively. Thermodynamic properties of the gas, such as viscosity and heat capacity, are set by functional dependencies on temperature, also described in Sec. Section 3.1. Heat transfer with the external environment is not taken into account, assuming adiabatic walls, which approximately corresponds to experimental conditions.

Table A1. Summary of the RANS cases to investigate the turbulent shock wave/boundary layer interaction problem: id – case index, Method – pressure-based, density-based or hybrid, Mesh – A0 or A1, TM – turbulence model, RKE or SA, respectively

id	Method	Mesh	TM
A0-PB-RKE	pressure-based	A0	RKE
A1-PB-RKE	pressure-based	A1	RKE
A1-PB-SA	pressure-based	A0	SA
A1-DB-RKE	density-based	A1	RKE
A1-DB-SA	density-based	A1	SA
A1-H-SA	hybrid	A1	SA

Table A1. Summary of the RANS cases to investigate the turbulent shock wave/boundary layer interaction problem: id – case index, Method – pressure-based, density-based or hybrid, Mesh – A0 or A1, TM – turbulence model, RKE or SA, respectively

id	Method	Mesh	TM
A0-PB-RKE	pressure-based	A0	RKE
A1-PB-RKE	pressure-based	A1	RKE
A1-PB-SA	pressure-based	A0	SA
A1-DB-RKE	density-based	A1	RKE
A1-DB-SA	density-based	A1	SA
A1-H-SA	hybrid	A1	SA

Figure A3. Stream-wise (a) and vertical (b) velocity profiles in the upstream boundary layer ($x=261$ and 291 mm) and along the interaction region ($x=330$, 350 and 370 mm)

Figure A3. Stream-wise (a) and vertical (b) velocity profiles in the upstream boundary layer ($x=261$ and 291 mm) and along the interaction region ($x=330$, 350 and 370 mm)

Figure A4. Visualization of the separated zone: streamlines obtained experimentally by LDV (a) and PIV (b) and numerically (c, case A1-H-SA) colored by the velocity magnitude

Figure A4. Visualization of the separated zone: streamlines obtained experimentally by LDV (a) and PIV (b) and numerically (c, case A1-H-SA) colored by the velocity magnitude

Appendix A.3. Results and Discussion

The main results are summarized in Figure A3, which presents comparison of the both LDV and PIV experiments and numerical results for the stream-vise ($U_x$) and vertical ($U_y$) velocities profiles at five different axial states, $x=261$ mm, 291 mm, 330 mm, 350 mm and 370 mm. One can see that satisfactory agreement both for $U_x$ and $U_y$ was achieved between RANS and measurements at the positions closed to the upstream and downstream location of the separation bubble. At the middle part of the reticulation zone (approx., $x=330$ mm), all simulations are failed to reproduce dimensions of the separated zone. Figure A4 demonstrates streamlines obtained experimentally (both LDV and PIV) and numerically (case A1-H-SA), where it can be observed some differences in the form of separated zones. There are also a few additional observations that can be made. The pressure-based and density-based approached revealed quite similar results, meanwhile the hybrid approach a little different from them. The results computed using the SA and RKE turbulence models are quite consistent as well. Of course, the low frequency unsteadiness like the fluctuations induced by the shock motion (frequency about 500 Hz) or the large broadband eddies in the separated zone induced by Kelvin-Helmholtz instability (about 6 kHz). was not captured by the RANS methods.

In practice, all methods capture the primary flow features, such as incident and reflected shock waves, expansion fans, and separation and relaxation zones, reasonably well within the mathematical frameworks defined by the imposed assumptions. The predicted mean velocity field in the upstream boundary layer, as well as the shock wave/boundary layer interaction downstream, also shows satisfactory agreement, indicating that the current numerical methods are capable of modeling this class of problems.