Numerical study of flow separation control over a circular hump using synthetic jet actuators

Shear layer flows are fundamental to many areas of fluid mechanics. They are found anywhere, where the flow experiences high velocity gradients, such as flows around aerodynamic or bluff bodies. Understanding their physics as well as trying to control such flows may have many practical applications, including drag reduction,^1,2 turbulence and mixing control, vectoring and increase of propulsion of aircraft, improvement of aerodynamic performance,^3,4 enhancing heat transfer in industrial devices, e.g., heat exchangers or avoiding undesirable structural frequencies in mechanical and civil engineering. Consequently, the research on reattachment of separated flow, partially or completely, to the surface, or delaying the onset of flow separation by controlling the natural instabilities has a potential for increasing the safety levels and system performance in a number of industries, especially in aeronautical and aerospace applications.

Synthetic Jet Actuators (SJAs) technology has been utilized as an excellent and promising means of active flow control in aeronautical and aerospace applications. The effectiveness of SJAs in controlling flow separation in the slat cutback region of a large passenger aircraft demonstrator was assessed in the Clean Sky 2 AFLoNext program,^5,6 leading to Technology Readiness Level (TRL) 3 systems ready for wind tunnel testing.⁷ They have also been tested for the increase in the axial turbulence intensity,⁸ control of dynamic stall,^9,10 vectoring,¹¹ mixing enhancement,¹² and thermal management.^13,14 The main purpose of synthetic jets utilization is generation of a zero-net mass flux with non-zero momentum flux by producing a train of vortex rings that can modify the characteristics of the cross flow, thus affecting the flow control. Using this technology has several advantages, including lack of complicated piping systems, a small size and mass, use of ambient fluid as the working fluid, high bandwidth, very low power consumption, and low cost. The use of this device is not frequently easy because of complex flow physics and huge computational expense of simulations of fluid flow essential for a successful utilization.

Figure 1 shows a typical synthetic jet actuator where an oscillating diaphragm is attached to a cavity equipped with an orifice or a slot. The oscillation can be most simply implemented by using disk-shaped piezoceramic transducers,^15,16 although the use of speakers¹⁷ or mechanical pistons¹⁸ has also been reported. Some configurations use two diaphragms.¹⁹ 

An oscillatory flow is generated by changing the cavity volume at a specific frequency. The fluid is drawn into the cavity from around the orifice during the diaphragm’s down-stroke, while during its up-stroke, it is ejected through the orifice, producing a vortex ring. A vortex ring is produced in each oscillation cycle, and hence a sequence of vortices is generated with time. The vortices spread away from the orifice under their own self-induced velocity fields—while “synthetizing” a jet-like velocity profile.

An accurate prediction of flow separation and its control is still a challenging issue for the Computational Fluid Dynamics (CFD) tools. Therefore, this field has remained a subject of ongoing research studies due to its relevance to many industrial and technological applications. A wide range of numerical approaches like Direct Numerical Simulation (DNS),²⁰ steady Reynolds-Averaged Navier–Stokes (RANS) or Unsteady RANS (URANS) simulations,^21–24 Large-Eddy Simulation (LES),^25,26 Detached Eddy Simulation (DES),²² and Implicit LES (ILES)²⁷ were applied to anticipate the turbulent flow separation and its control over a wall-mounted so-called “NASA hump” (developed by NASA as an experimental benchmark).²⁸ The Reynolds number is calculated based on the freestream velocity U and the hump chord length C.

The DNS simulations of Postl et al.²⁰ have been questioned on account of using unrealistic inflow velocity profiles, narrow span width, and inadequate grid resolution. Viken et al.²¹ carried out a CFD investigation to anticipate the hump model aerodynamic performance utilizing Full Unstructured Navier-Stokes 2-Dimensional (FUN2D), the structured TLNS3D Thin-Layer Navier-Stokes 3-Dimensional (TLNS3D), and structured Computational Fluids Laboratory - 3D (CFL3D) codes by applying the time-accurate RANS method. The numerical results were compared with experimental data of Seifert and Pack.²⁸ A significant decrease of separation bubble size was observed in actuated suction and oscillatory cases. LES results of You et al.²⁵ were in good agreement with the experimental data in predicting the flow quantities, such as mean velocity, reattachment length, and pressure coefficient as well as turbulence statistics. LES results of Saric et al.,²⁶ which employed Smagorinsky model, gave reasonable predictions in the unactuated and actuated cases with steady suction. However, the results show considerable differences in mean velocity profiles from the experimental results provided by Greenblatt et al.^29,30 in the oscillatory jet case. The DES results of Krishnan et al.²² are better in comparison with their RANS results, while there are considerable differences from the experimental data of Seifert and Pack.²⁸ ILES results of Morgan et al.²⁷ were in closer accordance with experimental results, though the applied Reynolds number in their study is around 20% of the experimental Reynolds number.

Hybrid RANS/LES methods have been used by Saric et al.³¹ and Jakirlic et al.³² and provided good comparisons with experiments. Also, Fadai-Ghotbi et al.³³ extended a RANS-LES model using the elliptic-blending approach to consider for the impact of wall kinematic blockage. Gritskevich et al.³⁴ suggested two modifications for hybrid CFD schemes, including Delayed Detached Eddy Simulation (DDES) and DDES with improved wall-modeling ability (IDDES) by calibrating these methods to the k–ω SST (Shear Stress Transport) background RANS model. Both modifications were assessed on a set of separated and attached separated flows, including wall-mounted hump, backward facing step developed channel, periodic hills, and hydrofoil with trailing edge separation.

