Numerical improvement to Glauert correction for the flow around a wind turbine

Abstract

Accurate and reliable computation of the aerodynamic characteristics of wind turbines is very important for the development of new efficient designs. The flow around a wind turbine is modeled by a permeable disc (PD), solved through the Unsteady Reynolds-Averaged Navier–Stokes equations (URANS), here named PD/URANS method. The finite volume method and a total variation diminishing (TVD) scheme solve numerically the flow governing equations. The turbulent flow in the wake of the wind turbine is simulated utilising a one-equation turbulence model. The Glauert correction calculation considers a uniform normal force distribution (C_T) on the virtual permeable disc applied to the flow, while the axial induction factor is obtained directly from the numerical solution of the URANS equations. The numerical axial induction factor obtained agrees fairly well with Glauert correction, except if the flow behind the turbine is highly unsteady and Reynolds number dependent.

Nomenclature

a

induction factor

$a_\infty$

free-stream isothermal speed of sound

$C_T$

thrust coefficient

D

turbine diameter

P

pressure

$M_\infty$

free-stream Mach number

$R_D$

Reynolds number

$S_z$

z-area vector component

$S_r$

r-area vector component

t

time

T

temperature

$v_r$

r-velocity component

$v_z$

z-velocity component

$v_\theta$

$\theta$ -velocity component

V

volume

$V_0$

wind speed

$\rho$

density

$\tau$

isothermal compressibility

$\upsilon_t$

turbulent kinematics viscosity

$\infty$

free-stream conditions

1.0 Introduction

Efficient and accurate computation of wind turbine aerodynamic characteristics is very important for the development of new efficient designs. Usually, the aerodynamic characteristics of a wind turbine are obtained either by experimental model measurements, direct numerical solution of the Navier–Stokes equations or iterative solution of the blade element and momentum method (BEM) equations. Experimental measurements are expensive and timely consuming, requiring special devices and measurement equipment, wind tunnel and small tolerance construction of the wind turbine model. On the other hand, the BEM method is based on a quite simplified theory that neglects viscous effects on the flow behind the wind turbine. Glauert [5] proposed a correction factor to overcome such a deficiency of the BEM method that has become popular in the aeronautical applications. In turn, the three-dimensional numerical simulation of the flow around a wind turbine is highly computationally expensive since it requires a quite well-refined and a complex rotational grid to be generated to discretise the three-dimensional domain around the rotor blades.

The BEM method does not take explicitly into account the viscous effect that induces circumferential/tangential velocities in the turbine wake. Such an influence tends to be smaller on longer blades operating at high tip speed ratios, as is the case of wind turbine blades. Thus, BEM tends to be a suitable method for wind turbines rather than for marine propellers. Glauert has evaluated his axial induction factor from experimental results. Therefore, Glauert correction incorporated the existing implicit viscous effect relating tangential and axial inducted velocities. Numerical results from Navier–Stokes equations are able to capture such viscous influence on the axial induction factor, confirming, improving and generalising Glauert correction.

Therefore, an alternative approach as the Permeable Disc/Navier–Stokes (PD/NS) method, utilised in Ref. (1), can be a very accurate and computationally efficient way to solve the flow around a wind turbine. The PD/NS method may greatly reduce the number of required grid cells since the rotor geometry is not resolved. A virtual disc on which the forces from the BE theory are distributed replaces the rotor, while the induction factors required from the blade element (BE) theory are obtained directly through the flow numerical solution. The BE theory equations and the flow field governing equations in cylindrical coordinates are solved simultaneously. The Navier–Stokes equations may be described based on two cylindrical coordinates only because the flow is assumed to be axisymmetric. Therefore, a much simpler two-dimensional computational grid will be sufficient to produce the numerical solution of the flow field governing equations, resulting in a computationally efficient but still accurate method.

Mikkelsen et al. [6] studied numerically the influence of a wind turbine rotor coning using the BEM method and the permeable disc model combined with PD/NS. The incompressible Navier–Stokes equations in cylindrical coordinates were solved for axisymmetric flow using the finite difference method for Reynolds number = 5000. Comparison between results obtained with both models has shown the limitation of the BEM method. The computation also demonstrated that the power coefficient based on the projection area of the actuator disc is invariant to coning.

Zhong et al. (2019) combined the effect of tip loss corrections into the PD/NS method. The study was conducted using an axisymmetric PD/NS solver to simulate the flow around the wind turbine. The numerical solution of the Reynolds Averaged Navier–Stokes equations in cylindrical coordinates was obtained through the finite volume method with second order upwind scheme and [7] turbulence model. The widely used [5] tip loss function F, used in the BEM, and a newly developed tip loss correction were discussed and evaluated. According to their conclusions, the new tip loss correction has shown superior performance in different flow conditions compared to the original [5] tip loss function.

Sharp [8] used a general momentum theory for an energy-extracting permeable disc, based on the PD/NS method, modeling the rotor with a multiple blades defined by a uniform circulation that includes the effects of wake rotation and wake expansion. According to his conclusions, the rotation of the wake has been accompanied by an additional fall in the static pressure if one compares it with across the permeable disc. The overall effect of wake rotation on energy extraction is very small for wind turbine operating at high tip speed ratios, but can significantly increase the predicted power output otherwise.

