A new class of actuator surface models for wind turbines

TL;DR

This paper introduces a novel actuator surface model for wind turbines that captures detailed blade and nacelle geometries, improving wake prediction accuracy on coarse meshes compared to traditional models.

Contribution

The paper develops a new actuator surface model for turbine blades and nacelle, incorporating detailed geometrical features and nacelle effects, enhancing simulation fidelity.

Findings

The model accurately predicts turbine wake characteristics at far wake locations.

It performs well on coarse meshes, reducing computational costs.

The approach is applicable to other bluff bodies in flow simulations.

Abstract

Actuator line model has been widely employed in wind turbine simulations. However, the standard actuator line model does not include a model for the turbine nacelle which can significantly impact turbine wake characteristics as shown in the literature. Another disadvantage of the standard actuator line model is that more geometrical features of turbine blades cannot be resolved on a finer mesh. To alleviate these disadvantages of the standard model, we develop a new class of actuator surface models for turbine blades and nacelle to take into account more geometrical details of turbine blades and include the effect of turbine nacelle. In the actuator surface model for blade, the aerodynamic forces calculated using the blade element method are distributed from the surface formed by the foil chords at different radial locations. In the actuator surface model for nacelle, the forces are distributed from the actual nacelle surface with the normal force component computed in the same way as in the direct forcing immersed boundary method and the tangential force component computed using a friction coefficient and a reference velocity of the incoming flow. The actuator surface model for nacelle is evaluated by simulating the flow over periodically placed nacelles. Both the actuator surface simulation and the wall-resolved large-eddy simulation are carried out. The comparison shows that the actuator surface model is able to give acceptable results especially at far wake locations on a very coarse mesh. It is noted that although this model is employed for the turbine nacelle in this work, it is also applicable to other bluff bodies. The capability of the actuator surface model in predicting turbine wakes is assessed by simulating the flow over the MEXICO (Model experiments in Controlled Conditions) turbine and a hydrokinetic turbine.

Equations

$$\mathbf{L}=\frac{1}{2}\rho C_{L}c|{\mathbf{V}}_{rel}{|}^2{\mathbf{e}}_L$$

$$\mathbf{D}=\frac{1}{2}\rho C_{D}c|{\mathbf{V}}_{rel}{|}^2{\mathbf{e}}_D,$$

$$\alpha =\varphi -\gamma ,$$

$${\mathbf{V}}_{rel}=u_x{\mathbf{e}}_x+(u_{\theta }-\Omega r){\mathbf{e}}_{\theta },$$

$$u_x=\frac{1}{c}{\int }_c\mathbf{u}\left(\mathbf{X}\right)\cdot {\mathbf{e}}_{x}ds$$

$$u_{\theta }=\frac{1}{c}{\int }_c\mathbf{u}\left(\mathbf{X}\right)\cdot {\mathbf{e}}_{\theta }ds,$$

$$\mathbf{u}\left(\mathbf{X}\right)=\sum _{\mathbf{x}\in g_{\mathbf{x}}}\mathbf{u}(\mathbf{x}){\delta }_h\left(\mathbf{x}-\mathbf{X}\right)V(\mathbf{x}),$$

