Unsteady Simulation of Transonic Buffet of a Supercritical Airfoil with Shock Control Bump

: The unsteady ﬂow characteristics of a supercritical OAT15A airfoil with a shock control bump were numerically studied by a wall-modeled large eddy simulation. The numerical method was ﬁrst validated by the buffet and nonbuffet cases of the baseline OAT15A airfoil. Both the pressure coefﬁcient and velocity ﬂuctuation coincided well with the experimental data. Then, four different shock control bumps were numerically tested. A bump of height h/c = 0.008 and location x B /c = 0.55 demonstrated a good buffet control effect. The lift-to-drag ratio of the buffet case was increased by 5.9%, and the root mean square of the lift coefﬁcient ﬂuctuation was decreased by 67.6%. Detailed time-averaged ﬂow quantities and instantaneous ﬂow ﬁelds were analyzed to demonstrate the ﬂow phenomenon of the shock control bumps. The results demonstrate that an appropriate “ λ ” shockwave pattern caused by the bump is important for the ﬂow control effect.


Introduction
Modern civil aircraft usually adopt a supercritical wing and operate at transonic speed [1]. Shockwaves appear at both the cruise condition and the buffet condition for a supercritical wing or airfoil [2]. The buffet, which is a self-sustained shock oscillation induced by periodic interactions between a shockwave and a separated boundary layer [3], is a key constraint of the flight envelope of a civil aircraft [4].
The transonic buffet phenomenon has been studied for decades. It was first observed by Hilton and Fowler [5] in a wind tunnel investigation. When a buffet occurs, the shockwave location and the lift coefficient, as well as the turbulent boundary layer thickness after the shockwave, change periodically [5]. Iovnovich and Raveh [6] classified transonic buffets into two types: one is where a shock oscillation appears on both the upper side and lower side of the wing, while the other is where an oscillation only exists on the upper side. The latter is more common for a supercritical wing of a civil aircraft [7]. Pearcey et al. [8] studied the relation between the flow separation near the trailing edge and the onset of a buffet, which is widely used in the determination of the buffet onset boundary [3]. When the shockwave strength is weak, which is the case for a supercritical airfoil, the Kutta wave is believed to play an important role in the shockwave oscillation [9,10]. A selfsustained feedback model of a shockwave buffet was proposed by Lee et al. [10]. The model illustrated that when the pressure fluctuation of the shockwave propagates downstream of the trailing edge, the sudden change in the disturbance forms a Kutta wave propagating upstream, and the shockwave is made to oscillate by the energy contained in the Kutta wave [10]. The Kutta wave, which propagates by the speed of sound [10], is also considered a sound wave generated at the trailing edge [11]. Chen et al. [12] further extended the feedback model to more types of shockwave oscillations. Crouch et al. [13] found that both the stationary mode and oscillatory mode of instability occur on an unswept buffet wing, and that the intermediate-wavelength mode has a flow feature of the airfoil buffeting mode.
Experimental [14,15] and numerical investigations [16][17][18] have been widely used in buffet studies in recent years. The transonic buffet characteristics of an OAT15A supercritical airfoil were experimentally measured [14] and widely adopted as a validation case for numerical study. The unsteady Reynolds-averaged Navier-Stokes (RANS) method [16] and detached eddy simulation [17] can resolve the large-scale pattern of a transonic buffet and determine the main frequency of the shockwave oscillation, while a large eddy simulation (LES) [18] can provide abundant details of the interaction between the shockwave and turbulence structures inside the boundary layer.
The buffet onset boundary determines the flight envelope of a civil aircraft [19]. Suppressing the shockwave oscillation and improving the buffet onset lift coefficient are vital for modern civil aircraft design. Vortex generators [20] and shock control bumps (SCBs) [21] are two common passive flow control devices used for transonic buffets. The vortex generator has been intensively studied and widely adopted on Boeing's civil aircraft to control the transonic buffet [22]. However, SCBs are still receiving considerable research attention and are not yet practical on civil aircraft. SCBs can be classified into two-dimensional (2D) bumps and three-dimensional (3D) bumps [23]. The 2D bump is identical along the spanwise direction, while the 3D bump can have different shapes in the spanwise direction.
The basic control mechanism of a 2D SCB is a combination of the compression wave and expansion wave near the shockwave [24]. A supercritical airfoil usually has a weak normal shock on the upper surface. If a bump is placed at the shockwave location, a compression wave can be formed when the flow moves upward from the bump and weakens the normal shock. Expansion waves downside of the bump are also beneficial for pressure recovery after the shockwave. With the compression and expansion waves, a "λ" shockwave structure is formed instead of a normal shock on the airfoil [25], meaning the wave drag can be reduced. As the benefit comes from the interaction between the shockwave and the bump, the drag reduction effect is relevant to the bump location and height. Sabater et al. [26] applied a robust optimization on the location and shape of a 2D SCB and achieved a drag reduction of more than 20% on an RAE2822 airfoil. Jia et al. [27] found that a 2D SCB can reduce the drag coefficient at a large lift coefficient; however, it increases the drag coefficient at a small lift coefficient. Adaptive SCBs [28,29], which change shapes according to the shockwave location, provide drag reduction effectiveness under a wide range of flow conditions. Lutz et al. [30] optimized adaptive SCBs and found the drag curve shows significant improvement compared to the initial airfoil.
The control mechanism of a 3D SCB includes not only the "λ" shockwave structure formed by the bump but also the vortical flow induced by the bump [31], which is similar to a small vortex generator [32]. The streamwise vortex induced by the edge of the 3D SCB is beneficial for the mixing of high-energy outer flow and low-energy boundary layer flow, consequently suppressing the flow separation. A 3D bump with a wedge shape [33,34] or hill shape [31,35,36] can effectively induce a vortical flow after the bump and change the shockwave behavior.
SCBs have shown great potential for buffet control. However, there are also some issues mentioned in the literature. One of the concerns is that the control efficiency is sensitive to the shock location [37]. This might decrease the performance at the off-design condition. An adaptive SCB can extend the effective range of shock control; however, an adaptive SCB might increase the structure complexity and maintenance cost [38].
From the numerical method perspective, the steady [39,40] or unsteady [32,41] RANS method has been widely used in the study of buffets or in the shape design of SCBs for supercritical airfoils. However, the control effect of SCBs is rarely studied by high-fidelity turbulence simulation methods such as the RANS-LES hybrid method or the LES method. As stated before, a buffet is a natural unsteady turbulent flow that includes periodic shockwaves and a thickened turbulent boundary layer under a high Reynolds number. A large eddy simulation can resolve the small-scale turbulent structures of the boundary and the generated acoustic waves of the tailing edge [18], which is a superior method to study the unsteady interaction mechanism between the shockwave and the SCB.
In this paper, the transonic buffet of an OAT15A supercritical airfoil with an SCB was numerically studied by a wall-modeled large eddy simulation (WMLES). The numerical method was validated by the experimental data of the OAT15A airfoil. Several 2D SCBs were adopted to test the control effect of the unsteady flow. The flow fields of both the nonbuffet case and buffet case were analyzed. To the authors' knowledge, this is the first study to use a large eddy simulation on the SCB of a supercritical airfoil. The time-averaged flow quantities, periodic flow patterns and pressure fluctuations were compared for the airfoil without or with an SCB.

