Investigation of Turbulence Intensity Effects on Tidal Turbine Wakes through the BEM–CFD Method

1. Introduction

Tidal energy offers a clean and efficient alternative to fossil fuels due to its reliability, high energy density, predictability, and durability[1]. Over the past two decades, key terms in tidal energy research, such as "tidal turbine" and "wake"[2], highlight the importance of understanding turbine hydrodynamics and wake behavior for optimizing turbine array designs. The wake from upstream turbines influence the inflow conditions of downstream turbines, reducing water kinetic energy and increasing unsteady flow fluctuations[3]. Consequently, the hydrodynamic impact of turbine wakes critically affects the array layout and the overall power performance[4].

Generally, there are two primary methods for exploring the characteristics of turbine wakes: experimentation and numerical modeling. Simmons et al. used experimental methods, and demonstrate that higher flow speeds and turbulence at the turbine tips accelerate wake recovery[3]. Chen conducted wake experiment and found that the wake turbulence was strong and anisotropic[5].

Maganga, Fabrice, and several other researchers have conducted many numerical studies, where the turbine was modeled in its actual shape within the numerical model[6,7,8,9]. Due to the high computational demands of modeling turbines in their real shape, this approach is not suitable for large turbine arrays. As a result, simplified turbine models, such as the porous disk and Blade Element Momentum Theory (BEMT) models, have been proposed. These models approximate the turbine's impact on the flow without detailing the complex geometry of actual blades. This simplification preserves essential wake characteristics while significantly lowering computational costs, making it suitable for wake analysis and array optimization.

Chen and Myers used the porous disk model in laboratory experiments to study wake deficit and shear stress variations, their findings showed that wake velocity rapidly decreases due to momentum loss. In a two-row array, the far wake region (>4D downstream) exhibited a higher velocity deficit compared to a single disk[4,10,11].Batten developed a numerical actuator disk model to investigate the wake characteristics [12]. Shives coupled an actuator disk approach to predict turbine power performance, and the results showed the arrangement of turbine arrays is crucial for effective energy capture[13]. Churchfield modelled tidal turbine using rotating actuator lines to study wake and power performance of different tidal arrays[14].

Ambient turbulence intensity plays a significant role in turbine wake dynamics, influencing the wake's shape, length. Edmunds demonstrated that increased turbulence intensity in the rotor region, particularly near the wake boundary, improves momentum transfer downstream[15]. Gaurier and Mycek revealed that turbine wake recovers faster at higher turbulence [16,17], Nuernberg found that the turbulence intensity in the wake field first increases and then decreases [18]. While high ambient turbulence intensity appears to aid in wake recovery, Mycek’s study showed that the mean power performance of a single turbine is hardly influenced by varying turbulence intensity[17], with the main impact occurring on the inflow of downstream turbine arrays. High turbulence intensity may pose a challenge to the turbine's fatigue and safe. Vinod and Maganga both found that high turbulent intensity caused greater fluctuations in turbine thrust[19,20,21].

To achieve the industrial goals of tidal energy, optimizing the layout of turbine arrays in marine environments is a key objective. Similar to the Jensen model commonly used in wind energy, analytical models for tidal turbines have been developed to predict wake characteristics. These models can be effectively integrated with large-scale marine models such as FVCOM, and can be applied in turbine array optimization of marine sites.

Lam developed an analytical turbine wake model by modifying the Gaussian probability distribution equation, demonstrating strong predictive accuracy for wake velocity [22,23]. In addition to the wake model for horizontal axis turbines, Lam also developed an analytical wake model for vertical axis turbines[24]. Jordan studied different turbine arrays performance using the analytical model FLORIS in ocean model Thetis[25], Brutto developed a wake equation that can be applied in conjunction with Jensen's equation to predict wake velocity[26]. While the wake characteristics and analytical models have been studied by many researchers, the effect of turbulence intensity on wake shape and how to incorporate this influence into the analytical model still requires further investigation.

The aim of this work is to investigate the influence of turbulence intensity on wake dynamics and to develop an analytical wake model that incorporates the effects of turbulence intensity. Section II outlines the theoretical framework of the BEM-CFD numerical model. Section III validates the numerical model by comparing its predictions with experimental data. The numerical results demonstrate strong agreement with the experimental data. In Section IV, we investigate the influence of varying background turbulence intensities, from 2% to 18%, on the wake's expansion and range. Section V proposes an analytical wake model to predict the wake and compares it with CFD data. Finally, summarizes the key conclusions of this study.

