Hybrid RANS/LES Turbulence Model Applied to a Transitional Unsteady Boundary Layer on Wind Turbine Airfoil

Abstract

A hybrid Reynolds-averaged Navier Stokes/large-eddy simulation (RANS/LES) turbulence model integrated with a transition formulation is developed and tested on a surrogate model problem through a joint experimental and computational fluid dynamic approach. The model problem consists of a circular cylinder for generating coherent unsteadiness and a downstream airfoil in the cylinder wake. The cylinder flow is subcritical, with a Reynolds number of 64,000 based upon the cylinder diameter. The quantitative dynamics of vortex shedding and Reynolds stresses in the cylinder near wake are well captured, owing to the turbulence-resolving large eddy simulation mode that was activated in the wake. The hybrid model switched between RANS and LES modes outside the boundary layers, as expected. According to the experimental and simulation results, the airfoil encountered local flow angle variations up to ±50°. Further analysis through a phase-averaging technique found phase lags in the airfoil boundary layer along the chordwise locations, and both the phase-averaged and mean velocity profiles collapsed into the Law-of-the-wall in the range of $0<y^{+}<50$. The features of high blade-loading fluctuations due to unsteadiness and transitional boundary layers are of interest in the aerodynamic studies of full-scale wind turbine blades, making the current model problem a comprehensive benchmark case for future model development and validation.

1. Introduction

Modern interest in wind energy emerged after the first oil crisis in the 1970s, and a concerted effort followed to develop viable utility-scale wind energy technologies. The megawatt wind turbine debuted in 1998, and by 2010, wind turbines in the 2 MW and 5 MW classes became the standard of onshore and offshore applications, respectively. An even larger 10 MW wind turbine has been designed and tested offshore [1]. The introduction of a 10 MW–20 MW wind turbine with rotor diameter over 200 m was also planned [2].

Wind turbines operate in the eddies and gusts generated by wakes from upstream turbines and atmospheric turbulence. The highly unsteady flow over the rotating blades is characterized by large transient changes in angle-of-attack, laminar-turbulent transition, flow separation and reattachment, and dynamic stall. Most experimental studies use model-scale wind turbines due to the high cost of full-scale testing, and past numerical studies have been conducted via multiple approaches of different fidelities. Rotor-scale large-eddy simulation (LES) of the blade boundary layer for a full-scale wind turbine still exceeds the state-of-the-art [3]. Therefore, lower-fidelity models such as the actuator-line model [4,5] are often used instead. Unfortunately, lower-fidelity models are incapable of predicting the aforementioned unsteady flow phenomena. The uncertainty of performance and loads predictions in the wind energy industry remain high compared with conventional energy sources.

While the LES of full-scale turbulent rotor flow is still too expensive for most practical problems, the next order of fidelity, Reynolds-averaged Navier Stokes (RANS) models cannot replicate many complex flow features. As such, the RANS approach currently falls short in precisely predicting the shaft torque at stall, as laminar/turbulent boundary layer transition and flow separation have both been shown to have first-order effects on blade aerodynamics predictions [5]. Full-scale computational fluid dynamics (CFD) studies have shown that increasing RANS fidelity for separated flows and wake regions by employing detached-eddy simulation (DES) can significantly improve predictions of blade torque near stall compared with RANS [6,7,8], and the application of the hybrid RANS/LES model yielded nearly identical solutions compared with RANS at low wind speed.

One of the major difficulties in full-scale CFD is the large diversity of both temporal and spatial scales in wind plant operation conditions. In the atmospheric boundary layer, the energy-containing motions are on the time scale of minutes and the length scale from 10 to 1000 m, while the attached blade surface viscous boundary layers have time scales on the order of milliseconds with length scales of importance down to 10 μm (e.g., see [3]). Thus, neither RANS nor LES models perform well across the entire time/space range of the full-scale problems.

The hybrid RANS/LES model in this study is based on a correlation-based four-equation RANS transition model, which integrates two local-variable transport equations with the shear stress transport (SST) $k-\omega$ turbulence model framework [9]. It is capable of capturing natural, bypass and separation-induced transitions [10]. An application of the transition model to a small, full-scale wind turbine in the literature showed improved agreement with experiments compared with RANS and engineering models [11]. For a wind turbine operating in the wake of an upstream turbine, blades experience large variations in angle of attack; the transition model alone helps capturing the load hysteresis loop in dynamic stall [12]. A recent study [13] assessed the effects of blade boundary layer transition inside atmospheric turbulence eddies through blade-resolved simulation and full-scale field experiments. The authors found that the boundary layer separation was closely correlated with the transition locations in off-design conditions, where the blades experienced massive separations and deep stall [13]. Despite the improvement, however, further development is needed for properly predicting the reattachment process.

In the current effort, a new surrogate model problem (Figure 1) is developed; the cylinder-airfoil design encompasses key flow physics such as large-scale, rotational unsteadiness and a transitional boundary layer over an aerodynamic body. A number of past experimental studies involving circular cylinders have focused on various flow problems, including aeroacoustics [14], spanwise correlation [15], vortex and wake dynamics [16,17], and drag reduction [14]. The 3D nature of turbulent flow causes phase shift in the spanwise direction, which accelerates the flow structure break-down and turbulent energy cascade [15]. A cylinder-plate configuration is one of the well-studied wake/body interaction problems, and experimental results have indicated strong correlations between the passage of wake vortices and flat-plate boundary layer dynamics, including the formation of secondary flow structures and transition onset [16]. Of interest to the wind turbine problem motivating the current work, the spacing between the upstream circular cylinder and downstream plate has a major effect on the dynamic loads on the plate. For instance, when plate is supported so that it is free to rotate around its centerline, the flow-driven oscillations of the plate at small spacings are dominated by the passage of shed vortices from the upstream cylinder, while larger spacings show dominance of self-forcing from wake features shedding off the plate [18]. It is interesting to note that equilibria, as well as major oscillations around these equilibria, are present for both bluff and aligned plate orientations.

