15. Wind Energy Modeling

Wind energy analysis is the primary application area for the Nalu-Wind development team. This section describes the theoretical basis of Nalu-Wind from a wind energy perspective, using nomenclature familiar to wind energy experts and mapping it to Nalu-Wind concepts and nomenclature described in previous sections. Hopefully, this will provide an easier transition for users familiar with WRF and SOWFA to Nalu-Wind.

In order to evaluate the energy output and the structural loading on wind turbines, the code must model: 1. the incoming turbulent wind field across the entire wind farm, and 2. the evolution of turbine wakes in turbulent inflow conditions and their interaction with the downstream turbines. First, the governing equations with all the terms necessary to model a wind farm are presented with links to implementation and verification details elsewhere in the theory and/or verification manuals. A brief description of Nalu-Wind’s numerical discretization schemes is presented next. This is followed by a brief discussion of the boundary conditions used to model atmospheric boundary layer (ABL) flows with or without wind turbines (currently modeled as actuator sources within the flow domain).

Currently Nalu-Wind supports two types of wind simulations:

Precursor simulations

Precursor simulations are used in wind applications to generate time histories of turbulent ABL inflow profiles that are used as inlet conditions in subsequent wind farm simulations. The primary purpose of these simulations are to trigger turbulence generation and obtain velocity and temperature profiles that have converged to a statistic equilibrium.

Wind farm simulation with turbines as actuator sources

In this case, the wind turbine blades and tower are modeled as actuator source terms by coupling to the OpenFAST libraries. Velocity fields are sampled at the blade and tower control points within the Nalu-Wind domain and the blade positions and blade/tower loading is provided by OpenFAST to be used as source terms within the momentum equation.

15.1. Governing Equations

We begin with a review of the momentum and enthalpy conservation equations within the context of wind farm modeling [CLM+12]. Equation (15.1) shows the Favre-filtered momentum conservation equation (Eq. (2.1)) reproduced here with all the terms required to model a wind farm.

(15.1)$$\underset{\mathbf{I}}{\underbrace{\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{u}}_i\right)}}+\underset{\mathbf{II}}{\underbrace{\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_i{\tilde{u}}_j\right)}}=-\underset{\mathbf{III}}{\underbrace{\frac{\partial p^{\prime }}{\partial x_j}{\delta }_{ij}}}-\underset{\mathbf{IV}}{\underbrace{\frac{\partial {\tau }_{ij}}{\partial x_j}}}-\underset{\mathbf{V}}{\underbrace{2\overline{\rho }{ϵ}_{ijk}{\mathrm{\Omega }}_j{u}_k}}+\underset{\mathbf{VI}}{\underbrace{(\overline{\rho }-{\rho }_{\circ })g_i}}+\underset{\mathbf{VII}}{\underbrace{S_i^u}}+\underset{\mathbf{VIII}}{\underbrace{f_i^T}}$$

Term $\mathbf{I}$ represents the time rate of change of momentum (inertia);

Term $\mathbf{II}$ represents advection;

Term $\mathbf{III}$ represents the pressure gradient forces (deviation from hydrostatic and horizontal mean gradient);

Term $\mathbf{IV}$ represents stresses (both viscous and sub-filter scale (SFS)/Reynolds stresses);

Term $\mathbf{V}$ describes the influence Coriolis forces due to earth’s rotation – see Sec. Section 2.2.1;

Term $\mathbf{VI}$ describes the effects of buoyancy using the Boussinesq approximation – see Section 2.2.2;

Term $\mathbf{VII}$ represents the source term used to drive the flow to a horizontal mean velocity at desired height(s) – see Section 2.7; and

Term $\mathbf{VIII}$ is an optional term representing body forces when modeling turbine with actuator disk or line representations – see Section 15.6.

In wind energy applications, the energy conservation equation is often written in terms of the Favre-filtered potential temperature, $\theta$, equation, as shown below

(15.2)$$\frac{\partial }{\partial t}\left(\overline{\rho }\tilde{\theta }\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_j\tilde{\theta }\right)=-\frac{\partial }{\partial x_j}{\hat{q}}_j$$

