Numerical Methods Comparison of Shock-Induced Separation on Transonic Axisymmetric Hump

Abstract

A transonic hump with a Mach number (M∞) of 0.875 was utilized to compare the prediction capabilities of RANS, URANS, IDDES, and IDDES-SLA for shock-induced separation. The results matched the experimental data concerning pressure, friction coefficients, velocity distribution, and other parameters. The improved delayed detached eddy simulation with a modified sub-grid scale (IDDES-SLA) method performed the best amongst tested numerical methods, demonstrating differences of 1.3% and 4.1% from the experimental results for separation and reattachment locations, respectively. The primary difference between IDDES and IDDES-SLA is reflected in detailed data, such as velocity and high-frequency characteristics, although both methods exhibit similar simulation capabilities for the main structural features of the flow field. The POD and PSD analyses of the flow field results demonstrated that the IDDES-SLA method was more capable of resolving higher modes of flow field.

1. Introduction

Shock wave/boundary layer interactions (SWBLIs) are prevalent in transonic, supersonic, and hypersonic flows. Since the initial observation of this phenomenon by Ferri [1] during experiments on supersonic airfoil flows, the investigation of SWBLIs has received significant emphasis. Through extensive research, it has been recognized that SWBLIs are a significant factor in the performance of high-speed vehicles [2,3,4]. Typically, as the gas traverses a shock wave, its pressure increases, exacerbating the adverse pressure gradient within the boundary layer, which facilitates boundary layer separation due to the increased thickness of the layer around the shock. Recent work has developed various supersonic models, such as sidewall [5,6], cylinder [7] or elliptical cylinder [8], and wedge [9], to uncover the underlying mechanisms of SWBLIs and to establish more accurate simulation methods for predicting aerodynamic characteristics.

In supersonic or hypersonic flow fields, the region before the shock wave is supersonic, while the region behind it is subsonic. At transonic speeds, different dynamic characteristics determine the shape, strength, and location of the shock wave in the interaction zone [10]. In general, a normal shock wave with “λ-structure” occurring at the foot is often observed in transonic flow, which is due to the thicker boundary layer, which will cause a series of wake compression waves to increase the local boundary-layer-edge pressure and decrease the Mach number before strong normal shock; the border of compression waves and normal shock will combine as “λ-structure”. Another important thing is the enhancement of adverse pressure gradient inside the boundary layer, which will significantly increase the possibility of flow separation and the flow field structure will be more complex, which makes a large challenge for simulation methods. In Davis’s study [11], the SWBLI of NACA 64A010 airfoil was investigated through experiment and CFD simulation at Mach 0.8. The results indicated that CFD methods were inadequate in accurately predicting the flow characteristics due to shock-induced boundary layer separation. In [12,13], DNS and WMLES methods were applied, respectively, to research the separating features about NACA 0012 airfoil. However, there are a few differences in the position of separation and attachment. Current evidence suggests that the accuracy of CFD methods employed in transonic separation conditions still requires further enhancement.

Shock-induced separation is one of the most important phenomena about SWBLI in transonic flow field. The boundary layer separation is generally displayed as a certain size of bubble on the object wall, and the separation bubble will mix with the original separation behind the shock wave in some cases, which could cause larger separation flow and more aerodynamic efficiency loss. Bachalo [14] developed a simple axisymmetric hump model for comparison with experiment in order to investigate the question about shock-induced separation and expected more CFD methods were verified. In 2019, NASA reconducted the experiment in the Sandia Transonic Wind Tunnel so as to supplement the measurement data [15]. The experiment used particle image velocimetry (PIV), pressure-sensitive paint (PSP), oil-flow visualization, and oil-film interferometry to measure properties of the incoming boundary layer, the shock/separation interaction region, the detached shear layer and reverse flow region, and reattachment. A blind CFD challenge was built for obtaining more precise CFD methods about the experiment, which would publish the results of measurements after the end of the challenge [16,17,18].

So many methods were applied in the challenge. Riley [19] used RANS and wall-modeled LES to simulate the flow field of the axisymmetric hump by Bachalo. In the investigation, one-equation Spalart–Allmaras and two-equation cubic $k-\epsilon$ turbulence models are employed and the results of RANS are used as inputs for WMLES. For the results of RANS, different turbulence models cause 16% difference in the separation length. Regarding WMLES, the separation point was advanced and the length of separation was shorter. In Gupta’s [20] work, the explicit filtering LES (EFLES) was applied and stationary ring-like structures were observed in the boundary layer at the apex of the hump, which was considered to associate with turbulence recovery. Meanwhile, the separation point was located at x/c = 0.73, which approximated 0.723 of Riley [19] in the methods of RANS and had 20% of difference by WMLES. And, in the paper, it was observed that the boundary of shock wave near the apex of the hump was away from the tunnel wall, which means there is no influence of shock wave reflection by the tunnel wall in the whole flow field. High-order method LES was applied in a quarter of axisymmetric hump models by Rahmani [21]. As a result, the separation location at x/c = 0.54 did not match the experiment (x/c = 0.68) well, but the reattachment point at x/c = 1.23 matched well with the experiment (x/c = 1.2). And the simulation resolved the main characteristics of flow field such as the shock with “λ -structure”. RANS with SST turbulence model [22] and WMLES [23] methods were used in simulations; however, there still exists a little difference in the results of flow field characteristics between experiment and simulations.