$$\phi(r)=\left\{\begin{array}[]{llr}\frac{1}{4}+\frac{\sin\left({\pi}\left(2|r|+1\right)/4\right)}{2\pi}-\frac{\sin\left({\pi}\left(2|r|-1\right)/4\right)}{2\pi},&&|r|\leq 1.5,\\ \frac{5}{8}-\frac{|r|}{4}-\frac{\sin\left({\pi}\left(2|r|-1\right)/4\right)}{2\pi},&&1.5\leq|r|\leq 2.5,\\ 0,&&2.5\leq|r|,\\ \end{array}\right.$$

$$C_{L,3D}=C_{L,2D}+f_L\left(C_{L,p}-C_{L,2D}\right),$$

$$C_{D,3D}=C_{D,2D}-f_D\left(C_{D,2D}-C_{D,0}\right),$$

$$f_L=\frac{1}{2\pi }\left(\frac{1.6(c/r)a-(c/r{)}^{\frac{d}{\Lambda }\frac{R}{r}}}{0.1267b+(c/r{)}^{\frac{d}{\Lambda }\frac{R}{r}}}-1\right),$$

$$f_D=\frac{1}{2\pi }\left(\frac{1.6(c/r)a-(c/r{)}^{\frac{d}{2\Lambda }\frac{R}{r}}}{0.1267b+(c/r{)}^{\frac{d}{2\Lambda }\frac{R}{r}}}-1\right),$$

$$C_L=F_1{C}_{L,3D},$$

$$C_D=F_1{C}_{D,3D},$$

$$F_1=\frac{2}{\pi }{\cos }^{-1}\left(\exp \left(-g\frac{B(R-r)}{2r\sin \varphi }\right)\right),$$

$$g=\exp \left(-c_1\left(B\Omega R/U_{\infty }-c_2\right)\right)+c_3,$$

$$\mathbf{f}(\mathbf{X})=\left(\mathbf{L}+\mathbf{D}\right)/c.$$

$$\mathbf{f}(\mathbf{x})=-\sum _{\mathbf{X}\in g_{\mathbf{X}}}\mathbf{f}(\mathbf{X}){\delta }_h\left(\mathbf{x}-\mathbf{X}\right)A(\mathbf{X}),$$

$${\mathbf{f}}_n(\mathbf{X})=\frac{h\left(-{\mathbf{u}}^d(\mathbf{X})+\tilde{\mathbf{u}}(\mathbf{X})\right)\cdot {\mathbf{e}}_n(\mathbf{X})}{\Delta t}{\mathbf{e}}_n(\mathbf{X}),$$

$$\tilde{\mathbf{u}}(\mathbf{x})={\mathbf{u}}^n(\mathbf{x})+{\mathbf{rhs}}^n(\mathbf{x})\Delta t,$$

$${\mathbf{f}}_{\tau }(\mathbf{X})=\frac{1}{2}{c}_f{U}^2{\mathbf{e}}_{\tau }(\mathbf{X})$$

$$c_f=0.37(\log R{e}_x{)}^{-2.584},$$

$${\mathbf{e}}_{\tau }(\mathbf{X})=\frac{\mathbf{u}\left(\mathbf{X}+h{\mathbf{e}}_n\left(\mathbf{X}\right)\right)}{|\mathbf{u}\left(\mathbf{X}+h{\mathbf{e}}_n\left(\mathbf{X}\right)\right)|}.$$

$$J\frac{\partial U^i}{\partial {\xi }^i}=0,$$

$$\frac{1}{J}\frac{\partial U^i}{\partial t}=$$

Full text

††articletype: RESEARCH ARTICLE

1 Introduction

Using the current state-of-the-art supercomputers, it is still impossible to simulate wind farms resolving all the relevant length scales from the thickness of the boundary layer over a wind turbine blade ($\approx$0.016 m at the trailing edge estimated using the equation ${\delta}/{c}=0.37\left({U_{\infty}c}/{\nu}\right)^{-0.2}$[1] with $\nu=1.5\times 10^{-5}$ $m^{2}/s$, and airfoil chord $c=1$ m and $U_{\infty}=100$ $m/s$ at the blade tip for a MW scale wind turbine with rotor diameter approximately 100 $m$, tip speed ratio 10 and incoming wind speed 10 $m/s$) to the thickness of an atmospheric boundary layer ($\approx$1000 $m$ [2]). Different rotor models including the actuator disk model [3, 4, 5], the actuator line model [6, 7] and the actuator surface model [8] have been developed in the literature to parametrize the interaction between the turbine blades and the incoming flow without directly resolving the boundary layer flow over turbine blades.

In the actuator disk model [9] the wind turbine rotor is modelled by a circular permeable disk. The impacts of incoming wind on wind turbines are taken into account by uniformly distributed forces on the disk. In the conventional actuator disk model, the forces are only applied along the turbine axial direction, which can be calculated either from one-dimensional momentum theory with a prescribed induction factor [3] or from the blade-element momentum theory [10]. In the work by Wu and Porté-Agel [11], the rotation was taken into account in the actuator disk model. It was found that the predictions from the actuator disk model with rotation show better agreement with the wind tunnel measurements reported in [12].

The actuator line model provides a more accurate approach to incorporate the rotational effect, finite blade number effect and the effect of non-uniform force distribution in the azimuthal direction. In this model, the turbine blade is represented by a rotating line with distributed forces, which are calculated from a blade element approach [13] combined with tabulated drag and lift coefficients of the foils. The employed drag and lift coefficients are usually from two-dimensional numerical simulations or experiments without considering the rotational effects. In order to take into account the three-dimensional rotational effects, corrections on the two-dimensional force coefficients, such as those proposed in [14], are usually employed. As discussed by Shen et al. [15], the lift needs to approach zero as the blade tip is approached. However, the lift on the tip predicted by the blade element approach is usually not zero especially for a pitched blade or a chambered foil. A tip-loss correction in addition to the correction for the three-dimensional effect is needed in order to have a correct lift distribution in the region near the blade tip as suggested by Shen et al. [16]. There are two major limitations of the standard actuator line model: i) The lack of an effective nacelle model. Nacelle induced coherent structures were shown to have a significant impact on turbine wake characteristics, such as velocity deficit in the near wake [17, 7, 18] and wake meandering at far wake locations for model wind turbines [19, 20, 21] and hydrokinetic turbines [17]. Both geometry-resolving, wall-modelled large-eddy simulations and actuator type large-eddy simulations were carried out in [17]. It is shown that the actuator line model without a nacelle model cannot accurately capture the wake meandering of the hydrokinetic turbine, and underpredicts the turbulence intensity at far wake locations. This was attributed to the fact that the presence of the nacelle induces spiral vortex breakdown of the hub vortex causing a large-scale energetic coherent vortex that precesses and expands radially outward to intercept and energize the tip shear layer at three to four rotor diameters downwind of the turbine [17]. The nacelle can be represented using a permeable disk with a specified drag coefficient. However, Yang et al. [7] showed that such a simple nacelle model is not able to take into account the nacelle effect accurately. The computed results from Tossas et al. [18] indicated that with their nacelle and tower model the velocity deficit in the near wake is predicted more accurately. The improvement on the prediction of the turbulence intensity, on the other hand, is very limited. ii) A finer mesh cannot resolve more geometrical features of the turbine blade. The results of Yang et al. [7] showed that grid-independent results cannot be obtained when the mesh is refined in the actuator line model simulations. This is because more physics cannot be resolved on a finer grid in the state-of-the-art actuator type models. The accuracy of the computed results mainly depends on the parameterizations employed in actuator type models.

The objective of this paper is to develop a new class of actuator surface models for turbine blade and nacelle to take into account nacelle geometry and better incorporating the geometrical effect of turbine blade. In the actuator surface model for blade, the forces are still computed using the blade element method but are distributed on the surface formed by the foil chords at different radial locations. In the actuator surface model for the nacelle, the normal forces are computed by satisfying the non-penetration boundary conditions as in the direct forcing immersed boundary method [22, 23, 24]. The tangential forces, on the other hand, are calculated using a friction coefficient and a reference velocity of the incoming flow. More details on the proposed actuator surface models are presented in Sec. 2 followed by a brief description of the governing equations and flow solver in Sec. 3. Some cases for validating the proposed actuator surface model are then reported in Sec. 4. Finally, conclusions are drawn in Sec. 5.

2 Actuator surface models

In this section, we present the proposed actuator surface models for the blade and nacelle. In both actuator surface models, we have two sets of independent meshes, i.e., the background Cartesian grid for the flow with its coordinate denoted by $\mathbf{x}$ ($x$, $y$, $z$ or $x_{1}$, $x_{2}$, $x_{3}$), and the Lagrangian grid following the actuator surfaces with its coordinate denoted by $\mathbf{X}$ ($X$, $Y$, $Z$ or $X_{1}$, $X_{2}$, $X_{3}$). For both models, the smoothed discrete delta function in [24] is employed as the kernel function for transferring information between the two meshes. The major difference between the actuator surface model for blade and the actuator surface model for nacelle is the way how the forces on the actuator surfaces are calculated. In Sec. 2.1 we present the force calculation and distribution methods of the actuator surface model for blade. And in Sec. 2.2 we introduce the forces calculation method for the actuator surface model for nacelle.