where, ${\hat{q}}_j$ represents the temperature transport due to molecular and SFS turbulence effects. Due to the high Reynolds number associated with ABL flows, the molecular effects are neglected everywhere except near the terrain. Potential temperature is related to absolute temperature by the following equation

$$\theta =T{\left(\frac{\overline{p}}{p_{\circ }}\right)}^{-\left(\frac{R}{c_p}\right)}$$

Under the assumption of ideal gas conditions and constant $c_p$, the gradients in potential temperature are proportional to the gradients in absolute temperature, i.e.,

$$[\frac{\partial T}{\partial t},\frac{\partial T}{\partial x},\frac{\partial T}{\partial y}]={\left(\frac{\overline{p}}{p_{\circ }}\right)}^{\left(\frac{R}{c_p}\right)}[\frac{\partial \theta }{\partial t},\frac{\partial \theta }{\partial x},\frac{\partial \theta }{\partial y}]$$

Furthermore, ignoring the pressure and viscous work terms in Eq. (2.12) and assuming constant density (incompressible flow), it can be shown that solving the enthalpy equation is equivalent to solving the potential temperature equation. The enthalpy equation solved in wind energy problems is shown below

(15.3)$$\frac{\partial }{\partial t}\left(\overline{\rho }\tilde{T}\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_j\tilde{T}\right)=-\frac{\partial }{\partial x_j}{q}_j$$

It is noted here that the terms ${\hat{q}}_j$ (Eq. (15.2)) and $q_j$ (Eq. (15.3)) are not equivalent and must be scaled appropriately. User can still provide the appropriate initial and boundary conditions in terms of potential temperature field. Under these assumptions and conditions, the resulting solution can then be interpreted as the variation of potential temperature field in the computational domain.

15.2. Turbulence Modeling

LES turbulence closure is provided by the Subgrid-Scale Kinetic Energy One-Equation LES Model or the standard Smagorinsky model for wind farm applications.

15.3. Numerical Discretization & Stabilization

Nalu-Wind provides two discretization approaches

Control Volume Finite Element Method (CVFEM)

Nalu-Wind uses a dual mesh approach (see Section 3.1) where the control volumes are constructed around the nodes of the finite elements within the mesh – see 15.1. The equations are solved at the integration points on the sub-control surfaces and/or the sub-control volumes.

Edge-Based Vertex Centered Scheme

The edge-based scheme is similar to the finite-volume approach used in SOWFA with the nodes at the cell center of the dual mesh.

15.1 Schematic of HEX-8 mesh showing the finite elements, nodes, and the associated control volume dual mesh.

The numerical discretization approach is covered in great detail in Section 3, the advection and pressure stabilization approaches are documented in Section 4 and Section 5 respectively. Users are strongly urged to read those sections to gain a thorough understanding of the discretization scheme and its impact on the simulations.

15.4. Time stepping scheme

The time stepping method in Nalu-Wind is described in the Fuego theory manual [Tea16] for the backward Euler time discretization. The implementation details of the BDF2 time stepping scheme used in Nalu-Wind is described here. The Navier-Stokes equations are written as

(15.4)$$\begin{array}{rl}{\mathbf{F}}_i({\rho }^{n+1},u_i^{n+1},P^{n+1})-\int {\frac{\partial \rho u_i}{\partial t}|}^{n+1}\mathrm{d}V& =0,\\ {\mathbf{F}}_i({\rho }^{n+1},u_i^{n+1},P^{n+1})-\frac{({\gamma }_1{\rho }^{n+1}{u_i}^{n+1}+{\gamma }_2{\rho }^n{u_i}^n+{\gamma }_3{\rho }^n{u_i}^{n-1})}{\mathrm{\Delta }t}\mathrm{\Delta }V& =0,\end{array}$$

where

$$\begin{array}{rl}{\mathbf{F}}_i({\rho }^{n+1}{u}_i^{n+1})& =-\int {\rho }^{n+1}{u}_i^{n+1}{u}_j^{n+1}{n}_j\mathrm{d}S+\int {\tau }_{ij}^{n+1}{n}_j\mathrm{d}S-\int P^{n+1}{n}_i\mathrm{d}S-\int ({\rho }^{n+1}-{\rho }_{\circ })g_i\mathrm{d}V,\\ & =-\int u_i^{n+1}{\dot{m}}^{n+1}+\int {\tau }_{ij}^{n+1}{n}_j\mathrm{d}S-\int P^{n+1}{n}_i\mathrm{d}S-\int ({\rho }^{n+1}-{\rho }_{\circ })g_i\mathrm{d}V.\end{array}$$