In the present work, the CFL3D solver (Version 6.7) was used for obtaining the aerodynamic characteristics of the transonic hump experimented by NASA. CFL3D is an open-source solver designed by NASA; the codes can be modified for adding new turbulence models. After version 6.3, CFL3D supported computing the full N-S equation, which was necessary to resolve the SWBILs accurately. The codes of IDDES-SLA were added in the present solver. In this paper, the half model was applied and RANS, URANS, IDDES, and IDDES-SLA methods with the SST turbulence model were used to discuss the accuracy of predictive ability. To validate the numerical methods and grid topology, three grids of different densities were tested, and the results were compared with experimental data. Ref. [20] has already shown that the tunnel wall had no reflective shock. Therefore, in this paper, the far field was used instead of the tunnel device, simplifying the simulation process without significantly reducing accuracy. This paper provides figures, such as $C_p{,C}_f,x-t$, to analyze the information of the flow field and compare the prediction capacity of four methods and obtained some significant conclusions written in the last section.

2. Methodology

2.1. Two-Equation Shear-Stress Transport (SST) Turbulence Model

The two-equation k-ω SST model is an eddy viscosity model formulated to close the Reynolds-Averaging Navier–Stokes equation. Essentially, it combines the $k-\epsilon$ model [24] and the Wilcox $k-\omega$ model [25,26]. Near the wall zone, the transferred $k-\epsilon$ model is combined with the Wilcox $k-\omega$ model. At the edge of the boundary layer and free shear layer, the transferred $k-\epsilon$ model is applied, and the standard $k-\epsilon$ model is used in the fully turbulent region far from the wall. The $k-\omega$ model contains $k$ equation and $\omega$ equation as follows:

$${\displaystyle \frac{\partial \left(\rho k\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho u_{j}k\right)}{\partial x_j}}=P-{\beta }^{\mathrm{*}}\rho \omega k+{\displaystyle \frac{\partial }{\partial x_j}}((\mu +{\sigma }_k{\mu }_t){\displaystyle \frac{\partial k}{\partial x_j}})$$

(1)

$$\begin{array}{cc}{\displaystyle \frac{\partial \left(\rho \omega \right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho u_j\omega \right)}{\partial x_j}}& ={\displaystyle \frac{\gamma }{{\nu }_t}}P-\beta \rho {\omega }^2+{\displaystyle \frac{\partial }{\partial x_j}}((\mu +{\sigma }_{\omega }{\mu }_t){\displaystyle \frac{\partial \omega }{\partial x_j}})\\ & +2(1-F_1){\displaystyle \frac{\rho {\sigma }_{\omega 2}}{\omega }}{\displaystyle \frac{\partial k}{\partial x_j}}{\displaystyle \frac{\partial \omega }{\partial x_j}}\end{array}$$

(2)

where $k$ means turbulent kinetic energy and $\omega$ means the rate of dissipation of turbulence per unit energy. $\beta , {\beta }^{*},\gamma ,{\sigma }_k and {\sigma }_{\omega }$ are the closure coefficients equal to 0.075, 0.09, 0.556, 0.85, and 0.5, respectively. The $\mu$ means the molecular viscosity, ${\nu }_t$ and ${\mu }_t$ mean the turbulent viscosity, and ${\nu }_t={\mu }_t/\rho$. $F_1$ relates to the mixing function of $k-\epsilon$ and $k-\omega$ model. The form of the $k$ equation will be changed in the following different methods.

2.2. DES Method and DDES Method

DES is a type of method blending RANS and LES. Based on the DES97 method with one-equation S-A turbulence model proposed by Spalart [27], Strelets [28] and Travin [29] modified the length contained in the two-equation SST turbulence model and developed the DES method with $k-\omega$ SST model. The hybrid length $l_{hyb}$ is defined as follows:

$$l_{hyb}=l_{DES}=min(l_{SST},C_{DES}\Delta )$$

(3)

where $\Delta$ represents the sub-grid scale; it can be determined by:

$$\Delta ={\Delta }_{max}=max(\Delta x,\Delta y,\Delta z,)$$

(4)

The ${\Delta }_{max}$ means the maximum scale of three directions of grid units. $C_{DES}$ is defined as:

$$C_{DES}=\left(1-F_1\right)C_{DES}^{k-\epsilon }+F_1{C}_{DES}^{k-\omega }$$

(5)

$C_{DES}^{k-\epsilon }$ and $C_{DES}^{k-\omega }$ are constants and, respectively, are equal to 0.61 and 0.78. The parameter $l_{SST}$ in Formula (3) represents the length scale of RANS, as follows:

$$l_{SST}={\displaystyle \frac{k^{1/2}}{\left({\beta }^{\mathrm{*}}\omega \right)}}$$

(6)

