2 The ALM in OpenFOAM

This section briefly describes the underlying ALM methodology as implemented in an atmospheric boundary-layer (ABL) solver in OpenFOAM (Churchfield et al. Reference Churchfield, Moriarty, Vijayakumar and Brasseur2010). The rotor blade or wing is discretized into a finite number (typically between 25 and 40) of actuator elements. Sectional lift and drag forces are computed at the centre of each of these actuator elements by interpolation of local flow-field data and using a table lookup for aerofoil properties. The computed blade forces at discrete actuator points are then projected onto the background Cartesian grid as spherical Gaussian body forces centred around each actuator element. The last term, $(F_{p})_{i}$ , in the following momentum equation corresponds to the body force:

(2.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{u}_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\overline{u}_{j}\overline{u}_{i})}{\unicode[STIX]{x2202}x_{j}}=-2\unicode[STIX]{x1D700}_{ijk}\unicode[STIX]{x1D6FA}_{j}\overline{u}_{k}-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}x_{i}}-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}-\frac{\unicode[STIX]{x2202}R_{ij}^{D}}{\unicode[STIX]{x2202}x_{j}}+\left(\frac{\unicode[STIX]{x1D70C}_{b}}{\unicode[STIX]{x1D70C}_{0}}-1\right)g_{i}+\frac{1}{\unicode[STIX]{x1D70C}_{0}}(F_{p})_{i}.\end{eqnarray}$$

In the present case, the momentum equation is for an ABL-LES, where $\unicode[STIX]{x2202}\overline{u}_{i}/\unicode[STIX]{x2202}t$ is the time rate of change of filtered velocity, $\unicode[STIX]{x2202}(\overline{u}_{j}\overline{u}_{i})/\unicode[STIX]{x2202}x_{j}$ is the convection term and $2\unicode[STIX]{x1D700}_{ijk}\unicode[STIX]{x1D6FA}_{j}\overline{u}_{k}$ is the Coriolis term due to planetary rotation. Furthermore, $\unicode[STIX]{x2202}\tilde{p}/\unicode[STIX]{x2202}x_{i}$ is the gradient of the modified pressure variable $\tilde{p}$ , which is the deviation in resolved-scale static pressure from its time-averaged value and lumped with one-third of the trace of the stress tensor $\unicode[STIX]{x1D70F}_{kk}/3$ . The driving pressure gradient is $\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x_{i}$ , while $\unicode[STIX]{x2202}R_{ij}^{D}/\unicode[STIX]{x2202}x_{j}$ is the divergence of the deviatoric part of the fluid stress tensor. The term $(\unicode[STIX]{x1D70C}_{b}/\unicode[STIX]{x1D70C}_{0}-1)g_{i}$ models buoyancy effects using the Boussinesq approximation, where $g_{i}$ represents gravity and $\unicode[STIX]{x1D70C}_{0}$ is the reference fluid density. OpenFOAM uses finite-volume discretization, with second-order-accurate spatial central differencing and a second-order-accurate Crank–Nicolson method for time integration. A pressure-implicit splitting operation (PISO) algorithm is used for the pressure, while the Smagorinsky model accounts for subgrid-scale effects.

The transformation of an elemental blade load $\boldsymbol{f}_{N,m}$ into a volumetric body force $\boldsymbol{F}_{p}$ is typically achieved by a Gaussian projection function, shown in vector form as

(2.2) $$\begin{eqnarray}\boldsymbol{F}_{p}(x_{P},y_{P},z_{P})=-\mathop{\sum }_{N}\mathop{\sum }_{m}\boldsymbol{f}_{N,m}(x_{N,m},y_{N,m},z_{N,m})\unicode[STIX]{x1D702}_{N,m},\end{eqnarray}$$

where the projection function is

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{N,m}=\frac{1}{\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x03C0}^{3/2}}\exp [-(|\boldsymbol{r}|/\unicode[STIX]{x1D716})^{2}].\end{eqnarray}$$

Here, $N$ is the blade index, $m$ is the actuator point index, $|\boldsymbol{r}|=|\boldsymbol{r}_{p}-\boldsymbol{r}_{N,m}|$ is the distance from grid cell $\boldsymbol{r}_{p}=(x_{P},y_{P},z_{P})^{\text{T}}$ to the actuator point $\boldsymbol{r}_{N,m}=(x_{N,m},y_{N,m},z_{N,m})^{\text{T}}$ and $\unicode[STIX]{x1D716}$ is a Gaussian projection radius. The elemental blade load is a function of the sectional lift and drag forces according to