Moens [9] studied the performance of tip-loss correction factor for advanced permeable disc (PD) method coupled with Large Eddy Simulation (LES) Navier–Stokes equations (PD/NS) for a wind farm configuration. The LES incompressible Navier–Stokes equations were solved using fourth-order finite difference scheme in a three-dimensional Cartesian grid. The classical [5] tip-loss factor was added to correct the tip and the root induced velocities at the rotor. The numerical results showed that the PD/NS performance clearly improves by the additional correction on the [5] tip-loss factor.

Simisiroglou [10] discussed a comparative analysis featuring two different permeable disc methods (PD/NS) [11, 12], and two alternative analytical wake models presented by Refs (13, 14), applied on the wind farm power production assessment. Simisiroglou [10] concluded that the PD/NS used by Refs (12, 14) outperforms all the other methods. In the two numerical approaches used in all analysis, the RANS equations and the k-ε turbulence model were solved using the commercial CFD code WindSim. The difference between [11, 12] results come from the equation used to calculate the thrust forces. In the Ref. (12) method, the thrust coefficient is updated at every time of iteration, whereas in the Ref. (11) method, it is kept constant during all the simulation.

The present work features the slightly compressible URANS equations solution in cylindrical coordinates using finite volume method and upwind total variation diminishing (TVD) scheme of Refs (2,3) using the [15] flux limiter. The [16] hypothesis and the [4] turbulence model simulate the turbulent flow in the wake generated by the disc. The governing equations depend only on two coordinates since the flow is assumed to be axisymmetric. Therefore, a two-dimensional computational grid is enough to solve the flow around the permeable disc. A uniform normal force distribution (C _T ) on the virtual permeable disc is applied to the flow field, and the axial induction factor at the permeable disc is obtained directly from the Unsteady Reynolds Averaged Navier-Stokes (URANS) equations numerical solution. Finally, the numerical results are compared with the traditional momentum theory results based on the [5] correction and [17] correction to evaluate accuracy.

2.0 Mathmatical formulation

The Reynolds-Averaged Navier–Stokes (Reynolds Averaged Navier-Stokes) equations in cylindrical coordinates are obtained by decomposing the Navier–Stokes equations dependent variables into time averaged and fluctuating components. The time averaged governing equations obtained by using the mass-weighted procedure suggested by Ref. (18) are shown below:

(1) \begin{equation}\frac{{\partial \bar \rho }}{{\partial t}} + \frac{1}{r}\frac{\partial }{{\partial r}}\big( {\bar \rho r{{\tilde v}_r}} \big) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\big( {\bar \rho {{\tilde v}_\theta }} \big) + \frac{\partial }{{\partial z}}\big( {\bar \rho {{\tilde v}_z}} \big) = 0\end{equation}

(2) \begin{align}\frac{{\partial \bar \rho {{\tilde v}_r}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \bar \rho r\tilde v_r^2}}{{\partial r}} + \frac{1}{r}\frac{{\partial \bar \rho {{\tilde v}_r}{{\tilde v}_\theta }}}{{\partial \theta }} + \frac{{\partial \bar \rho {{\tilde v}_r}{{\tilde v}_z}}}{{\partial z}} - \frac{{\bar \rho \tilde v_\theta ^2}}{r} + \frac{{\partial \bar p}}{{\partial r}} &= \frac{1}{r}\frac{\partial }{{\partial r}}\Big( {r{{\bar \tau }_{rr}} - r{{\overline {\rho {{v''}_r}v''} }_r}} \Big) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\Big( {{{\bar \tau }_{\theta r}} - {{\overline {\rho {{v''}_r}v''} }_\theta }} \Big) \nonumber \\ & \quad+ \frac{\partial }{{\partial z}}\Big( {{{\bar \tau }_{zr}} - {{\overline {\rho {{v''}_r}v''} }_z}} \Big) - \frac{{{{\bar \tau }_{\theta \theta }} - {{\overline {\rho {{v''}_\theta }v''} }_\theta }}}{r}\end{align}

(3) \begin{align}&\frac{{\partial \bar \rho {{\tilde v}_\theta }}}{{\partial t}} + \frac{1}{r}\frac{{\partial \bar \rho r{{\tilde v}_r}{{\tilde v}_\theta }}}{{\partial r}} + \frac{1}{r}\frac{{\partial \bar \rho \tilde v_\theta ^2}}{{\partial \theta }} + \frac{{\partial \bar \rho {{\tilde v}_\theta }{{\tilde v}_z}}}{{\partial z}} + \frac{{\bar \rho {{\tilde v}_r}{{\tilde v}_\theta }}}{r} + \frac{1}{r}\frac{{\partial \bar p}}{{\partial \theta }} \nonumber \\ & \quad= \frac{1}{r}\frac{\partial }{{\partial r}}\Big( {r{{\bar \tau }_{r\theta }} - r{{\overline {\rho {{v''}_r}v''} }_\theta }} \Big) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\Big( {{{\bar \tau }_{\theta \theta }} - {{\overline {\rho {{v''}_\theta }v''} }_\theta }} \Big) + \frac{\partial }{{\partial z}}\Big( {{{\bar \tau }_{z\theta }} - {{\overline {\rho {{v''}_\theta }v''} }_z}} \Big) + \frac{{{{\bar \tau }_{\theta r}} - {{\overline {\rho {{v''}_\theta }v''} }_r}}}{r}\end{align}