Changing $l_{SST}$ to $l_{hyb}$ and bringing it to Formula (6) and $k$ equation, Formula (2) can be rewritten as DES form:

$${\displaystyle \frac{\partial \left(\rho k\right)}{\partial t}}+{\displaystyle \frac{\partial (\rho u_{j}k)}{\partial x_j}}=P-{\displaystyle \frac{\rho k^{3/2}}{l_{hyb}}}+{\displaystyle \frac{\partial }{\partial x_j}}((\mu +{\sigma }_k{\mu }_t){\displaystyle \frac{\partial k}{\partial x_j}})$$

(7)

According to the above formulas, the specific form of $k$ equation is determined by the scale of turbulence length.

The hybrid length defined in the original DES method is overly simplistic; as a result, the LES mode can intrude into the boundary layer, reducing vortex viscosity and causing nonphysical phenomena such as grid-induced separation (GIS). To address this issue, Menter and Kuntz [30] developed the Delayed DES (DDES) model based on the two-equation $k-\omega$ SST model. The model optimized the hybrid length $l_{hyb}$ as:

$$l_{hyb}=l_{DDES}=min[l_{SST},{\displaystyle \frac{C_{DES}{\Delta }_{max}}{1-F_{SST}}}]$$

(8)

The primary difference lies in the introduction of the blending formulation $F_{SST}$, which significantly determines the transition between RANS and LES. It is easy to know, when $F_{SST}$ is equal to 0, Formula (8) will change to the standard DES model (Formula (3)). In this paper, the $F_{SST}$ is described as follows:

$$\begin{gathered} F_2=tanh\left(ar{g}_2^2\right) \\ ar{g}_2=max\left(2{\displaystyle \frac{\sqrt{k}}{{\beta }^{\mathrm{*}}\omega d}},{\displaystyle \frac{500\nu }{d^2\omega }}\right) \end{gathered}$$

(9)

2.3. IDDES Method

The Improved Delayed DES (IDDES) method was proposed by Shur [31] for resolving the relatively strong logarithmic layer mismatch (LLM) between the inner RANS and the outer LES regions. The IDDES method combined the features between DDES and WMLES and would accelerate the RANS transfer to LES at the separation region. The hybrid length and blending formulation were obviously changed in the IDDES method. $l_{hyb}$ was determined as:

$$l_{lyb}=l_{DDES}={\tilde{f}}_d\cdot (1+f_e)\cdot l_{SST}+(1-{\tilde{f}}_d)\cdot C_{DES}\Delta$$

(10)

Parameters ${\tilde{f}}_d$ and $f_e$ were functions defining the implementation of DDES and the wall-modelled LES modes. The definition is in ref. [31]. The sub-grid scale is defined as:

$$\Delta =min\left[max(C_{w}d,C_w{\Delta }_{max},\Delta y),{\Delta }_{max}\right]$$

(11)

In Formula (11), $d$ was the distance from the grid to the wall, $\Delta y$ was the wall-normal scale of the grid, and $C_w$ was a constant and equal to 0.15. Compared to ${\Delta }_{max}$, the sub-grid scale used in IDDES could further decrease the viscosity of vortex, which was beneficial for increasing the space resolution ratio of the flow field.

2.4. Improved Sub-Grid Scale (IDDES-SLA) IDDES Method

As for frequent sub-grid scale, ${\Delta }_{max}$ was commonly defined as Formula (4). If the grid length in three directions is significantly different, the local vortex viscosity predicted by ${\Delta }_{max}$ would be larger and obviously decreased the analyzed ability regarding turbulence of DES [32]. Tang [33] applied an improved sub-grid scale model proposed by Shur [32] and used ${\Delta }_{\mathrm{S}\mathrm{L}\mathrm{A}}$ instead of ${\Delta }_{max}$. The specific expression is as follows:

$${\Delta }_{SLA}={\tilde{\Delta }}_{\omega }{F}_{KH}\left(<VTM>\right)$$

(12)

In this paper, an improved sub-grid scale model was used (IDDES-SLA). The details of model construction can be found in [33].

2.5. Proper Orthogonal Decomposition (POD)

Proper Orthogonal Decomposition (POD) is an analysis method that can obtain the primary characteristics of a flow field and identify the structures of different energy levels in the turbulence field. The basic concept is establishing a set of standard orthogonal basis functions from the series spaces of a high-order nonlinear flow system. The linear combination of established standard orthogonal basis functions would express the whole-time flow field approximately.

The information $u_i$ of transient flow could be described as:

$$u_i={\displaystyle \frac{1}{N}}{\sum }_{j=1}^N{u}_j+v_i$$

(13)

in which ${\displaystyle \frac{1}{\mathrm{N}}}\sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{u}}_{\mathrm{j}}$ means the average flow and ${\mathrm{v}}_{\mathrm{i}}$ means the transient value and $N$ is the total number of the field snapshot. For the transient part, the POD method decomposed it into linear combinations of a series of low-order optimal orthogonal bases as:

$$v_i={\sum }_{j=1}^N{a}_j\left(i\right)P_j$$