and ${\gamma }_i$ are factors for BDF2 time discretization scheme (see Section 13). As is typical of incompressible flow solvers, the mass flow rate through the sub-control surfaces is tracked independent of the velocity to maintain conservation of mass. The following conventions are used:

$$\begin{array}{rl}{\varphi }^{\ast }& =\text{ Predicted value of }\varphi \text{ at }n+1\text{ time step before linear solve}\\ \hat{\varphi }={\varphi }^{\ast \ast }& =\text{ Predicted value of }\varphi \text{ at }n+1\text{ time step after linear solve}\end{array}$$

The Newton’s method is used along with a linearization procedure to predict a solution to the Navier-Stokes equations at time step $n+1$ as

(15.5)$$\begin{array}{rl}{\mathbf{A}}_{ij}\delta u_j& ={\mathbf{F}}_i^{\ast }-\frac{({\gamma }_1{\rho }^{\ast }{u_i}^{\ast }+{\gamma }_2{\rho }^n{u_i}^n+{\gamma }_3{\rho }^n{u_i}^{n-1})}{\mathrm{\Delta }t}\mathrm{\Delta }V,\\ \text{where }\delta u_j& =u_i^{\ast \ast }-u_i^{\ast },\\ {\mathbf{A}}_{ij}& =(\frac{{\gamma }_1{\rho }^{\ast }}{\mathrm{\Delta }t}\mathrm{\Delta }V{\delta }_{ij}-{\frac{\partial F_i}{\partial u_j}|}^{\ast }),\\ \text{and }{\mathbf{F}}_i^{\ast }& =-\int u_i^{\ast }{\dot{m}}^{\ast }+\int {\tau }_{ij}^{\ast }{n}_j\mathrm{d}S-\int P^{\ast }{n}_i\mathrm{d}S-\int ({\rho }^{\ast }-{\rho }_{\circ })g_i\mathrm{d}V.\end{array}$$

After each Newton or outer iteration, ${\varphi }^{\ast \ast }$ is a better approximation to ${\varphi }^{n+1}$ compared to ${\varphi }^{\ast }$. $\rho \ast$ and ${\dot{m}}^{\ast }$ are retained constant through each outer iteration. $\mathbf{F}({\rho }^{\ast }{u}_i^{\ast \ast })$ is linear in $u_i^{\ast }$ and hence

(15.6)$${\mathbf{F}}_i^{\ast }={\frac{\partial F_i}{\partial u_j}|}^{\ast }{u}_j^{\ast }-\int P^{\ast }{n}_i\mathrm{d}S-\int ({\rho }^{\ast }-{\rho }_{\circ })g_i\mathrm{d}V$$

Applying Eq. (15.6) to Eq. (15.5), we get the linearized momentum predictor equation solved in Nalu-Wind.

(15.7)$$\begin{array}{rl}{\mathbf{A}}_{ij}\delta u_j& ={\frac{\partial F_i}{\partial u_j}|}^{\ast }{u}_j^{\ast }-\int P^{\ast }{n}_i\mathrm{d}S-\int ({\rho }^{\ast }-{\rho }_{\circ })g_i\mathrm{d}V\\ & -\frac{({\gamma }_1{\rho }^{\ast }{u_i}^{\ast }+{\gamma }_2{\rho }^n{u_i}^n+{\gamma }_3{\rho }^{n-1}{u_i}^{n-1})}{\mathrm{\Delta }t}\mathrm{\Delta }V\\ {\mathbf{A}}_{ij}\delta u_j& =(\frac{{\gamma }_1{\rho }^{\ast }}{\mathrm{\Delta }t}\mathrm{\Delta }V{\delta }_{ij}-{\frac{\partial F_i}{\partial u_j}|}^{\ast }){u_j}^{\ast }-\int P^{\ast }{n}_i\mathrm{d}S-\int ({\rho }^{\ast }-{\rho }_{\circ })g_i\mathrm{d}V\\ & -\frac{({\gamma }_2{\rho }^n{u_i}^n+{\gamma }_3{\rho }^{n-1}{u_i}^{n-1})}{\mathrm{\Delta }t}\mathrm{\Delta }V\\ {\mathbf{A}}_{ij}\delta u_j& ={\mathbf{A}}_{ij}{u}_j^{\ast }-\int P^{\ast }{n}_i\mathrm{d}S-\int ({\rho }^{\ast }-{\rho }_{\circ })g_i\mathrm{d}V\\ & -\frac{({\gamma }_2{\rho }^n{u_i}^n+{\gamma }_3{\rho }^{n-1}{u_i}^{n-1})}{\mathrm{\Delta }t}\mathrm{\Delta }V\end{array}$$