Cappelli and Mansour³⁵ studied the flow separation over a hump model. Benefits and drawbacks of RANS models, including Spalart–Allmaras, k–ε, k–ω, and k–ω SST, were estimated by the OpenFOAM software. However, only the unactuated case was considered, and the slot flow control cases were not studied. Wang et al.³⁶ carried out a bidimensional study on the effect of SJs on the control of flow separation over a hump model. The vortex dynamics was investigated using flow visualization and Particle Image Velocimetry (PIV) methods. The full elimination of the separation region was reported for optimized actuating frequencies.

Yagiz et al.³⁷ and Duda and Fares³⁸ implemented the lattice Boltzmann method to study the flow separation over a wall mounted hump model. Feng and Wang³⁹ studied the active control of the wakes of a circular cylinder utilizing synthetic jets. The measurement of flow field was carried out using the PIV system. Six types of vortex shedding modes are classified under synthetic jet control. The vortex dynamics analysis showed the regular variations of wake vortex trajectory and the vortex circulation for these typical shedding modes. Fisher et al.⁴⁰ performed 2D URANS simulations for NASA 2D hump design to investigate the potential advantages of varying the inflow and outflow directions of a SJ, referred to as the Bi-Directional SJ (BDSJ), compared to the traditional SJ. Kara et al.⁴¹ studied the possible usage of a Sweeping Jet (SWJ) actuator to delay the flow separation by performing a range of two-dimensional URANS simulations for the two-dimensional NASA hump model with a built-in SWJ actuator design. This was to evaluate the function of the SWJ actuator in comparison with the unactuated and experimental measurements. Koklu^42,43 studied a passive and active separation control with Vortex Generators (VGs) and a steady discrete blowing. The findings revealed that the separation region was significantly contracted using VGs. The results showed that a better control was achieved by the implementation of SWJ rather than suction and zero-net-mass-flow actuators, especially for high amplitude excitations. In a numerical and experimental study, Koklu⁴⁴ analyzed the flow separation control over the NASA hump model by investigating the impact of frequency and excitation amplitude in steady and unsteady modes. Despite the over prediction of separation bubble size by numerical analysis, he reported that the unsteady excitation had a better performance than the steady excitation and slightly outperformed the SWJ actuators. The steady suction was reported as the most influential method.

In another numerical study, Tang and Agarwal⁴⁵ investigated the flow separation control over a NASA hump by considering a uniform blowing jet and a SJ at the Reynolds number of 1 × 10⁶ based on the hump chord and Mach number of 0.09, and by employing k–ω SST and Spalart–Allmaras turbulence models. They found that for the uniform blowing jet and jet velocities in excess of 85 m/s, the separated flow is fully reattached to the surface. However, in the synthetic jet case, the flow becomes fully reattached when the jet velocity exceeds 49 m/s. Kim and Kim⁴⁶ conduced a numerical study of the installation conditions for the fluidic oscillators embedded in a hump surface using three-dimensional RANS equations with a modified shear stress transport model to improve the flow separation control over the hump. Aram and Shan⁴⁷ investigated the synchronization of an array of sweeping exiting jets to improve the flow control on a wall mounted hump utilizing the Improved Delayed Detached Eddy Simulation (IDDES) model. Xu et al.⁴⁸ conducted a CFD analysis of the impact of a Co-flow Jet (CFJ) active flow control on the NASA hump. They used the high fidelity in-house CFD code Flow Acoustics Structure Interaction Package (FASIP) with the two-dimensional URANS equations with the one-equation Spalart–Allmaras model. The flow was fully reattached with the blowing location at 50% of the hump chord and suction location at 70% of the hump chord at C_µ = 0.0077 with the CFJ power coefficient (P_c) of 0.0032 and the energy coefficient (C_E) of 0.0034. The study revealed that the optimum location for CFJ suction, with the lowest energy consumption, is at the location where the slope of the hump surface reaches the minimum value (67% of the hump chord).

The full simulation of the interaction between the external flow field and the flow inside the cavity of SJA requires the numerical solution of unsteady, turbulent flow in time-dependent domains, with sufficiently precise computational schemes in both space and time. The main problem affiliated with the full simulation approaches is the computational cost. The numerical calculation of the cavity flow needs considerable computational resources; occasionally comparable with those are required to resolve the external flow. For geometries fitted with multiple actuators, the mesh requirements for the actuators could significantly exceed those of the external flow and would considerably come up with the computational cost. The space–time accurate calculation needs very small space/time scales, e.g., orifice diameters of the order of a millimeter and actuator vibration of the order of kHz. In addition, it must be emphasized that the reduced-order SJA models not considering sufficiently the mutual interaction with the external flow field are not appropriate for cases where a strong coupling exists, e.g., in numerical simulation of flow control with multiple actuators in close vicinity. Another issue is the three-dimensional effects, which are ignored in most two-dimensional numerical studies due to the computational costs. However, in two-dimensional simulations, the interaction of flow structures in regions between adjacent orifices is ignored. In addition, the three-dimensional effects that are more noticeable in the orifice with decreasing jet-to-cross-flow momentum ratio are ignored, which leads to less accurate computations.