2.1 Actuator surface model for blade

In the actuator surface model for blade, the blade geometry is represented by a surface formed by the chord lines at different radial locations of a blade as shown in Fig. 1(a). The forces are calculated in the same way as in the actuator line model. The lift ($\mathbf{L}$) and drag ($\mathbf{D}$) at each radial location are calculated by

[TABLE]

and

[TABLE]

where $C_{L}$ and $C_{D}$ are the lift and drag coefficients, $c$ is the chord length, $\mathbf{V}_{rel}$ is the relative incoming velocity, and $\mathbf{e}_{L}$ and $\mathbf{e}_{D}$ are the unit vectors in the directions of lift and drag, respectively. The $C_{L}$ and $C_{D}$ are functions of the Reynolds number and the angle of attack. The angle of attack $\alpha$ shown in Fig. 1(b) is computed at each radial location by

[TABLE]

where $\phi=-\tan^{-1}(u_{x}/(u_{\theta}-\Omega r))$, and $\gamma$ is the angle including the blade twist and the blade pitch, the latter of which is specified to a given value instead of from a turbine control algorithm. The relative incoming velocity $\mathbf{V}_{rel}$ at each radial location is computed by

[TABLE]

where $\Omega$ is the rotational speed of the rotor, $\mathbf{e}_{x}$ and $\mathbf{e}_{\theta}$ are the unit vectors in the axial and azimuthal directions, respectively. The $u_{x}$ and $u_{\theta}$ are the axial and azimuthal components of the flow velocity averaged over the chord for each radial locations, which are computed using

[TABLE]

and

[TABLE]

where $\mathbf{X}$ denote the coordinates of the grid points on the actuator surfaces. Generally the grid points on the actuator surfaces do not coincide with any background nodes. We employ a smoothed discrete delta function (i.e. the smoothed four-point cosine function) proposed by Yang et al. [24] to interpolate $\mathbf{u}\left(\mathbf{X}\right)$ from the values on the background grid nodes as follows:

[TABLE]

where $\mathbf{x}$ are the coordinates of the background grid nodes, $g_{\mathbf{x}}$ is the set of the background grid cells, $V=h_{x}h_{y}h_{z}$ ($h_{x}$, $h_{y}$ and $h_{z}$ are the grid spacings in the $x$, $y$ and $z$ directions, respectively) is the volume of the background grid cell, $\delta_{h}\left(\mathbf{x}-\mathbf{X}\right)=\frac{1}{V}\phi\left(\frac{x-X}{h_{x}}\right)\phi\left(\frac{y-Y}{h_{y}}\right)\phi\left(\frac{z-Z}{h_{z}}\right)$ is the discrete delta function, and $\phi$ is the smoothed four-point cosine function [24], which is expressed as

[TABLE]

in which $r_{i}=(x_{i}-X_{i})/h_{i}$ ($i$=1, 2, 3).

The blade rotation causes the stall delay phenomenon at the inboard sections of the blade, which increases the lift coefficients and decreases the drag coefficients as compared with the corresponding two-dimensional airfoil data. To account for such three-dimensional rotational effect, the stall delay model developed by Du and Selig [25] is employed to correct the lift and drag coefficients from two-dimensional experiments or numerical simulations. In Du and Selig’s model, the corrected lift and drag coefficients ($C_{L,3D}$ and $C_{D,3D}$) are calculated as follows:

[TABLE]

and

[TABLE]

where $C_{L,p}=2\pi\left(\alpha-\alpha_{0}\right)$, $C_{D,0}$ the two-dimensional drag coefficient at zero angle of attack, and the correction functions $f_{L}$ and $f_{D}$ are determined by

[TABLE]

and

[TABLE]

respectively, where $\Lambda=\Omega R/\sqrt{U^{2}+(\Omega R)^{2}}$, $U$ is the incoming wind speed, $R$ is the rotor radius, $a$, $b$ and $d$ are the empirical correction factors. In this work, a, b and d are equal to 1 as in Du and Selig’s paper [25].

Non-zero force can exist at the blade tip when the pitch angle is nonzero or a chambered foil is used [15]. This is in contradiction with the physical understanding that the force should tend to zero at the tip due to pressure equalization as discussed by Shen et al. [15]. To correct this non-physical force behaviour, the tip-loss correction proposed by Shen et al. [16, 15] is applied to the drag and lift coefficients computed from Eqs. (9) and (10). With the tip loss correction, the employed $C_{D}$ and $C_{L}$ are calculated as

[TABLE]

and

[TABLE]

where

[TABLE]

in which $B$ is the number of blades, $g$ is computed by

[TABLE]

where $c_{1}$, $c_{2}$ and $c_{3}$ are the correction coefficients, which are equal to 0.125, 21 and 0.1 as in [16], respectively.

After calculating the lift ($\mathbf{L}$) and drag ($\mathbf{D}$), the force per unit area on the actuator surface at each radial location is then calculated as

[TABLE]

Note that the above expression essentially means that the lift and drag on the blade are uniformly distributed in the chordwise direction. To calculate the forces on the background mesh for the flowfield, the computed forces on the actuator surfaces are then distributed to the background grid nodes as follows:

[TABLE]

where $g_{\mathbf{X}}$ is the set of the actuator surface grid cells and $A$ is the area of the actuator surface grid cell. The same discrete delta function as in Eq. (7) is employed. It is noted the negative sign is because that $\mathbf{f}(\mathbf{x})$ represents the forces of the actuator surfaces acting on the flow, while $\mathbf{f}({\mathbf{X}})$ denotes the forces of the flow acting on the actuator surfaces.

2.2 Actuator surface model for nacelle

Resolving the boundary layer over a nacelle of an utility scale wind turbine, which becomes almost impossible for wind farm scale simulations, requires a much finer mesh as compared with the mesh for resolving the thickness of an atmospheric boundary layer. To take into account the nacelle effects on a relatively coarse mesh, in this section we present an actuator surface model for parametrizing the nacelle effect without directly simulating the boundary layer flows over the nacelle.