$u_i^{\ast \ast }$ will not satisfy the continuity equation. A correction step is performed at the end of each outer iteration to make $u_i^{\ast \ast }$ satisfy the continuity equation as

$$\begin{array}{rl}{u}_i^{n+1}& =u_i^{\ast \ast }-\frac{{\tau }_3}{\rho }\mathbf{G}\mathrm{\Delta }{P}^{\ast \ast }\\ \text{where }\mathrm{\Delta }{P}^{\ast \ast }& =P^{\ast \ast }-P^{\ast }\end{array}$$

As described in Section 12.1, the continuity equation to be satisfied along with the splitting and stabilization errors is

(15.8)$$\mathbf{D}\rho u^{\ast \ast }=b+({\mathbf{L}}_{\mathbf{1}}-\mathbf{D}{\tau }_3\mathbf{G})\mathrm{\Delta }{P}^{\ast \ast }+({\mathbf{L}}_{\mathbf{2}}-\mathbf{D}{\tau }_2\mathbf{G})P^{\ast }$$

where $b$ contains any source terms when the velocity field is not divergence free and the other terms are the errors due to pressure stabilization as shown by Domino [Dom06]. The final pressure Poisson equation solved to enforce continuity at each outer iteration is

(15.9)$$\begin{array}{rl}{u}^{n+1}& =u^{\ast \ast }-\frac{{\tau }_3}{\rho }\mathbf{G}\mathrm{\Delta }{P}^{\ast \ast }\\ b+({\mathbf{L}}_{\mathbf{1}}-\mathbf{D}{\tau }_3\mathbf{G})\mathrm{\Delta }{P}^{\ast \ast }& +({\mathbf{L}}_{\mathbf{2}}-\mathbf{D}{\tau }_2\mathbf{G})P^{\ast }\\ & =\mathbf{D}(\rho u^{n+1})=\mathbf{D}(\rho \hat{u})-\mathbf{D}({\tau }_3\mathbf{G}\mathrm{\Delta }{P}^{\ast \ast })\\ b+{\mathbf{L}}_{\mathbf{1}}\mathrm{\Delta }{P}^{\ast \ast }& =\mathbf{D}(\rho \hat{u})-({\mathbf{L}}_{\mathbf{2}}-\mathbf{D}{\tau }_2\mathbf{G})P^{\ast }\\ -{\mathbf{L}}_{\mathbf{1}}\mathrm{\Delta }{P}^{\ast \ast }& =\mathbf{D}\rho \hat{u}+\mathbf{D}{\tau }_2\mathbf{G}{P}^{\ast }-{\mathbf{L}}_{\mathbf{2}}{P}^{\ast }\\ -{\mathbf{L}}_{\mathbf{1}}\mathrm{\Delta }{P}^{\ast \ast }& =-\mathbf{D}\rho \hat{u}-\mathbf{D}{\tau }_2\mathbf{G}{P}^{\ast }+{\mathbf{L}}_{\mathbf{2}}{P}^{\ast }+b\end{array}$$

Thus, the final set of equations solved at each outer iteration is