In this paper, three dimensional wall-resolved URANS simulations are performed to predict the flow behavior over a circular hump model equipped with synthetic jets. A different configuration was dictated by the need to easily change the “azimuthal angle” of the position of the SJA array in the current study in order to find its optimum—cf. Ja’fari et al.⁴⁹ However, here the focus is on CFD solutions, for both unactuated and actuated cases, conducted using the OpenFOAM package, and this is also used to generate computational grids. In contrast to many previous research works, in the current model, the SJA cavity has been fully considered in the numerical simulation, which adds to the novelty of the approach. Therefore, the coupling between the external field and the internal flow of SJA is considered that leads to a better prediction of SJA performance for flow control. Also, by performing the three-dimensional simulations, the three-dimensional features of SJAs and interaction of flow structures in regions between adjacent orifices are captured in this study. To facilitate these approaches, new “merging” and “stitching” techniques are implemented using the OpenFOAM package that considerably reduces the computational cost. The numerical results are compared with the experimental data obtained by the Hot Wire Anemometry (HWA) and Particle Image Velocimetry (PIV) measurement methods. This article is arranged as follows: The experimental setup and the computational methodology are discussed in Sec. III. The results are illustrated in Sec. IV. The analysis of the numerical simulation results as well as the comparisons with the experimental data is provided in Sec. V. Finally, a summary of findings and concluding comments are provided in Sec. VI.

This section gives a very brief explanation of the experimental setup (Sec. III A), an illustration of the geometry and the computational grid of the problem (Sec. III B), and the numerical method used in the study (Sec. III C).

The experimental investigations were conducted in the University of Huddersfield closed loop low-speed wind tunnel described already by Ja’fari et al.⁴⁹ This has the top speed of 25 m/s, the test section cross section of 500 × 500 mm² and the length of 1000 mm, and the free-stream turbulence level (at the wind speeds utilized in the current study) of ∼0.66%. The circular hump model has a chord length of 200 mm with the apex height of 30 mm and the radius curvature of 181.7 mm. Its span is 500 mm to fill the width of the test section (cf. Fig. 2). Twelve SJAs are embedded in the model in the spanwise direction; each actuator cavity (with the depth of 5.5 mm and diameter of 32 mm) is equipped with three orifices to inject the vortical structures into the boundary layer. Thus, altogether, there are 36 orifices uniformly distributed across the model span with the orifice-to-orifice distance of 13.33 mm. The orifices are 1 mm in diameter and 1.5 mm in depth (length).

In this research, a Constant Temperature Anemometry (CTA) system (MiniCTA 54 T42) provided by the manufacturer Dantec has been utilized to perform the measurements. The system is composed of a 55P11 single-sensor miniature probe with a tungsten wire of 5 µm diameter, a probe support, a BNC cable, and 50 Ω probe cable to join the bridge unit to the hot wire probe, an analog/digital board, and the anemometer with signal conditioning circuit. The probe can be traversed within the tunnel domain using a 2D in-house-built traversing system based on “oriental motor” (with stepping motors of the AZM46AK model and stepping motor drivers of the AZD-KD model) with the nominal spatial (X–Y) resolution of 0.005 mm. PIV measurements were conducted using a LaVision system consisting of a pulsed Nd:Yag laser (model Nano-L-200-15), a CCD camera (model ImagePro4M) with the Macro-Nikon lens of Sigma 105 mm f/2.8, and a CNC-IMS8 controller (isel Germany AG) used for traversing the camera. LaVison seeding particle producer (model Sv-23113) was implemented to seed the flow field with Di-Ethyl-Hexyl-Sebacat (DEHS) oil. More details of the experimental setup and calibration procedures for both techniques can be found in Jafari.⁵⁰ The experimental setup is shown in Fig. 3.

Clearly, the computational cost of solving a CFD problem for the entire wind tunnel cross section would be prohibitively expensive, especially as the modeling success in this type of problem depends strongly on the ability to resolve the flow within the boundary layer using a fine mesh. Therefore, in the current study, only a 40 mm wide “slice” of the physical domain was considered, with the remaining domain dimensions (1000 mm length and 500 mm height of the test section) being preserved (cf. Fig. 4). In essence, the 40 mm width covers one SJA cavity with three neighboring orifices opening out to the cross flow and allows setting up the symmetry boundary conditions on the side walls of the “slice.” The choice of the width is dictated by an attempt to model possible 3D interactions of the flow structures, emerging from the neighboring orifices, further downstream.

The OpenFOAM has been implemented to generate the computational grid and later for pre-processing and post-processing of the simulation results. Two features of the OpenFOAM, namely blockMesh and snappyHexMesh tools, were utilized to generate the computational grids in both unactuated and actuated cases.

To reduce the computational cost in the unactuated case, the blockMesh tool was used to generate a structured mesh without considering the actuator. To size up the near-wall mesh in the computation, standard correlations for a fully turbulent boundary layer over a smooth flat plate under a streamwise zero pressure gradient were used,⁵¹ 

where δ, Re_x, and C_f are the thickness of the boundary layer, the Reynolds number, and skin friction coefficient, respectively. y⁺, u_τ, and ν are dimensionless wall distance, friction velocity, and kinematic viscosity of the fluid, respectively. For a wall-resolved URANS simulation, the first cell height of the computational domain has to be of the same order of y⁺ = 1. Then, by considering the free stream velocity of 7 m/s at the inlet of the test section, the first cell height at the inlet should be set to 0.04 mm. The generated structured mesh in the X–Y plane is shown in Fig. 5, and the mesh specifications are provided in Table I.

In the actuated case, the approach had to be different due to the complexities of the actuator geometry and the small orifice size. It was preferable to use the best features of both the blockMesh and snappyHexMesh tools by merging and stitching two separate meshes (so-called the masterMesh and slaveMesh) in order to produce the computational grid. The main advantage of this technique is the reduction of computational cost by avoiding the refinement of cells in all directions during the refinement phase using the snappyHexMesh tool, as well as avoiding the refinement of the computational blocks in the neighborhood of orifice blocks using the blockMesh tool.

