Evaluation of the Effects of Actuator Line Force Smearing on Wind Turbines Near-Wake Development

In the present study, the effects of the actuator line force smearing on the predicted near-wake development of a model wind turbine are evaluated. To exclude the uncertainties related to the blade loads computation, the loads of a corresponding blade-resolving simulation are imposed. Three force projection functions are applied, which all radially scale with the local chord length: an isotropic Gauss force distribution, an anisotropic Gauss distribution with distinct projection widths in the chord and thickness directions, and a Gauss-Gumbel projection function which creates an airfoil-like force distribution in the chord direction. Comparisons with corresponding blade-resolving simulations and experimental results show a good agreement of the tip vortex core size, position and circulation. Comparison of the mid-wake instantaneous streamwise velocity with blade-resolving simulations shows a thinner vortex sheet and stronger effects of the trailed vortices in the actuator line simulation. Although the bound vortex shape scales with the projection function, no significant difference was observed in the resulting tip vortices and near-wake flow field between the three projection functions used. It is attributed to the similar projection width in the tip vortex region for all cases considered. With the chosen grid refinement in the tip vortex area and a smearing width that scales with the chord length, tip vortices that compare with high-resolution experiments and high-fidelity blade-resolving simulation could be generated with the actuator line approach.


Introduction
The actuator line (AL) modeling of wind turbine rotor blades is typically used for numerical wake studies, where satisfyingly fine wake structures are produced at reduced computational costs when compared to simulations with resolved turbine boundary layers. Originally introduced by Sørensen and Shen [1], the AL approach combines the blade element theory to estimate the blade aerodynamic loads with high-fidelity methods to simulate the inflow and wake. While the unsteady aerodynamics at airfoil scale is per se not captured by this approach, the rotor-and blade-scale aerodynamics can be captured by the blade element theory. The latter provides the radial circulation distribution which determines the trailed vorticity and thus the wake. The resulting flow field features the major wake characteristics that control wake dynamics, like the blade tip vortices, which play a major role in wake development and breakdown mechanisms.
However, due to the modeling of the blade aerodynamic forces, differences are observed in the wake when compared to fully resolved simulations. Troldborg et al. [2] compared the wake velocity deficit and the turbulence level in the wake of a turbine modeled as a fully-resolved rotor, an AL rotor and an actuator disc rotor. The results indicated in uniform inflow a similar velocity deficit between the approaches but a much lower turbulence level in the AL and AD simulations. These differences were significantly reduced in turbulent inflow. In a wake comparison between fully-resolved and AL simulations for both an offshore case and a case in complex terrain, Weihing et al. [3] observed a similar wake structure for both approaches but the AL simulations overpredicted the velocity fluctuations in the mid part of the wake. Moreover, the tip vortices were more smeared and more elliptic with the actuator lines.
The two key elements of the AL method are i) the accurate computation of the blade loads from the local freestream velocity vector at each time step and ii) the projection of the aerodynamic forces on the flow field as body-forces.
Step i) demands a robust velocity sampling method along with high quality lift and drag data or, alternatively, generalized blade load correction models. To address i), Churchfield et al. [4] introduced a new integral velocity sampling method. An approach based on the lifting-line theory was chosen by Meyer Forsting et al. [5], [6] and similarly by Dag and Sørensen [7], that corrects the extraction of the velocity vector at the line position for low AL nodes resolution.
Step ii) is necessary to apply the discrete aerodynamic forces, computed for each actuator element, in form of a continuous blade loading onto the three-dimensional flow field.
The form of the projection of the body-forces greatly influences the near-wake velocity distribution around an actuator line. In the original AL formulation, the smearing function takes the form of a three-dimensional Gaussian [1] and has likely been introduced in order to avoid numerical instabilities. Recently, a more physical meaning of the force projection function has been striven for. Shives and Crawford [8] proposed a variation of the isotropic Gaussian that varies in the radial direction, where the Gaussian width is proportional to the airfoil chord length c, and thus varies along the blade span. Built upon this, Jha et al [9] introduced an equivalent elliptic planform chord to determine the Gaussian width in order to reduce overloading at the blade tip. Churchfield et al. [4] implemented an anisotropic Gaussian force projection, which width is proportional to the airfoil sectional geometry (chord length and airfoil thickness). It turned out that the anisotropic force projection generated a more realistic and compact blade wake when compared to a constant projection width along the blade. The tip vortices and their resultant downwash is better captured than with an isotropic force projection. Martínez-Tossas et al. estimated using a 2D approach that the optimum width for an isotropic and an anisotropic force projection lies in the range 0.25c and 0.2 − 0.4c, respectively, when compared to an analytical airfoil solution [10]. The optimal anisotropic projection width was applied to a wind turbine rotor and comparison with BEM calculations showed a very good agreement of the blade loads, attributed to a good prediction of the tip vortex core size [11]. Schollenberger et al. proposed for aircraft propellers an anisotropic Gaussian-Gumbel force projection, which mimics the lift distribution of an airfoil in chord direction [12].
The present study focuses on the effects of the force smearing on the near-wake flow field of a model wind turbine. The studied case corresponds to a low tip speed ratio (TSR), so that interaction effects between tip vortices are unlikely. The blade forces extracted from a high-fidelity Detached-Eddy Simulation (DES) are imposed, so that uncertainties related to the blade loads computation from the local inflow velocity and airfoil data are removed. A similar approach was used by Troldborg et al. [2] for uniform inflow conditions, but no parametric study of the AL projection function was performed. Although axial uniform inflow conditions are considered, flow separation occurred on the blade which was captured by the DES simulations. The same background grids were used in both the AL and the DES simulation, so that differences in the flow field mainly result from the modeling approach of the blade and its parameters. The present study aims at evaluating the effects of the actuator line force smearing on the near-wake structure, using width parameters based on recent progresses on AL modeling.
First, the methodology of the present approach is presented and the numerical setup and AL parameters are introduced. Then, the results for various projection functions are compared to experimental results and blade-resolving simulations and discussed. Finally, a short summary