$$\begin{array}{rl}{\mathbf{A}}_{ij}\delta u_j& ={\mathbf{A}}_{ij}{u}_j^{\ast }-\int P^{\ast }{n}_i\mathrm{d}S-\int ({\rho }^{\ast }-{\rho }_{\circ })g_i\mathrm{d}V\\ & -\frac{({\gamma }_2{\rho }^n{u_i}^n+{\gamma }_3{\rho }^{n-1}{u_i}^{n-1})}{\mathrm{\Delta }t}\mathrm{\Delta }V\\ -{\mathbf{L}}_{\mathbf{1}}\mathrm{\Delta }{P}^{\ast \ast }& =-\mathbf{D}\rho \hat{u}-\mathbf{D}{\tau }_2\mathbf{G}{P}^{\ast }+{\mathbf{L}}_{\mathbf{2}}{P}^{\ast }+b\\ u_i^{n+1}& =u_i^{\ast \ast }-\frac{{\tau }_3}{\rho }\mathbf{G}\mathrm{\Delta }{P}^{\ast \ast }\end{array}$$

15.4.1. Approximations for the Schur complement

Nalu-Wind implements two options for approximating the Schur complement for the split velocity-pressure solution of the incompressible momentum and continuity equation. The two options are:

$$\begin{array}{rl}\tau ={\tau }_1={\tau }_2={\tau }_3& =\mathrm{\Delta }t\mathrm{Original}\mathrm{implementation}\\ & =(A_{ii}{)}^{-1}\mathrm{Alternate}\mathrm{algorithm}\end{array}$$

where $A_{ii}$ is the diagonal entry of the momentum linear system. The latter option is similar to the SIMPLE and PIMPLE implementations in OpenFOAM and is used for simulations with RANS and hybrid RANS-LES models with large Courant numbers.

15.4.2. Underrelaxation for momentum and scalar transport

By default, Nalu-Wind applies no underrelaxation during the solution of the Navier-Stokes equations. However, in RANS simulations at large timesteps some underrelaxation might be necessary to restore the diagonal dominance of the transport equations. User has the option to specify underrelaxation through the input files. When underrelaxation is applied, the advection and diffusion contributions to the diagonal term are modifed by dividing these terms by the underrelaxation factor. It must be noted that the underrelaxation is only applied to the advective and viscous contributions in the diagonal term and not the time derivative term.

$$A_{ii}=-\frac{\sum _{i\ne j}{A}_{ij}}{\omega }+\frac{{\gamma }_1\rho \mathrm{\Delta }V}{\mathrm{\Delta }T}$$

The pressure update can also be underrelaxed by specifying the appropriate relaxation factor in the input file. When this option is activated, the full pressure update, in a given Picard iteration step, is used to project the velocity and mass flow rate and the relaxation is applied to the pressure solution at the end of the Picard iteration.

15.5. Initial & Boundary Conditions

This section briefly describes the boundary conditions available in Nalu-Wind for modeling wind farm problems. The terrain and top boundary conditions are described first as they are common to precusor and wind farm simulations.

15.5.1. Initial conditions

Nalu-Wind has the ability to initialize the internal flow fields to uniform conditions for all pressure, velocity, temperature, and TKE ($k$) in the input file <initial_conditions.constant>. Nalu-Wind also provides a user function to add perturbations to the velocity field to trigger turbulence generation during precursor simulations. To specify more complex flow field conditions, a temperature profile with a capping inversion for example, users are referred to pre-processing utilities available in NaluWindUtils library.

15.5.2. Terrain (Wall) boundary condition

Users are referred to Section 9.3 for the treatment of the terrain BC using roughness models. For enthalpy, users can provide a surface heat flux for modeling stratified flows.

15.5.3. Top boundary condition

For problems with minimal streamline curvature near the upper boundary (e.g. nearly flat terrain, negligible turbine blockage), a symmetry BC (slip wall) can be when modeling wind farm problems. By default a zero vertical temperature gradient will be imposed for the enthalpy equation when the symmetry boundary condition is used. If a non-zero normal temperature gradient is required to drive the flow to a desired temperaure profile, e.g., a capping inversion, then the abltop BC can be used. In this case the user_data input normal_temperature_gradient: value will set the normal temperature gradient to value at the top boundary.

For cases with significant terrain features or significant turbine blockage, the abltop BC can also be used to achieve an open boundary that allows for both inflows and outflows at the domain top. See the abltop BC documentation for details.