2. CFD-BEM Coupled Model

The Blade Element Momentum (BEM) theory is a two-dimensional hydrofoil method widely used in the design of blades for wind and tidal turbines. Viterna and Myers applied BEM theory to the design of wind turbine and tidal turbine, respectively [27,28].

To improve computational efficiency, some researchers have successfully applied the BEM method to simplify tidal turbines in CFD numerical models. This approach simplifies the turbine by dividing the rotating plane into different annular sections, with each section calculated individually.

Guo and Masters both utilized a BEM-CFD numerical model to simplify the horizontal-axis tidal turbine and investigate the wake recovery and other associated characteristics[29,30]. In addition to horizontal-axis turbines, vertical-axis turbines can also be modeled using the BEM-CFD method [31]. Given the extensive research on BEM-CFD modeling of tidal turbine wakes[32,33,34], the theory has been thoroughly introduced in the literatures. This paper provides only a brief overview of the coupled BEM-CFD model.

The coupled BEM-CFD model represents the full turbine rotor geometry using the BEM method in the numerical domain, while the surrounding flow is governed by the Reynolds-averaged Navier–Stokes equations. The numerical simulations were conducted using ANSYS Fluent. The governing equations for momentum and mass conservation are presented below.

Mass conservation:

$${\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial (\rho U_i)}{\partial x_i}}=0$$

(1)

Momentum conservation:

$${\displaystyle \frac{\partial (\rho U_i)}{\partial t}}+{\displaystyle \frac{\partial (\rho U_i{U}_j)}{\partial x_j}}=-{\displaystyle \frac{\partial P}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}\left[\mu ({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}})\right]+{\displaystyle \frac{\partial }{\partial x_j}}(-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j})+\rho g_i+S_i$$

(2)

Here, $\rho$ denotes the density of water, t represents time, $U_i$ is the velocity of the water averaged over time, $x_i$ is the spatial distance, $P$ is the mean pressure, $-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}$ is the Reynolds stress, which must be resolved using a turbulence model, $\mu$ is the viscosity, $g_i$ is the component of the gravitational acceleration, and $S_i$ is a source term added to the $i=x$, y, or z momentum equations. The user-defined function of Fluent was used to define the added source term.

The resultant velocity of the blade element, $W=\sqrt{{\left(U_{\infty }(1-a)\right)}^2+{\left(\Omega r(1+b)\right)}^2}$, is based on BEM theory and a and b are the axial and tangential flow induction factors, respectively. When coupled with CFD, $U_{\infty }(1-a)$ in this equation can be replaced by the local velocity $U_u$ and the tangential flow induction factor b can be calculated from the local velocities $U_v$ and $U_w$, $b={\displaystyle \frac{\sqrt{U_v^2+U_w^2}}{\Omega r}}$. The axial and tangential momentum loss terms of the disk region, which represent the turbine in the CFD domain are $S_a$ and $S_t$, respectively.

$$S_a={\displaystyle \frac{B\rho c{W}^2(C_l\cos \varphi +C_d\sin \varphi )}{4\pi re}}{F}_t$$

(3)

$$S_t={\displaystyle \frac{B\rho c{W}^2(C_l\sin \varphi -C_d\cos \varphi )}{4\pi re}}{F}_t$$

(4)

Here, c is the length of the blade chord at radius r, r is the radius of the annular ring, e is the disk thickness, which equals 0.1–0.2 of the blade chord length[33], B is the number of blades, φ is the flow angle in the plane of rotation, C_l and C_d are the lift and drag coefficients of the hydrofoil, respectively, and F_t is the blade tip loss factor, which is expressed as:

$$F_t={\displaystyle \frac{2}{\pi }}\mathrm{arc} \cos \left[\exp {\displaystyle \frac{-B(R-r)}{2r\sin \varphi }}\right]$$

(5)

where R is the turbine diameter.

A flowchart of the coupled model is shown in Figure 1.

2.1. Validation of the Numerical Model

The validation experiment was conducted in the circulating flume at the National Ocean Technology Center, Tianjin, China. A two-blade horizontal axis turbine was positioned at the center of the flume as shown in Figure 2. The dimensions of the flume were 70 m in length, 1.6 m in width, and 2 m in height.