(4) \begin{align}\frac{{\partial \bar \rho {{\tilde v}_z}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \bar \rho r{{\tilde v}_r}{{\tilde v}_z}}}{{\partial r}} + \frac{1}{r}\frac{{\partial \bar \rho {{\tilde v}_z}{{\tilde v}_\theta }}}{{\partial \theta }} + \frac{{\partial \bar \rho \tilde v_z^2}}{{\partial z}} + \frac{{\partial \bar p}}{{\partial z}} &= \frac{1}{r}\frac{\partial }{{\partial r}}\Big( {r{{\bar \tau }_{zr}} - r{{\overline {\rho {{v''}_z}v''} }_r}} \Big) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\Big( {{{\bar \tau }_{\theta z}} - {{\overline {\rho {{v''}_z}v''} }_\theta }} \Big) \nonumber \\ & \quad + \frac{\partial }{{\partial z}}\Big( {{{\bar \tau }_{zz}} - {{\overline {\rho {{v''}_z}v''} }_z}} \Big)\end{align}

where the viscous stresses are:

\begin{equation*}\begin{array}{l}{\tau _{rr}} = \mu \!\left[ {2\dfrac{{\partial {v_r}}}{{\partial r}} - \dfrac{2}{3}\left( {\vec \nabla \cdot \vec v} \right)} \right]\\[1em]{\tau _{\theta \theta }} = \mu \!\left[ {2\!\left( {\dfrac{1}{r}\dfrac{{\partial {v_\theta }}}{{\partial \theta }} + \dfrac{{{v_r}}}{r}} \right) - \dfrac{2}{3}\left( {\vec \nabla \cdot \vec v} \right)} \right]\\[1em]{\tau _{zz}} = \mu \!\left[ {2\dfrac{{\partial {v_z}}}{{\partial z}} - \dfrac{2}{3}\left( {\vec \nabla \cdot \vec v} \right)} \right]\end{array}\end{equation*}

(5) \begin{align}{\tau _{r\theta }} &= \mu \!\left[ {r\frac{\partial }{{\partial r}}\left( {\frac{{{v_\theta }}}{r}} \right) + \frac{1}{r}\frac{{\partial {v_r}}}{{\partial \theta }}} \right]\\{\tau _{\theta z}} &= \mu \!\left( {\frac{{\partial {v_\theta }}}{{\partial z}} + \frac{1}{r}\frac{{\partial {v_z}}}{{\partial \theta }}} \right)\nonumber \\{\tau _{zr}} &= \mu \!\left( {\frac{{\partial {v_z}}}{{\partial r}} + \frac{{\partial {v_r}}}{{\partial z}}} \right)\nonumber\end{align}

and

(6) \begin{equation}\tilde f = \frac{{\overline {\rho f} }}{{\bar \rho }}\end{equation}

From basic thermodynamics, the state of a homogeneous substance can be defined by only two properties. The mathematical closure of the URANS set of equations requires the extra physical equation:

(7) \begin{align}\rho = \rho \!\left( {p,T} \right)\end{align}

Taylor expansion with respect to the free-stream conditions yields:

(8) \begin{align}\rho &= {\rho _\infty } + \left( {p - {p_\infty }} \right){\left. {\frac{{\partial \rho }}{{\partial p}}} \right)_\infty } + \left( {T - {T_\infty }} \right){\left. {\frac{{\partial \rho }}{{\partial T}}} \right)_\infty } + \frac{{{{\left( {p - {p_\infty }} \right)}^2}}}{{2!}}{\left. {\frac{{{\partial ^2}\rho }}{{\partial {p^2}}}} \right)_\infty } \nonumber \\ & + \left( {p - {p_\infty }} \right)\left( {T - {T_\infty }} \right){\left. {\frac{{{\partial ^2}\rho }}{{\partial p\partial T}}} \right)_\infty } + \frac{{{{\left( {T - {T_\infty }} \right)}^2}}}{{2!}}{\left. {\frac{{{\partial ^2}\rho }}{{\partial {T^2}}}} \right)_\infty } + \cdots \end{align}

For homogeneous temperature field, T=T _∞ , Taylor expansion simplifies to:

(9) \begin{equation}\rho = {\rho _\infty } + \left( {p - {p_\infty }} \right){\left. {\frac{{\partial \rho }}{{\partial p}}} \right)_\infty } + \frac{{{{\left( {p - {p_\infty }} \right)}^2}}}{{2!}}{\left. {\frac{{{\partial ^2}\rho }}{{\partial {p^2}}}} \right)_\infty } + \cdots \end{equation}