Numerical Method
In this paper, a wall-modeled large eddy simulation was adopted to compute the flow field of an OAT15A supercritical airfoil. The solver was developed based on a finite-volume CFD code [42,43]. The inviscid flux was computed based on a blending of an upwind scheme and a central difference scheme [44,45]. The time advancing method is a threestep Runge-Kutta method. Vreman's eddy viscosity model [46] was used to model the subgrid-scale turbulence motion. Based on the estimation of Choi and Moin [47], the grid requirement of a wall-resolving large eddy simulation is proportional to Re 13/7 L , and L is the character length of the flow field; in contrast, the grid requirement of a wall-modeled simulation is proportional to Re L . In this study, the computational cost of a wall-resolving LES was too high because the Reynolds number is 3 × 10 6 . Consequently, an equilibrium stress-balanced equation was solved to model the energy-containing eddies in the inner turbulent boundary layer [44,45]. The wall-modeled equation was solved on an additional one-dimensional grid with 50 grid points to obtain the velocity profile of the inner layer of the boundary layer. The shear stress provided by the wall-modeled equation was applied as the boundary condition of the Navier-Stokes equation. The velocity of the upper boundary of the wall-modeled equation was interpolated from the third layer of the LES grid. The numerical method is not the focus of this paper. Consequently, detailed numerical schemes and formulas are not presented here. Readers can refer to our previous work on high-Reynolds separated flows [44,45].

Computational Grid
The computational grid used in this paper was a C-type grid, which was generated by an in-house grid generation code that solves Poisson's equation. The computational domain was [−18c, 20c] in the x direction and [−20c to 20c] in the y direction (c is the chord length). Although the bump shape was two-dimensional, we carried out a threedimensional simulation because of the three-dimensional nature of the turbulence. The spanwise length was 0.075c, which is larger than the WMLES in the study by Fukushima et al. [18]. The grid had 1945 points in the circumferential direction and 261 points in the normal direction. It had 51 points evenly located in the spanwise direction. The trailing edge thickness was 0.005c, and a grid block with 21 × 51 × 249 points was located on the trailing edge region. The total grid number was 26.2 million. Figure 1a shows the grid of the x-y plane of the baseline OAT15A airfoil. The figure is plotted by every two points of the grid. Figure 1b shows the wall grid near the leading edge of the airfoil. The first layer height of the grid was 2.0 × 10 −4 c. The increasing ratio in the normal direction was less than 1.1. The first grid layer can be located in the logarithmic layer of the turbulent boundary layer because the wall-modeled equation can provide the shear stress boundary condition of the wall. When the SCB was installed on the airfoil surface, the wall surface was deformed to simulate the effect of the bump while keeping the grid topology and grid number the same.
When the SCB was installed on the airfoil surface, the wall surface was deformed to simulate the effect of the bump while keeping the grid topology and grid number the same.
(a) Surface grid in the x-y plane (b) Wall grid near the leading edge

Validation of the Baseline OAT15A Airfoil
The baseline OAT15A airfoil was adopted as the validation case. The flow condition was Ma = 0.73 and Re = 3.0 × 10 6 . [14] The nondimensional time step was ΔtU∞/c = 1.37 × 10 −4 . Approximately 180 time units (tU∞/c) were computed for each case, and the last 100 time units were applied for flow field statistics. It cost approximately 2.0 × 10 5 core hours for each flow case on an Intel cluster with a 2.0 GHz CPU.
The Δx + and Δy + of the first grid layer were collected to demonstrate the grid resolution of the present computation. Figure 2a shows the Δx + and Δy + of the LES grid collected by an angle of attack α = 2.5°. Δy + was approximately 50 on the upper surface and 80 on the lower surface, while Δx + was approximately twice the Δy + . The first grid point in the wall normal direction was located at the logarithmic layer of a turbulent boundary layer, which satisfies the requirement for a WMLES. The Δy + of the first grid layer of the wallmodeled equation (Δy + wm) was always less than 1.0, as shown in Figure 2b. This is sufficient to describe the viscous sublayer of the turbulent boundary layer. The spanwise correlation [43] of the pressure fluctuation was computed to test the spanwise domain. As shown in Figure 3, the correlations of both x/c = 0.4 and 0.8 damp fast in the spanwise direction, which demonstrates that the present spanwise domain size (0.075 c) is adequate.
(a) Δx + and Δy + of the LES grid (b) Δy + of the wall-modeled grid