The water depth was maintained at 1.2 m, and the incoming flow was set to 0.489 m/s. The turbine had a diameter of 0.8 m, with a geometric scale of 1:15. The incident flow velocity corresponded to a real-life velocity of 1.9 m/s.

The time-varying velocity components were measured using a NORTEC acoustic doppler velocimeter, which captures velocity data with high temporal resolution. The sampling frequency was set to 200 Hz. Four curves, labeled 01, 02, 03, and 04, are shown in Figure 2, with curve 01 representing the centerline. The adjacent measuring points for each line were spaced 20 cm apart, and the four lines were separated by 10 cm. Further experimental details can be referenced in the study by[35].

The numerical model was validated against experimental results. The domain of the numerical model matched that of the experimental setup, with the incoming velocity set to 0.489 m/s. The computational domain dimensions were 25 m × 1.6 m × 1.2 m. The CFD parameters are detailed in Table 1.

The results from the coupled numerical and experimental models are presented in Figure 3. The dots represent the experimental data, while the line depicts the simulation data. The equation for the deficit coefficient is provided in Eq. (7).

$$Deficit=1-{\displaystyle \frac{\overline{U_x}}{\overline{U_0}}}$$

(7)

Here, $\overline{U_0}$ is the free-stream flow velocity and $\overline{U_x}$ is the wake velocity.

Table 2 presents the mesh independence results, showing the wake deficit at 8D downstream of centerline. Four cases were calculated in total, with δ representing the relative error, Δthe absolute error, and the experimental wake deficit at 8D being 0.20. The results indicate that as the number of mesh nodes increases, the wake deficit approaches the experimental results. The minimum deficit occurs in Case01, with 5.1 × 10⁶ mesh nodes, while the second closest result is Case02, with 6.3 × 10⁵ mesh nodes. For meshes coarser than case02, the error between simulation and experiment becomes relatively large.

Since the deficit between Case01 and Case02 is very small, the focus shifts to computational efficiency. The computational time for Case 02 is under 1 hour, whereas for Case 01, it is nearly 6 hours. Therefore, for both accuracy and efficiency, the mesh with 6.3 × 105 nodes (Case 02) was selected for the subsequent analysis.

The experimental data indicate that the wake deficit decreases as the axial distance increases. The deficit coefficient around the centerline is larger than that behind the turbine tip. Specifically, the deficit of the centerline is 0.4 at 4 X/D (axial distance/diameter of turbine), 0.29 at 6 X/D, and 0.09 at 12 X/D; For curve 04, the deficit is 0.21 at 4 X/D, 0.18 at 6 X/D, and 0.08 at 12 X/D. The error analysis between the simulation data and the experimental data is presented in Table 3, where $\delta$ represents the relative error and $\Delta$ represents the absolute error.

As the axial distance increases, both $\delta$ and $\Delta$ decrease. The discrepancy between the numerical and experimental results is slightly larger at shorter axial distances. This difference may be attributed to the absence of swirling in the coupled model, which contrasts with the rotating blades of a real turbine. Consequently, the near-wake region exhibits a greater variance. However, as the axial distance increases, the impact of swirling diminishes, leading to improved agreement between the numerical and experimental models. Given that this coupled numerical model is intended primarily for large array studies, this level of agreement is considered acceptable.

2.2. Analysis of the Turbulence Intensity

We calculated five cases with varying background turbulence intensity (TI) values—2%, 6%, 10%, 14%, and 18%—using the numerical model. The background TI was defined in the velocity inlet boundary conditions of ANSYS Fluent. To study the influence region of the wake, the domain of the numerical model was set to 20 m × 4 m × 4 m, with an incident velocity of 2 m/s. All other parameters remained consistent with those of the validated numerical model.

Figure 4 illustrates the velocity contours for X/D ranging from −4 to 20 and Z/D ranging from −4 to 4, for turbulence intensities (TI) of 2% and 18%. In both cases, the wake regions expand. The low-speed area for TI = 2% is larger than that for TI = 18%, indicating that the water flow in the wake with TI = 18% recovers faster than that with TI = 2%. The extent of the contour lines under the same velocity conditions initially expands and then contracts.At X/D = 20, the 2.0 m/s contour line for TI = 2% is still expanding, while the 2.0 m/s contour line for TI = 18% is already contracting. This phenomenon is discussed in detail in the following section.