From the definition of isothermal compressibility, see Ref. (19):

(10) \begin{equation}\tau = \frac{1}{\rho }{\left. {\frac{{\partial \rho }}{{\partial p}}} \right)_T}\end{equation}

By assuming constant isothermal compressibility, it follows:

(11) \begin{align} &{\left. {\frac{{\partial \rho }}{{\partial p}}} \right)_T} = \tau \rho \nonumber\\ &{\left. {\frac{{{\partial ^2}\rho }}{{\partial {p^2}}}} \right)_T} = \tau {\left. {\frac{{\partial \rho }}{{\partial p}}} \right)_T} = {\tau ^2}\rho \nonumber \\ &{\left. {\frac{{{\partial ^3}\rho }}{{\partial {p^3}}}} \right)_T} = {\tau ^2}{\left. {\frac{{\partial \rho }}{{\partial p}}} \right)_T} = {\tau ^3}\rho \\ &\vdots \nonumber \end{align}

Substituting Equation (11) into Taylor expansion, yields:

(12) \begin{equation}\rho = {\rho _\infty } + {\rho _\infty }\left( {p - {p_\infty }} \right)\tau + {\rho _\infty }\frac{{{{\left( {p - {p_\infty }} \right)}^2}}}{{2!}}{\tau ^2} + \cdots \end{equation}

Since the isothermal compressibility can be thought as a very small physical parameter (10^–5 m²/N for air, and 10^–10 m²/N for water, at p = 1 atm), then the second and higher order terms can be neglected from Equation (12), yielding:

(13) \begin{equation}\rho \approx {\rho _\infty } + {\rho _\infty }\left( {p - {p_\infty }} \right)\tau \end{equation}

Replacing Equation (13) into the continuity equation, shown in Equation (1), it follows:

(14) \begin{align}&\tau {\rho _\infty }\left[ {\frac{{\partial \bar p}}{{\partial t}} + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {\bar pr{{\tilde v}_r}} \right) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\left( {\bar p{{\tilde v}_\theta }} \right) + \frac{\partial }{{\partial z}}\left( {\bar p{{\tilde v}_z}} \right)} \right] + {\rho _\infty }\left( {1 - {p_\infty }\tau } \right) \nonumber \\ & \quad \left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{{\tilde v}_r}} \right) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\left( {{{\tilde v}_\theta }} \right) + \frac{\partial }{{\partial z}}\left( {{{\tilde v}_z}} \right)} \right] = 0\end{align}

Equation (14) is the continuity equation now expressed in terms of pressure, which can be satisfied if:

(15) \begin{align}&\frac{{\partial \bar p}}{{\partial t}} + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {\bar pr{{\tilde v}_r}} \right) + \frac{1}{r}\frac{\partial }{{\partial \theta }}\left( {\bar p{{\tilde v}_\theta }} \right) = 0\nonumber\\&{\rm{and}}\nonumber \\&{p_\infty } = \frac{1}{\tau } \end{align}

In non-dimensional form, Equation (15) may be written as:

(16) \begin{align}&\frac{{\partial {{\bar p}^*}}}{{\partial {t^*}}} + \frac{1}{{{r^*}}}\frac{\partial }{{\partial {r^*}}}\left( {{{\bar p}^*}{r^*}\tilde v_r^*} \right) + \frac{1}{{{r^*}}}\frac{\partial }{{\partial \theta }}\left( {{{\bar p}^*}\tilde v_\theta ^*} \right) = 0 \\&{\hspace{-110pt}\rm{and}}\nonumber \\ &p_\infty ^* = 1\nonumber \end{align}

where

(17) \begin{equation}v_r^* = \frac{{{v_r}}}{{{a_\infty }}},\quad v_\theta ^* = \frac{{{v_\theta }}}{{{a_\infty }}},\quad {r^*} = \frac{r}{D},\quad {p^*} = \frac{p}{{{\rho _\infty }a_\infty ^2}},\quad {t^*} = \frac{{t{a_\infty }}}{D},\quad {a_\infty } = \sqrt {{{\left( {\frac{{\partial p}}{{\partial \rho }}} \right)}_T}} = \sqrt {\frac{1}{{{\rho _\infty }\tau }}} \end{equation}

Hereafter, the star symbol of the dimensionless variables in Equation (16) will be dropped for sake of simplicity. Equation (16) together with the incompressible momentum equations gives a convenient system of equations for slightly compressible flows. The system of equations below stands for axisymmetric flows in the non-dimensional form, including the [16] hypothesis in Equations (2), (3), and (4), to substitute the Reynolds stresses terms.

(18) \begin{equation}\frac{{\partial \bar p}}{{\partial t}} + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {\bar pr{{\tilde v}_r}} \right) + \frac{\partial }{{\partial z}}\left( {\bar p{{\tilde v}_z}} \right) = 0\end{equation}