15.5.4. Inlet conditions

Time histories of inflow velocity and temperaure profiles can be provided as inputs (via I/O transfer) to drive the wind farm simulation with the desired flow conditions. See Section 8.3 for more details on this capability. Driving a wind farm simulation using velocity and temperature fields from a mesoscale (WRF) simulation would require an additional pre-processing steps with the wrftonalu utility.

15.5.5. Outlet conditions

See the description of open BC for detailed description of the outlet BC implementation. For wind energy problems, it is necessary to activate the global mass correction as a single value of pressure across the boundary layer is not apprpriate in the presence of buoyancy effects. It might also be necessary to fix the reference pressure at an interior node in order to ensure that the Pressure Poisson solver is well conditioned.

15.6. Wind Turbine Modeling

Wind turbine rotor and tower aerodynamic effects are modeled using actuator source representations. Compared to resolving the geometry of the turbine, actuator modeling alleviates the need for a complex body-fitted meshes, can relax time step restrictions, and eliminates the need for turbulence modeling at the turbine surfaces. This comes at the expense of a loss of fine-scale detail, for example, the boundary layers of the wind turbine surfaces are not resolved. However, actuator methods well represent wind turbine wakes in the mid to far downstream regions where wake interactions are important.

Actuator methods usually fall within the classes of disks, lines, surface, or some blend between the disk and line (i.e., the swept actuator line). Most commonly, the force over the actuator is computed, and then applied as a body-force source term, $f_i$ (Term $\mathbf{VIII}$), to the Favre-filtered momentum equation (Eq. (15.1)).

The body-force term $f_i$ is volumetric and is a force per unit volume. The actuator forces, $F_i^{\prime }$, are not volumetric. They exist along lines or on surfaces and are force per unit length or area. Therefore, a projection function, $g$, is used to project the actuator forces into the fluid volume as volumetric forces. A simple and commonly used projection function is a uniform Gaussian as proposed by Sorensen and Shen [SorensenS02],

$$g(\overrightarrow{r})=\frac{1}{{\pi }^{3/2}{ϵ}^3}{e}^{-{\left(\left|\overrightarrow{r}\right|/ϵ\right)}^2},$$

where $\overrightarrow{r}$ is the position vector between the fluid point of interest to a particular point on the actuator, and $ϵ$ is the width of the Gaussian, that determines how diluted the body force become. As an example, for an actuator line extending from $l=0$ to $L$, the body force at point $(x,y,z)$ due to the line is given by

(15.10)$$f_i(x,y,z)={\int }_0^{L}g\left(\overrightarrow{r}\left(l\right)\right)F_i^{\prime }\left(l\right)\text{d}l.$$

Here, the projection function’s position vector is a function of position on the actuator line. The part of the line nearest to the point in the fluid at $(x,y,z)$ has most weight.

The force along an actuator line or over an actuator disk is often computed using blade element theory, where it is convenient to discretize the actuator into a set of elements. For example, with the actuator line, the line is broken into discrete line segments, and the force at the center of each element, $F_i^k$, is computed. Here, $k$ is the actuator element index. These actuator points are independent of the fluid mesh. The point forces are then projected onto the fluid mesh using the Gaussian projection function, $g(\overrightarrow{r})$, as described above. This is convenient because the integral given in Equation (15.10) can become the summation

(15.11)$$f_i(x,y,z)=\sum _{k=0}^{N}g({\overrightarrow{r}}^k)F_i^k.$$

This summation well approximates the integral given in Equation (15.10) so long as the ratio of actuator element size to projection function width $ϵ$ does not exceed a certain threshold.

Presently, Nalu-Wind uses an actuator line representation to model the effects of turbine on the flow field; however, the class hierarchy is designed with the potential to add other actuator source terms such as actuator disk, swept actuator line and actuator surface capability in the future. The ActuatorLineFAST class couples Nalu-Wind with NREL’s OpenFAST for actuator line simulations of wind turbines. OpenFAST is a aero-hydro-servo-elastic tool to model wind turbine developed by the National Renewable Energy Laboratory (NREL). The ActuatorLineFAST class allows Nalu-Wind to interface as an inflow module to OpenFAST by supplying the velocity field information.