Methodology
Two right-handed Cartesian coordinate systems are considered: i) a global system (x, y, z), centered on the rotor center with the x axis positive in the axial inflow direction and z axis positive away from the ground and ii) a local airfoil system (x c , x t , x r ), centered on a segment of the blade with x c in the chord direction and pointing towards the trailing edge, x t in the airfoil thickness direction and x r oriented along the blade line and pointing towards the blade tip.

Numerical Solver
The simulations have been conducted with the flow solver FLOWer. FLOWer is a compressible, finite-volume and block-structured solver [13], originally developed by the German Aerospace Center (DLR). It is the central link in the process chain for wind turbine simulations that has been established at the Institute of Aerodynamics and Gas Dynamics (IAG) of the University of Stuttgart [14]. In the present study, Delay-Detached Eddy (DDES) simulations are performed with a hybrid SST turbulence model. In the background and tip vortex meshes, the 5 th order WENO scheme is applied in order to reduce numerical dissipation of the vortical structures. The grids of the different components are overset with use of the Chimera grid overlapping method. The blade loads of the actuator lines are added as a body-force field to the Navier-Stokes equations. More details about the projection functions used are given in the following.

The frozen Actuator Line approach
As mentioned in the introduction, for the present investigations, the blade loads are not computed from the local flow field as it is traditionally the case for the AL. Instead, steady blade loads extracted from a corresponding high-resolution DES simulation are imposed in order to evaluate exclusively the effects of the force projection onto the flow field. It is worth mentioning that some studies have shown an influence of the smearing length scale on the computed loads [15]. Here, the assumption of frozen blade loads without interplay between the form of the projection function and the aerodynamic forces is made. Three variations of the force distribution function are used in the present study, that all consider the geometric properties of the blade. Accordingly, the projection width varies along the blade radius as following: • isotropic Gaussian: is a linear function of the chord length c, according to Shives and Crawford [8]. • anisotropic Gaussian: is a linear function of the chord length c whereby different proportion parameters are chosen in chordwise and spanwise direction of the airfoil in order to mimic the aspect ratio of the airfoil, according to Churchfield et al. [4]. • anisotropic Gaussian-Gumbel: is a linear function of the chord length c with different proportion parameters in chordwise and spanwise direction, similar to the anisotropic Gaussian projection. Moreover, a Gumbel function is used in chordwise direction, allowing a force distribution more like that of the airfoil, according to Schollenberger et al. [12].
The body forces compute respectively as follows: Here, ∆x c , ∆x t and ∆x r are the distance of a grid cell from the center of the AL element i in the chord, thickness and radial direction of the local frame of reference of the airfoil, c , t and r the distribution width in the respective directions for anisotropic distributions, and f i is the aerodynamic loading of the blade element considered. In previous studies, Shives and Crawford [8] recommended ≈ 0.25c for an isotropic force distribution, Churchfield et al. [4] used c = 0.85c and t = 0.85t, with t the airfoil thickness, and Martínez-Tossas et al. [10] have shown c = 0.4c and t ≤ 0.2c to be optimal parameters for a 2D flow around a flat plate and a Joukowski airfoil.
According to previous investigations, numerical stability is ensured if /∆x > 2 − 4, with ∆x the local grid size [16]. Subsequently, a minimal user-defined width min is implemented here, which traduces in a filter function for the width parameter j in any spatial direction j The MEXICO model rotor [17] is considered here, which features a radius R = 2.25m and has been widely used for code comparison studies in the field of rotor and wake aerodynamics. In this context, the geometry of the rotor has already been simulated with FLOWer in the past [18], [19]. In the present study, the blade is modeled with an actuator line and to account for the flow around the nacelle, its exact geometry is considered. The background and nacelle grids from Weihing [19] are used. A low tip speed ratio is considered and the parameters of the operating point, that was taken from the experiments of the Mexnext campaign [20], are summarized in Table 1.