(14)

$P_j$ is the mode basis function and $a_j\left(i\right)$ is the corresponding mode coefficients at the time $i$. Through the analysis of different modes, the main flow structures would be revealed and the energy accumulation in the field would be elucidated.

2.6. Simulation Setup

The transonic axisymmetric hump model was designed by Bachalo [14], and the specific data are in [16]. In this paper, the model was simplified to half. According to the research of the quarter model [19,21], for this kind of model with circumferential symmetry and a small cross-flow effect, using half mode for calculation was able to ensure the accuracy of the flow field and reduce the difficulty of calculation. In addition, from the simplified model, it could obtain how the model design influences the IDDES simulation results. Due to no impact of reflective shock wave, the far field boundary was used instead of the tunnel wall.

Figure 1 displays the axisymmetric hump. The chord length of the hump was 2.53 inches. The original point was located at the center of the starting section of the hump, and the apex was located on the middle. The freestream was along the x-axis in this simulation. Figure 2 shows the far field of the model, which was a half portion of a cylinder with a radius of 20$c$ and a length of 200$c$. The center lines of the hump and cylinder were at the same position. The back section of the hump was expended to the boundary of the far field along the x-axis and maintained the original radius. The whole hump model was beyond the middle position a bit in the far field for obtaining the clean freestream.

A structured grid method was applied and the fine surface grids of the far field and hump are displayed in Figure 3. “O-Block” topology strategy was adopted for creating a high-quality grid of the boundary layer. Figure 4 shows the node distribution along the x-axis and circumference of boundary layer grids in the far field out of the hump wall side. From the experiment and other references, both shock wave position and separation point were located after the apex of the hump, and the reattachment point was at about $x/c=1.2$. In this paper, the high-density grid operation began slightly before the apex and extended to the wake region at the position of $x/c=1.8$. The node quantity in the wake region was in consideration of the resolution of the wake characteristic. The quantity of total grid and node distribution is in

Table 1. Node information of different grids.

Type	Total Grid	Circumferential Nodes of Hump	Wake Region Nodes Along x-Axis
Fine	46 million	165	905
Medium	16 million	113	637
Coarse	5 million	77	449

. Three types of grids with regards to space resolution were applied in order to verify grid independence in this paper, and the height of the first layer for all three grids is set smaller than 10⁻⁵ times of $c$, which would ensure the value of Y-Plus is approximate to 1 order.

For boundary conditions, the hump was set as a no-slip wall, and the outer surfaces of the far field were set to far-field conditions with specified air parameters. The plane cutting the hump was set as a symmetry plane. The Mach number was 0.875 and the Reynolds number was 1 million, in which the requirements were the same as the experiment [17]. In CFL3D solver, the dimensionless time step “$t$” was applied, which was determined as $t=\tilde{t}\widetilde{a_{\infty }}/\widetilde{L_R}$. $\tilde{t}$ was real dimension time step, $\widetilde{a_{\infty }}$ was the sound speed of the far field, and $\widetilde{L_R}$ was the reference length used by CFL3D code and it was determined as $\widetilde{L_R}=\tilde{L}/L_{ref}$. $\tilde{L}$ was the characteristic length and the chord length was selected as $\tilde{L}$. $L_{ref}$ was the corresponding length in the grid (dimensionless). In the present work, it meant the length value of chord length in the grid. In this paper, the dimensionless time step was set to 0.007, while the real time step was about 5 × 10⁻⁶ s.

3. Results and Discussion

3.1. Grid Resolution and Field Validation

The Q value was utilized to characterize the flow field details in the wake region [34], which was determined as follows:

$$Q_c=-{\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial z}}\right)}^2\right]-\left({\displaystyle \frac{\partial u}{\partial y}}\right)\left({\displaystyle \frac{\partial v}{\partial x}}\right)-\left({\displaystyle \frac{\partial u}{\partial z}}\right)\left({\displaystyle \frac{\partial w}{\partial x}}\right)-\left({\displaystyle \frac{\partial v}{\partial z}}\right)\left({\displaystyle \frac{\partial w}{\partial y}}\right)$$

(15)

where x, y, and z represent the co-ordinates of the flow in the streamwise, lateral, and wall-normal directions, and u, v, and w represent the velocity along the x, y, and z directions, respectively. For visualizing the vortex intensity, the Q value was colored by Mach number for three kinds of grid, as shown in Figure 5.

The results of the IDDES-SLA method were visualized to validate the flow field characteristics. It is evident that the capability for vortex analysis improved with increasing grid resolution. A greater number of small-scale vortices were captured in the medium and fine grids, while the vortex structure was less defined in the coarse grid, potentially resulting in the loss of turbulent details and inaccuracies in aerodynamic characteristics. From the Q isosurface figures, the scales of vortex caught by the cases of medium grid and fine grid were similar, which qualitatively proved that the medium grid is enough to analyze the flow field. As for the Mach number, few differences could be found in the fine grid and medium grid. In the coarse grid, the distribution of the Mach number after shock wave was not symmetric and significantly different to the other two grids.