In this model, the nacelle geometry is represented by the actual surface of the nacelle with distributed forces. As shown in Fig. 2, the force acting on the actuator surface is decomposed into the normal component and the tangential component. From a physical point of view, the actuator surface of the nacelle cannot act as a sink or source of mass not matter how the boundary layer over the nacelle is modelled. Based on this consideration, the normal component of the force on the actuator surface is computed by satisfying the non-penetration constraint. The non-penetration constraint can be met by directly reconstructing the wall-normal velocity near the actuator surface as in the sharp interface immersed boundary method [22, 26, 27, 28]. However, a physically reasonable interpolation scheme is difficult to implement because of the very coarse gird employed in the actuator type simulations. Moreover direct velocity reconstruction cannot ensure a smooth representation of the nacelle geometry because of the very coarse grid. In this work we represent the nacelle geometry as a diffused interface the same as in the direct forcing immersed boundary method [23, 24, 29]. In this actuator surface model for nacelle, the normal component of the force acting on the surface per unit area is calculated in the same way as in the direct forcing immersed boundary method, which is expressed as follows:

[TABLE]

where $\mathbf{u}^{d}(\mathbf{X})$ is the desired velocity on the nacelle surface, $\mathbf{e}_{n}(\mathbf{X})$ is the unit vector in the normal direction of the nacelle surface, $h=(h_{x}h_{y}h_{z})^{1/3}$ is the length scale of the local background grid spacing (different values of $h$ have been tested without noticing any significant differences), $\mathbf{\tilde{u}}(\mathbf{X})$ is the estimated velocity on the actuator surface, which is interpolated from the corresponding estimated velocity on the background mesh using Eq. (7) with the $\mathbf{\tilde{u}}(\mathbf{x})$ on the background mesh computed as

[TABLE]

where $\Delta t$ is the size of the time step, the right-hand-side term $\mathbf{rhs}^{n}$ includes the convection, pressure gradient and diffusion terms computed from the quantities of previous time step $n$.

The tangential force acting on the surface depends on the incoming velocity, the surface geometry and the complex near-wall turbulence. Neither of the last two can be directly captured in the actuator type simulations by using the no-slip boundary conditions or the shear stress boundary conditions for wall models [30, 31]. In the present actuator surface model, we assume that the tangential force is proportional to the local incoming velocity, and the effects of surface geometry and near-wall turbulence can be parameterized using a single parameter, the friction coefficient $c_{f}$. Based on this assumption, the tangential force acting on the surface per unit area is computed as

[TABLE]

where $c_{f}$ is calculated from the empirical relation proposed by F. Schultz-Grunow [1] for turbulent boundary layers with zero pressure gradient, which is expressed as follows:

[TABLE]

where $Re_{x}$ is the Reynolds number based on the incoming velocity and the distance from the upstream edge of the immersed body. It is noticed that this expression is invalid for the region with a large pressure gradient. For the present nacelle simulation, however, zero pressure gradient can be regarded as a reasonable assumption for the most part of the turbine nacelle. The framework presented in this work is applicable to more complex geometries if the corresponding distribution of $c_{f}$ is available from experiments or high fidelity simulations. In Eq. (21), $U$ is the reference velocity of the incoming flow. In the present work, it is selected as the magnitude of the incoming streamwise velocity. The direction of the tangential force $\mathbf{e}_{\tau}(\mathbf{X})$, on the other hand, is determined by the local tangential velocity relative to the velocity at the nacelle surface at points located at $h$ away from the wall, i.e., point B in Fig. 2, as follows:

[TABLE]

After computing the forces on the actuator surface of the nacelle, the forces on the background mesh are distributed from the actuator surface mesh in the same way as in Eq. (18). It is noted that the forces are distributed over the closest five cells defined by the width of the employed smoothed cosine discrete delta function. This force distribution width will be very small if we have a sufficiently high grid resolution to resolve the boundary layer. For the grids using in the actuator type simulations, however, this force distribution width is so large that it may cause non-physical velocity field near the nacelle and affect the calculation of the relative incoming velocity for the blade. To remedy this problem, we simply let the relative incoming velocity ($\mathbf{V}_{rel}$ in Eq. (4)) of the two closest radial locations to the nacelle surface be equal to that of the third closest location away from the nacelle.

3 Flow solver

The governing equations are the three-dimensional, unsteady, filtered continuity and Navier-Stokes equations in non-orthogonal, generalized, curvilinear coordinates, which read in compact tensor notation (repeated indices imply summation) as follows ($i,j=1,2,3$):

[TABLE]

[TABLE]

where $x_{i}$ and $\xi^{i}$ are the Cartesian and curvilinear coordinates, respectively, $\xi_{l}^{i}={\partial\xi^{i}}/{\partial x_{l}}$ are the transformation metrics, $J$ is the Jacobian of the geometric transformation, $u_{i}$ is the $i^{th}$ component of the velocity vector in Cartesian coordinates, $U^{i}$=${(\xi_{m}^{i}/J)u_{m}}$ is the contravariant volume flux, $g^{jk}=\xi_{l}^{j}\xi_{l}^{k}$ are the components of the contravariant metric tensor, $\rho$ is the density, $\mu$ is the dynamic viscosity, $p$ is the pressure, $f_{l}(l=1,2,3)$ are the body forces introduced by the turbine blade and nacelle computed using the actuator surface models presented in Sec. 2, and $\tau_{ij}$ represents the anisotropic part of the subgrid-scale stress tensor, which is modelled by the dynamic eddy viscosity subgrid-scale stress model [32]. The governing equations are discretized in space using a second-order accurate central differencing scheme, and integrated in time using the fractional step method [28]. An algebraic multigrid acceleration along with GMRES solver is used to solve the pressure Poisson equation. A matrix-free Newton-Krylov method is used for solving the discretized momentum equation.

4 Test cases