(19) \begin{align}\frac{{\partial \bar \rho {{\tilde v}_r}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \bar \rho r\tilde v_r^2}}{{\partial r}} + \frac{{\partial \bar \rho {{\tilde v}_r}{{\tilde v}_z}}}{{\partial z}} - \frac{{\bar \rho \tilde v_\theta ^2}}{r} + \frac{{\partial \bar p}}{{\partial r}} & = \frac{1}{r}\frac{\partial }{{\partial r}}\left[ {r\!\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\left( {\frac{{\partial {{\tilde v}_z}}}{{\partial r}} + \frac{{\partial {{\tilde v}_r}}}{{\partial z}}} \right)} \right] \nonumber \\ & \quad+ \frac{\partial }{{\partial z}}\left( {2\!\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{{\partial {{\tilde v}_z}}}{{\partial z}}} \right) - \frac{2}{{{r^2}}}\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right){\tilde v_r}\end{align}

(20) \begin{align}\frac{{\partial \bar \rho {{\tilde v}_\theta }}}{{\partial t}} + \frac{1}{r}\frac{{\partial \bar \rho r{{\tilde v}_r}{{\tilde v}_\theta }}}{{\partial r}} + \frac{{\partial \bar \rho {{\tilde v}_\theta }{{\tilde v}_z}}}{{\partial z}} + \frac{{\bar \rho {{\tilde v}_r}{{\tilde v}_\theta }}}{r} &= \frac{1}{r}\frac{\partial }{{\partial r}}\left[ {{r^2}\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{\partial }{{\partial r}}\left( {\frac{{{{\tilde v}_\theta }}}{r}} \right)} \right] \nonumber \\ & \quad + \frac{\partial }{{\partial z}}\left( {\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{{\partial {{\tilde v}_\theta }}}{{\partial z}}} \right) + \left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{\partial }{{\partial r}}\left( {\frac{{{{\tilde v}_\theta }}}{r}} \right)\end{align}

(21) \begin{equation}\frac{{\partial \bar \rho {{\tilde v}_z}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \bar \rho r{{\tilde v}_r}{{\tilde v}_z}}}{{\partial r}} + \frac{{\partial \bar \rho \tilde v_z^2}}{{\partial z}} + \frac{{\partial \bar p}}{{\partial z}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left[ {r\!\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\left( {\frac{{\partial {{\tilde v}_z}}}{{\partial r}} + \frac{{\partial {{\tilde v}_r}}}{{\partial z}}} \right)} \right] + \frac{\partial }{{\partial z}}\left( {2\!\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{{\partial {{\tilde v}_z}}}{{\partial z}}} \right)\end{equation}

where υ _t is the turbulent kinematic viscosity obtained by utilising [4] turbulence model, shown below in cylindrical coordinates – see Ref. (4) for further details:

(22) \begin{align}\frac{{\partial {\upsilon _t}}}{{\partial t}} + \frac{1}{r}\frac{{\partial r{\upsilon _t}{v_r}}}{{\partial r}} + \frac{{\partial {\upsilon _t}{v_z}}}{{\partial z}} &= \frac{1}{r}\frac{\partial }{{\partial r}}\left[ {\frac{r}{\sigma }\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{{\partial {\upsilon _t}}}{{\partial r}}} \right] + \frac{\partial }{{\partial z}}\left[ {\frac{1}{\sigma }\left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)\frac{{\partial {\upsilon _t}}}{{\partial z}}} \right] \nonumber \\ & \quad + {c_{b1}}\left( {1 - {f_{t2}}} \right)\hat S{\upsilon _t} - \left[ {{c_{w1}}{f_w} - \frac{{{c_{b1}}}}{{{k^2}}}{f_{t2}}} \right]{\left( {\frac{{{\upsilon _t}}}{d}} \right)^2} + \frac{{{c_{b2}}}}{\sigma }\left[ {{{\left( {\frac{{\partial {\upsilon _t}}}{{\partial r}}} \right)}^2} + {{\left( {\frac{{\partial {\upsilon _t}}}{{\partial z}}} \right)}^2}} \right]\end{align}

3.0 Numerical formulation

The governing equations for three-dimensional and axisymmetric flows, written in the integral, vector and conservative forms, are solved using the finite volume method and the upwind TVD scheme, see Refs (2, 3), using the [15] flux limiter. Equation (23) shows the approximated URANS equations using finite differences and finite volumes for the time and space discretisation, respectively. Figure 1 illustrates the finite volume (in red) utilised in the spatial discretisation of the governing equations.

(23) \begin{equation}Q_{i,j}^{n + 1} = Q_{i,j}^n - \frac{{\Delta t}}{{{V_{i,j}}}}\left[ {{{\left( {\vec P \cdot \vec S} \right)}_{i + 1/2,j}} + {{\left( {\vec P \cdot \vec S} \right)}_{i - 1/2,j}} + {{\left( {\vec P \cdot \vec S} \right)}_{i,j + 1/2}} + {{\left( {\vec P \cdot \vec S} \right)}_{i,j - 1/2}}} \right] + \Delta t{H_{i,j}}\end{equation}