15.6.1. Nalu-Wind – OpenFAST Coupling Algorithm

A nacelle model is implemented using a Gaussian drag body force. The model implements a drag force in a direction opposite to velocity field at the center of the Gaussian. The width of the Gaussian kernel is determined using the reference area and drag coefficient of the nacelle as shown by Martinez-Tossas. [MartinezT17]

$${ϵ}_d=\sqrt{2c_{d}A/\pi }$$

where $c_d$ is the drag coefficient, and $A$ is the reference area. This value of ${ϵ}_d$ ensures that the momentum thickness of the generated wake is of the right size. The velocity sampled at the center of the Gaussian is corrected to obtain the upstream velocity.

$${{\tilde{u}}_i}_{\mathrm{\infty }}=\frac{1}{1-c_{d}A/(4\pi {ϵ}_d^2)}{{\tilde{u}}_i}_c$$

where ${\widetilde{u_i}}_c$ is the velocity at the center of the Gaussian and ${{\tilde{u}}_i}_{\mathrm{\infty }}$ is the free-stream velocity used to compute the drag force. The drag body force is then

$$f_d(x,y,z)=\frac{1}{2}{c}_{d}A{{\tilde{u}}_i}_{\mathrm{\infty }}^2\frac{1}{{\pi }^{3/2}{ϵ}_d^3}{e}^{-|\overrightarrow{r}{|}^2/{ϵ}_d^2}$$

where $\overrightarrow{r}$ is the position vector between the fluid point of interest and the center of the Gaussian force.

The actuator line implementation allows for flexible blades that are not necessarily straight (pre-bend and sweep). The current implementation requires a fixed time step when coupled to OpenFAST, but allows the time step in Nalu-Wind to be an integral multiple of the OpenFAST time step. At present, a simple time lagged FSI model is used to interface Nalu-Wind with the turbine model in OpenFAST:

• The velocity at time step at time step $n$ is sampled at the actuator points and sent to OpenFAST,

• OpenFAST advances the turbines up-to the next Nalu-Wind time step $n+1$,

• The body forces at the actuator points are converted to the source terms of the momentum equation to advance Nalu-Wind to the next time step $n+1$.

This FSI algorithm is expected to be only first order accurate in time. We are currently working on improving the FSI coupling scheme to be second order accurate in time.

15.6.2. Nalu-Wind – Actuator Disk Model via OpenFAST

An actuator disk model is implemented in Nalu-Wind by using an OpenFAST actuator line to sample the flow and compute the forcing. The actuator line is held stationary which leads to computational savings during execution because there is only 1 search operation in the initial setup.

The forces are gathered at each actuator line point, and the total force at each discrete radial location ($r_j$ where $j\in [1,N_R]$) is computed using (15.12).

(15.12)$${\mathbf{F}}_{total}(r_j)=\sum _{i=1}^{N_B}\mathbf{F}(r_j,{\theta }_i)$$

where $N_B$ and $N_R$ are the number of blades and number of radial points respectively.

${\mathbf{F}}_{total}(r_j)$ is then spread evenly across the original actuator line points and additional ‘swept-points’ that are added in between the actuator lines. The swept-points are always uniformly distributed azimuthally, but the number of swept points can either be non-uniformly or uniformly distributed along the radial direction (left and right images in figure 15.2). The non-uniform distribution uses the distance between points along the embedded actuator line blades as the arc-length between points in the azimuthal direction. This is the default behavior. If uniform spacing is desired then num_swept_pts must be specified in the input deck. This is the number of points between the actuator lines, so in figure 15.2 the num_swept_pts is 3.

15.2 Actuator Disk with non-uniform (left) and uniform (right) sampling in the azimuthal direction.

The force that is spread across all the points at a given radius is then calculated as (15.13).

(15.13)$$\mathbf{f}(r_j)=\frac{{\mathbf{F}}_{total}(r_j)}{N_B\ast (N_{S,j}+1)}$$

where $N_{S,j}$ is the number of swept points for a given radius. The index j is used because this value varies between radii when non-uniform sampling is applied.