(2.4) $$\begin{eqnarray}\boldsymbol{f}_{N,m}(x_{N,m},y_{N,m},z_{N,m})={\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{0}V_{rel}^{2}c~(c_{l}~\boldsymbol{e}_{l}+c_{d}\boldsymbol{e}_{d})\unicode[STIX]{x1D6E5}_{b},\end{eqnarray}$$

where $c$ is the local chord length, $c_{l}$ and $c_{d}$ are sectional aerofoil lift and drag coefficients, $\boldsymbol{e}_{l}$ and $\boldsymbol{e}_{d}$ are unit vectors perpendicular and parallel to the sectional relative velocity vector $\boldsymbol{V}_{rel}$ and $\unicode[STIX]{x1D6E5}_{b}$ is the width of the sectional actuator element. It should be noted that $c_{l}$ and $c_{d}$ are determined from table lookup using the local angle of attack $\unicode[STIX]{x1D6FC}$ , which is the difference between the directions of $\boldsymbol{V}_{rel}$ and the sectional chord line (Sørensen & Shen Reference Sørensen and Shen2002; Calaf, Meneveau & Meyers Reference Calaf, Meneveau and Meyers2010; Lu & Porté-Agel Reference Lu and Porté-Agel2011; Meyers & Meneveau Reference Meyers and Meneveau2012; Chatelain et al. Reference Chatelain, Backaert, Winckelmans and Kern2013). It is recommended that the Gaussian projection radius $\unicode[STIX]{x1D716}$ be based on the grid spacing $\unicode[STIX]{x1D6E5}_{grid}$ (Troldborg et al. Reference Troldborg, Sørensen and Mikkelsen2010) as

(2.5) $$\begin{eqnarray}\text{ALM:}\quad \unicode[STIX]{x1D716}/\unicode[STIX]{x1D6E5}_{grid}=2.\end{eqnarray}$$

Blade tip loads are typically overpredicted, if no appropriate tip-loss correction (Shen et al. Reference Shen, Mikkelsen, Sørensen and Bak2005) is used; see § 4. Recent work by the authors was driven by the hypothesis that the Gaussian projection radius $\unicode[STIX]{x1D716}$ is related to the primary mode of the blade loading within an ALM concept. These modifications to the ALM, in the following referred to as ALM*, led to improved prediction of integrated power and sectional blade loads (Jha et al. Reference Jha, Churchfield, Moriarty and Schmitz2014) and also captured smaller-scale turbulence-induced unsteadiness of sectional blade loads (Jha & Schmitz Reference Jha and Schmitz2016). This was accomplished using a variable Gaussian radius $\unicode[STIX]{x1D716}(r)$ , which was determined using an equivalent elliptic planform $c^{\ast }(r)$ of the same aspect ratio, $AR$ , as the actual blade as

(2.6) $$\begin{eqnarray}\text{ALM*:}\quad \unicode[STIX]{x1D716}/c^{\ast }=C_{1}(\unicode[STIX]{x03C0}AR)=\unicode[STIX]{x1D716}_{0}/c_{0},\end{eqnarray}$$

where $C_{1}=0.25(n_{max}\unicode[STIX]{x1D6E5}_{grid}/R)$ . Furthermore, $\unicode[STIX]{x1D716}_{0}$ and $c_{0}$ are respective values at the midspan of the blade or wing, with $\unicode[STIX]{x1D716}_{0}=n_{max}\unicode[STIX]{x1D6E5}_{grid}$ and $n_{max}$ being determined for a given $\unicode[STIX]{x1D6E5}_{grid}$ through $C_{1}\in [0.025,0.030]$ . While this was a step towards more accurate modelling of rotor blades keeping $\unicode[STIX]{x1D6E5}_{b}/\unicode[STIX]{x1D6E5}_{grid}\geqslant 1.5$ , the projection function (equation (2.3)) remained a spherical Gaussian that overlapped with adjacent body-force projections; see figure 1. While the elliptic $\unicode[STIX]{x1D716}$ -based ALM* is more accurate compared with the standard ALM and is consistent with grid refinements, the overlap is not fully eradicated. A different approach was recently taken by Churchfield et al. (Reference Churchfield, Schreck, Martinez-Tossas, Meneveau and Spalart2017) and Martinez-Tossas et al. (Reference Martinez-Tossas, Churchfield and Meneveau2017), where an optimized Gaussian spreading width becomes a function of the blade chord distribution, an approach that can be considered to be a generalization of earlier works (Watters & Masson Reference Watters and Masson2010; Shives & Crawford Reference Shives and Crawford2013), although it requires considerably finer grids and associated higher computational cost.

Figure 1.

Volumetric body-force projection (ALM, ALM $^{\ast }$ versus ACE).

Figure 2.

Three-dimensional (ALM $^{\ast }$ ) versus two-dimensional (ACE) Gaussian kernel function.

3 Actuator curve embedding in OpenFOAM

The conceptual idea of ACE is to remove the inconsistency in force overlap inherent to standard ALM and ASM approaches. This can be achieved if (i) the kernel function is chosen in such a way that the projection of blade forces is circumferential only in the region normal to the local blade section, (ii) each grid point has a unique body force associated with a specific section of an actuator element and (iii) the volume integral of the total projected blade force is conserved. It should be noted that the force overlap is unphysical and can only be controlled by an appropriate spacing between actuator points (Jha et al. Reference Jha, Churchfield, Moriarty and Schmitz2014), and is inconsistent for curved lifting lines (or aeroelastically deforming blades) due to a dynamic overlap region in convex or concave curvature.

The ACE concept is inspired by the vorticity-embedding methodology due to Steinhoff & Ramachandran (Reference Steinhoff and Ramachandran1988) and subsequent implementation into a Eulerian/Lagrangian free-wake solver for hovering rotorcraft (Schmitz et al. Reference Schmitz, Bhagwat, Moulton, Caradonna and Chattot2009). As opposed to sequentially finding a sphere of influence around each actuator point, a reverse approach is followed in ACE. Here, the flow field experiences the influence of a blade/wing as one entity. This is realized by associating each grid cell with a unique actuator element along the lifting curve, rather than being influenced by multiple actuator points; see figure 2.

Figure 3.

Actuator curve embedding – computation of $p_{n}$ and $p_{s}$ at each $(x_{P},y_{P},z_{P})$ .

3.1 Actuator curve embedding – geometric parameterization

For a given grid point $\boldsymbol{r}_{p}=(x_{P},y_{P},z_{P})^{\text{T}}$ , each actuator-element pair $m$ and $m+1$ has unique values for the tangential coordinate $p_{s}$ along the respective actuator element and normal distance $p_{n}$ to the actuator element, as shown in figure 3. The tangential coordinate $p_{s}$ is normalized by the width of the actuator element, $\unicode[STIX]{x1D6E5}_{b}=\Vert \boldsymbol{r}_{N,m+1}-\boldsymbol{r}_{N,m}\Vert$ , between $\boldsymbol{r}_{N,m}$ and $\boldsymbol{r}_{N,m+1}$ , and is determined by

(3.1) $$\begin{eqnarray}p_{s}=\frac{(\boldsymbol{r}_{p}-\boldsymbol{r}_{N,m}).(\boldsymbol{r}_{N,m+1}-\boldsymbol{r}_{N,m})}{\unicode[STIX]{x1D6E5}_{b}^{2}}.\end{eqnarray}$$

A given grid point $\boldsymbol{r}_{p}$ is associated with the pair of actuator points $m$ and $m+1$ for which $p_{s}\in [0,1]$ . This simple criterion prevents any body-force spreading beyond the edges of the lifting line, see (G3) and (G4) in figure 3, and associates a given grid cell $\boldsymbol{r}_{p}$ with a specific actuator segment, or results in a grid cell not being affected by body-force projection. The corresponding normal distance, $p_{n}$ , is then expressed as

(3.2) $$\begin{eqnarray}p_{n}=\frac{\Vert (\boldsymbol{r}_{p}-\boldsymbol{r}_{N,m})\times (\boldsymbol{r}_{p}-\boldsymbol{r}_{N,m+1})\Vert }{\unicode[STIX]{x1D6E5}_{b}}.\end{eqnarray}$$

It can be noted from figure 3 that two different grid cells, (G1) and (G2), can be associated with the same pair of actuator points, $\boldsymbol{r}_{N,m}$ and $\boldsymbol{r}_{N,m+1}$ , but still have unique values for $p_{s}$ and $p_{n}$ . For a lifting line $N$ with curvature, special cases may arise for grid point $\boldsymbol{r}_{p}$ facing either the convex (G5) or concave (G6) side,

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \text{convex:}\quad (p_{s_{m}}>1)\wedge (p_{s_{m+1}}<0)\rightarrow p_{n}=\Vert \boldsymbol{r}_{p}-\boldsymbol{r}_{N,m+1}\Vert ,\quad p_{s}=1; & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \text{concave:}\quad p_{s_{m}},~p_{s_{m+1}}\in [0,1]~\rightarrow ~p_{n}=\min (p_{n_{m}},p_{n_{m+1}}),\quad p_{s}=p_{s_{m}}\vee p_{s_{m+1}}. & \displaystyle\end{eqnarray}$$

It should be noted that in the concave case, $p_{s}$ is chosen from either actuator element $m$ or $m+1$ depending on the minimum normal distance $p_{n}$ , possibly resulting in an index shift $m\rightarrow m+1$ . The resulting body force experienced at grid point $\boldsymbol{r}_{p}$ becomes $\boldsymbol{f}_{p_{s}}$ , which is determined as an interpolated value between $\boldsymbol{f}_{N,m}$ and $\boldsymbol{f}_{N,m+1}$ through

(3.5) $$\begin{eqnarray}\boldsymbol{f}_{p_{s}}=\boldsymbol{f}_{N,m}+p_{s}(\boldsymbol{f}_{N,m+1}-\boldsymbol{f}_{N,m}).\end{eqnarray}$$

Consequently, a continuous, rather than discrete, body-force distribution is obtained along the actuator line, thus allowing quasi-two-dimensional kernel integration normal to the associated actuator line element. This stands in contrast to ALM and ALM*, where each grid cell can be influenced by multiple discrete actuator points and experience cumulative effects. In ACE, on the other hand, each grid cell is influenced by a continuous and unique body force $\boldsymbol{f}_{p_{s}}$ , which need not be at a physical actuator point.

3.2 Actuator curve embedding – quasi-two-dimensional kernel integration

A quasi-two-dimensional (2D) Gaussian kernel function is the key element to eliminate both of the non-physical body-force overlap projections beyond the blade tips as seen in the standard ALM or ALM*. It should be noted that this kernel function had been suggested earlier by Mikkelsen (Reference Mikkelsen2003) and more recently by Shives & Crawford (Reference Shives and Crawford2013) for finite span wings. Here, the blade force, including lift and drag forces, is distributed in a disc normal to the local actuator element. Conservation of blade and body forces is satisfied through normalization with the local actuator-element width $\unicode[STIX]{x1D6E5}_{b}$ . In this respect, the volumetric Gaussian distribution of the body force in ACE is fundamentally different from that used in the ALM and ALM* implementations. The respective body force per unit volume is computed as

(3.6) $$\begin{eqnarray}\boldsymbol{F}_{p}(x_{P},y_{P},z_{P})=-\mathop{\sum }_{N}\frac{1}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6E5}_{b}}\boldsymbol{f}_{p_{s}}\exp [-(p_{n}/\unicode[STIX]{x1D716})^{2}],\end{eqnarray}$$