In this section, we evaluate the performance of the proposed actuator surface models for turbine blade and nacelle. The nacelle body case in Sec. 4.1 is to evaluate the capability of the actuator surface model for nacelle. The MEXICO (Model experiments in Controlled Conditions) turbine [33, 34, 35] case in Sec. 4.2 is to show the capability of the actuator surface models in predicting the tip vortices and velocity distribution in the near wake. The capability of the new actuator surface models on the wake meandering prediction is demonstrated by the hydrokinetic turbine case in Sec. 4.3. In turbine simulations the proposed actuator surface model is abbreviated as ASB-ASN when both the actuator surface blade and actuator surface nacelle models are employed, and ASB when only the blade model is employed. The conventional actuator line simulations are also carried out for comparison. The actuator line model for turbine blade is abbreviated as ALB. The actuator line model for blades with the actuator surface model for nacelle is abbreviated as ALB-ASN. Throughout this paper, $x$, $y$ and $z$ represent the streamwise, spanwise and vertical directions, respectively.

4.1 Flow over periodically placed nacelles

In this section, the actuator surface model for nacelle is applied to simulate the flow over periodically placed nacelles to mimic wind farm scenarios. The nacelle body is composed of a hemisphere facing upstream and a circular cylinder in its downstream as shown in Fig. 3. Free-slip boundary condition is applied in the crosswise directions. Periodic boundary condition is employed in the streamwise direction. The Reynolds number based on the cylinder radius $R$ and the freestream velocity is 1000. The size of the computational domain is $L_{x}\times L_{y}\times L_{z}=30R\times 20R\times 20R$. The downstream end of the nacelle body is located at the origin of the streamwise coordinate. Three grids are employed with the number of grid nodes $N_{x}\times N_{y}\times N_{z}=502\times 348\times 348$, $N_{x}\times N_{y}\times N_{z}=115\times 151\times 151$ and $N_{x}\times N_{y}\times N_{z}=153\times 80\times 80$, grid spacings near the nacelle $\Delta x=R/8,R/3.75$ and $R/2.5$, $\Delta y=\Delta z=R/33,R/7.5$ and $R/5$, and the corresponding time steps $0.01R/U$, $0.1R/U$ and $0.1R/U$ for the fine, medium and coarse grids, respectively. Large-eddy simulation with the dynamic subgrid-scale model [32] is employed for the flow simulation [36]. On the coarse grid, two actuator type large-eddy simulations are carried out: one uses the proposed actuator surface model for nacelle; the other one represents the nacelle by a permeable disk located at the downstream end of the nacelle. In the permeable disk simulation, the mean drag force computed in the coarse grid actuator surface simulation is uniformly distributed on the disk, which gives a drag coefficient of $C_{D}=0.48$ defined using the freestream velocity and the area $\pi R^{2}$. On the medium grid, only the actuator surface simulation is carried out. On the fine grid, the geometry-resolving, wall-resolved large-eddy simulation is carried out with the sharp-interface immersed boundary method developed by Ge et al. [28] for resolving the boundary layer and near wall turbulence, from which the computed results is employed to gauge the performance of the proposed actuator surface model. In both actuator surface and immersed boundary simulations, the number of surface triangular cells is 2652. In the permeable disk simulation, the number of surface triangular cells is 234.

The contours of the instantaneous streamwise velocity computed from the coarse grid actuator surface simulation is compared with the wall-resolved LES in Fig. 4. As seen, the large-scale flow structure in the nacelle wake is well predicted by the actuator surface model on a very coarse grid.

The mean velocity profiles computed from the actuator type models are compared with that from the wall-resolved LES in Fig. 5. As seen, the velocity deficit computed from actuator surface model agrees well with that from the wall-resolved LES at $1R$ downstream. At $3R$, $5R$ and $7R$ downstream locations, however, the velocity deficits computed from the actuator surface model on the coarse grid are significantly larger. For further downstream locations, on the other hand, the velocity deficits predicted by the actuator surface model show good agreement with that computed by the wall-resolved LES. The velocity deficits computed by the actuator surface model on the medium grid show good agreements with that from the wall-resolved LES at all downstream locations. The permeable disk model significantly underpredicts the velocity deficits at all downstream locations, although the drag force with the same magnitude is applied on the permeable disk.

The turbulence kinetic energy profiles computed from the actuator type models are compared with that from the wall-resolved LES in Fig. 5. At $1R$ downstream location, the actuator surface model gives acceptable predictions on both coarse and medium grids. At $3R$ downstream location, the actuator surface model significantly underpredicts the turbulence kinetic energy on the coarse grid. On the medium grid, on the other hand, a significant improvement is observed. At $7R$ to further downstream locations, the turbulence kinetic energy profiles computed on coarse and medium grids are similar with each other and are larger that that from the wall-resolved LES. It is also noticed that the peaks of turbulence kinetic energy predicted by the actuator surface simulations are different from that predicted by wall-resolved LES at $3R$, $5R$ and $7R$ downstream locations. The actuator disk model, on the other hand, predicts nearly negligible turbulence kinetic energy at all downstream locations.

To further evaluate the performance of the actuator surface model, the power spectral density (PSD) predicted by the actuator surface simulation is compared with that from the wall-resolved LES at different downstream locations in Fig. 7. It is seen the low frequency, energy-containing part of the PSD predicted by the actuator surface model shows an overall good agreement with the wall-resolved LES. Both the actuator surface and wall-resolved simulations predict the inertial range of the PSD with the slope of $-5/3$. However, the inertial range predicted in the actuator surface simulation is narrower than that from the wall-resolved LES at $1R$, $3R$ and $5R$ downstream locations. This implies that the low spatial resolution is the main reason for the discrepancy observed in the turbulence kinetic energy profiles in the near wake locations (Fig. 6), where the small scale energetic turbulent flow structures are not resolved by the grid. At far wake locations, both energy-containing part and inertial part of the PSD predicted by the actuator surface model agree well with that computed from the wall-resolved LES, even though the very high frequency, dissipation range is not resolved in the actuator surface simulation.

4.2 Flow over the MEXICO turbine