Figure 1. Finite volume (in red) utilised in the spatial discretisation of the governing equations.

where

(24) \begin{equation}\vec P \cdot \vec S = \left( {{E_e} - {E_v}} \right){S_z} + \left( {{F_e} - {F_v}} \right){S_r}\end{equation}

(25) \begin{equation}Q = r\!\left\{ {\begin{array}{*{20}{c}}p\\{{v_r}}\\{{v_z}}\\c\\{{\upsilon _t}}\end{array}} \right\}\quad \quad {E_e} = r\!\left\{ {\begin{array}{*{20}{c}}{p{v_r}}\\[1pt]{v_r^2 + \dfrac{p}{\rho }}\\[1pt]{{v_z}{v_r}}\\{c{v_r}}\\{{\upsilon _t}{v_r}}\end{array}} \right\}\quad \quad {F_e} = r\!\left\{ {\begin{array}{*{20}{c}}{p{v_z}}\\{{v_z}{v_r}}\\[1pt]{v_z^2 + \dfrac{p}{\rho }}\\{c{v_z}}\\[1pt]{{\upsilon _t}{v_z}}\end{array}} \right\}\end{equation}

(26) \begin{equation}{E_v} = \left( {{\upsilon _t} + \frac{{{M_\infty }}}{{{R_D}}}} \right)r\!\left\{ {\begin{array}{*{20}{c}}0\\[3pt]{2\dfrac{{\partial {v_r}}}{{\partial r}}}\\[1em]{\left( {\dfrac{{\partial {v_r}}}{{\partial z}} + \dfrac{{\partial {v_z}}}{{\partial r}}} \right)}\\[1em]0\\[5pt]{\dfrac{1}{\sigma }\dfrac{{\partial {\upsilon _t}}}{{\partial r}}}\end{array}} \right\}\quad {F_v} = \left( {{\upsilon _t} + \dfrac{{{M_\infty }}}{{{R_D}}}} \right)r\!\left\{ {\begin{array}{*{20}{c}}0\\[3pt]{\left( {\dfrac{{\partial {v_r}}}{{\partial z}} + \dfrac{{\partial {v_z}}}{{\partial r}}} \right)}\\[1em]{2\dfrac{{\partial {v_z}}}{{\partial z}}}\\[1em]0\\[3pt]{\dfrac{1}{\sigma }\dfrac{{\partial {\upsilon _t}}}{{\partial z}}}\end{array}} \right\}\quad H = r\!\left\{ {\begin{array}{*{20}{c}}{{H_p}}\\{{H_{qr}}}\\{{H_{qz}}}\\{{H_c}}\\{{H_{\upsilon t}}}\end{array}} \right\}\end{equation}

(27) \begin{equation}{H_p} = 0\end{equation}

(28) \begin{equation}{H_{qr}} = - 2{\upsilon _e}\dfrac{{{v_r}}}{{{r^2}}} + \dfrac{p}{{\rho r}}\end{equation}

(29) \begin{equation}{H_{qz}} = F\end{equation}

(30) \begin{equation}{H_c} = 0\end{equation}

(31) \begin{equation}{H_{\upsilon t}} = {c_{b1}}\left( {1 - {f_{t2}}} \right)\hat S{\upsilon _t} - \left[ {{c_{w1}}{f_w} - \frac{{{c_{b1}}}}{{{k^2}}}{f_{t2}}} \right]{\left( {\frac{{{\upsilon _t}}}{d}} \right)^2} + \frac{{{c_{b2}}}}{\sigma }\left[ {{{\left( {\frac{{\partial {\upsilon _t}}}{{\partial r}}} \right)}^2} + {{\left( {\frac{{\partial {\upsilon _t}}}{{\partial z}}} \right)}^2}} \right]\end{equation}

In Equation (29), F is the permeable disc force applied on the flow. Figure 2 illustrates the upper part of the permeable disc (profile view), the stream lines around it and the force applied on the flow, defined by Equation (32), where V is the volume of the computational cell, C _T is the thrust coefficient, and S _z is the z component of the vector of area.