Validation of the Baseline OAT15A Airfoil
The baseline OAT15A airfoil was adopted as the validation case. The flow condition was Ma = 0.73 and Re = 3.0 × 10 6 . [14] The nondimensional time step was ∆tU ∞ /c = 1.37 × 10 −4 . Approximately 180 time units (tU ∞ /c) were computed for each case, and the last 100 time units were applied for flow field statistics. It cost approximately 2.0 × 10 5 core hours for each flow case on an Intel cluster with a 2.0 GHz CPU.
The ∆x + and ∆y + of the first grid layer were collected to demonstrate the grid resolution of the present computation. Figure 2a shows the ∆x + and ∆y + of the LES grid collected by an angle of attack α = 2.5 • . ∆y + was approximately 50 on the upper surface and 80 on the lower surface, while ∆x + was approximately twice the ∆y + . The first grid point in the wall normal direction was located at the logarithmic layer of a turbulent boundary layer, which satisfies the requirement for a WMLES. The ∆y + of the first grid layer of the wall-modeled equation (∆y + wm ) was always less than 1.0, as shown in Figure 2b. This is sufficient to describe the viscous sublayer of the turbulent boundary layer. The spanwise correlation [43] of the pressure fluctuation was computed to test the spanwise domain. As shown in Figure 3, the correlations of both x/c = 0.4 and 0.8 damp fast in the spanwise direction, which demonstrates that the present spanwise domain size (0.075 c) is adequate. When the SCB was installed on the airfoil surface, the wall surface was deformed to simulate the effect of the bump while keeping the grid topology and grid number the same.
(a) Surface grid in the x-y plane (b) Wall grid near the leading edge