The study on the effect of cylinder wake on a downstream airfoil is of particular interest for our research. The highly-turbulent wake of the cylinder induces transition of the airfoil boundary layer and postpones separation on the airfoil. Past studies found that the attached boundary layer further reduces the vorticity magnitude in the airfoil wake, which improves its aerodynamic characteristics [14]. Detailed experimental results have been obtained [17], providing quantitative data for comparison and insights that guide model developments. Leveraging the new insights, we developed a novel hybrid RANS/LES approach which includes transition modeling, and, thus, has the advantages of both resolving the energy-containing turbulence structures generated upstream of the blade and also improving the fidelity of physics modeled in the blade boundary layer region. Some preliminary computational results of circular cylinder flow and airfoil boundary layer were published previously [19].

The objectives of this study include: (1) design of a small-scale surrogate wind-tunnel model using a collaborative CFD-EFD (experimental fluid dynamics) approach, capturing key flow physics; (2) perform high-fidelity experiments which can provide validation-quality data and guide the development of turbulence models; (3) incorporate the Langtry-Menter 4-equation transition model [20,21] into a hybrid RANS/LES framework and validate resulting computational approach; and (4) provide analysis of the results to further our understanding of unsteady transitional boundary layers on wind-turbine airfoils.

Section 2 introduces the design of the surrogate model and its model parameters. Section 3 presents the turbulence model used in the current study and the corresponding simulation setups. Section 4 contains simulation results of cylinder flow using the hybrid transition model, also including the comparison with experimental data. Section 5 is an analysis of the unsteady airfoil boundary layer and comparison between the simulation and experimental results. Finally, Section 6 concludes the paper.

2. Model Formulation and CFD-based Experimental Design

The aerodynamics of wind plants consists of eddies, boundary-layers, and wakes which fall in a vast range of length and time scales (

Table 1. Range of scales and modelling techniques for full-scale wind turbines.

Time Scale (s)	Length Scale (m)	Blade Re	Reduced Frequency	Modeling Technique
$\sim \mathcal{O}\left({10}^{-3}-{10}^1\right)$	$\sim \mathcal{O}\left({10}^{-6}-{10}^2\right)$	$\sim \mathcal{O}\left({10}^7\right)$	$\sim \mathcal{O}\left({10}^{-2}-{10}^0\right)$	CFD, actuator methods, BEM

). The full-scale problem is too complicated and computationally expensive for fundamental model development and verification, validation, and uncertainty quantification (VVUQ). A typical study of a 5 MW full-scale single rotating blade in the atmosphere with well-resolved boundary layers requires over 50 million cells and thousands of processors, and each intermediate output at a single time step from the CFD solver occupies over 20 GB of storage [3]. Therefore, we are motivated to develop a surrogate model.

The current study adopts a complementary CFD-EFD approach, and preliminary simulation data are used to guide the experimental design. The experimental data is used for turbulence model validation purposes, while the comparison and analysis of numerical and experimental data lead to the fully-described flow physics of the unsteady boundary layer.

An overview of the experimental and numerical procedure is shown in Figure 2. The flow chart consists of red blocks for numerical work, blue blocks for experimental work and black blocks for joint efforts. All the experiments and corresponding analysis were conducted by co-author D.R. Cadel and the details can be found in the article by Cadel et al. [22].

The surrogate model in our study consists of a modified NACA 6-series airfoil [23] embedded in the wake of a circular cylinder, which is used as the surrogate atmospheric turbulence generator. Figure 3 is a sketch of the experiment setup in the test section, the freestream velocity is at $U\infty = 26 \mathrm{m}/\mathrm{s}$ and turbulence intensity $Ti = 1\%.$ The diameter of the circular cylinder is $D = 3.8 \mathrm{cm}$, and the airfoil chord length c = 10.2 cm, which corresponds to $R{e}_D = 6.4 \times {10}^4$ and $R{e}_c=1.7\times {10}^5$.

The key parameter of this surrogate model is the spacing L from cylinder center to airfoil leading edge, in order to minimize the direct effects on unsteady boundary layer physics from the upstream cylinder, the airfoil is placed outside the immediate wake of the cylinder. Various spacings from $1.5D$ up to $10D$ were test using the potential flow solver in OpenFOAM (Version 3.0), and the pressure recovery curves are plotted from cylinder rear stagnation point to airfoil leading edge as shown in Figure 4, as the airfoil moves further than 9$D$ downstream the pressure recovers to near freestream level ($C_p = 0$) at $\sim 4D$. The furthest airfoil position is also limited by the experimental equipments, e.g., the total length of the test section and particle image velocimetry (PIV) camera blockage from test section supporting brackets, the spacing is finally decided at $L = 40.6 \mathrm{cm}$, i.e., $10.67D$. Flow conditions and model parameters are summarized in

Table 2. Test section flow conditions and surrogate model parameters.

Freestream Conditions	Test Section Dimensions	Circular Cylinder	NACA63215b Airfoil
$U_{\infty }=26 \mathrm{m}/\mathrm{s}$ $Ti=1\%$	Streamwise (x): $121.9\mathrm{cm}$	Diameter $D=3.81\mathrm{cm}$	Chord $c=10.16 \mathrm{cm}$
	Transverse (y): $71.1 \mathrm{cm}$	$R{e}_D=6.4\times {10}^4$	$R{e}_c=1.7\times {10}^5$
	Spanwise(z): $71.1 \mathrm{cm}$	Spacing $L=10.7D$	Reduced frequency $k=1.5$

.