A monitor point was set on the surface of the hump, as shown in Figure 6a. The point was located directly above the original point in the $z-x$ plane, and the fluctuations of dimensionless density are plotted in Figure 6b. The dimensionless density was determined as $\rho ={\rho }^{*}/{\rho }_{\infty }^{*}$, where the ${\rho }^{\mathrm{*}}$ was local dimensional density and ${\rho }_{\infty }^{\mathrm{*}}$ was dimensional density of inlet flow at the far field. It can be discovered visibly that density vibration at the monitor point has clear periodicity. The amplitude about the coarse grid was a bit lower than the other two grids in the whole period; in opposite, the fine grid and medium grid possessed a similar amplitude. Figure 6c shows the velocity distribution along the z-axis. The start position was across the monitor point and at the location $z/c=0.4$, where $z/c$ is the ratio of the z co-ordinate value and chord length $c$. The experimental data were published by Sandia Wind Tunnel [15]. The results simulated by the three grid types exhibited good agreement with the experimental data, particularly when $z/c$ > 0.5, in which $u/u_{\infty }$ > 0.9 and the flow was beyond the boundary layer. In the boundary layer, there were minor discrepancies compared to the experimental results in all three grids. In particular, the results produced by coarse grids fit the experimental data poorly, which would obtain SWBLIs inaccurately and deeply influence the flow field in the wake region. For both the fine and medium grids, the most noticeable and acceptable difference with the experiment was in the area with the highest rate of velocity changes. In other regions, both grids simulated the flow accurately. Overall, the monitor point was situated in front of the shock wave, and the flow characteristics at this location significantly influenced the simulation results due to the interaction between the shock and the boundary layer. The cases of fine grid and medium grid performed well at this point.

Pressure coefficient and friction coefficient were provided by the experiment and the comparisons are displayed in Figure 7. The error bars of the experiment were plotted in the meanwhile. In Figure 7a, before the peak of negative pressure, all three grid cases fit the experimental data well. At the peak, the negative $C_p$ was much higher in the coarse grid condition, and the results were getting worse after encountering the shock wave; the condition did not become better until $x/c$ greater than 1.6, at which point the position was far away from reattachment point (refer to Figure 7b) and the flow became relatively steady. Looking at the $C_p$ value in the cases of fine grid and medium grid, it performed well whether at the peak or after shock wave. Figure 7b shows the surface friction coefficient. Overall, the simulation results had certain deviations and fluctuations compared to the experiment, a condition observed in many other references [19,21]. These proved that the current CFD methods still have some weakness in the prediction of $C_f$. In contrast, the errors produced by the coarse grid were more obvious than the fine grid and medium grid. Compared to the specific data in

Table 2. The specific location of separation and reattachment for experiment and three types of grids.

Type	Separation $\mathbf{Location} (\mathit{x}/\mathit{c}$)	Reattachment $\mathbf{Location} (\mathit{x}/\mathit{c}$)
Experiment	0.68	1.2
Fine	0.675	1.214
Medium	0.671	1.249
Coarse	0.664	1.324

, the separation points were very close to the experiment in the three cases; however, the significant deviation about reattachment point appeared in the coarse grid.

In this part, different flow filed characteristics were compared in three quantities of grids, and the results approved that the coarse grid could not simulate well the details of the field. The comparison between a fine grid and medium grid indicated that improvement in accuracy was very limited, although the grid quantity increased a lot. In the following analysis, the results of the medium grid are applied to discuss the predictive ability for different methods.

3.2. Discussion of Different Methods and Flow Field

The analysis of grid resolution confirmed the accuracy of the medium grid. In this section, different methods were applied for further investigation based on the medium grid. The RANS, URANS, IDDES, and IDDES-SLA were named from Case1 to Case4, respectively.

Table 3. Separation and reattachment location.

CaseID	Separation Location (x/c)	Difference From Experiment (%)	Reattachment Location (x/c)	Difference From Experiment (%)
Experiment	0.68	/	1.2	/
CASE1: RANS	0.61	10.2	1.137	5.3
CASE2: URANS	0.611	10.1	1.253	4.4
CASE3: IDDES	0.712	4.7	1.275	6.3
CASE4: IDDES-SLA	0.671	1.3	1.249	4.1
HpMusic (p = 2)	0.54	22.9	1.23	2.5
Riley et al., SA	0.72	5.7	/	/
Riley et al., KE	0.7	2.9	/	/
Riley et al., WMLES	0.9	27.8	/	/
Gupta et al., Medium	0.76	11.1	1.14	5.1