Numerical setup
A one-third model is considered here, with one rotor blade embedded in a 120 • section of a cylindrical background mesh, as represented in Figure 1. Periodic boundary conditions are applied laterally, resulting in the simulation of a three-bladed rotor. Farfield boundary conditions are applied at the outer boundaries of the 120 • background mesh, except the inner boundary that is modeled as an Euler wall. The background, tip vortex and geometry-conform nacelle grids from the study of Weihing [19] are used. They feature a near-wake and a tip-vortex refinement area extending 1.5R into the wake with a respective grid resolution of 0.01x0.01x0.01m 3 and 0.004×0.005×0.006m 3 in axial, radial and tangential directions. Finally, a cylindrical mesh is created for the actuator line with a resolution of 0.004×0.004×0.004m 3 . The whole setup results in 118 Mio cells and a time-step equivalent to a blade azimuthal displacement of 1/32 • is used. The actuator line is modeled with 124 unevenly distributed nodes with a higher density at the blade root and blade tip, where higher gradients in radial direction are expected. The imposed thrust and driving force, extracted from a corresponding high-resolution Bernoulli-based DES (BDES) simulation from Weihing [19] -a DES variation with advanced shielding functions -, are represented in Figure 2. The forces show characteristic changes in the gradient, particularly in the mid-blade section where different airfoil types are blended. As the low TSR operating point features flow separation, due to an approximate angle of attack of 15 • , high temporal fluctuations are observed in the loads of the fully turbulent BDES simulation. Only the mean distribution is considered here. Figure 4 shows the isosurface of the resulting body-force component in streamwise direction for the three projection functions. As mentioned in the previous section, the projection width depends linearly on the local blade chord length, resulting in a decreasing j value with increasing radius. The thrust shows the opposite trend, as observed in Figure 2,  Figure 1. Numerical setup with background (grey) mesh and its refinement area (turquoise), nacelle mesh(red), tip vortex refined mesh (yellow) and blade refinement mesh (green). For more clarity, only every 8 th line is represented.

Results
To analyse the effects of the different approaches on the near-wake structure, three domains of interest are separately analysed. First, the bound vortex is inspected, then the tip vortices properties are evaluated and compared to fully-resolved BDES simulations and experimental data, and finally the mid-section of the wake is considered. Figure 5 shows the vortex structure of the near-wake flow field for the three force distribution functions. All cases feature a similar structure comprising a bound vortex along the AL, a vortex helix emerging from the blade-tip area and near-wall vortical structures that are created at the vicinity of the nacelle. A thin trailing vortex sheet is visible in the vicinity of the trailing edge within the refined blade domain, which is more persistent in case of the Gauss-Gumbel projection function. The isosurfaces show a wider bound vortex for the anisotropic case ( c = 0.4c and t = 0.2c) than for the isotropic case with = 0.23c or the anisotropic Gauss-Gumbel force distribution, which presents the most compact vortex structure. The scaling of the bound vortex with the force distribution size is in line with the observations of Churchfield et al. [4].

Bound vortex
To allow a closer investigation of the local influence of the force distribution function, Figure 6 shows the local vorticity distribution in a cut through the actuator line. The vorticity reflects the asymmetric force distribution and a more blurred bound vortex is observed in the cases with larger projection width . At the close vicinity of the actuator element, the streamlines are influenced by the form of the projection function. Farther away, about 1 − 2c away from the actuator element, the streamlines of all cases are on top of each other. The thin vortex layer developing in the wake of the blade is due to the narrow distribution width in thickness direction and is more compact with the Gauss-Gumbel distribution, as t is lowest in this case.  Figure 5. Isosurface of the λ 2 -criterion for the isotropic Gaussian (left), anisotropic Gaussian (center) and Gauss-Gumbel (right) force distribution. Figure 6. Vorticity of the bound vortex and streamlines in a cut through the AL at r/R=0.7 for the isotropic (left), anisotropic (center) and Gauss-Gumbel (right) force distribution. Figure 7 compares the evolution of the AL vortex properties downstream the rotor with the BDES vortex helix and the high-resolution PIV measurements of the MexNext measurement campaign [20]. The vortex cores have been identified through the vorticity peak at their center, and their radius evaluated as the half distance between extremas of the circumferential velocity. The circulation was obtained via integration of the vorticity over a constant circular area of 0.2 m centered on the core [21]. Figure 7 shows a relatively good agreement of the vortex position, core radius and maximal vorticity at the core center for all approaches.