Validation of the Baseline OAT15A Airfoil
The baseline OAT15A airfoil was adopted as the validation case. The flow condition was Ma = 0.73 and Re = 3.0 × 10 6 . [14] The nondimensional time step was ΔtU∞/c = 1.37 × 10 −4 . Approximately 180 time units (tU∞/c) were computed for each case, and the last 100 time units were applied for flow field statistics. It cost approximately 2.0 × 10 5 core hours for each flow case on an Intel cluster with a 2.0 GHz CPU.
The Δx + and Δy + of the first grid layer were collected to demonstrate the grid resolution of the present computation. Figure 2a shows the Δx + and Δy + of the LES grid collected by an angle of attack α = 2.5°. Δy + was approximately 50 on the upper surface and 80 on the lower surface, while Δx + was approximately twice the Δy + . The first grid point in the wall normal direction was located at the logarithmic layer of a turbulent boundary layer, which satisfies the requirement for a WMLES. The Δy + of the first grid layer of the wallmodeled equation (Δy + wm) was always less than 1.0, as shown in Figure 2b. This is sufficient to describe the viscous sublayer of the turbulent boundary layer. The spanwise correlation [43] of the pressure fluctuation was computed to test the spanwise domain. As shown in Figure 3, the correlations of both x/c = 0.4 and 0.8 damp fast in the spanwise direction, which demonstrates that the present spanwise domain size (0.075 c) is adequate.
(a) Δx + and Δy + of the LES grid (b) Δy + of the wall-modeled grid     Figure 4 shows the lift and drag coefficient histories at three different angles of attack. When a buffet does not occur (α = 2.5°), the lift and drag coefficients have oscillations with small amplitudes. Then, the amplitude of oscillation is increased at α = 3.0°. When α = 3.5°, a clear periodic oscillation pattern appears on both the lift and drag coefficients, which demonstrates that a buffet occurs. The lift coefficient was time-averaged for comparison with the experimental data, as shown in Figure 5. The experimental data were extracted from [48]. There were two series of experimental data in the reference. One was from the S3MA wind tunnel, and the other was from the T2 wind tunnel. The Reynolds number in that experiment was 6.0 × 10 6 , which is higher than the present computation. Although the flow condition was slightly different from the present computation, the time-averaged lift coefficients of the present computation were located between the two series of experimental data. The mean pressure coefficient was compared with the experimental data of [14]. As shown in Figure 6, the pressure coefficients at three different angles of attack coincide well with the experimental data, except that the shockwave location is slightly downstream of the experimental data.   Figure 4 shows the lift and drag coefficient histories at three different angles of attack. When a buffet does not occur (α = 2.5°), the lift and drag coefficients have oscillations with small amplitudes. Then, the amplitude of oscillation is increased at α = 3.0°. When α = 3.5°, a clear periodic oscillation pattern appears on both the lift and drag coefficients, which demonstrates that a buffet occurs. The lift coefficient was time-averaged for comparison with the experimental data, as shown in Figure 5. The experimental data were extracted from [48]. There were two series of experimental data in the reference. One was from the S3MA wind tunnel, and the other was from the T2 wind tunnel. The Reynolds number in that experiment was 6.0 × 10 6 , which is higher than the present computation. Although the flow condition was slightly different from the present computation, the time-averaged lift coefficients of the present computation were located between the two series of experimental data. The mean pressure coefficient was compared with the experimental data of [14]. As shown in Figure 6, the pressure coefficients at three different angles of attack coincide well with the experimental data, except that the shockwave location is slightly downstream of the experimental data. The lift coefficient was time-averaged for comparison with the experimental data, as shown in Figure 5. The experimental data were extracted from [48]. There were two series of experimental data in the reference. One was from the S3MA wind tunnel, and the other was from the T2 wind tunnel. The Reynolds number in that experiment was 6.0 × 10 6 , which is higher than the present computation. Although the flow condition was slightly different from the present computation, the time-averaged lift coefficients of the present computation were located between the two series of experimental data.   Figure 4 shows the lift and drag coefficient histories at three different angles of attack. When a buffet does not occur (α = 2.5°), the lift and drag coefficients have oscillations with small amplitudes. Then, the amplitude of oscillation is increased at α = 3.0°. When α = 3.5°, a clear periodic oscillation pattern appears on both the lift and drag coefficients, which demonstrates that a buffet occurs. The lift coefficient was time-averaged for comparison with the experimental data, as shown in Figure 5. The experimental data were extracted from [48]. There were two series of experimental data in the reference. One was from the S3MA wind tunnel, and the other was from the T2 wind tunnel. The Reynolds number in that experiment was 6.0 × 10 6 , which is higher than the present computation. Although the flow condition was slightly different from the present computation, the time-averaged lift coefficients of the present computation were located between the two series of experimental data. The mean pressure coefficient was compared with the experimental data of [14]. As shown in Figure 6, the pressure coefficients at three different angles of attack coincide well with the experimental data, except that the shockwave location is slightly downstream of the experimental data. The mean pressure coefficient was compared with the experimental data of [14]. As shown in Figure 6, the pressure coefficients at three different angles of attack coincide well with the experimental data, except that the shockwave location is slightly downstream of the experimental data.   Figure 7a shows the mean Mach contour at α = 3.5°. It is clear that a low-speed region appears near the trailing edge of the airfoil. The mean shockwave is also smeared because of the periodic oscillation of the shock location. Figure 7b demonstrates the root mean square (RMS) of the pressure fluctuations. It is normalized by the dynamic pressure of the freestream condition. The maximum pressure fluctuation appears near the shockwave location. Figure 8 compares the RMS of the streamwise velocity fluctuation with the experimental data of the S3Ch wind tunnel [49]. The values in Figure 8 are dimensional, and the contours have the same legend. The contours of the computation and experiment are quite consistent in both shape and value, which validates that the present computational method is reliable.    Figure 7a shows the mean Mach contour at α = 3.5 • . It is clear that a low-speed region appears near the trailing edge of the airfoil. The mean shockwave is also smeared because of the periodic oscillation of the shock location. Figure 7b demonstrates the root mean square (RMS) of the pressure fluctuations. It is normalized by the dynamic pressure of the freestream condition. The maximum pressure fluctuation appears near the shockwave location. Figure 8 compares the RMS of the streamwise velocity fluctuation with the experimental data of the S3Ch wind tunnel [49]. The values in Figure 8 are dimensional, and the contours have the same legend. The contours of the computation and experiment are quite consistent in both shape and value, which validates that the present computational method is reliable.  Figure 7a shows the mean Mach contour at α = 3.5°. It is clear that a low-speed region appears near the trailing edge of the airfoil. The mean shockwave is also smeared because of the periodic oscillation of the shock location. Figure 7b demonstrates the root mean square (RMS) of the pressure fluctuations. It is normalized by the dynamic pressure of the freestream condition. The maximum pressure fluctuation appears near the shockwave location. Figure 8 compares the RMS of the streamwise velocity fluctuation with the experimental data of the S3Ch wind tunnel [49]. The values in Figure 8 are dimensional, and the contours have the same legend. The contours of the computation and experiment are quite consistent in both shape and value, which validates that the present computational method is reliable.    Figure 9 further compares the wall pressure fluctuation at α = 3.5°. The experimental data were extracted from [14]. The computed location of the peak fluctuation is slightly downstream, which is consistent with the pressure coefficient distribution (Figure 6c). However, the peak value of the present computation is close to the experiment, which means that the primary shock buffet phenomenon is well captured by the computation. The motivation of the present paper was to study the effect of the shock control bump. The variation caused by the bump, rather than the accurate shockwave location, is our focus in the following sections. Consequently, the numerical method is validated to be practicable.

Shape of the Shock Control Bump
In this section, a two-dimensional SCB is numerically tested and compared with the baseline configuration. Several bumps with different heights and locations are analyzed. The geometries of the bumps are controlled by the Hicks-Henne function [27], as shown in Equation (1). The bump function f(xB) is superposed on the baseline airfoil. h/c is the relative height of the bump, xB/c is the central location and lB is the bump length. Figure  10 shows the geometries of bumps with different parameters. Four bumps are computed in this paper, as shown in Figure 10a. The lengths of the bumps are 0.2c. The first three bumps are all located at xB/c = 0.55, which is the shockwave location of the baseline OAT15A airfoil. The relative heights h/c of the first three bumps are 0.004, 0.008 and 0.012. The fourth bump is located slightly upstream xB/c = 0.50 and has the same height as Bump 2. Figure 10b shows the geometries of the airfoil superposed with the bumps.   [14]. The computed location of the peak fluctuation is slightly downstream, which is consistent with the pressure coefficient distribution (Figure 6c). However, the peak value of the present computation is close to the experiment, which means that the primary shock buffet phenomenon is well captured by the computation. The motivation of the present paper was to study the effect of the shock control bump. The variation caused by the bump, rather than the accurate shockwave location, is our focus in the following sections. Consequently, the numerical method is validated to be practicable.  Figure 9 further compares the wall pressure fluctuation at α = 3.5°. The experimental data were extracted from [14]. The computed location of the peak fluctuation is slightly downstream, which is consistent with the pressure coefficient distribution (Figure 6c). However, the peak value of the present computation is close to the experiment, which means that the primary shock buffet phenomenon is well captured by the computation. The motivation of the present paper was to study the effect of the shock control bump. The variation caused by the bump, rather than the accurate shockwave location, is our focus in the following sections. Consequently, the numerical method is validated to be practicable.