The mesh is designed with a cell size of 0.04 × 0.04 mm² (X–Y plane) at the bottom left corner of the computational domain. This is expanded by exponential stretching with a ratio of 1.03 in the Y direction and 1.08 in the X direction. The cell size in the spanwise direction for the slaveMesh is 0.04 mm, while for the masterMesh, it is 0.8 mm to reduce computational costs. The final computational grid of the actuated case as well as a detailed view of the slaveMesh that is stitched to the masterMesh is shown in the X–Y plane in Fig. 6.

The specifications of the final mesh for the actuated case are summarized in Table II. It is worth mentioning that the overall number of generated cells, either merely by the snappyHexMesh tool or blockMesh tool (with a minimum number of blocks for actuators to be able to define the structured mesh), are 26 393 697 and 24 546 455, respectively. Therefore, ∼50% reduction in mesh size has been achieved by using merging and stitching techniques.

The flow equations obtained from the conservation of mass and momentum for an isothermal incompressible fluid are as follows:

The flow can be divided into two mean and fluctuating parts for the sake of analyzing the effect of turbulence. This process is called Reynolds decomposition and is the starting point of the Reynolds-Averaged Navier–Stokes (RANS) method. The Reynolds-decomposition would be

where φ represents an arbitrary property with mean part of $φ̄$ and the fluctuating part of φ′. The application of the Reynolds-decomposition to Eqs. (4)–(7) yields the so-called Reynolds-Averaged Navier–Stokes equations as follows:

The Reynolds-decomposition introduces the specific source terms on the right-hand side of Eqs. (4)–(6). These terms represent the so-called Reynolds-stresses, which are unknown. In addition of the above equations, two additional equations from a turbulence model are used to determine the unknown Reynolds-stresses.

Many turbulence models have been proposed to predict the flow behavior; each model has its own benefits and drawbacks depending on the flow characteristics, geometry, and other problem parameters. In the current numerical analysis, the k–ω Shear Stress Transport (SST) model was used to predict the flow behavior in both un-actuated and actuated cases. It was implemented to overcome the drawbacks of the standard k–ω model due to the dependency on the freestream values on both k and ω. The SST model is a combination of k–ε and k–ω models in such a manner that the k–ε model is implemented for regions away from the walls and the k–ω model is used for the near wall regions by a blending function. It can predict flow behavior for the adverse pressure gradient conditions and flow separation. The turbulence specific dissipation rate equation is given by

and the turbulent kinetic energy by

The turbulent viscosity is obtained using

The values of constants in the equations are given in Table III.

The turbulent kinetic energy, k, and the turbulence intensity, I, are related to each other through Eq. (16) as follows:

where u_avg is the mean flow velocity.

For internal flows, such as the wind tunnel flow, the turbulence intensity at the inlet generally depends on the upstream history of the flow. The measured turbulence intensity of the wind tunnel is around 0.66%. The turbulence length scale that represents the size of eddies having big energy in a turbulent flow is set equal to the pitch of the honeycomb positioned in the settling chamber upstream of the contraction part (9.5 mm).

At the boundary representing the wind tunnel test section inlet, the turbulence quantities are defined, while the zero gradient boundary condition is assumed at the test section outlet. The gradient of all quantities in the spanwise direction are set to zero at the side wall of the “slice” in Fig. 3. The fixed velocity value of 7 m/s is prescribed at the inlet, and the atmospheric pressure is assumed as at the outlet. At all walls with no slip boundary conditions, all turbulent quantities, except ω, are set to zero while ω satisfies Eq. (17) proposed by Menter,⁵² 

Here, Δy is the distance to the next cell away from the wall. In the actuated case, the time varying boundary condition is assumed at the inlet of the cavity as follows:

It should be noted that the piezoelectric diaphragm is supposed to move in a piston-like motion. In the experiment, a function generator with an amplifier is used to excite the piezoceramic transducers at the bottom of the cavity, with the Helmholtz frequency of the cavity, f = 960 Hz, and a known voltage Vpp of 23 V. The amplitude U₀ in Eq. (9) is strictly unknown. However, the desired peak exit jet velocity at the orifice outlet is known from experiments. Therefore, to find the amplitude U₀, a separate case is solved where the slaveMesh is used in quiescent flow conditions and the amplitude U₀ is varied until the exit jet velocity at the orifice outlet is matched to the experimental values in an equivalent quiescent condition. For the sake of the presented numerical analysis, a specific experimental case is chosen where the ratio of peak exit jet velocity to the free stream is 1.5, and the angular location of the SJA array is 9.5° measured from the hump apex (cf. Fig. 2).

In the current study, the pisoFoam solver is implemented to solve governing equations of isothermal flow [Eqs. (9)–(14)]; this is a transient solver for incompressible, turbulent flow using Pressure-Implicit with Splitting of Operators (PISO) algorithm. Time steps of 1.5 × 10⁻⁵ and 1 × 10⁻⁶ s are used for unactuated and actuated cases, respectively. The maximum Courant number for the unactuated case is 0.80, while for the actuated case, it fluctuated between 0.27 and 0.75 due to the time varying boundary condition of the velocity at the cavity inlet.

The calculations were carried out utilizing one of the nodes of 3M Buckley Innovation Centre (3M BIC) HPC system at the University of Huddersfield, with 24 cores. The average time required for each iteration for the actuated case was 40 s. It took 2400 node-h = 2400 × 24 CPU h (1 node-h = 1 node × 1 h) to reach a flow time of 0.216 s for the actuated case. Thus, the computations in the actuated case were performed for about 3.5 months.