The validation experiments are conducted in Virginia Tech’s 0.7 m blower type subsonic wind tunnel with a closed $71.1\mathrm{cm} \times 71.1 \mathrm{cm}\times 127 \mathrm{cm}$ acrylic-walled test section. Two laser/camera systems were used to acquire PIV data. For the cylinder wake measurements and time resolved boundary layer measurements, a diode-pumped Nd-YLF frequency doubled 527 nm laser (Photronics Industries International, Inc., Ronkonkoma, NY, USA) was used, with a maximum pulse energy of 70 mJ at 2920 Hz. Photron Fastcam SA-1 cameras with 1024 × 1024 pixel sensors were used with this laser. For the non-time resolved airfoil cross-section measurements, an EverGreen 200 Nd:YAG laser (Quantel, Bozeman, MT, USA) was used with a maximum pulse energy of 200 mJ at 4 Hz. Two Imager pro X 4M double shutter charge-coupled device (CCD) cameras (LaVision, Ypsilanti, MI, USA) with 2048 × 2048 sensor size were used with this system. Due to low stereoscopic sensitivity of the spanwise component in the boundary layer measurements only the in-plane velocity data are analyzed. DEHS oil (mean diameter 1 μm) was used for cylinder wake measurements, and a smoke consisting of 50% pharmaceutical grade glycerin and 50% deionized water (0.2–0.3 μm diameter) was used for all other PIV data. Seed was introduced immediately downstream of the fan blower, well upstream of the honeycomb grids, and was seen to be dispersed and uniform in the test section. Higher seed density is introduced upstream for the airfoil unsteady boundary layer measurements due to high magnification boundary layer planes.

Owing to the limited spatial resolution, the uncertainty in PIV data is estimated at 1% of the velocity magnitude [17]. Both the airfoil and cylinder models are chrome-plated in order to mitigate laser flare during PIV testing. Flow field data is acquired in both the cylinder near wake and on the suction side of the airfoil. Extra effort is made to obtain airfoil near-surface measurements (see the airfoil near-surface planes depicted in Figure 5). For the cylinder wake measurements, 10,658 total realizations were acquired over four runs, and were subsequently normalized using the measured plenum pressure in order to combine sets for analysis. Airfoil measurements over the full chord were acquired with 2000 realizations for each inflow case. For the boundary layer profile measurements, 2728 images were used for the steady inflow case at each station, and 5456 images were used for the periodic inflow case at each station. A 40 mm × 40 mm field of view with grid spacing of 0.6 mm was acquired for the cylinder wake measurements. Airfoil measurements over the full chord had a grid spacing of 0.5 mm. For the boundary layer profile measurements, 15 mm × 15 mm fields of view with grid spacing of 0.23 mm. Many further details of the experiment are discussed in past published works [22,24].

Although the airfoil Reynolds number is two orders of magnitude lower than full-scale problem, and the imposed unsteadiness from upstream cylinder is much higher in frequency, our surrogate model still captures key features of the unsteady physics as seen in full-scale problems. A future experiment at near-full-scale conditions is planned in the Virginia Tech Stability Wind Tunnel (183 cm × 183 cm cross section).

3. Turbulence Model and Simulation Setups

3.1. Turbulence Model

3.1.1. Delayed-Detached-Eddy Simulation

Our hybrid RANS/LES transition model uses the $k-\omega$ SST as the base model and integrates the Delayed Detached-Eddy Simulation (DDES) formulation with correlation-based four-equation transition model of Langtry and Menter, which was implemented in OpenFOAM by Nandi [13] as part of the PSU CyberWind Facility.

The DDES length scale $\tilde{d}$ blends local RANS lengths scale with the LES length scale via a blending function $f_d$ that depends on local flow solution [25]:

$$\tilde{d}\equiv l_{RANS}-f_d\max \left(0, l_{RANS}-l_{LES}\right)$$

(1)

where $l_{RANS}$ is the RANS length scale. The DES coefficient $C_{DES}=0.65$, and $\mathsf{\Delta }\equiv \sqrt[3]{\mathsf{\Delta }x\cdot \mathsf{\Delta }y\cdot \mathsf{\Delta }z}$ or $\mathsf{\Delta }\equiv \max \left(\mathsf{\Delta }x, \mathsf{\Delta }y,\mathsf{\Delta }z\right)$. The latter definition of $\mathsf{\Delta }$ is usually preferred as the measurement of local grid size for its improved performance on high-aspect-ratio cells. $f_d$ is the blending function:

$$f_d\equiv 1-\tanh {\left(20{\mathrm{r}}_{\mathrm{d}}\right)}^3$$

(2)

From Equation (2), it is clear that the hybrid model operates in RANS mode as $f_d=0$. It switches to LES mode when $f_d=1$ and $\tilde{d}=C_{DES}\mathsf{\Delta }$, which is equivalent to the original DES97 formulation [26]. ${\mathrm{r}}_{\mathrm{d}}$ is designed to work with the eddy-viscosity model:

$$r_d\equiv \frac{{\nu }_t+\nu }{\sqrt{U_{i,j}{U}_{i,j }}{\kappa }^2{d}^2 }$$

(3)

where ${\nu }_t$ and $\nu$ are eddy and molecular kinematic viscosity respectively. $\mathsf{\kappa }=0.41$ is the Karman constant and $d$ is the nearest wall distance.

The hybrid length scale $\tilde{d}$ replaces $d$ in the destruction term of turbulence kinetic energy (TKE) equation. As the model switches to LES mode, TKE is resolved down to the local grid size. Compared with original DES97 [20], the DDES model delays the mode transition outside the boundary layer by applying $f_d$. Several major issues, e.g., modeled stress depletion, that limited the application of DES97 are significantly alleviated.

Since the overall effect of adopting the hybrid length scale $\tilde{d}$ is reducing the turbulence length scale, leading to higher dissipation of TKE, the above formulation is transformed into a function $F_{DDES}$:

$$F_{DDES}\equiv \frac{{\mathrm{D}}_{\mathrm{k}}^{\prime }}{{\mathrm{D}}_{\mathrm{k}}}$$

(4)

where ${\mathrm{D}}_{\mathrm{k}}^{\prime }$ is the hybrid destruction term and ${\mathrm{D}}_{\mathrm{k}}$ is the RANS destruction term, both are inversely proportional to the turbulence length scale. F_DDES is directly multiplied with the dissipation term in TKE equation, i.e., for $F_{DDES}=1$, the original RANS formulation is recovered; for $F_{DDES}>1$, the turbulence model switches to LES mode.