Figure 5 illustrates the TI contours for boundary conditions set at (a) 2% and (b) 18%. In Figure 5a, TI increases notably as water passes through the tidal turbine, decreasing as axial distance increases.Figure 5b demonstrates that the TI across the entire domain gradually diminishes with increasing axial distance. As water flows through the tidal turbine, TI experiences a significant increase; however, due to the higher background TI, the elevated TI induced by the wake gradually blends with it. The area exhibiting increased TI due to the wake is relatively limited, less distinct compared to the wake effects observed with a background TI of 2%.

The axial results along the centerline for each case were extracted from the numerical model and are depicted in Figure 6. Figure 6a displays the velocity deficit curve, while Figure 6b shows the TI curve. The deficit coefficient decreases with increasing axial distance and turbulence intensity. The velocity deficit curves for TI = 2% and TI = 6% exhibit close proximity. At X/D=2, the deficit ranges from 0.4–0.6, whereas at X/D = 20, it ranges from 0.15–0.02. This indicates that the turbine wake recovers more swiftly with higher TI.

Figure 6b depicts the TI curve along with the axial distance, showing a notable increase as water passes through the turbine, followed by a decrease with axial distance. The values of TI across all cases are quite similar initially, but distinctions become clearer with greater axial distances. Specifically, TI is approximately 27% at X/D = 2, ranging between 13% and 16% at X/D = 8, and between 8%–12% at X/D = 20.

Figure 7 displays the radial velocity results, illustrating how the influence region varies with changes in axial distance and TI. As axial distance increases, the velocity deficit diminishes under consistent TI conditions. The shapes of the radial deficit curves remain similar, with the primary distinction being the peak intensity around Z/D = 0. This peak decreases in magnitude as axial distance increases.

To discuss this phenomenon in detail, Figure 8 illustrates the influence region for each case where the deficit exceeds 0.05, showing variations with different TI and axial distances.

For background TI of 2%, the influence region ranges from -1.18 to 1.18 (Z/D) at X/D=4, expands to -1.20 to 1.20 (Z/D) at X/D=6, further to -1.25 to 1.25 (Z/D) at X/D=10, and reaches -1.28 to 1.28 (Z/D) at X/D=20, continuing to increase with axial distance.

At 6% background TI, the influence region spans from -1.18 to 1.18 (Z/D) at X/D=4, widens to -1.24 to 1.24 (Z/D) at X/D=10 and X/D=14, and narrows to -1.19 to 1.19 (Z/D) at X/D=20, initially increasing and then decreasing beyond X/D=14.

With 10% background TI, the influence region expands from X/D=4 to X/D=6 and contracts for X/D>8. For 18% background TI, the influence region contracts for X/D>4.

These observations highlight how the extent of the influence region varies with both axial distance and background turbulence intensity.

The range of the turbine wake initially expands and then contracts with increasing axial distance. This contraction starts at different locations depending on the TI values. As TI increases, the onset of wake contraction moves closer to the turbine.

Figure 9 presents the radial results of turbulent intensity (TI) for various cases. In Figure 9a, which depicts TI at an axial distance of X/D=4, the TI values at Z/D=0 (the centerline) across all cases are approximately 20%. As the wake region approaches the centerline, the TI values increase, with varying degrees depending on the background TI values. The magnitude of TI increase due to the turbine is more pronounced when the background TI is low.

Figure 9b displays the radial curve of turbulent intensity (TI) at an axial distance of X/D=10. The TI values at Z/D=0 for different cases range between 12% and 15%. At higher background TI values such as 18% or 14%, TI increases only slightly. When the background TI values are 2%, 6%, or 10%, the TI values around the centerline are similar. The heightened TI induced by the turbine gradually blends with the background TI in these cases.

Figure 9c illustrates the radial curve of turbulent intensity (TI) at X/D=20. Significant changes are observed primarily when the background TI is 2% or 6%. When the background TI values are 18%, 14%, or 10%, the TI values remain nearly unchanged. This suggests that the increase in TI caused by the turbine diminishes, and the TI of the wake becomes well mixed with the background TI at this location. As the background TI increases, the rate at which the wake-field TI mixes with the background TI accelerates.

2.2. Analytical Wake Model