Tip vortex properties
It is noticeable that the core size is slightly overpredicted -and reciprocally the maximum vorticity slightly underpredicted -in all AL simulations when compared to both experiments  Figure 7. Comparison of the vortex properties : vortex core position y/R, vortex core size c/R, vortex circulation Γ and the normalized vorticity at the center of the vortex core Ω y c tip /U tip and blade resolving simulations. The higher slope of the vortex aging for AL simulations than for BDES simulation suggests a higher numerical diffusion, although the same grid is used in the tip vortex region. This is attributed to the choice of a 5 th order WENO spatial scheme for the AL simulation compared to the higher-order WENO-SYMBO spatial scheme for the BDES simulations. The vortex circulation is slightly overpredicted in the AL simulations, and more generally overpredicted in the CFD simulations when compared to the experiments. This is attributed to the presence of a thin vortex sheet rolling around the tip vortex in the simulations, which is not present in the integration area for the experimental data.
Interestingly, no clear distinction can be made in the tip vortex properties for the different force projection functions. In particular, differences could be expected in the core size as the bound vortex shape scales with the force distribution. Considering the stability condition in Equation (4) and the chosen proportionality constants, all the variations of the force projection investigated here result in a similar projection width in the outer blade tip region. It is believed that the size of the tip vortex core is significantly influenced by the bound vortex width in the blade tip area, which explains why the tip vortex properties are hardly influenced by the smearing function chosen here. Figure 8 shows the streamwise velocity distribution in an axial plane through the wake. As the resulting flow fields are very similar for all AL simulations, only the results with use of the Gauss-Gumbel projection function are compared here to the blade-resolving simulation. Both the tip vortices and the vortex sheet in the wake of the blade are observed in this plane. In the nacelle region, the flow field is very similar in both approaches, and the hub vortex is well captured by the DDES method in the AL setup. The velocity fluctuations observed in the midpart of the wake are comparable in strength for both approaches. The vortex sheet in the AL wake is however thinner than the one predicted by the BDES simulations. This is attributed to the viscous velocity deficit in the blade wake, which is particularly thick due to the massive flow separation prevailing at the considered inflow conditions and which is not captured in the AL simulation. The BDES simulation moreover features flow separation, which cause loads fluctuations and add unsteadiness into the flow field. Thus, the higher complexity of the blade wake in the BDES simulation leads to the breakdown into small-scale turbulent structures, that are not observed in the AL simulations. The trace of the trailed vortex in the mid-part of the blade, due to the loads gradient related to the airfoil change in this section, is observed in the AL simulation. Its effects on the velocity distribution are more intense when compared to the BDES simulation, as in the AL setup the trailed vortex does not overlap with small-scale turbulence.

Summary and conclusions
In the present study, simulations of the MEXICO rotor were performed with the actuator line approach in order to evaluate the effects of the force smearing on the near-wake development. Three force projection functions were applied, which all radially scale with the local chord length: an isotropic Gauss force distribution, an anisotropic Gauss distribution with distinct projection widths in the chord and thickness directions, and a Gauss-Gumbel projection function which creates an airfoil-like force distribution in the chord direction. For the isotropic and anisotropic Gaussian cases, the projection width was chosen according to recommendations from the literature. In order to asses only the effects of the force smearing of the AL approach, without the uncertainties related to the blade load computation, the mean blade loads extracted from the BDES simulation were imposed at the actuator line.
It was observed that the bound vortex shape scales with the force projection function, being hence more compact and featuring a higher vorticity with smaller distribution width. A thin vortex layer was observed in the wake of the actuator line, which evolves in a helical shape through the wake. Comparison of the resulting small-scale velocity fluctuations in the mid-part of the wake with blade-resolving BDES simulations show comparable fluctuation amplitude. However, the vortical sheet is wider in the BDES simulation due to the viscous blade shear layer and the massive flow separation. Moreover, the trace of the trailed vortices from the mid-part of the blade is more intense in the actuator line simulation. Finally, no significant difference was observed in the tip vortex properties for the three force projection function used here. It is attributed to the similar distribution width in the tip vortex area, which conditions the tip vortex form.
To conclude, tip vortices that compare well with high-resolution experiments and high-fidelity blade-resolving simulation were generated with the AL approach under two conditions: the projection width scales with the chord length and sufficiently fine grid cells are used in the blade tip area. The grid size in the tip region corresponds here to 0.1c tip , with c tip a representative chord length for the blade tip region. The quality of the tip vortex prediction is a major requirement in wake breakdown and wake interaction studies, for which the present findings could find application. Further investigations using blade loads computation from airfoil data instead of specifying them are planned to validate those results for future applications.

Acknowledgements
The authors gratefully acknowledge the HLRS Stuttgart for providing computational resources.