3.1.2. Transition Model Formulation

The transition model [9] introduces two extra transport equations for turbulence intermittency $\gamma$ and transition momentum thickness Reynolds number $\widetilde{R{e}_{\theta }}$. The latter is used to calculate the critical momentum thickness Reynolds number ${\widetilde{Re}}_c$, and it marks the location that intermittency first starts to increase in the boundary layer. The use of empirical correlations in relating ${\widetilde{Re}}_{\theta }$ to ${\widetilde{Re}}_c$ eliminates non-local dependence, which is particularly important in today’s large-scale parallel computing. The production of TKE is multiplied by the turbulence intermittency $\gamma$, which varies from 0 (laminar) to 1 (fully turbulent).

Statistically, the turbulence intermittency is the ratio of the total amount of TKE produced over time to the TKE predicted by the turbulence model (fully-turbulent). Therefore, the intermittency filed is multiplied with the TKE production of fully-turbulent prediction. However, in current model, the intermittency field is calculated by a transport equation as other turbulence fields. Instead of representing a statistical state of flow, the turbulence intermittency field $\gamma$ only serves as a means to reduce TKE production. Two modifications deviate $\gamma$ from its physical definitions: 1) $\gamma$ = 1 in freestream where flow is laminar 2) $\gamma$ > 1 in separated flow. The first modification enables the model to predict the effects of large freestream turbulence level on the laminar boundary layer and the second modification helps capturing the separation-induced transition by making TKE grows rapidly after laminar separation. In order to distinguish $\gamma$ from the physical definition of turbulence intermittency, it is called the Langtry-Menter (LM) turbulence intermittency in the current study.

In order to prevent unphysical damping and dissipation of turbulence kinetic energy in the freestream upstream of the cylinder, a weak source term is added to the turbulence model equation which cancels the dissipation term in the freestream [27]. Inside boundary layer, the production and dissipation of turbulence fields are several orders of magnitude higher than the added source term and, thus, its effect is negligible.

3.2. Simulation Setups

All simulations in our study are conducted using OpenFOAM, an open-source CFD toolbox written in C++ which includes solvers, utilities, and libraries for finite-volume discretization of the field equations of computational mechanics. Transient, incompressible flow is solved using a hybrid PISO-SIMPLE algorithm, RANS and hybrid RANS/LES turbulence models. Block-structured meshing is accomplished using Pointwise. Simulations were conducted on the Virginia Tech Advanced Research Computing supercomputers BlueRidge (408-node Cray CS-300 cluster) [28] and Cascade (196-node 32-core Broadwell processors) [29].

The computational domain is designed to match the experimental facility and, therefore, has the same width as the test section ($71.1 \mathrm{cm}$). In the spanwise direction, periodic boundary condition is used, and the spanwise dimension is set to one chord length. Figure 6 shows both an overall and zoomed view of the mesh. In order to capture transition, the near-wall cells on the cylinder, airfoil and tunnel walls are kept at height of $y^{+}\approx 1$, where $y^{+}=y{u}_{\tau }/\nu$ and $u_{\tau }$ is the wall friction velocity. We use 40 cells in the spanwise direction in order to keep the cell size between the cylinder and airfoil nearly isotropic to better capture the scales of the energy-containing eddies in the wake of the cylinder. The total cell count is approximately 10 million, and the simulation runs on 320 cores at fixed time step $\mathsf{\Delta }t\approx 3\times {10}^{-6} \mathrm{s}$. The time step size is chosen to keep the maximum Courant number less than 1, which is essential for resolving the boundary layers. Further reducing the time step significantly increases the computation time without discernible change in solution.

4. Circular Cylinder Flow

Flow over circular cylinder is a well-studied problem, abundant experimental and numerical data is available from low-Reynolds number laminar flow up to supercritical high-Reynolds number turbulent flow. In our test, the circular cylinder Reynolds number ${\mathrm{Re}}_{\mathrm{D}}=6.4\times {10}^4$, which is in the subcritical flow regime. Cylinder experiences laminar separation and flow bursts into turbulence inside the shear layer. Alternative vortices are formed from shear layer roll-up and shed downstream with constant Strouhal number. Beyond the subcritical regime, the cylinder boundary layer transition initiates before separation and delays the separation points, which significantly reduces the pressure drag and total drag, thus the named “drag crisis”.

The objective of studying circular cylinder flow is to test the hybrid transition model performance and confirm that we can accurately predict unsteady flow off the cylinder. The spanwise dimension is fixed at $2.67D$ (one chord length), and a grid independence study was conducted using three systematically refined mesh. Grid statistics of the coarse, medium and fine mesh is summarized in

Table 3. Mesh statistics and force coefficients of the grid independence study.

Grid	Average ${\mathit{y}}^{+}$ on Surface	Cells around Cylinder	Total Cell Count	$\overline{{\mathit{C}}_{\mathit{D}}}$	${\mathit{C}}_{\mathit{L}}^{\prime }$	$\mathit{S}\mathit{t}$
Coarse	$\approx 2\left(y=2\times {10}^{-5}\mathrm{m}\right)$	$120$	328,680	0.64	0.16	0.256
Medium	$\approx 1\left(y=1\times {10}^{-5}\mathrm{m}\right)$	240	2,768,566	1.31	0.67	0.193
Fine	$\approx 0.5\left(y=5\times {10}^{-6}\mathrm{m}\right)$	480	21,790,080	1.35	0.71	0.190

.

All three simulations use the hybrid transition model as introduced above, same numerical schemes are used to solve the linear system with fixed time step of $3\times {10}^{-6} \mathrm{s}$. The results from coarse mesh predict different flow physics than the other two cases, the boundary layer separates later, which results in higher shedding frequency as seen in supercritical cases. The differences in force and shedding predictions between medium and fine mesh are within $5\%$. Thus, the medium mesh is adopted in the current study for its performance and moderate computational cost.