Shape of the Shock Control Bump
In this section, a two-dimensional SCB is numerically tested and compared with the baseline configuration. Several bumps with different heights and locations are analyzed. The geometries of the bumps are controlled by the Hicks-Henne function [27], as shown in Equation (1). The bump function f(xB) is superposed on the baseline airfoil. h/c is the relative height of the bump, xB/c is the central location and lB is the bump length. Figure  10 shows the geometries of bumps with different parameters. Four bumps are computed in this paper, as shown in Figure 10a. The lengths of the bumps are 0.2c. The first three bumps are all located at xB/c = 0.55, which is the shockwave location of the baseline OAT15A airfoil. The relative heights h/c of the first three bumps are 0.004, 0.008 and 0.012. The fourth bump is located slightly upstream xB/c = 0.50 and has the same height as Bump 2. Figure 10b shows the geometries of the airfoil superposed with the bumps.

Shape of the Shock Control Bump
In this section, a two-dimensional SCB is numerically tested and compared with the baseline configuration. Several bumps with different heights and locations are analyzed. The geometries of the bumps are controlled by the Hicks-Henne function [27], as shown in Equation (1). The bump function f (x B ) is superposed on the baseline airfoil. h/c is the relative height of the bump, x B /c is the central location and l B is the bump length. Figure 10 shows the geometries of bumps with different parameters. Four bumps are computed in this paper, as shown in Figure 10a. The lengths of the bumps are 0.2c. The first three bumps are all located at x B /c = 0.55, which is the shockwave location of the baseline OAT15A airfoil. The relative heights h/c of the first three bumps are 0.004, 0.008 and 0.012. The fourth bump is located slightly upstream x B /c = 0.50 and has the same height as Bump 2. Figure 10b shows the geometries of the airfoil superposed with the bumps.
The computational settings are the same as the baseline OAT15A airfoil in Section 2. The computational grids of the bump configurations are deformed from the baseline grid, which keeps the computation appropriate to compare with the baseline OAT15A airfoil.

Time-Averaged Characteristics of the Airfoil with Bumps
The aerodynamic performance of the airfoil with bumps was computed and compared with the baseline airfoil. The nondimensional time step is ΔtU∞/c = 1.37 × 10 −4 . The flow field of the baseline airfoil is used as the initial condition of the bump cases. More than 130 time units (tU∞/c) are computed for each case. The last 80 time units are used for time averaging to obtain the aerodynamic coefficients. Six cases were computed, which are α = 3.5° for all four bumps and α = 2.5° for Bump 2 and Bump 4. Figure 11 shows the histories of the lift and drag coefficients in the computation. It is obvious that the periodic oscillation of the lift and drag coefficients at α = 3.5° is suppressed by Bumps 2, 3 and 4, while Bump 1 has very little effect on the lift oscillation. It can be seen that the bumps may have a negative influence on the lift coefficient. The lift coefficient is more or less decreased at α = 3.5°, while it is significantly decreased at the nonbuffet condition α = 2.5°. Moreover, the drag coefficient is also influenced by the bumps. The drag coefficient at α = 2.5° is notably increased, as shown in Figure 11d.
The time-averaged values of the aerodynamic coefficients are listed in Table 1. The relative variations (in percentage) in the coefficients compared with the baseline configuration are also listed in the parentheses of the table. A tendency of the bumps is that the lift coefficient is consistently decreased by the four bumps at both α = 2.5° and α = 3.5°. The lift coefficient at α = 2.5° is decreased by approximately 5%, and it is decreased by approximately 1% to 5% at α = 3.5°. With increasing bump height (Bumps 1, 2 and 3), the lift drop becomes increasingly significant. If the location of the bump is moved upstream (Bumps 2 and 4), the lift decrement caused by the bump becomes larger. The influence of the bump on the drag coefficient is distinctly different under the nonbuffet and buffet conditions. The drag coefficient is increased at α = 2.5°, and it is reduced at α = 3.5°. Consequently, the lift-to-drag ratio is increased at the buffet case α = 3.5°, except for the highest bump (Bump 3), while the lift-to-drag ratio is decreased at the nonbuffet case α = 2.5°. Consequently, the bump is beneficial for the buffet case and might be harmful for the nonbuffet case, which coincides with our previous study based on Reynolds-averaged Navier-Stokes computation [27]. The pitching moment is also listed in the table. The reference points of the moment are x/c = 0.25 and y = 0. The negative value means a nosedown moment. It can be seen that the absolute values of the moment are all decreased by the SCBs; however, the variations are not significant. The computational settings are the same as the baseline OAT15A airfoil in Section 2. The computational grids of the bump configurations are deformed from the baseline grid, which keeps the computation appropriate to compare with the baseline OAT15A airfoil.