was inherited from Rahmani [21] and the results of this paper were added. The table displays the specific separation and reattachment location and the difference from the experiment in different CFD methods or turbulent models. The difference in separation between the experimental and simulation results from Case1 to Case4 progressively decreased. The RANS method applied in Case1 produced steady results, different from this; from Case2 to Case4, the unsteady variation was considered and the averaged results in a whole period were used. Thus, Case1 to Case4 was closer to real flow. Theoretically, the accuracy of IDDES and IDDES-SLA was higher than URANS, and the separation results approved this opinion; even the difference was lower to 1.3% through the IDDES-SLA. Not similar to separation, the difference for reattachment was about 5% in Case1 to Case4; these differences were acceptable and close to other references. There were no significant differences in reattachment locations among the four methods. And Case4 was the closest to the experiment. Compared to other references, even though LES methods like HpMusic and WMLES were applied, the results did not meet expectations, particularly at the separation location. But the reattachment locations displayed in Table 3 were similar. These may approve that, although LES has strong analytical ability for flow fields, it is necessarily suitable for the current transonic hump structure. In contrast, the RANS methods (Case1, KE and SA turbulent model applied by Riley) with weaker parsing ability can accurately obtain the separation position. The entire flow field did not separate before the shock wave and, due to the high viscosity of the RANS method, the simulation of this scenario was more accurate. However, the LES method was easily affected by the grid, so it may not be able to simulate well. In the separated flow field after the shock wave, the LES method becomes more accurate in calculating the fine structures of the field, thus obtaining better simulation results. The results of HpMusic model well illustrate this opinion. As for the IDDES method, especially after improving the sub-grid discriminant scale (IDDES-SLA), it could reduce the influence of the grid near the wall and capture sufficiently accurate turbulent structures in the wake region; the simulation results for the entire flow field are highly accurate.

Pressure and friction coefficients on the hump surface intersected by $z-x$ plane are shown in Figure 8a,b. From the perspective of $C_p$, there were no significant differences in the distribution between Case1 and Case4 before the separation occurred. However, in the separation zone after the shock wave, the four methods began to show some differences. The results of Case1 and Case2 were relatively close, while the results of Case3 and Case4 were relatively close. After $x/c$ exceeded 0.8, the trend of changes in Case3 and Case4 was basically consistent with the experiment, and the $C_p$ at most positions correspond well with the experiment. Although Case1 and Case2 were close to the experimental trend, there was still a significant difference. It was not until x/c exceeded 1.6 that the experimental results gradually became consistent. This may be due to the fact that Case1 uses the RANS method, which cannot simulate the separated flow well. Case2 added transient factors on the basis of Case1, which could improve the accuracy of simulating the flow in the wake region to a certain extent. Case3 and Case4 used RANS method near the wall and LES method far away from the wall in the separation zone, so the simulation of the flow field structure in the wake zone was more accurate and closer to real flow. Therefore, the average results matched the experimental results more closely. For $C_f$, similar to the $C_p$ results, Case1 and Case2 based on RANS and Case3 and Case4 based on IDDES methods remained consistent but, overall, there were significant differences between these two methods. In the range of $x/c$ from 0.6 to 0.8, Case1 and Case2 showed a clear trend of first decreasing and then increasing, while Case3 and Case4 did not exhibit this characteristic, which may be due to the method itself. As $x/c$ increased, although the calculation results of the four methods were consistent with the experimental trend, the errors became progressively larger. Among them, the errors of Case3 and Case4 are smaller than those of Case1 and Case2, indicating that the IDDES method has stronger simulation ability for separated flow than the RANS method. Based on other references, current simulations of $C_f$, regardless of the method used, have significant differences from experiments, and this is not an isolated case. The distribution of $C_f$ on the surface of Case1 was described in Figure 8c, and the separation position and reattachment position were clearly observed. These two positions had certain changes along the circumference, especially the reattachment position. This may be mainly due to the unsteady characteristics of the entire flow along the circumferential direction reflected in Case1, and the division of grid blocks may also have a certain impact on the results.

The velocity distribution along the z-axis on the hump surface from $x/c=0.05$ to 1.5 is shown in Figure 9. The $r_{Surf}$ was the distance from the $z$ co-ordinate to the hump surface in the $z-x$ plane, where the $r_{Surf}$ of the hump surface was 0. From the whole view of the velocity distribution, Case1 to Case4 matched all the experimental values well when the $r_{Surf}/c$ was greater than 0.08, but there was a certain difference in the range of 0 to 0.08. Especially between the separation point and the reattachment point, the difference was very obvious. Within the range of $x/c$ from 0.055 to 0.6, the velocity distribution of Case1 and Case2 remained almost unchanged but, as $x/c$ increased, there was a significant difference between the velocity distribution near the wall and the experimental values. Case3 and Case4 were in good agreement with experimental values within this range. When $x/c$ was greater than 0.65, it was basically near the shock wave and separation region. At this point, all cases from Case1 to Case4 showed significant differences from the experimental values, but the overall trend of velocity distribution remained similar. As $x/c$ increased, in the wake region, there was a significant difference between Case1 and Case2 and the experiment. The prediction of the reverse velocity distribution was significantly greater than the experiment itself, and there was a clear bulge in the figure. However, Case3 and Case4 can well match the experiment. Case4, which specifically used the IDDES-SLA method, had a better fit with the experiment than Case3, which also demonstrated the advantages of improving the sub-grid discriminant scale. At the flow reattachment, the difference between Case1 and Case2 in the experiment gradually decreased. After guiding $x/c$ to greater than 1.4, Case1 to Case4 remained basically consistent with the experimental values. It can be seen that the four methods of RANS, URANS, IDDES, and IDDES-SLA had good agreement with experimental values outside the wake region. However, in the wake region, RANS and URANS methods cannot accurately simulate the surface velocity distribution of the hump. This was mainly due to the high modeling of the RANS method, which cannot accurately analyze the complex flow in the wake region. Although there was a certain gap between the IDDES method and experimental values, the overall resolution of the wake region is significantly better than the RANS method.