This section presents selected numerical results. Section IV A is concerned with the flow over the hump model in the unactuated case. Section IV B is relevant to the simulation of SJAs in quiescent conditions (Sec. IV B), which explores the physics of vortex ring formation. Finally, Sec. IV C presents the numerical results for the flow over the hump in the actuated case.

In this case, the solution was stabilized after 12 000 iterations while the time step was considered 0.000015 s, which gives the maximum Courant number of 0.80. Keeping the Courant number below 1 helps to maintain both accuracy and stability of the solution. Figure 7 shows the turbulent kinetic energy profile and the velocity boundary layer profile vs experimental data at X = 0.3 m (before the flow reaches the hump), respectively. The experimental data were obtained using Hot Wire Anemometry (HWA) technique.

Figure 8 shows the time averaged contours of turbulent kinetic energy, velocity, and pressure fields. The turbulent kinetic energy contour is depicted in Fig. 8(a). The instantaneous turbulent kinetic energy was calculated through Eq. (16). Figure 8(b) illustrates the contour of kinematic pressure. It should be noted that the kinematic pressure is defined as the static pressure of the fluid divided by its density, as done in the OpenFOAM for the simplification of the system of equations. Figure 8(c) shows the contour of velocity field. Figure 9(a) shows the time-averaged velocity field, while Fig. 9(b) demonstrates the streamline profiles in the wake region. The profile of streamlines has been obtained by employing an adaptive integrator (Runge–Kutta 4–5). The seeds for the streamlines are generated by using a high resolution line source in the Y-direction, which is passed through the wake region. For further comparisons and analyses, Fig. 10 shows the time averaged fields of velocity vector and velocity magnitude, as well as the streamline profiles from the experimental measurements. For PIV measurements, 100 images were captured in a double frame mode in the unactuated case to calculate the average velocity field.

For the sake of generation of the train of vortex rings in the geometrical configuration, such as that shown in Fig. 6(b), it is important to simulate the movement of the diaphragm at the bottom of the SJA cavity. This is implemented by setting a time varying boundary condition as stated by Eq. (9) at the diaphragm patch (treated as the cavity “inlet”). The amplitude of the excitation wave to maintain the desired velocity magnitude of 10.5 m/s at the orifice outlet is found by using the slaveMesh for quiescent flow conditions. After a few trial-and-error solution cases, the value of 0.0365 m/s (at a frequency of 960 Hz) was found for the velocity amplitude to satisfy the desired jet velocity at the orifices’ outlet. This exit jet velocity maintains the selected Velocity Ratio (VR), i.e., the ratio of exit jet velocity to the free stream velocity, of 1.5 which is used for simulation of SJA with cross flow. It is worth mentioning that the continuity equation does not apply for the SJAs and one cannot simply relate the velocity values at the diaphragm patch to the orifice outlet by continuity equation due to the existence of compressibility effects in a resonant cavity.

Figure 11 depicts the history of the normalwise component of the velocity in the center of the outlet of the middle orifice of SJA. Figures 12(a) and 12(b) depict the pressure and velocity contours in the Y–Z plane cutting through all three orifices. Figure 12(c) depicts the velocity field in the X–Y plane cutting through the middle orifice of SJA. All graphs in Fig. 12 are for the “blowing phase” of the SJA cycle.

Figure 13(a) shows the profile of the vorticity contour in the Z–Y plane that is calculated by Eq. (19) as follows:

Figure 13(b) depicts the vorticity magnitude along the middle orifice wall on the right side of the orifice center (dashed line shown in Fig. 13(a).

In this section, the numerical results are presented for the case where the flow over the hump model is subjected to actuation due to synthetic jets. The location of SJA implementation is important to have a successful flow control operation to influence the boundary layer. The best location to implement the SJA is advised to be somewhere upstream and near to the separation point. The estimated separation point for the unactuated case by the numerical analysis was at X = 0.58325 m. Therefore, the SJA was implemented for a selected distance of 18 mm upstream of the separation point. In the experimental environment, this is obtained by rotating the hump model around the central axis (aligned to the spanwise direction) by the angle of 9.5°. The simulation of the actuated case was performed for only one selected value of SJA position due to the huge computational cost of the problem as mentioned in Sec. III B.

Figure 14 shows the instantaneous velocity fields in the wake region for four selected times of t = 0.05, 0.075, 0.16, and 0.21 s. Figure 15 depicts four samples of instantaneous velocity field obtained from the PIV measurements.

Figures 16 and 17 depict the time averaged velocity field profiles of the actuated case obtained from the numerical simulations and PIV measurements, respectively. The time averaged velocity fields for the numerical solution are computed for the time interval of 0–0.21 s in three planes in the spanwise direction, including planes Z = −0.01333, Z = 0, and Z = 0.01333 m. The planes contain the centers of the orifices.

Figure 18 represents the profiles of the time averaged vorticity magnitude and its standard deviation in adjacent of the wall regions and in the vicinity of SJAs. Figure 19 depicts the profile of time averaged wall shear stress on the hump surface and in the vicinity of SJA. The wall shear stress is computed at the patch containing the hump surface and the test section floor by Eq. (20) as follows:

where R is the shear stress symmetric tensor retrieved from the turbulence model and n is the patch normal vector (into the domain).

Figure 20 shows the profiles of time averaged velocity components above the hump trailing edge in plane Z = −0.01333 m.

In this section, first, the results of the solutions for the unactuated case are discussed (Sec. V A). This is followed by the analysis of the simulation results for the SJAs in quiescent conditions (Sec. V B). Finally, the results of the simulation in the actuated case are discussed in Sec. V C.

Figure 7(a) depicts the turbulent kinetic energy profile obtained from the numerical solution, while Fig. 7(b) shows the comparison between numerical and experimental results in terms of the velocity profile within the boundary layer at X = 0.30 m. The velocity fluctuations reach the maximum in near wall region (Y = 2 mm) but vanish toward the solid wall and the free stream flow on either side of the maximum. It should be noted the turbulence length scale is considered equal to the honeycomb pitch (9.5 mm). This quantity represents the size of eddies having big energy in a turbulent flow in relation to the physical size of the problem. Because the turbulence eddies are restricted by the geometry of the problem, e.g., side walls of the test section, this length scale should logically be smaller than the dimensions of the model and test section. Here, the honeycomb cell size provides an additional dimensional constraint on the turbulence length scale.

The time-averaged results predict the thickness of the boundary layer around 7 mm at a distance of 0.30 m from the test section inlet (X = 0.30 m). The predicted value is in very good agreement with the HWA measurements results as shown in Fig. 7(b)—both indicate the thickness of the boundary layer of 7 mm. It is worth mentioning that that, however, the traverse unit resolution was 0.005 mm, and the first measurement position has to be at the height of 0.25 mm above the floor due to the risk of damaging the hot wire probe in close proximity of the solid wall.

Investigation of the time averaged profiles of the turbulent kinetic energy in Fig. 8(a) reveals that the most energy dissipation occurs in the wake region downstream of the hump apex due to the high velocity fluctuations in this region. Indeed, the conversion of potential energy of the flow to kinetic energy of the eddies occurs in the wake region when the boundary layer is “peeled away” and separated from the surface and subsequently takes the form of vortices and eddies in the wake region. Figure 8(b) shows that the flow experiences the maximum pressure at the hump leading edge due to the sudden impingement onto the hump body. As it is seen, the pressure becomes minimum around the hump apex and the velocity is maximum because of the Venturi effect. The contours of velocity magnitude [Fig. 8(c)] and the velocity profile at X = 0.30 m [as shown in Fig. 7(b)] indicate that the flow is fully developed upstream of the hump leading edge.

Investigation of the velocity field [Fig. 9(a)] and streamline profiles [Fig. 9(b)] shows that the predicted reattachment point with zero value of velocity (which is between the forward and reverse flow downstream of the separation point), is at the location of X = 0.62469 m. Similarly, the investigation of velocity vectors shows that the flow is separated approximately at X = 0.58325 m. This would be obtained by exploring for the position where the velocity vectors just above the wall become normal to the surface of the hump.

An analogous investigation of the flow field from the experimental PIV data [e.g., Fig. 10(a)] shows that the flow separation and reattachment occur at X = 0.5820 m and X = 0.62105 m, respectively. The predicted location of the core of the separation bubble from numerical analysis is in very good agreement with the streamline patterns achieved from postprocessing of PIV measurements as shown in Fig. 10(c). In both cases, the predicted location is around X = 0.6070 m. A comparison of predicted location of critical points of flow in the wake region by experiments (Fig. 10) with the predicted values by numerical analysis (Fig. 9) shows a very good agreement between numerical and experimental data—which also includes the velocity magnitude distribution shown in Figs. 9(a) and 10(b).

The history of the velocity normalwise component in the center of the outlet of the middle orifice of the SJA is depicted in Fig. 11. As it is seen, the fluid velocity during the blowing phase is bigger than the fluid velocity during the suction phase. During the blowing phase, the fluid leaves the actuator through a relatively narrow orifice, hence having a relatively high velocity, characteristic of a jet of fluid. On the other hand, during the suction phase, the fluid is drawn into the cavity from the whole region that exists outside of the orifice. This leads to a difference between cycle halves, with the blowing velocity being bigger than suction velocity. This feature of SJA can be helpful to ascertain the velocity direction and consequently the diaphragm displacement for measurements with Hot Wire Anemometry technique without a need for a displacement transducer alongside limitations of space to implement such facilities.

Investigation of the pressure and velocity contours as well as velocity field in the vicinity of the middle orifice of the SJA in the blowing phase [Figs. 12(a)–12(c)] shows that the fluid reaches the maximum velocity (or minimum pressure) before the orifice exit and the velocity is maximum on the orifice axis where there is no frictional contact. At the orifice outlet, some kinetic energy loss occurs due to the area expansion, and consequently, the exit jet experiences a velocity magnitude drop away from the orifice. The flow is drawn back into the cavity from the environment around the orifice during the diaphragm’s down-stroke. During the diaphragm up-stroke, the fluid is ejected through the orifice, producing a vortex ring as shown in Fig. 12(c). A repetition of such vortex ring generation in consecutive diaphragm oscillation cycles leads to a sequence of vortices that propagate away from the orifice under their own self-induced velocity (thus creating a “synthetic” jet). The vortex ring is also formed inside the cavity, the coupling of external field with the internal flow is captured as it is seen in Fig. 12(c), while this feature has been ignored in SJAs simulations without considering the cavity. Figure 13(a) shows the distribution of vorticity magnitude in the Y–Z plane. The evolution of vortex rings can also be observed from this distribution. Figure 13(b) depicts the vorticity magnitude along the middle orifice wall on the right side of the orifice center. As it is seen, the maximum of vorticity strength occurs at the orifice inlet/outlet and at the orifice wall where the shear rate is high. The vorticity strength fades away at a distance of around 9 mm above the orifice as it is seen in the vorticity magnitude distribution. The main criterion for the formation of a synthetic jet flow is that the vortex rings generated during the blowing phase can travel a sufficient distance away from the orifice exit, and thus cannot be drawn back into the cavity during the suction phase. A non-dimensional parameter, referred to as Stokes number, can be described as the ratio of the unsteady force to the viscous force. It is typically associated with the vortex rings roll-up from SJAs with a circular orifice [Eq. (21)]. Zhou et al.⁵³ reported that that the minimum Stokes number value to allow for the happening of the roll-up of vortex rings is about 10,

In this research, the value of Stokes number is 20.2, showing that the SJAs utilized here can make possible the happening of the roll-up of vortex rings.

In this section, the results of the simulation of the interaction of SJA with uniform cross flow is discussed. The samples of the instantaneous velocity fields shown in Fig. 14 demonstrate the unsteadiness of the flow in the wake region, in the form of wave-like patterns with superimposed eddy structures. Four samples of instantaneous velocity field obtained from PIV measurements are depicted in Fig. 15. A qualitative comparison of the flow features leads to a conclusion that the flow physics is well represented by the URANS solver. Clearly, these unsteady features should be “filtered out” when time averaging of results is conducted, while at the same time revealing the recirculating separation bubble.

The time averaged velocity fields for the actuated case using URANS simulations are shown in Fig. 16 for the three X–Y planes coinciding with the Z coordinates of the three orifices of the SJA. Interestingly, the comparison of velocity fields in the spanwise direction [i.e., between Figs. 16(a)–16(c)] shows that the size of the bubble in not same in the three planes. Therefore, the location of the reattachment point is not constant in the spanwise direction, and the flow is not symmetric with respect to the mid-plane of the computational domain (Z = 0). Figure 16(a) shows that the size of bubble is bigger in plane Z = −0.01333 m, while it is smaller downstream of the middle orifice in plane Z = 0.

Table IV summarizes the comparison between CFD and experimental results for the unactuated case by giving the X-coordinates of separation and reattachment points and the length of the recirculation region. Given the inherent uncertainties in defining the locations of the characteristic points in the flow from PIV, the agreement between CFD predictions and experimental results is indeed very good.

Table V gives similar comparisons between CFD and PIV for the actuated case. In this instance, the comparisons have also been made between different Z planes in the CFD model, corresponding to the locations of the three orifices connected to a single cavity. These were subsequently averaged. Comparing these averaged features from the CFD model with PIV data, it appears that the separation and reattachment points are predicted slightly more upstream by CFD (by about 5–6 mm) while the length of the recirculation zone is predicted relatively well. Also, importantly the reduction of the recirculation zone in response to the actuation is in very good agreement; for CFD, the length is reduced from 0.04144 to 0.02865 m, while for PIV, it appears to be reduced from 0.03905 to 0.03063 m. All these features can be identified by comparing the flow field data presented in Figs. 16 and 17.

The possible reasons for the difference in CFD and PIV results can be relevant to limited size of the interrogation windows used to calculate the velocity vectors based on cross-correlating the intensity distributions. In addition, the PIV method is not able to measure the velocity component along the Z-axis (i.e., toward or away from the camera). Therefore, this component might not only be missed, but could also introduce an interference in the data for the X/Y-components. As already discussed, the bubble size in the spanwise direction is not constant. Some of these issues could be, in principle, resolved by using stereoscopic PIV or 3D “tomographic” PIV methods, but these were not available here. On the other hand, the drawbacks of URANS technique itself should not be forgotten. The URANS models are probably less accurate in capturing transient turbulent structures in comparison to other simulation techniques, such as Large Eddy Simulation (LES). Also, URANS accuracy varies from one place to another place (e.g., in the near wall region with a high velocity gradient and out of the wake region). Finally, the structured mesh is used in the unactuated case while a combination of structured and unstructured meshes is used in the actuated case to reduce the computational cost. However, the accuracy of URANS simulation with structured mesh is higher than in simulations with unstructured mesh.

The profiles of time averaged vorticity magnitude and its standard deviation in near wall regions and in the vicinity of the SJA are depicted in Fig. 18. Figure 18(a) depicts the trace of the pair of counter rotating vortices generated on the hump surface due to the interaction of vortex rings with the cross flow. Indeed, the synthetic jet structures from Fig. 13(a) get torn by the cross flow to create pairs of counter rotating vortices. Figure 18(b) shows the difference between the magnitude of the vorticity relevant to the pairs of the counter rotating vortices and the adjacent flow field in the near wall regions.

Investigation of the profile of time averaged wall shear stress on the hump surface and in the vicinity of SJA reveals that the pair of counter rotating vortices have responded to the surface shear stress and has formed a pattern as it is seen in Fig. 19. It appears that if the distance between the orifices was slightly smaller, the flow actuation generated by the SJA would probably be able to cover the regions with lower shear stress to avoid the inception of instabilities of the flow and the onset of separation. This observation can be used to modify the experimental design and improve the performance of the flow control using SJA. A comparison of the time averaged vorticity and wall shear stress profiles also shows the existence of the relationship between the near wall vortices and wall shear stress. This issue is in agreement with the DNS results by Guo and Li⁵⁴ where they showed that the wall shear stress is associated with the near-wall streamwise vortices. Also, Kravchenko et al.⁵⁵ reported that the higher skin‐friction values were correlated with streamwise vortices positioned nearer to the wall utilizing a DNS database.

Investigation of the profiles of time averaged velocity components above the hump trailing edge in plane Z = −0.01333 m (Fig. 20) shows that the value of streamwise component of the velocity is much bigger than the spanwise and normal components of the velocity in the near wall region. Therefore, a bigger contribution of spanwise vorticity in the vorticity magnitude is expected. The inflection points of velocity profiles are seen in the near wall region that directly influences the shear rate. The magnitude of the normalwise and streamwise vorticities adjacent to the wall is about 2% and 10% of the magnitude of the spanwise vorticity, respectively. Consequently, the spanwise and streamwise vortices as the main characteristic of the turbulent boundary layer in the wake region are the important factors affecting the wall shear stress. This observation is in agreement with the experimental study by Shen et al.⁵⁶ where they showed the relation between the wall shear stress (streamwise and spanwise) and near-wall flow structures (streamwise, spanwise, and outside structures) using digital holographic microscopy.

The delay of boundary layer flow separation is caused by the introduction of vortex rings and, consequently, generation of pairs of counter rotating vortices into the boundary layer by the implementation of SJA. The generated vortical structures interact with the cross flow and transfer the high momentum flow from the outer flow into the near wall regions and the retarded boundary layer is re-energized. The momentum coefficient [Eq. (22)] as a key factor in the SJAs performance is directly proportional to the square of the Velocity Ratio (VR) and the entire area of throats of actuators (A_SJAs),

It should be noted that the A_ref in Eq. (22) is the reference area of the hump model that is the product of the hump span and the hump chord. The momentum coefficient value for one SJA at VR = 1.50 is 3.48 × 10⁻⁵, which is in accordance with the values described in several research (Ciuryla et al.,⁵⁷ Farnsworth et al.,⁵⁸ and Tang et al.⁵⁹). Although the momentum convected to the retarded boundary layer is not considerable compared to the freestream characteristic momentum, the increased momentum is yet effective in alleviating the momentum inadequacy in the turbulent boundary layer and changing the local flow field in such a way that the flow separation is delayed.

Wall resolved URANS simulations of fluid flow over the hump model are performed by employing the k–ω SST model using OpenFOAM software. Merging and stitching techniques are utilized to get the best features of the blockMesh and snappyHexMesh grid generation tools of OpenFOAM software that was very helpful to save computational cost. The full simulation of SJAs was performed by considering the cavity, hence capturing the coupling between the external field and internal flow of the SJAs. Also, performing the three-dimensional simulations made it possible to capture the three-dimensional features of SJAs and interaction of flow structures in regions between adjacent orifices. Results show the effectiveness of SJA on flow control by delaying the flow separation and pushing back the reattachment point toward the hump trailing edge due to changing of flow structures by the interaction of generated vortical structures with separated shear layers. The prediction of flow behavior in the wake region by URANS simulations is in reasonable agreement with PIV results. CFD simulations revealed that the location of the separation and reattachment points is not constant in the spanwise direction and the flow is asymmetric with respect to the mid-plane of the test section (Z = 0). The simulations suggest that a smaller distance between orifices may provide a more effective flow separation control by delaying the onset of instabilities and as a result the delay of flow separation phenomenon.

The numerical analysis showed that in addition of spanwise vortices, the streamwise vortices are influential on the wall shear stress. Although the predictions of separated flows in fluid dynamics field is always challenging, the comparison of current numerical results with experimental data shows that the URANS approach employing the k–ω SST model can be used as a reasonable strategy for prediction of characteristics of separated turbulent flows and their control.

Mohammad Ja’fari acknowledges the EPSRC, UK) (Grant No. EP/R023328/1) for the financial support. The authors are grateful to the 3M BIC Center of University of Huddersfield for access to the High Performance Computing facility, and we are also thankful to Dr Steve Andrews for his support.