Time-Averaged Characteristics of the Airfoil with Bumps
The aerodynamic performance of the airfoil with bumps was computed and compared with the baseline airfoil. The nondimensional time step is ∆tU ∞ /c = 1.37 × 10 −4 . The flow field of the baseline airfoil is used as the initial condition of the bump cases. More than 130 time units (tU ∞ /c) are computed for each case. The last 80 time units are used for time averaging to obtain the aerodynamic coefficients. Six cases were computed, which are α = 3.5 • for all four bumps and α = 2.5 • for Bump 2 and Bump 4. Figure 11 shows the histories of the lift and drag coefficients in the computation. It is obvious that the periodic oscillation of the lift and drag coefficients at α = 3.5 • is suppressed by Bumps 2, 3 and 4, while Bump 1 has very little effect on the lift oscillation. It can be seen that the bumps may have a negative influence on the lift coefficient. The lift coefficient is more or less decreased at α = 3.5 • , while it is significantly decreased at the nonbuffet condition α = 2.5 • . Moreover, the drag coefficient is also influenced by the bumps. The drag coefficient at α = 2.5 • is notably increased, as shown in Figure 11d.
The time-averaged values of the aerodynamic coefficients are listed in Table 1. The relative variations (in percentage) in the coefficients compared with the baseline configuration are also listed in the parentheses of the table. A tendency of the bumps is that the lift coefficient is consistently decreased by the four bumps at both α = 2.5 • and α = 3.5 • . The lift coefficient at α = 2.5 • is decreased by approximately 5%, and it is decreased by approximately 1% to 5% at α = 3.5 • . With increasing bump height (Bumps 1, 2 and 3), the lift drop becomes increasingly significant. If the location of the bump is moved upstream (Bumps 2 and 4), the lift decrement caused by the bump becomes larger. The influence of the bump on the drag coefficient is distinctly different under the nonbuffet and buffet conditions. The drag coefficient is increased at α = 2.5 • , and it is reduced at α = 3.5 • . Consequently, the lift-to-drag ratio is increased at the buffet case α = 3.5 • , except for the highest bump (Bump 3), while the lift-to-drag ratio is decreased at the nonbuffet case α = 2.5 • . Consequently, the bump is beneficial for the buffet case and might be harmful for the nonbuffet case, which coincides with our previous study based on Reynolds-averaged Navier-Stokes computation [27]. The pitching moment is also listed in the table. The reference points of the moment are x/c = 0.25 and y = 0. The negative value means a nose-down moment. It can be seen that the absolute values of the moment are all decreased by the SCBs; however, the variations are not significant. A more interesting phenomenon shown in Table 1 is that regardless of whether the lift and drag are increased or decreased, the RMSs of the coefficients are all suppressed by the bumps. With increasing bump height, the suppression effect of the lift oscillation is more significant. In the present computation, Bump 2 has the highest increment of L/D α = 3.5° and a good suppression effect on the shockwave oscillation. The L/D of the buffet case is increased by 5.9%, and the RMS of the lift coefficient fluctuation is decreased by 67.6%; consequently, it is considered a good SCB design, with an appropriate height and location.   A more interesting phenomenon shown in Table 1 is that regardless of whether the lift and drag are increased or decreased, the RMSs of the coefficients are all suppressed by the bumps. With increasing bump height, the suppression effect of the lift oscillation is more significant. In the present computation, Bump 2 has the highest increment of L/D α = 3.5 • and a good suppression effect on the shockwave oscillation. The L/D of the buffet case is increased by 5.9%, and the RMS of the lift coefficient fluctuation is decreased by 67.6%; consequently, it is considered a good SCB design, with an appropriate height and location. Figure 12 shows the time-averaged Mach contour of the four bumps. Bump 1 has the lowest height, meaning the shock pattern is similar to the baseline configuration in Figure 7a. A clear "λ"-type shock is formed near the shockwave foot of Bump 2 and Bump 4. However, the "λ" shock of Bump 4 is stronger than that of Bump 2, which results in more lift loss and less drag reduction. Bump 3 has the highest bump height. As shown in Figure 12c, the bump pushes the first shock upstream and forms a second shock on the bump; consequently, the aerodynamic performance of this bump is worse than the others.
Aerospace 2021, 8, x FOR PEER REVIEW 10 of 18 Figure 12 shows the time-averaged Mach contour of the four bumps. Bump 1 has the lowest height, meaning the shock pattern is similar to the baseline configuration in Figure  7a. A clear "λ"-type shock is formed near the shockwave foot of Bump 2 and Bump 4. However, the "λ" shock of Bump 4 is stronger than that of Bump 2, which results in more lift loss and less drag reduction. Bump 3 has the highest bump height. As shown in Figure  12c, the bump pushes the first shock upstream and forms a second shock on the bump; consequently, the aerodynamic performance of this bump is worse than the others.  Figure 13 presents the pressure coefficients of the airfoil with bumps. For the drag reduction cases (α = 3.5°), the strong shockwave of the baseline configuration is divided into two weak shockwaves. In contrast, for the drag increment cases (α = 2.5°), the shockwave is also divided, but the strength of the second shockwave is even stronger than that of the baseline airfoil.   Figure 13 presents the pressure coefficients of the airfoil with bumps. For the drag reduction cases (α = 3.5 • ), the strong shockwave of the baseline configuration is divided into two weak shockwaves. In contrast, for the drag increment cases (α = 2.5 • ), the shockwave is also divided, but the strength of the second shockwave is even stronger than that of the baseline airfoil.  Figure 12 shows the time-averaged Mach contour of the four bumps. Bump 1 has the lowest height, meaning the shock pattern is similar to the baseline configuration in Figure  7a. A clear "λ"-type shock is formed near the shockwave foot of Bump 2 and Bump 4. However, the "λ" shock of Bump 4 is stronger than that of Bump 2, which results in more lift loss and less drag reduction. Bump 3 has the highest bump height. As shown in Figure  12c, the bump pushes the first shock upstream and forms a second shock on the bump; consequently, the aerodynamic performance of this bump is worse than the others.  Figure 13 presents the pressure coefficients of the airfoil with bumps. For the drag reduction cases (α = 3.5°), the strong shockwave of the baseline configuration is divided into two weak shockwaves. In contrast, for the drag increment cases (α = 2.5°), the shockwave is also divided, but the strength of the second shockwave is even stronger than that of the baseline airfoil.   Figure 14 shows the time-averaged wall friction coefficient of the airfoil. It is worth noting that the wall friction is computed by the wall-modeled equation. The lower surface is hardly influenced by the upper surface bumps, meaning only the lower surface c f of the baseline airfoil is plotted in the figure. The wall friction coefficient of the upper surface is strongly affected by the bumps. The friction first declines when approaching the bump and then increases on the bump, which coincides with the tendency of the flow acceleration on the bump (Figure 12).  Figure 14 shows the time-averaged wall friction coefficient of the airfoil. It is worth noting that the wall friction is computed by the wall-modeled equation. The lower surface is hardly influenced by the upper surface bumps, meaning only the lower surface cf of the baseline airfoil is plotted in the figure. The wall friction coefficient of the upper surface is strongly affected by the bumps. The friction first declines when approaching the bump and then increases on the bump, which coincides with the tendency of the flow acceleration on the bump (Figure 12).     Figure 14 shows the time-averaged wall friction coefficient of the airfoil. It is worth noting that the wall friction is computed by the wall-modeled equation. The lower surface is hardly influenced by the upper surface bumps, meaning only the lower surface cf of the baseline airfoil is plotted in the figure. The wall friction coefficient of the upper surface is strongly affected by the bumps. The friction first declines when approaching the bump and then increases on the bump, which coincides with the tendency of the flow acceleration on the bump (Figure 12).