In this section, we simulate the flow over the MEXICO turbine to evaluate the performance of the proposed actuator surface models on predicting the velocity and tip vortices at turbine near wake locations. The size of the computational domain is $L_{x}\times L_{y}\times L_{z}=4.2D\times 4D\times 4D$. The number of grid nodes is $N_{x}\times N_{y}\times N_{z}=673\times 502\times 502$. The mesh is uniform across the turbine and in the streamwise direction with grid spacing $D/160$, and is stretched in the spanwise and vertical directions at regions away from the turbine. The size of time step is $T/750$, where $T$ is rotational period of the turbine rotor. In the actuator surface simulations, the number of triangular cells for actuator blades and nacelle is 1461 and 5464, respectively. In the actuator line simulation, the number of segments is 141. Free-slip boundary condition is applied at the boundaries in the spanwise and vertical directions. The flow at the inlet is uniform. Three cases with three different tip speed ratios ($\text{TSR}=\Omega R/U$ where $\Omega$ is the rotational speed of the rotor, $R$ is the radius of the rotor and $U$ is the incoming wind speed), i.e. 4.2, 6.7 and 10, are carried out. In the experiment, different tip speed ratios were realized by adjusting the incoming wind speed. In the LES, on the other hand, different rotational speeds are employed for different tip speed ratios. The Reynolds number based on the incoming wind speed and the rotor diameter is $Re=6.9\times 10^{6}$ the same for the cases with different tip speed ratios. Both ASB-ASN and ALB simulations are carried out. In the ASB-ASN simulation, the nacelle is simply represented by a circular cylinder with diameter $0.2D$ and length $0.2D$.

Before presenting the computed results in the turbine near wake, the bound circulations along the blade predicted by the actuator surface model and the conventional actuator line model are compared with measurements in Fig. 8. The bound circulation is computed by $\Gamma=\frac{L}{\rho|\mathbf{V}_{rel}|}$. As seen, the location with maximum bound circulation is located near the blade root for TSR=4.2, in the middle for TSR=6.7, while close to the tip for TSR=10. Such dependence of the bound circulation distributions on the values of TSR is well captured by the proposed actuator surface model. The amplitudes of the predicted bound circulation also show an acceptable agreement with the measurements except for the TSR=10 case, where the discrepancies between the computed values and measured ones are somewhat larger than the other two cases. The actuator line model predicts with good accuracy the bound circulation for the TSR=6.7 and 10 cases, while underpredicts the bound circulation for the TSR=4.2 case.

The capability of the proposed actuator surface model in predicting the tip vortices is examined in Fig. 9, in which we show the comparison of the computed tip vortex circulations with the measurements for the three cases with different tip speed ratios. To calculate the circulation of tip vortices, the center of a tip vortex core is selected as the position with the maximum vorticity magnitude. The circulation is computed by integrating the velocity along a circle centered at the vortex core center with radius of $0.05D$. The computed tip vortex circulations show overall good agreements with the measurements for the actuator surface model. Some discrepancies are observed for the TSR=10 case for the third and the fourth tip vortices downwind from the rotor for the actuator surface model. While the circulation of tip vortices predicted by the actuator line model agrees well the measurements for the TSR=10 case. From the comparison of the tip vortex circulations among the three TSR cases, it is interesting to notice that they are in the same order when normalized using $U$ and $D$, although the dimensional circulations are significantly different for the three cases.

The contours of the instantaneous streamwise velocity are shown in Fig. 10 for different tip-speed ratios for both ASB-ASN and ALB simulations. As seen, a jet exists at the center in the ALB results. On the other hand, a wake from the nacelle forms in the ASB-ASN results. This nacelle wake substantially interacts with the outer part of the wake and causes intensive velocity fluctuations in the turbine wake. Although significantly differences in turbulence statistics between ASB-ASN and ALB results are expected, this work will only focus on the mean velocity profiles for which the wind tunnel measurements are available for validation purpose.

The computed axial ($\langle u_{a}\rangle$) and radial ($\langle u_{r}\rangle$) velocities are compared with measurements at different radial locations in Fig. 11 and Fig. 12, respectively. As seen, the radial velocity $\langle u_{r}\rangle$ increases as the flow approaches turbine, and decreases rapidly in the immediate downwind of the turbine for all the four considered radial locations for all the three cases. It is also seen that such increase of the radial velocity is somewhat larger when TSR is higher. From these comparisons, we can see that the developed actuator surface model can accurately predict the downwind variations of radial velocity. For the axial velocity, a sudden decrease and increase in the near wake of the turbine are observed at $r/D=0.52$ for all the three tip-speed ratios for the measured data. This is, however, because of the corrupted PIV photos at these locations caused by the light reflected on the turbine hub as noticed in the paper by Nilsson et al. [35]. The experimental data from the corrupted PIV images will be ignored for the comparisons. It is seen that overall good agreements with measurements are obtained for the axial velocities although some discrepancies are observed for the TSR=10 case. The actuator line model, on the other hand, underpredicts the velocity deficits for the TSR=6.7 and 10 cases.

4.3 Flow over a hydrokinetic turbine

The capability of the developed actuator surface models in predicting the wake meandering at far wake locations is assessed by simulating the turbulent flow over the axial-flow hydrokinetic turbine at Saint Anthony Falls Laboratory, University of Minnesota [37]. Kang, Yang and Sotiropoulos [17] carried out the geometry-resolving, actuator disk and actuator line simulations of this turbine. It was found that the conventional actuator line model without a nacelle model underpredicts the turbulence intensity in the far wake because it cannot accurately capture the inner-outer wake interaction and thus the wake meandering in the far wake.