4.1. Instantaneous Flow Field

The time history of cylinder lift and drag coefficients are shown in Figure 7, the mean drag coefficient $C_d=1.3$ and the lift coefficient $C_l$ averages at 0. Both force coefficients are plotted against the non-dimensionalized time $T=t{U}_0/D$, i.e., $T=700$ equals to $1\mathrm{s}$ simulation time. The first $0.2 \mathrm{s}$ history is neglected to exclude any transient effects. The spectral analysis of the oscillatory force on circular cylinder shows a dominant lift force frequency at $132 \mathrm{Hz}$, which corresponds to Strouhal number $St=0.193$. The wake shedding frequency from our PIV experiment is approximately $126 \mathrm{Hz}$ and $St=0.185$. Complied experimental data of $St$ for circular cylinder at Reynolds number around $\mathrm{64,000}$ is in the range of $0.186-0.196$ [30], the simulation data lies in the middle of the range.

The instantaneous pressure distribution on the circular cylinder shows the front stagnation line oscillates with eddy shedding, which agrees with hot-wire data from experiments [31]. The separation points on both sides oscillates with amplitude of approximately $\pm 4°$ with the mean separation angle of $\theta =85°$ from the front stagnation point. The instantaneous LM intermittency distribution on cylinder is plotted at four phases over one shedding cycle as shown in Figure 8. The peak near the front stagnation point is a numerical artifact which is caused by the diffusion of the freestream LM intermittency( ${\gamma }_{\infty }=1$). Between $\theta =\pm 85°$ the distribution of $\gamma$ remains low and steady over the entire shedding cycle and the boundary layer remains laminar. Transition initiates after separation. Various experiments were conducted in the range of subcritical Reynolds number, studies found the separation angle oscillates between $78°$ and $90°$ from the front stagnation point for circular cylinder of $R{e}_D=\mathrm{65,000}$ [32].

The stability of the shear layer is correlated with the intensity of eddy viscosity. Compared with the experimental data, simulation results predict less stochastic structures in cylinder wake than experiment and the shear layer is more coherent in the streamwise direction (Figure 9), which leads to longer a recirculation zone in the cylinder wake.

Figure 10 is an example of instantaneous iso-surface of Q-criterion colored by vorticity magnitude; it outlines the shedding vortices as they travel downstream and diffuse outwards at the same time. The spanwise band on the circular cylinder right before separation indicates a uniform separation position. Extensive experiments found significant effects of cylinder aspect ratio on the correlation length in low Reynold number regimes ($R{e}_D<{10}^4$). This dependency diminishes beyond $R{e}_D>\mathrm{20,000}$ and correlation length $L_c\sim 3D$ [33].

The blending of the hybrid model is fulfilled by the blending function $F_{DDES}$, which is directly multiplied to the destruction term in the TKE equation (Equation (4)). Compared with the RANS formulation, modeled turbulence is reduced in all regions of $F_{DDES}>1$.

Inside cylinder boundary layer, $F_{DDES}$ remains value of 1 (Figure 11), which preserves the boundary layer to RANS, it is essential for the hybrid model to work as desired since the boundary layer mesh cannot support resolved turbulence structures. Outside the cylinder wake, turbulence is close to freestream level, and the turbulence model operates in RANS mode, hence $F_{DDES}=1$.

4.2. Mean Flow Field and Flow Statistics

The comparison of mean velocity components $\overline{u}$ and $\overline{v}$ in cylinder wake up to $1D$ downstream is shown in Figure 12. The shape of PIV interrogation widow limits the acquired data in a rectangular region that is tangent to the cylinder surface, all data is normalized by the freestream velocity $U_{\infty }$. The numerical prediction of the length of recirculation zone and width of the wake as seen in Figure 12a are within $5\%$ difference compared with experimental data. A quantitative comparison between EFD and CFD results for the stream-wise and transverse velocities is presented in Figure 12e. The experimental results exhibit a slight asymmetry in the stream-wise velocity profile, which deteriorates the comparison; however, otherwise, the quantitative agreement of the results is satisfactory.

The streamlines of mean streamwise velocity from simulation in Figure 13 shows the recirculating flow extends to $0.725D$ downstream of the rear stagnation point, and the width of the wake, which determines the eddy shedding frequency, is approximately $1.1D$. The transverse velocity component $\overline{v}$ is antisymmetric about the center wake line (Figure 12b), as expected for transverse mean flow in a plane wake, and the main vortex formation region is centered at $1.3D$ downstream of the cylinder center. Another notable feature in Figure 12b is the small counter-rotating region near cylinder rear surface, which was also captured by experiments as separation bubbles.

Reynolds stress distributions from simulations and experiments are presented in Figure 14. The streamwise normal component $\overline{u^{\prime }{u}^{\prime }}$ peaks at approximately $1D$ downstream in both experiments and simulation results. The vertical normal stress $\overline{v^{\prime }{v}^{\prime }}$ reaches a maximum at $1.3D$ and $1.4D$ in the simulation and experiment, respectively, which indicates that the shed vortex detaches from the cylinder slightly further downstream in the experiments. The distribution of the shear stress component $\overline{u^{\prime }{v}^{\prime }}$ gives average length of eddy formation.

Similar to $\overline{v^{\prime }{v}^{\prime }}$, the streamwise location of its peak from simulation is located $0.1D$ upstream of experimental data. Generally, the contours indicate that the simulation predicts peak turbulence levels that are located slightly nearer the cylinder than the experiment indicates, with maximum and minimum values similar across the experiments and simulations.

The local flow fluctuations in cylinder wake determines the boundary layer dynamics on the downstream objects. Instantaneous velocity components are sampled at frequency $3\times {10}^6 \mathrm{Hz}$ and the flow angle is calculated at four locations downstream of the cylinder center: $2D$, $4D$, $8D$ and $10.7D$.