Building upon the previous research, we studied the wake distribution formula for tidal stream turbines. Based on Lam’s formula, which is derived from a Gaussian probability distribution[22], we made improvements. The original study only considered specific turbulence intensities of 3%, 5%, 8%, and 15%, and developed separate wake distribution formulas for each turbulence intensity. In this study, we extended Lam’s approach by incorporating the turbulence intensity parameter into the formula. Compared to the original formula, the new formula includes three additional coefficients, namely a, b, and c, allowing the formula to reflect the influence of turbulence intensity in marine environments. The formula is expressed as follows:

$$U_{x,r}=U_{\infty }-\mathrm{c}(U_{\infty }-U_{\min })e^A$$

(8)

$$A=-a{(r/(b(x-R)))}^2$$

(9)

Empirical equations incorporating turbulent intensity are proposed to predict the recovery of the minimum velocity as follows:

$${\displaystyle \frac{V_{\min }}{V_0}}=0.9903+0.0997\left({\displaystyle \frac{x}{D}}\right)+2.5917T_i-0.0031{\left({\displaystyle \frac{x}{D}}\right)}^2-0.0436T_i\left({\displaystyle \frac{x}{D}}\right)$$

(10)

Considering that the axial distance between consecutive rows of tidal turbine arrays in the ocean is unlikely to be less than 4D, the above formula is applicable for axial distances greater than or equal to 4D.

The efflux velocity is:

$$V_0=V_{\infty }\sqrt{1-C_T}$$

(11)

The additional coefficient, incorporating turbulent intensity and axial distance, is:

$$a=-0.0227+0.0149x-0.5523T_i+0.0020x^2+0.3233x{T}_i$$

(12)

$$b=0.2523T_i+0.0565$$

(13)

$$c=1.0730-0.0088x-0.0478T_i$$

(14)

In this study, the turbine’s C_T=0.8 and V_∞=2m/s，The comparison between the above formula and numerical simulations is shown in Figure 10.

Table 4 presents the comparison of errors between the analytical wake model and CFD data. Data within the radial distance range of -2R to 2R, with a 2D interval, and the axial distance range of 4D to 12D are selected. It can be observed that the R² value ranges from 0.7861 to 0.9944, and the RMSE ranges from 0.0022 to 0.0355. The proposed predictive formula in this study effectively predicts the variation in wake velocity field with different turbulence intensities.

3. Conclusions

This study investigates the influence of ambient turbulence intensity on the wake dynamics of tidal turbines using a coupled Blade Element Momentum and Computational Fluid Dynamics model. The primary objective was to understand how different levels of background TI affect wake shape, recovery, and interaction with the surrounding flow. The numerical model, validated against experimental data, demonstrates good agreement, thereby confirming its reliability.

1. Wake Influence Region and Turbulence Intensity: The wake influence region initially expands and then contracts as the background TI increases. The point at which contraction begins moves closer to the turbine with higher TI values. As background TI increases, the turbulent intensity within the wake merges more quickly with the surrounding flow, leading to a less pronounced increase in TI within the wake region and indicating faster homogenization with the background TI.

2. Wake Recovery and Velocity Deficit: Higher background TI significantly accelerates the wake recovery process, reducing the extent of low-speed areas. The velocity deficit decreases with increasing axial distance from the turbine, with the deficit around the centerline being higher compared to areas behind the turbine tip. This highlights the non-uniformity of wake recovery and emphasizes the role of background TI in wake dynamics.

3. Implications for Turbine Array Layouts: This study highlights the critical role of considering ambient turbulence intensity (TI) in the design and optimization of tidal turbine arrays. The proposed formula, which reflects the influence of varying turbulence intensities, can be conveniently applied to the study of wake interactions in turbine arrays. Higher TI affects wake dynamics, potentially influencing the overall efficiency and performance of the arrays. However, the impact of high turbulence intensity on the structural integrity of the turbines was not addressed in this study and warrants further investigation.

In summary, besides demonstrating that higher turbulence intensity enhances wake recovery, our study also reveals notable variations in the wake field structure under different turbulence intensities. As turbulence intensity increases, the expansion trend and the affected region of the wake become more pronounced. Beyond the widely studied wake velocity variations, we also investigated the overall turbulence intensity changes within the wake field. Furthermore, we propose a practical analytical wake model that can rapidly predict the wake velocity, providing a useful tool for future studies and applications. This paper provides new insights and tools for optimizing turbine array configurations and elucidating wake evolution mechanisms.