Figure 10 displays the Q isosurface at the $z-x$ section. In four cases, the shock wave was obvious and the $\lambda$-structure appeared at the root of the shock wave, which was fit to other research. The main reason for the appearance of structure may be due to the appearance of secondary waves after the shock wave, but fusion occurred above the two shock waves, ultimately forming this structure. After the shock wave, all four cases generated separated vortices. The difference is that Case1 and Case2 had larger vortex structures, but Case2 had a slightly smaller volume than Case1. On the other hand, Case3 and Case4 had similar vortex structures at the same location but with significantly smaller volumes. In the separation zone, there were a large number of finely divided vortex structures, which were not present in Case1 and Case2. Comparing Case3 and Case4, there were no significant differences between the two methods in capturing vortex structures in the entire wake region, and both demonstrated the ability to simulate fine flow field structures. Through spatial slicing, it can be seen that RANS and URANS methods have a high degree of modeling for separated flow fields and cannot accurately capture small vortex structures, while IDDES and IDDES-SLA methods had a higher degree of refinement. According to the vortex details in four methods, the RANS and URANS were worse in predicting the field details and the instability state of the shear layer. One of the most obvious phenomena was that the size of the rigid vortex after shock wave in RANS and URANS was bigger than the other two DES methods. Although there still was a small degree of modeling in predicting the $K-H$ instability phenomena of the shear layer in IDDES and IDDES-SLA, the description of field details was better. Theoretically, the vortex viscosity of the sub-grid model ${\Delta }_{SLA}$ was smaller than ${\Delta }_{max}$; thus, the detailed structures of IDDES-SLA in the field were characterized well.

The results from Case2 to 4 exhibited transient characteristics. Monitoring areas were set up at 5%, 10%, and 20% of chord height positions on the surface of the hump in $z-x$ section to monitor the changes in Mach number over time in the $x/c$ range from −0.8 to 2.0. Figure 11a shows the monitor position and Figure 11b shows the contour of Mach number of Case2 to Case4. It can be seen that, in the wake region at the 5% C monitoring position, the Mach number of Case2 fluctuated less with time, while both Case3 and Case4 exhibited strong unsteady fluctuations. At the 10%C monitoring location, Case2 showed no significant fluctuations over time, while the fluctuations in the wake regions of Case3 and Case4 decreased. Near the 20%C, the Mach number fluctuations over time for all three cases were not significant. It can be seen that, at least from the 20%C height outward, the entire flow field mainly exhibited steady characteristics. At the monitoring positions of 10%C and 20%C, there was a high Mach number region near $x/c=$ 0.8 to 1.0 in Case2, but the region did not appear at the same position in Case3 and Case4. This region was compressed within a very small range near $x/c$ = 0.8, which may be mainly due to the fact that, although the URANS method has transient characteristics, it cannot accurately simulate the fine flow field structure in the wake region. This also indicated that the IDDES method performs well in simulating the temporal and spatial characteristics of the flow field.

For further comparison, the contour of time-averaged streamwise velocity (u/u_ref), radial velocity(v/u_ref), and Reynolds shear stress (u′v′) for Case3 and Case4 are displayed in Figure 12. The difference between Case3 and Case4 in terms of streamwise velocity was mainly reflected in details, and there was not a significant difference in macro aspects. The shock structure shown in Case4 is more pronounced and enlarged, and the overall streamwise velocity in the wake region is slightly lower than that in Case3. In terms of radial velocity, there were significant differences between Case3 and Case4, but the range of the wake region remained basically the same. In terms of Reynolds stress, the macroscopic distribution of Case3 and Case4 was basically the same, but Case4 had a clearer boundary at the shock wave position. From the above comparison, it can be seen that there was no significant difference in the macroscopic distribution simulation results of flow field characteristics between IDDES and IDDES-SLA methods. The difference mainly lay in some details of the flow field, which was mainly due to the improvement of the original IDDES sub-grid discrimination method.