(32) \begin{equation}F = \left\{\kern-4pt {\begin{array}{*{20}{l}}&0\qquad{{\textrm{away from the disk}}}\\[5pt]&{\dfrac{{{S_z}{C_T}}}{{2V}}} \quad{{\textrm{on the disk}}}\end{array}} \right.\end{equation}

Figure 2. Flow field around the upper part of the permeable disc (profile view).

The initial and boundary conditions required to solve the problem of the flow around the wind turbine are presented below:

(33) \begin{equation}{\textrm{Symmetry boundary: }}\left\{\kern-1pt {\begin{array}{*{20}{l}}{{v_r} = 0}\\[3pt]{\dfrac{{\partial {v_z}}}{{\partial r}} = 0}\\[7pt]{\dfrac{{\partial p}}{{\partial r}} = 0}\\[7pt]{\dfrac{{\partial {\upsilon _t}}}{{\partial r}} = 0}\end{array}} \right.\end{equation}

(34) \begin{equation}{\textrm{Inflow boundary: }}\left\{\kern-1pt {\begin{array}{*{20}{l}}{{v_r} = 0}\\{{v_z} = {U_\infty }}\\{p = {p_\infty }}\\[3pt]{{\upsilon _t} = \dfrac{3}{{{R_D}}}}\end{array}} \right.\end{equation}

(35) \begin{equation}{\textrm{Upper boundary: }}\left\{\kern-1pt {\begin{array}{*{20}{l}}{{v_r} = 0}\\{{v_z} = {U_\infty }}\\{p = {p_\infty }}\\[3pt]{{\upsilon _t} = \dfrac{3}{{{R_D}}}}\end{array}} \right.\end{equation}

(36) \begin{equation}{\textrm{Outflow boundary: }}\left\{\kern-1pt {\begin{array}{*{20}{c}}{\dfrac{{\partial {v_r}}}{{\partial z}} = 0}\\[5pt]{\dfrac{{\partial {v_z}}}{{\partial z}} = 0}\\[3pt]{p = {p_\infty }}\\[3pt]{\dfrac{{\partial {\upsilon _t}}}{{\partial z}} = 0}\end{array}} \right.\end{equation}

(37) \begin{equation}{\textrm{Initial conditions: }}\left\{\kern-1pt {\begin{array}{*{20}{c}}{{v_r} = 0}\\{{v_z} = {U_\infty } = 0.2}\\{p = {p_\infty } = 1}\\{{\upsilon _t} = \dfrac{3}{{{R_D}}}}\end{array}} \right.\end{equation}

4.0 Results and discussion

In the classical momentum theory, the wind turbine is replaced by a permeable disc. The disc is considered frictionless, without rotational velocity component in the wake [20]. Equation (38) shows the thrust coefficient as a function of the axial induction factor a, obtained from the momentum theory, while Equation (39) shows the definition of the axial induction factor, where u is the flow velocity at the permeable disc and V _o is the wind velocity [20].

(38) \begin{equation}{C_T} = 4a\!\left( {1 - a} \right)\end{equation}

where

(39) \begin{equation}u = \left( {1 - a)} \right){V_0}\end{equation}

However, experiments have shown that the assumption of the ideal wind turbine holds only for axial induction factor less than 0.4. Figure 3 shows experimental data found in the literature, basic momentum theory and (5, 17) corrections for C _T as a function of a.

Figure 3. Thrust coefficient C _T as a function of the axial induction factor a.

From Fig. 3, one can see that as C _T increases, the expansion of the wake also increases. Thus the velocity jump from V _o to u ₁ in the wake also increases (see Fig. 4). The reason that the basic momentum theory is not valid for values of a greater than 0.4 is that the free shear layer at the edge of the wake becomes unstable if the velocity jump is too high and eddies are formed, transporting momentum from the outer flow into the wake (known as Kelvin & Helmholtz instability).

Figure 4. The expansion of the wake and the velocity jump in the wake.

When the axial induction factor becomes larger than 0.4, the basic momentum theory fails, as can be seen in Fig. 3. Different empirical relations between the thrust coefficient C _T and the axial induction factor a can be constructed to fit better the measurements, see Refs (5, 17), given by Equations (40) and (41), respectively.

(40) \begin{equation}{C_T} = \left\{ {\begin{array}{*{20}{c}}{4a\!\left( {1 - a} \right)} \qquad\qquad\qquad\,\,\,\,{a \le 1/3}\\[3pt]{4a\!\left( {1 - \dfrac{1}{4}\left( {5 - 3a} \right)a} \right)}\qquad{a \gt 1/3}\end{array}} \right.\end{equation}

(41) \begin{equation}{C_T} = \left\{ {\begin{array}{*{20}{c}}{4a\!\left( {1 - a} \right)} \qquad\qquad\quad\,\,\,{a \le {a_c}}\\[3pt]{4\Big( {a_c^2 + \left( {1 - 2{a_c}} \right)a} \Big)}\qquad{a \gt {a_c}}\end{array}} \right.\end{equation}

In the present work, the thrust coefficient defined as a function of the axial induction factor is obtained numerically in order to be compared with the [5] correction factor. Results were obtained for three different Reynolds numbers, R _D =1.0 × 10⁴, 1.0 × 10⁶, and 1.0 × 10⁸, the numerical results will be discussed hereafter.

Figures 5, 6 and 7 present the axial induction factor field and stream lines around the wind turbine, for Reynolds number equal to 1.0 × 10⁸ and three different thrust coefficients, C _T =0.47, 1.08, and 1,57. High and low axial induction factors are printed in red and blue, respectively. As the thrust coefficient increases, the blockage of the wind turbine increases and higher values of the axial induction factor are observed behind the wind turbine, with higher deflections of the stream lines. In Fig. 7, where the thrust coefficient is equal to 1.56, the blockage of the wind turbine is too high forming a vortex behind it.

Figure 5. Axial induction factor around the wind turbine, C _T = 0.47 and R _D =1.0 × 10⁸.

Figure 6. Axial induction factor around the wind turbine, C _T = 1.08 and R _D = 1.0 × 10⁸.

Figure 7. Axial induction factor around the wind turbine, C _T = 1.57 and R _D = 1.0 × 10⁸.

Figure 8 presents the time series of the mean axial induction factor at the permeable disc, for Reynolds number equal to 1.0 × 10⁸ and three different thrust coefficients: C _T = 0.47, 1.08 and 1.57. For thrust coefficients equal to 0.47 and 1.57, steady state solutions are obtained and the axial induction factors reach constant values. On the other hand, for thrust coefficient equal to 1.08, an unsteady solution is obtained and the axial induction factor oscillates, as can be seen in Fig. 8. This unsteady solution for the thrust coefficient equal to 1.08 is probably caused by Kelvin & Helmhotz instability and is not predicted either by Refs. [5] or [17] corrections, but it was clearly captured here in the present work.

Figure 8. Time series of the axial induction factor for three different thrust coefficients and R _D = 1.0 × 10⁸.

Figure 9 presents a comparison of Refs. [5, 17], basic momentum theory, experimental data and present work numerical results for the thrust coefficient as a function of the axial induction factor. Numerical results were obtained for three different Reynolds numbers equal to 1.0 × 10⁴; 1.0 × 10⁶; 1.0 × 10⁸. The agreement with [5] correction is quite remarkable.

Figure 9. Comparison of the thrust coefficient obtained numerically for different Reynolds numbers and [5] correction.

For Reynolds number equal to 1.0 × 10⁴, the numerical results agree well with [5] correction for the entire range of the axial induction factor, except for the range 0.75 < a< 1.0. On the other hand, the numerical results for Reynolds numbers equal to 1.0 × 10⁶ and 1.0 × 10⁸ show a fairly good agreement with [5] correction, except for the ranges 0.75 < a < 1.0 and 0.4 < a < 0.7, where the flow around the permeable disc is unsteady, as discussed before in Fig. 8. Figure 10 shows the range of axial induction factor 0.4 < a < 0.7, where the flow around the permeable disc is unsteady for Reynolds numbers equal to 1.0 × 10⁶ and 1.0 × 10⁸.

Figure 10. Thrust coefficient for three different Reynolds number showing the unsteady behaviour.

Figure 11 presents authors’ suggestion (in red) for the thrust coefficient as function of the induction factor obtained by curve fitting of the numerical results to be used with the BEM. Equation (42) presents the mathematical representation of the suggested correction.

(42) \begin{equation}\begin{array}{l}{C_T} = \left\{ {\begin{array}{*{20}{c}}{4a\!\left( {1 - a} \right)} &\,\,\,\,\,{a \le 0.3}\\[3pt]{ - 1.035{a^3} + 1.807{a^2} + 0.3058a + 0.6237} &\,\,\,\,\,{a \gt 0.3}\end{array}} \right.\\{R_D} \le 1.0 \times {10^4}\end{array}\end{equation}

Figure 11. Present study correction based in numerical simulation.

Figure 12 presents the grid sensitive test for two different grid refinements, 100 × 150 and 150 × 150, and for Reynolds number equal to 1.0 × 10⁸. The change in grid refinement did not affect the thrust coefficient as a function of the axial induction factor. Therefore, the grid with 100 × 150 grid points could be utilised in the entire investigation.

Figure 12. Grid sensitivity test for two different grid refinements and R_D = 1.0 × 10⁸.

5.0 Conclusions

The PD/URANS method utilised to study the viscous and incompressible flow around a wind turbine for three different Reynolds numbers generated quite accurate and reliable results.

Calculation of the axial induction factor field around the permeable disc, time series of the averaged axial induction factor on the permeable disc surface, and the thrust coefficient as a function of the axial induction factor could be performed in a quite efficient way.

For the three Reynolds numbers studied, the agreement between the numerical results and original [5] correction for the thrust coefficient as a function of the axial induction factor is remarkable for the entire range of the axial induction factor, except for the range where the flow around the disc is very unsteady. But it is well known that [5] correction is not able to capture the unsteady state behaviour and the Reynolds number dependency of the flow around a wind turbine. In aeronautical practical applications [5], correction is taken as a reasonable approximation for the thrust coefficient as a function of the axial induction factor, but, unfortunately, it was obtained empirically from a curve fitting from quite scattered experimental data set, as seen in Fig. 9. Therefore, the results obtained numerically using the PD/URANS method can be considered to be much more reliable than those obtained using [5] correction.

The PD/URANS method solves three dimensional and axisymmetric flows around wind turbines in cylindrical coordinates. In that case, the flow description depends only on two coordinates (r, z), allowing to be solved in a two-dimensional computational grid. Therefore, the PD/URANS method can be very efficient, accurate and relatively small computational time, compatible with standard desktop/laptop computers.

The BEM is usually applied on a first approximation to study the aerodynamic characteristics of a wind turbine, but instead of using the traditional [5] correction, the suggested correction in Equation (42) improves the accuracy of the BEM method. However, the PD/URANS method may be used even for more accurate studies.