Fluctuation Characteristics of the Airfoil with Bumps
The fluctuations of the baseline airfoil and airfoil with SCBs are compared in this section. The RMS of the pressure fluctuation on the upper surface is shown in Figure 16. It is clear that the pressure fluctuation of the baseline airfoil is significantly suppressed by the bumps. Bump 2 has the smallest RMS value of the pressure fluctuation among the five cases, which is coincident with the fluctuation of the lift and drag coefficients. Not only the fluctuation near the shockwave but also the fluctuation near the trailing edge is suppressed by Bump 2. Bump 3 and Bump 4 also have favorable suppression effects on the pressure fluctuation.

Fluctuation Characteristics of the Airfoil with Bumps
The fluctuations of the baseline airfoil and airfoil with SCBs are compared in this section. The RMS of the pressure fluctuation on the upper surface is shown in Figure 16. It is clear that the pressure fluctuation of the baseline airfoil is significantly suppressed by the bumps. Bump 2 has the smallest RMS value of the pressure fluctuation among the five cases, which is coincident with the fluctuation of the lift and drag coefficients. Not only the fluctuation near the shockwave but also the fluctuation near the trailing edge is suppressed by Bump 2. Bump 3 and Bump 4 also have favorable suppression effects on the pressure fluctuation.  Figure 17 shows the RMS contours of the pressure fluctuations for Bump 2 and Bump 4. It is clear that the high-fluctuation region near the shockwave is decreased when compared with the baseline configuration (Figure 7b). The "λ" shock region has a high pressure fluctuation, while the shockwave above the "λ" shock is quite stable. It is interesting that the shockwave far away from the airfoil fluctuates greatly. This demonstrates that the expansion wave forms near the leading edge region, weakening the strength of the shockwave and inducing unsteadiness, which is a basic principle of the design of a supercritical airfoil [50]. Three streamwise locations on the airfoil are sampled to compare the power spectral density (PSD) of the pressure fluctuation, as shown in Figure 18. The locations x/c = 0.5 and 0.6 are near the shockwave, and x/c = 0.8 is located in the low-speed region near the trailing edge. The dimensional quantities in Figure 18 are scaled based on the length of the experiment [14]. It is clear that a peak frequency appears on the baseline airfoil, which is the buffet frequency of the baseline airfoil. With the SCBs, the PSD of the pressure fluctuation is significantly reduced at x/c = 0.5. The buffet frequency of the baseline configuration is clear even at the downstream x/c = 0.8. Bump 1 has the lowest height, meaning  Figure 17 shows the RMS contours of the pressure fluctuations for Bump 2 and Bump 4. It is clear that the high-fluctuation region near the shockwave is decreased when compared with the baseline configuration (Figure 7b). The "λ" shock region has a high pressure fluctuation, while the shockwave above the "λ" shock is quite stable. It is interesting that the shockwave far away from the airfoil fluctuates greatly. This demonstrates that the expansion wave forms near the leading edge region, weakening the strength of the shockwave and inducing unsteadiness, which is a basic principle of the design of a supercritical airfoil [50].