The diffusion and break down of turbulence structures lead to more evenly local flow angles distribution. At $10.7D$ downstream, where the airfoil sits in the surrogate model, the distribution peaks at $\pm 25°$ with extreme values up to $\pm 50°$ (Figure 15). Flow angle at airfoil position is also extracted from PIV data, compared with simulation data at $10.7D$, both peak at $\pm 25°$ with probability density value around 0.23. The experiment shows slightly shallower trough, which indicates more low amplitude fluctuations that are not captured in simulations.

5. Unsteady Boundary Layer on Airfoil

The cylinder flow analysis presented in the last section demonstrated the ability of the Langtry-Menter transition model implemented in a DDES formulation in predicting unsteady flow features. This section presents the unsteady flow dynamics on the downstream airfoil, the features of high blade loading fluctuations due to unsteadiness and transitional boundary layers are of interest in the fluid dynamics of full-scale wind turbine blades, making the surrogate problem a comprehensive test case for model development and validation. Along with the cylinder-airfoil setup case, the same airfoil in clean flow is also tested. The only difference is the removal of the upstream cylinder (9 million total cells).

5.1. Instantaneous Flow Field

Due to the periodic nature of the cylinder wake, the boundary layer response on the airfoil also shows periodic dynamics. The periods of the boundary layer dynamics is obtained from the lift and drag forces action on the airfoil, several other flow fields are analyzed at multiple phases within one cycle. With the presence of downstream airfoil, the cylinder shedding frequency increases from $132 \mathrm{Hz}$ to $137 \mathrm{Hz}$, the existence of downstream airfoil alters the vortex shedding by $3.6\%.$ Our wind tunnel experiments did not test the cylinder flow in the full-model setup.

The mean and RMS of $C_l$ and $C_d$ for both cases and XFOIL result are summarized in

Table 4. Mean and RMS force coefficient of airfoil.

	${\mathit{C}}_{\mathit{l}}$ Mean	${\mathit{C}}_{\mathit{l}}$ RMS	${\mathit{C}}_{\mathit{d}}$ Mean	${\mathit{C}}_{\mathit{d}}$ RMS
Airfoil in Clean flow	0.16	$0.16$	0.016	0.016
XFOIL	0.17	/	0.015	/
Airfoil in Cylinder wake	0.10	0.43	−0.026	0.042

. Airfoil in clean flow has much smaller fluctuations in the force coefficient, mean value matches the RMS value. In contrast, airfoil in cylinder flow experiences significant oscillations, the RMS values of both lift and drag are significantly higher than the corresponding mean value.

The spectral analysis shows the lift fluctuation coincides with the cylinder vortex shedding at $137\mathrm{Hz}$ and the drag force fluctuation is exactly twice the rate of the lift force. We conclude that the periodic vortex shed from cylinder is the main driver of the unsteady flow around airfoil.

The blending function $F_{DDES}$ is used to regulate TKE destruction as mentioned in last section. Inside the boundary layer, where the specific dissipation $\omega$ is high, $F_{DDES}$ remains value of 1 and the original RANS formulation is preserved (see an example distribution in Figure 16). Outside the boundary layer $F_{DDES}>1$ in the region of fine mesh (i.e., small $\mathsf{\Delta }$) and high TKE (inside eddies).

From Figure 16, the distribution of $F_{DDES}$ is concentrated around the cylinder and airfoil where the mesh is refined for transition prediction. As the grid size increases from cylinder surface to $\sim 1.5D$ downstream of the rear stagnation point, $F_{DDES}$ decreases to $\sim 1.4$ in the outer regions of the eddies, which indicates extra $40\%$ modeled TKE is destructed compared with RANS model. In the core of the eddies, where local TKE is high, $F_{DDES}$ remains $\sim 2$. As the vortex interacts with the airfoil, where mesh is refined again, the eddies are stretched around the leading edge and $F_{DDES}$ resumes to the same level as in the cylinder near wake.

Figure 17 is the instantaneous eddy viscosity distribution around airfoil in the cylinder wake. The alternating vortex structures exist on either side of the airfoil, notice there is a shear layer leaving the airfoil trailing edge which interacts with the cylinder vortices.

There is a strong correlation between $F_{DDES}$ and ${\nu }_t$ outside the boundary layer. $F_{DDES}$ reaches maximum in the center of the vortices and the eddy viscosity, which is proportional to TKE, is attenuated (e.g., see suction side near the leading edge). The rotational flow interacts with the airfoil boundary layer and cause local increase in TKE and eddy viscosity, the blending function $F_{DDES}$ follows the same trend. By the time the eddy reaches trailing edge, the attenuation effect of $F_{DDES}$ decreases the eddy viscosity significantly.

Inspecting the fluid force acting on the airfoil, we found the cylinder vortices are the main driver of the airfoil unsteadiness. The normalized vorticity field ${\tilde{\mathsf{\omega }}}_{\mathrm{z}}={\mathsf{\omega }}_{\mathrm{z}}/\left({\mathrm{U}}_{\infty }\mathrm{D}\right)$ is concentrated inside the Karman vortex and the width of the wake grows as the they travel downstream as shown in Figure 18, the coherent turbulence structures break up and stretched as they interact with the airfoil boundary layer (Figure 19).

The LM turbulent intermittency field gives a clear indication of laminar boundary layer separation, the airfoil boundary layer in clean flow remains laminar before separating and the transition location after separation is clearly visible from LM intermittency distribution on airfoil surface as shown in Figure 20a. Boundary layer transition initiates on both sides at approximately $x/c=0.6$, which is confirmed by the skin friction distribution (Figure 20b). The separation bubble extends $0.1c$ downstream of the transition point and the boundary layer flow reattaches at $x/c=0.7$ on both sides.

Airfoil boundary layer inside cylinder wake is fully turbulent and attached at all time. The surface LM intermittency remains low ($<0.1$ over $90\%$ of the airfoil surface), which indicates the existence of laminar viscous sublayer in the attached boundary layer. Experiments looked at boundary layer parameters from $x/c=0.3$ to $x/c=0.8$ at $0.1c$ interval and found it remained attached and turbulent at all locations, which confirms our observation from simulation.