The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.

Mohammad Ja’fari: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Artur J. Jaworski: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Aldo Rona: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

A_ref

reference area [Eq. (22)]

A_SJAs

area of SJAs [Eq. (22)]

C_f

skin friction coefficient [Eq. (2)]

C_μ

momentum coefficient [Eq. (22)]

D_o

orifice diameter [Eq. (21)]

f

frequency (Hz) [Eq. (21)]

I

turbulence intensity [Eq. (16)]

k

turbulence kinetic energy (m/s)² [Eqs. (14) and (16)]

n, s

normal and tangential vectors

p

pressure (Pa) [Eqs. (5)–(7) and (10)–(12)]

R

shear stress tensor (Pa) [Eq. (20)]

Re

Reynolds number [Eqs. (1) and (2)]

St

Stokes number [Eq. (21)]

t

time (s)

U

velocity (m/s) [Eqs. (4)–(7) and (9)–(14)]

u, v, w

velocity components (m/s) [Eqs. (5)–(7) and (10)–(12)]

VR

velocity ratio [Eq. (22)]

X, Y, Z

Streamwise, normalwise, and spanwise coordinates

Greek symbols

δ

boundary layer thickness (m) [Eq. (1)]

μ

dynamic viscosity (kg/m s) [Eqs. (5)–(7) and (10)–(12)]

ν

kinematic viscosity (m²/s) [Eqs. (3) and (21)]

ρ

density (kg/m³) [Eqs. (4)–(7) and (10)–(14)]

τ

wall shear stress (Pa) [Eq. (20)]

ω

turbulence specific dissipation rate [Eqs. (13) and (17)], vorticity (1/s) [Eq. (19)]

Subscripts

avg

referring to average

stddev

referring to standard deviation