Fluctuation Characteristics of the Airfoil with Bumps
The fluctuations of the baseline airfoil and airfoil with SCBs are compared in this section. The RMS of the pressure fluctuation on the upper surface is shown in Figure 16. It is clear that the pressure fluctuation of the baseline airfoil is significantly suppressed by the bumps. Bump 2 has the smallest RMS value of the pressure fluctuation among the five cases, which is coincident with the fluctuation of the lift and drag coefficients. Not only the fluctuation near the shockwave but also the fluctuation near the trailing edge is suppressed by Bump 2. Bump 3 and Bump 4 also have favorable suppression effects on the pressure fluctuation.  Figure 17 shows the RMS contours of the pressure fluctuations for Bump 2 and Bump 4. It is clear that the high-fluctuation region near the shockwave is decreased when compared with the baseline configuration (Figure 7b). The "λ" shock region has a high pressure fluctuation, while the shockwave above the "λ" shock is quite stable. It is interesting that the shockwave far away from the airfoil fluctuates greatly. This demonstrates that the expansion wave forms near the leading edge region, weakening the strength of the shockwave and inducing unsteadiness, which is a basic principle of the design of a supercritical airfoil [50]. Three streamwise locations on the airfoil are sampled to compare the power spectral density (PSD) of the pressure fluctuation, as shown in Figure 18. The locations x/c = 0.5 and 0.6 are near the shockwave, and x/c = 0.8 is located in the low-speed region near the trailing edge. The dimensional quantities in Figure 18 are scaled based on the length of the experiment [14]. It is clear that a peak frequency appears on the baseline airfoil, which is the buffet frequency of the baseline airfoil. With the SCBs, the PSD of the pressure fluctuation is significantly reduced at x/c = 0.5. The buffet frequency of the baseline configuration is clear even at the downstream x/c = 0.8. Bump 1 has the lowest height, meaning Three streamwise locations on the airfoil are sampled to compare the power spectral density (PSD) of the pressure fluctuation, as shown in Figure 18. The locations x/c = 0.5 and 0.6 are near the shockwave, and x/c = 0.8 is located in the low-speed region near the trailing edge. The dimensional quantities in Figure 18 are scaled based on the length of the experiment [14]. It is clear that a peak frequency appears on the baseline airfoil, which is the buffet frequency of the baseline airfoil. With the SCBs, the PSD of the pressure fluctuation is significantly reduced at x/c = 0.5. The buffet frequency of the baseline configuration is clear even at the downstream x/c = 0.8. Bump 1 has the lowest height, meaning the buffet is slightly suppressed by Bump 1, and the buffet frequency is slightly increased. Bump 2 is the most effective SCB because there is almost no frequency peak for the Bump 2 configuration. the buffet is slightly suppressed by Bump 1, and the buffet frequency is slightly increased. Bump 2 is the most effective SCB because there is almost no frequency peak for the Bump 2 configuration.  Figure 19 shows the streamwise velocity contours of the first grid layers of the baseline airfoil and airfoil with bumps. As the present computation adopted a wall-modeled simulation, the grid spacing of the first grid layer is located at Δy+ of approximately 50, as shown in Figure 2a; consequently, the first grid layer of the LES has a considerably large flow speed. The white regions in the contours of Figure 19 are blanked when u < 0.0. Fifteen instantaneous snapshots with Δt =1.37c/U∞ are plotted along the vertical direction in each figure. It is clear that the high-speed region of the baseline airfoil (Figure 19a) oscillates periodically. When Bump 1 is installed on the upper surface (Figure 19b), the oscillation is weakened, but periodical oscillation still exists. The high-speed regions of Bumps 2, 3 and 4 are almost constant over time. However, strong flow acceleration regions appear at the bump locations of Bumps 3 and 4. Moreover, an instantaneous reverse flow (u < 0.0) can be seen before Bumps 3 and 4 at 0.4 < x/c < 0.5, which means that an inappropriate set of SCBs may induce an unexpectedly low-speed region before the bump. Although Bump 2 has a good control effect for the shock oscillation, an instantaneous reverse flow (u < 0.0) region also exists after the bump, which coincides with the mean flow contour shown in Figure 12b.  Figure 19 shows the streamwise velocity contours of the first grid layers of the baseline airfoil and airfoil with bumps. As the present computation adopted a wall-modeled simulation, the grid spacing of the first grid layer is located at ∆y+ of approximately 50, as shown in Figure 2a; consequently, the first grid layer of the LES has a considerably large flow speed. The white regions in the contours of Figure 19 are blanked when u < 0.0. Fifteen instantaneous snapshots with ∆t = 1.37c/U ∞ are plotted along the vertical direction in each figure. It is clear that the high-speed region of the baseline airfoil ( Figure 19a) oscillates periodically. When Bump 1 is installed on the upper surface (Figure 19b), the oscillation is weakened, but periodical oscillation still exists. The high-speed regions of Bumps 2, 3 and 4 are almost constant over time. However, strong flow acceleration regions appear at the bump locations of Bumps 3 and 4. Moreover, an instantaneous reverse flow (u < 0.0) can be seen before Bumps 3 and 4 at 0.4 < x/c < 0.5, which means that an inappropriate set of SCBs may induce an unexpectedly low-speed region before the bump. Although Bump 2 has a good control effect for the shock oscillation, an instantaneous reverse flow (u < 0.0) region also exists after the bump, which coincides with the mean flow contour shown in Figure 12b.  Figure 20 shows the instantaneous streamwise velocity contours of the airfoils. The unsteady flow separation near the trailing edge is strong for the baseline airfoil. With bump installation, trailing edge separation is suppressed. A "λ"-type shockwave can be seen near Bump 2. A secondary shockwave is formed for Bumps 3 and 4.