To demonstrate the variations of the velocity magnitude inside boundary layer, the instantaneous velocity component that is parallel to the airfoil surface at $x/c=0.3$, $0.5$ and $0.8$ within one cycle are plotted in Figure 21. The instantaneous profile is averaged within 8 phase bins using phase averaging based upon the shedding phase of the cylinder vortex. No reverse flow is observed in the vicinity of the wall at the all positions, therefore there is no boundary layer separation at all time. Compare velocity profile at $x/c=0.3$ with $x/c=0.8$, it is noticeable that the boundary layer has thickened as flow travels downstream. There is phase lag across three locations, which will be inspected further in phase-averaged profiles.

5.2. Mean Flow Field and Flow Statistics

5.2.1. Mean Pressure and Velocity Field

Figure 22 is the predictions of mean pressure distribution from simulation results of airfoil in clean flow and cylinder wake. The pressure distribution in clean flow agrees with XFOIL data from leading edge to the reattachment locations. The CFD data predicts lower pressure gradient than XFOIL data in the turbulent boundary layer after reattachment.

Mean flow velocity around the airfoil in clean flow and cylinder wake is extracted from both experimental and simulation data (Figure 23). Note that the PIV data are only available on the suction side of the airfoil due to limitations of the optical arrangement. Flow reattaches on both sides of the airfoil in clean flow at $x/c=0.7,$ which is marked by the sudden change in mean velocity field (Figure 23a,b). Comparing the EFD and CFD results for the clean flow, the pressure gradient induced by the airfoil appears to have a slightly different far-field behavior; however, based upon the comparison of Figure 23e at $x/c$ = 0.5, small differences between the experimental and computation boundary conditions are likely to drive this difference: the discrepancy between results is nearly constant at all distances from the airfoil.

The presence and influence of the mean wake flow is clearly evident in the case of the airfoil in the cylinder wake (Figure 23c,d). A comparison between the EFD and CFD results at $x/c$ = 0.5 (Figure 23e) lead to some notable observations. In contrast to the clean flow case, even the boundary layer is highly pronounced and of a scale comparable to the large-scale mean flow gradients. Both the experiment and simulation show a similar thickness for the boundary layer this station. As is evident in the contours, but highlighted by the line plots, the mean flow due to the circular cylinder wake has some notable differences between the experiment and simulation. The gradient of the wake is greater for the CFD results, indicating that turbulent diffusion in the experiment was greater. The question of small differences in boundary conditions must be reconsidered in the assessment of these differences. Further, the small differences in the near-cylinder turbulence seen in Figure 14 will result in changes in the evolution of the turbulent wake development over the 10.67 diameters of spacing between the cylinder and airfoil. These results, while encouraging for the application of transition models in the DDES framework, highlight the challenges involved with providing experimental boundary conditions of sufficient fidelity for high confidence turbulence model assessment and improvement.

5.2.2. Reynolds Stress

Figure 24 is a comparison of three normalized Reynolds stress components around airfoil in cylinder wake. The experimental data is only available on the suction side and the PIV system resolution limited the measurement to the outer region of boundary layer.

The magnitude and distribution of all three Reynolds stress components from simulation agree with experimental data, however the experimental results contain more stochastic structures from the decay of cylinder wake eddies. The shear stress term $\overline{u^{\prime }{v}^{\prime }}$ peaks at leading edge due to high distortion of turbulence structures in cylinder wake, same feature is absent in the case of clean inflow.

The most noticeable feature of Reynolds stress field is the redistribution of fluctuation energy between the normal stress components. The streamwise component $\overline{u^{\prime }{u}^{\prime }}$ increases significantly from the leading edge while the transverse component $\overline{v^{\prime }{v}^{\prime }}$ reduces at the similar rate. The redistribution of the normal stress as approaching the solid surface is an invisicid effect described by rapid distortion theory (RDT, Figure 25). The Hunt and Graham (HG)’s results for an instantaneously appearing flat-plate [34] are rescaled with $L_{11}=0.030 \mathrm{m}$ and $L_{22}=0.060 \mathrm{m}$ to fit the measured data in experiments and simulations. The experimental and simulation results that are represented by dot and dash lines follow the same trends as the theoretical prediction with slightly higher gradient in $\overline{u^{\prime }{u}^{\prime }}$ close to the wall. According to Pope [35], the Reynolds stress component is proportional to the correspond strain rate and the pressure fluctuation $p^{\prime }$, as the strain rate component in transverse direction reduces on airfoil surface, the transverse normal stress $\overline{v^{\prime }{v}^{\prime }}$ is redistributed to the streamwise and spanwise components. The redistribution features are absent from the clean flow case since the pressure fluctuation $p^{\prime }$ must present to initiate the mechanism [17]. The redistribution of Reynolds stress also has significant effects on turbulence energy of the streamwise and transverse velocity components $u$ and $v$.

5.2.3. Mean Turbulence Kinetic Energy

The mean grid-resolved TKE around airfoil which sits in clean inflow and cylinder wake are shown in Figure 26, in both cases the TKE field is normalized by $U_{\infty }^2$. TKE starts to grow downstream of the laminar separation point on the airfoil in the clean flow case, the distribution is concentrated inside the separated boundary layer, which confirms that RANS model is preserved.

For the airfoil in cylinder wake, TKE is concentrated in the thin boundary layer at leading edge and it diffuses outwards as traveling downstream. The peak at the trailing edge is due to the free shear layer leaves the airfoil.

5.2.4. Phase-Averaged Boundary Layer Profiles

Probes are placed from airfoil surface up to $0.1c$ in the wall-normal direction. The phase of a band limited signal $s\left(\xi \right)$ is determined by Hilbert transform [36] $h\left(t\right)$:

$$h\left(t\right)=\frac{1}{\pi } {{\displaystyle \int }}_{-\infty }^{+\infty }\frac{s\left(\xi \right)}{\xi -t}d\xi$$