where the summation over the blade index $N$ is included for the sole purpose of completeness, for grid points $\boldsymbol{r}_{p}$ being located on symmetry planes between lifting lines. It should be noted that the total blade force $\boldsymbol{f}_{p_{s}}$ is conserved for the general case of $\unicode[STIX]{x1D6E5}_{grid}\neq \unicode[STIX]{x1D6E5}_{b}$ as the force per unit volume $\boldsymbol{F}_{p}$ is multiplied by a scaled cell volume $\unicode[STIX]{x0394}V^{2/3}\cdot \unicode[STIX]{x1D6E5}_{b}$ when injected into the momentum equation (2.1). Verification tests were conducted for straight as well as curved lifting lines approximated as second- and fourth-order polynomials (Jha Reference Jha2015). A variable Gaussian projection radius $\unicode[STIX]{x1D716}(r)$ is determined using

(3.7) $$\begin{eqnarray}\text{ACE:}\quad \unicode[STIX]{x1D716}/c^{\ast }=C_{1}(\unicode[STIX]{x03C0}AR)[C_{2}\cdot TSR+C_{3}]\end{eqnarray}$$

as a generalized version of (2.6), where $TSR$ is the ratio of lifting-line tip speed to free-stream speed, commonly referred to as the tip-speed ratio. Parametric studies conducted primarily on the NREL Phase VI rotor and analysed by Jha (Reference Jha2015) showed that a weaker $TSR$ dependence is necessary due to, as opposed to the parked case, the lifting line not traversing parallel to the sectional relative velocity vector $\boldsymbol{V}_{rel}$ . From ALM*, $C_{1}\in [0.025,0.030]$ , and $C_{2}=0.1410$ , $C_{3}=-0.3609$ should be used in conjunction with (3.7) for typical operating conditions of $TSR\geqslant 4$ , with $C_{2}=0$ and $C_{3}=1$ in the case of a parked blade or wing.