POD method was applied to further compare the flow field resolution capability of IDDES and IDDES-SLA. The dimensionless density was used to perform the POD analysis of 500 flow field snapshots when the flow variation periods were stable. The gap of snapshots was 0.0001 s. Figure 13 showed the POD energy accumulation curve from mode 0 to mode 250 of two methods. It was obvious that the accumulation rate was faster for the IDDES method, where even the energy percent was close to 50% in the first 50 modes, while the IDDES-SLA was about 30%. It could be seen that the curves grow faster in the low modes of IDDES; however, the growth was moderate in the whole 250 modes of IDDES-SLA. This suggested that the change in the sub-grid determination method further improved the ability of the IDDES method to portray the details of the flow field and attenuated the degree of method modeling. The characteristics of low modes would be displayed and recognized more accurately in the IDDES-SLA. Figure 14 shows the spatial distribution of POD mode1, mode3, and mode5 of the IDDES and IDDES-SLA flow field. It should to be noted that the mode coefficients of all POD modes were 1, so the specific value of contours could not explain the meaning of field snapshots. The contours ensured the consistent upper and lower limits. In general, the primary characteristics of three low-level POD modes were similar, which improved the basic features of the flow field consistently. The shape of separation bubbles caused by SWBLIs were marked by a red dotted line in mode 1; the two bubbles was basically the same. This also demonstrated that the two IDDES methods are able to accurately catch the obvious SWBLI field features. The shock wave was marked in mode1 of the IDDES method. Compared with IDDES-SLA, the area was bigger and the period of vibrations was more apparent. This may be due to the deep modeling of the IDDES method, to the extent that the accuracy was reduced in the shock region. The separation line was drawn in the contours. From all contours, the separation position was in a good correlation with the location of the shock wave root in the boundary layer. The separation location of IDDES ($x/c=0.712$) was slightly backward compared to IDDES-SLA ($x/c=0.671$); it may be caused by the lager area of shock wave. The reattachment lines were also painted. In mode 5, the vortex structures after the reattachment line of IDDES results were disordered; there was a high probability that the features in this region were deeply modeled. Again, in IDDES-SLA, the regular vibration was shown in the whole weak region. The different results in two methods intuitively demonstrated that the IDDE-SLA method was more accurate in predicting the current transonic SWBLIs flow field. The power spectral density (PSD) of dimensionless density of Strouhal number ($S_t$) from 0.01 to 2 at the monitor point ($x/c=1.0$, $r/c=0.44$) is displayed in Figure 15. The $S_t$ was determined as $S_t=fL/V_{\infty }$, in which $f$ means the frequency of dimensionless density after Fast Fourier Transform (FFT) of 500 field snapshots. The length $L$ was the chord length of the hump and the $V_{\infty }$ was the velocity of the far field. It could be seen that two methods have similar characteristics at $S_t=0.1$, from 0.1 to 1.0; the differences gradually appeared but the slope of all PSD vibration was close to −5/3, which suggested that the resolution of turbulence structures of two methods was credible. After $S_t=1.0$, the PSD of IDDES diverged from −5/3 obviously and the IDDES-SLA still kept the original slope. It was supposed that IDDES-SLA has more credible high-frequency turbulent features, and the method had better resolution capacity for small turbulence structures of the flow field.

4. Conclusions

In the present study, the transonic hump model designed by Bachalo [14] was simulated and compared using RANS, URANS, IDDES, and IDDES-SLA methods, all within the $k-\omega$ SST turbulent model, and the results were compared with the experimental data from Lynch [15]. The half model was used in the present work. And the tunnel wall was replaced by a far field boundary, which, due to the boundary of the shock wave, was far from the tunnel wall and the interaction was so weak that it can be ignored. A 16 million structured grid was applied in CFL3D solver at the condition that the Mach number was 0.875 and Reynolds number was 1 million. The grid independence was proved by comparing the field details and some specific data among three resolutions of grids.

For method comparison, the separation location simulated by RANS and URANS had 10% difference from the experiment, while the IDDES and IDDES-SLA was 4.7% and 1.3%, respectively. For the reattachment location, there was approximately a 5% difference among the four methods. The most accurate results predicted by the IDDES-SLA method indicated that the flow separated at $x/c$ = 0.61 and reattached at $x/c$ = 1.249.

Pressure and friction coefficients were analyzed, and the IDDES methods showed better accuracy in the wake region, possibly due to the high modeling of RANS methods. For velocity distribution, the results simulated by RANS and URANS matched the experiment well, except the separated region ($x/c$ from 0.65 to 1.4). IDDES and IDDES-SLA were basically consistent with the experiment in the whole zone, especially the IDDES-SLA. It can be observed from the Q isosurface at the $z-x$ section that IDDES methods had better resolution capacity and can capture the fine vortex structures, and the RANS methods only can obtain the macroscopic aerodynamic characteristics accurately.

The temporal characteristics of URANS, IDDES, and IDDES-SLA were analyzed. It was shown that the flow field displayed time-independence beyond the height of 20% chord length. In the transient region, the variation was not obvious for URANS and the IDDES methods were opposite. The difference between IDDES and IDDES-SLA was basically in the detail data of the flow field. As for macroscopic characteristics or the resolution of vortex, the two methods performed with great prediction capacity.

The deeper comparison work was conducted for IDDES and IDDES-SLA applying the POD analysis method. The results supposed that the original IDDES method was not sufficient in the resolution of small turbulence structures and would lose some period characteristics in the weak region. After changing the sub-grid method in IDDES-SLA, the high-frequency features would be more credible.

This study primarily focused on comparing four CFD methods but conducted limited analyses of phenomena such as shock wave vibrations and the influence of periodic characteristics. These aspects will be explored further in future work.