(5)

The phase is calculated as:

$$\varphi =co{s}^{-1} \left(\frac{real\left(h\left(t\right)\right)}{\left|h\left(t\right)\right|} \right)$$

(6)

The simulation data contains over 140 periods of the oscillatory flow around the airfoil, each cycle is subdivided into 16 phases with a width of $22.5°$ phase angles. The velocity profile is then averaged within each phase of all cycles. In order to plot the profile in terms of non-dimensional parameters $u^{+}$ and $y^{+}$, least square method is adopted to find the value of friction velocity $u_{\tau }$ that best fit the profile in viscous sublayer ($y^{+}<5$) and log layer ($30<y^{+}<100$). Figure 27 is a plot of the friction velocity distribution for each phase at every chordwise positions. There is a consistent phase-shift and decrease in $u_{\tau }$ as moving downstream, from the amount of phase-shift from station $x/c=0.3$ to $x/c=0.8$ together with the flow period, the wave speed is $22 \mathrm{m}/\mathrm{s}$.

Phase-averaged boundary layer profiles from simulation and experiments are plotted in Figure 28 together with the Spalding’s single formula for the “Law of the Wall” [37], viz.:

$$y^{+}=u^{+}+0.1108\left\{ e^{0.4u^{+}}-1-0.4u^{+}-\frac{{\left(0.4u^{+}\right)}^2}{2!}-\frac{{\left(0.4u^{+}\right)}^3}{3!} \right\}$$

(7)

which is valid through the viscous sublayer and the log layer.

The Log Law curve follows the standard “Law of the Wall” formulation:

$$u^{+}=\frac{1}{\kappa }\ln y^{+}+5$$

(8)

The experimental data is available from $y^{+}\sim 20$ due to the limited PIV system resolution near surface. Two formulations of “Law of the Wall” converge from approximately $y^{+}=30$ and the mean velocity profile is also plotted on the same figure.

There are several notable observations: although the least square method that fits the curves takes the data between $y^{+}<5$ and $30<y^{+}<100$ into account, every curve at each chordwise position converges at approximately $y^{+}=50$ and diverges above it, which indicates the cylinder wake dynamics dominate the wake region above $y^{+}>50$. The simulation data shows the magnitude of the oscillation in the wake region gradually grows in chordwise progression, while the experimental data predicted nearly constant variations $(\sim u^{+\prime }=5$) of the boundary layer profile at $y^{+}\sim 300$. The flow in the boundary layer wake region carries more energy in the simulation results.

Even in the highly-unsteady cylinder wake, the Law of the Wall is still valid in the low-momentum viscous inner regions, which confirms the validity of its adoption in many wall-modeled RANS and LES applications. Finally, the Law of the Wall agrees better with the phase and mean profiles of downstream locations, where the boundary layers are significantly thicker.

6. Conclusions and Future Work

The main objective of this study is to find an improved turbulence modeling approach that can capture the flow physics presented in the boundary layer of a wind turbine blade inside the atmospheric boundary layer and upstream turbine wake. A surrogate model was designed to keep the computational and experimental costs manageable while preserving the key flow physics present in the full-scale problem. We realized the LES-type turbulence model is beyond the computing power available currently even just for the surrogate model, while also noting promise from previous studies that show the value of using turbulence models that include boundary layer transition formulation into the RANS model for improved blade aerodynamics predictions.

The delayed detached eddy simulation (DDES) was chosen as the hybrid modeling methodology for its superior performance over RANS in massively-separated flow regions at affordable computational cost. The implementation of the transition model was verified against the simulation data by its original authors and validated using well-established experimental data. The integrated hybrid transition turbulence model was then applied to the circular cylinder at subcritical Reynolds number. Compared with the complementary experimental data, the mean flow dynamics, vortex shedding frequency and Reynolds stresses in the cylinder near wake were well captured by transition hybrid model.

Applying the transition hybrid model to the full-model led to some major discoveries in unsteady flow physics. The passage of the eddies dominates the pressure variations on the airfoil. Despite large local flow angle variations, there was no separation on the airfoil in the cylinder wake. Examining the Reynolds stress distribution, we found redistribution of the normal stress components, which is an inviscid effect predicted by rapid distortion theory. Inspecting the blending function $F_{DDES}$ distribution, the hybrid model operated in RANS mode inside the cylinder and airfoil boundary layer, where the local mesh resolution cannot afford resolving turbulence structure. The phase-averaged boundary layer profiles on airfoil followed the law-of-the-wall in the inner region.

The major downside of our current surrogate model is the scale mismatch compared with the full-scale problem as summarized in Table 1, and the authors have interest in addressing this through future work. The inflow unsteadiness is the key parameter in determining the blade dynamics, other than the blade root, and the reduced frequency of the full-scale blade is approximately one or two orders of magnitude lower than our model. By manipulating the definition of Strouhal number $St=fD/U$ and reduced frequency $k=\omega b/U$ we have $k=\pi St c/D$, $St$ maintains $\sim 0.2$ for ${10}^3<R{e}_D<{10}^6$, such that the reduced frequency $k\sim 0.63c/D$ is determined by the ratio of the chord length to the cylinder diameter—a scale-up model design will focus on increasing the cylinder size. Another benefit of the larger cylinder is increasing the turbulence length scale that can engulf the whole airfoil as the atmospheric turbulence does to the full-scale turbine blade. The flow field containing the larger diameter cylinder would contain key physics similar to the harmonic gust problem described by the Sears function.

In conclusion, the transition hybrid RANS/LES model showed its potential in calculating unsteady flow physics of a model wind turbine blade. Compared with RANS model, it resolves the turbulence structures down to the local grid size, which are averaged and smeared in unsteady RANS simulations; compared with LES models, the computational cost is affordable for most practical problems, and the meshing process is less stringent since the new model can, at worst, operate in RANS mode. The integration of the transition model enabled the accurate prediction of laminar separation, transition and flow reattachment, which were proven to be first-order factors in calculating blade aerodynamics.