The computational domain is $L_{x}\times L_{y}\times L_{z}=16D\times 5.5D\times 2.3D$, where $D$=0.5 m is the rotor diameter. Two different grids are employed: a coarse grid with the grid spacing $D/50$ and a fine grid with grid spacing $D/100$ around the turbine. The numbers of grid nodes are $N_{x}\times N_{y}\times N_{z}=802\times 277\times 116$ and $556\times 341\times 201$ with the corresponding time step $0.001D/U$ and $0.002D/U$ for the fine and coarse grids, respectively. For the coarse grid actuator surface cases, the number of triangular cells for actuator blades and nacelle is 912 and 1108, respectively. In the coarse grid actuator line cases, the number of segments is 78. For the find grid case, the number of triangular cells for actuator blades and nacelle is 1983 and 1636, respectively. On the coarse grid, four simulations are carried out using ASB-ASN, ASB, ALB and ALB-ASN models, respectively. On the fine grid only the ASB-ASN simulation is carried out. The smoothed discrete delta function is employed for the force distribution. In the coarse grid simulation, the $h_{x}$, $h_{y}$ and $h_{z}$ in Eq. (8) are equal to the grid spacing $D/50$. In order to ensure the same width of force distribution for the two different grids, the $h_{x}$, $h_{y}$ and $h_{z}$ are set as $D/50$ in the fine grid simulation. However, no significant difference is observed when $h_{x}$, $h_{y}$ and $h_{z}$ are equal to the grid spacing ($D/100$) for the fine grid simulation. The Reynolds number based on the bulk velocity and the rotor diameter is $2\times 10^{5}$. In the simulations, all the cases are first carried out until the total kinetic energy reaches a quasi-steady state, and subsequently the flowfields are averaged for $370T$ and $185T$ for the coarse grid and fine grid simulations, respectively, where $T=2\pi/\Omega$ is the rotor revolution period. The rotational speed of the rotor is $9.28$ $rad/s$, which gives a tip-speed ratio of 5.8 based on the bulk velocity (0.4 m/s).

We first examine the vertical profiles of the streamwise velocity and turbulence kinetic energy at different downstream locations. Comparison of the computed streamwise velocity profiles with the measurements is shown in Figs. 13 and 14. As expected, the ASB and ALB simulations without a nacelle model predict higher streamwise velocity near the hub at near wake locations ($1D$, $2D$ and $3D$). With the actuator surface model for nacelle, significant improvements are obtained on predicting the streamwise velocity near hub in the near wake for both actuator surface and actuator line blade models. For further downstream locations from $4D$ to $10D$, the predicted streamwise velocities, on the other hand, are nearly the same for all the simulations.

Turbulence kinetic energy at turbine’s far wake locations can be considered as the footprint of wake meandering. In order to evaluate the capability of the developed actuator surface model in capturing wake meandering in turbine’s far wake, we compare the computed turbulent kinetic energy ($k$) profiles with measurements in Figs. 15 and  16. At $1D$ turbine downstream, in the region near the bottom tip, the computed $k$ agrees well with the measurements for both ASB and ASB-ASN simulations. Near the top tip region, the computed $k$ from all the simulations is lower than that from the experiment, while the $k$ computed from the ASB-ASN on the fine grid, the ALB and ALB-ASN simulations shows a better agreement. The largest differences among different simulation are observed in the region near the hub: the $k$ from the simulation with the actuator surface nacelle model is lower than the measurements; the $k$ from the simulations without the nacelle model, on the other hand, is larger than the measurements. At $2D$ turbine downstream, the $k$ profiles computed from all the ASB-ASN and ASB simulations are in good agreement with measurements in the region near the top tip, which are overpredicted by the ALB-ASN and ALB models. The $k$ profiles predicted by the simulations with the actuator surface nacelle model are in good agreement with measurements in the near hub region, which are overpredicted by the simulations without the nacelle model. At $3D$ downstream from the turbine, the $k$ profiles computed from all the simulations are in good agreement with measurements. At $4D$, $5D$ and $6D$ downstream from the turbine, the $k$ profiles computed from the simulations with the actuator surface nacelle model are in good agreement with the measured values. The $k$ computed from the simulations without the actuator surface nacelle, on the other hand, is smaller than that from experiments at these locations. It is noticed that $3D$, $4D$ and $5D$ downstream from the turbine are the locations where maximum $k$ in the streamwise direction is observed. For further turbine downstream locations at $7D$, $8D$ and $9D$, the $k$ profiles computed from all the simulations show good agreements with measurements, although the $k$ computed from the simulations without a nacelle model is still slightly smaller than measurements.

In order to show the different wake meandering behavior for the simulations with and without a nacelle model more intuitively, we show in Fig. 17 the contours of instantaneous streamwise velocity at different time instants. It is seen that the jet from the hub of the turbine remains columnar until $3D$ turbine downstream for the ASB simulation. The inner wake from the ASB-ASN simulation, on the other hand, expands radially from $2D$ turbine downstream and interacts with the outer tip shear layer for about $1D$ to $2D$ distance.

To further examine the differences between the wake meandering characteristics predicted by the ASB and ASB-ASN simulations, the power spectral densities at different downstream locations are compared in Fig. 18. At $0.5D$ downstream from the turbine higher power spectral density is observed at low frequencies for the ASB-ASN simulation in comparison with ASB simulation (which is approximately ten times larger at the dominant frequency $fD/U=0.41$). The rotor frequency and the blade passing frequency are observed clearly on the PSD profile computed from the ASB simulation, which are not observed for the ASB-ASN simulation for which the inner wake behind the turbine is dominated by the wake from the nacelle. At $1D$ turbine downstream, only the rotor frequency is noticeable from the PSD computed in the ASB simulation. The peak frequency at around $fD/U=0.4$ observed at $0.5D$ turbine downstream in the ASB-ASN simulation does not exist anymore at $1D$ turbine downstream. While the PSD at low frequency from the ASB-ASN simulation is still significantly larger than that from the ASB simulation. At $2D$ turbine downstream, the differences of the PSD magnitudes at low frequencies between the ASB-ASN and ASB simulations becomes even larger in comparison with that at $0.5D$ and $1D$ turbine downstream. The wake meandering frequency of this turbine is $fD/U=0.28$ as reported in Chamorro et al.’s paper [37]. At $2D$, $3D$, $4D$, $5D$ and $6D$ turbine downstream locations, similar meandering frequency around 0.2 is well predicted by the ASB-ASN simulation. Such frequency, on the other hand, is observed at much further turbine downstream locations at $4D$, $5D$ and $6D$ for the ASB simulation. At $3D$ turbine downstream location, the PSD at low frequency is still significantly larger for the ASB-ASN simulation in comparison with the ASB simulation. It is also noticed that PSD at high frequency, on the other hand, is larger for the ASB simulation. At $4D$ and further downstream locations, the differences between ASB-ASN and ASB simulations become very small for the whole spectrum.

Fig. 17 has shown that the ASB-ASN model is able to predict the intensive interaction between the inner wake and outer shear layers, which plays an important role in triggering the wake meandering at far wake locations. To more clearly demonstrate the interaction between the inner wake and the outer wake, we plot in Fig. 19 the time-averaged spanwise vorticity for (a) and (b) from the actuator surface simulations, and (c) from the geometry-resolving simulation of Kang, Yang and Sotiropoulos [17]. As seen, the sign of the vorticity within the outer tip shear layers is opposite from that within the corresponding inner shear layers for the ASB simulation as shown in Fig. 17(a). The vortices within the inner region remain strong until $3D$ $\sim$ $4D$ turbine downstream, which prohibits the radially expansion of the outer wake to the wake center. On the other hand, the sign of the vorticity within the outer tip shear layers is the same as that within the corresponding inner shear layers for the ASB-ASN simulation and the geometry-resolving simulation as shown in Fig. 17(b) and (c), respectively. The vortices within the inner region gradually expands and joins with the vortices from the outer layers at approximately $2D$ turbine downstream, for which it happens somewhat earlier for the geometry-resolving simulation. Such interactions result in wider high-magnitude vorticity regions from approximately $2.5D$ to $6D$ turbine downstream for the simulation with the nacelle model, where the turbulent kinetic energy is underpredicted by the ASB simulation as shown in Fig. 15. Because of the different flows within the inner region, the radial expansions of the outer shear layers arrive the wake center at different downstream locations, which are approximately $3D$ (which is nearly the same for the geometry-resolving simulation as shown in (c)) and $4D$ for the simulations with and without a nacelle model. It is noticed that these two locations approximately coincide with the locations where the wake meandering frequency starts to be picked up for the ASB-ASN and ASB simulations, respectively.

5 Conclusions

In this paper, a new class of actuator surface models for wind and hydrokinetic turbines is proposed, which can take into account more geometrical features of turbine blade and nacelle as compared to the standard actuator line model. In the actuator surface model for blade, the blade is represented by the surface formed by foil chord lines at different radial locations. The forces on the blade are calculated using the blade element approach and are uniformly distributed in the chordwise direction. In the actuator surface model for nacelle, the nacelle is represented by the actual geometry of the nacelle. The normal component of the forces on the nacelle surface is calculated in the same way as in the direct forcing immersed boundary method. The magnitude of the tangential force is calculated using a friction coefficient and a reference velocity, in which the friction coefficient is computed using Schultz and Grunow’s formula [1], and the magnitude of the reference velocity is equal to the instantaneous streamwise incoming velocity upstream of the nacelle. The direction of the tangential force is taken the same as the local tangential velocity relative to the nacelle surface. It is noticed that the proposed actuator surface model for nacelle can be applied to other bluff bodies. Instead of using Schultz and Grunow’s formula, the friction coefficient may also be set to a constant value when the Reynolds number effect is insignificant or the friction drag is negligible in comparison with the form drag. For complex geometries when Schultz and Grunow’s relation fails, more accurate distributions of the friction coefficient need to be employed, which can be obtained from wind tunnel experiments or wall-resolved large-eddy simulations.

The actuator surface model for nacelle is evaluated by simulating the flow over periodic placed nacelles. The model predicted results are compared with that from the wall-resolved large-eddy simulation with the nacelle geometry represented by the sharp-interface immersed boundary method [26, 28]. The comparison shows relatively large discrepancies at near wake locations, where the wake recovery rate and the turbulence intensity computed from the model are lower than the wall-resolved large-eddy simulation predictions. At far wake locations, on the other hand, overall good agreement is obtained.

The capability of the actuator surface models in turbine wake predictions is evaluated by simulating the MEXICO turbine and the hydrokinetic turbine employed in [17]. Three different tip-speed ratios are considered in the MEXICO turbine simulations. The comparison of the computed results with the wind tunnel measurements shows that the proposed actuator surface models can accurately predict the strength of the tip vortices and the velocity in the wake of the MEXICO turbine. The capability of the actuator surface models for wake meandering predictions is examined using the hydrokinetic turbine case. The comparison of the computed results with measurements shows that the actuator surface models can accurately predict the velocity deficits in the near wake, the turbulence kinetic energy and the characteristic wake meandering frequency in the far wake. The importance of the actuator surface model for nacelle in the wake meandering prediction is examined by another simulation with only the actuator surface model for blade. It is shown that the simulation with only the actuator surface model for blade underpredicts the turbulence kinetic energy at approximately $4D\sim 7D$ turbine downstream locations and fails to predict the wake meandering frequency at $2D$ and $3D$ turbine downstream locations.

It is noticed that the present actuator surface model for blade is different from the actuator surface model developed by Shen, Zhang and Sørensen [8], in which the forces on the blade are computed from the pressure coefficients on blade surfaces and are non-uniformly distributed in the foil chordwise direction. Shen et al.’s model is more accurate in representing chordwise force distribution. To achieve this more accurate representation, it requires considerably higher spatial and temporal resolution, which may not be suitable for wind farm simulations. For the present actuator surface model, although the computational cost related to velocity interpolation (Eq. (7)) and force distribution (Eq. (18)) is significantly higher than that of actuator line simulations, the computational time added by actuator surface representation (Eqs. (7) and (18)) still only occupies a very small percentage of the total computational time. Using a background grid with resolution comparable to that in actuator line simulations, the present model provides an affordable way for field-scale wind farm simulations taking into account some geometrical effects of nacelle and blades.

\ack

This work was supported by Department of Energy DOE (DE-EE0002980, DE-EE0005482 and DE-AC04-94AL85000), Sandia National Laboratories and Xcel Energy through the Renewable Development Fund (grant RD4-13). Computational resources were provided by Sandia National Laboratories, National Renewable Energy Laboratory, and the University of Minnesota Supercomputing Institute.