Comparison of Free Vortex Wake and BEM Structural Results Against Large Eddy Simulations Results for Highly Flexible Turbines Under Challenging Inﬂow Conditions

. Throughout wind energy development, there has been a push to increase wind turbine size due to the substantial economic beneﬁts. However, increasing turbine size presents several challenges, both physically and computationally. Modeling large, highly ﬂexible wind turbines requires highly accurate models to capture the complicated aerodynamic response due to large deﬂections and nonstraight blade geometries. Additionally, development of ﬂoating offshore wind turbines requires modeling techniques that can predict large rotor and tower motion. Free vortex wake (FVW) methods model such complex 5 physics while remaining computationally tractable to perform the many simulations necessary for the turbine design pro-cess. Recently, a FVW model—cOnvecting LAgrangian Filaments (OLAF)—was added to the National Renewable Energy Laboratory engineering tool OpenFAST to allow for the aerodynamic modeling of highly ﬂexible turbines along with the aero-hydro-servo-elastic response capabilities of OpenFAST. In this work, FVW and low-ﬁdelity blade-element momentum (BEM) structural results are compared to high-ﬁdelity simulation results for a highly-ﬂexibly downwind turbine for varying TI, shear 10 exponent, and yaw misalignment conditions. Through these comparisons, it was found that for all considered quantities of interest, SOWFA, OLAF, and BEM results compare well for steady inﬂow conditions with no yaw misalignment. For OLAF results, this strong agreement was consistent for all yaw misalignment values. The BEM results, however, deviated signiﬁcantly more from SOWFA results with increasing absolute yaw misalignment. Differences between OLAF and BEM results were dominated by yaw misalignment angle, with varying shear exponent and TI leading to more subtle differences. Overall, 15 OLAF results were more consistent than BEM results when compared to SOWFA results under challenging inﬂow conditions. Percent differences of time-averaged means of BEM results increase with increasing absolute yaw angle, averaging 5 . 46% for no yaw misalignment and increasing to 15 . 9% for ± 30 ◦ yaw misalignment. OLAF results show no dependence on yaw standard deviations the results. These results show consistent percent differences for all shear exponents across 260 all QoIs. For mean results, percent differences deviate at the highest shear exponent only, with the median percent difference increasing. For standard deviation results, the median percent difference remains consistent across all shear exponents, but the these differences are not for most other QoIs, it is important to note that changing shear exponent can have a signiﬁcantly different effect on OLAF compared to BEM results for certain important QoIs.

OLAF solves for the turbine wake in a time-accurate manner, which allows the vortices to convect, stretch, and diffuse. The OLAF model is based on a Lagrangian approach, in which the turbine wake is discretized into Lagrangian 70 markers. There are many methods of representing the wake with Lagrangian markers (Branlard (2017)). In this work, a hybrid lattice/filament method is used, as depicted in Figure 1. Here, the position of the Lagrangian markers is shown in terms of wake age, ζ, and azimuthal position, ψ, though 75 in the code they are defined in terms of Cartesian coordinates. A lattice method is used in the near wake of the blade. The near wake spans over a user-specified angle. Though past research has indicated that a near-wake region of 30 • is sufficient ; Ananthan et al. (2002)), it has been shown that a larger near wake is required for high thrust and other challenging conditions. This is further investigated in this work. After the near wake region, the wake is assumed to instantaneously roll up into a tip and root vortex, which are assumed to be the most dominant features for the remainder of the 80 wake ). Each Lagrangian marker is connected to adjacent markers by straight-line vortex filaments. The wake is discretized based on the spanwise location of the blade sections and a specified time step (dt), which may be different from the time step of AeroDyn15. The wake is allowed to move and distort with time, thus changing the wake structure as the markers are convected downstream. To limit computational expense, the root and tip vortices are truncated after a specific distance downstream from the turbine. The wake truncation violates Helholtz's first law and hence introduces an erroneous 85 boundary condition. To alleviate this, the wake is "frozen" in a buffer zone between a specified buffer distance and the wake length. In this buffer zone, the markers convect at the average ambient velocity. In this way, truncation error is minimized ). The buffer zone is typically chosen as the convected distance over one rotor revolution.
As part of OpenFAST, induced velocities at the lifting line/blade are transferred from OLAF to AeroDyn15 and used to compute the effective blade angle of attack at each blade station, which is then used to compute the aerodynamic forces on 90 3 https://doi.org/10.5194/wes-2021-130 Preprint. Discussion started: 12 January 2022 c Author(s) 2022. CC BY 4.0 License. the blades. The OLAF method returns the same information as the BEM method, but allows for more accurate calculations in areas where BEM assumptions are violated, such as those discussed above.
The governing equation of motion for a vortex filament is given by where r is the position vector of a Lagrangian marker and V is the velocity. Vortex filament velocity is a nonlinear function 95 of the vortex position, representing a combination of the freestream and induced velocities (Hansen (2008)). The induced velocities at each marker, caused by each straight-line filament, are computed using the Biot-Savart law, which considers the locations of the Lagrangian markers and the intensity of the vortex elements ): Here, Γ is the circulation strength of the filament, dl is an elementary length along the filament, and r is the vector between a 100 point on the filament and the control point x. The regularization factor F ν is introduced because of the singularity that occurs in the Biot-Savart law at the filament location (Vatistas et al. (1991)). Viscous effects prevent this singularity from occurring and diffuse the vortex strength with time. OLAF uses a different regularization value for the blades and the wake, and the wake regularization factor increases with time according to a parameter referred to as the core spreading eddy viscosity.
The circulation along the blade span is computed using a lift coefficient-based method. This method determines the circula-105 tion within a nonlinear iterative solver that utilizes the polar data at each control point on the lifting line. The algorithm ensures that the lift obtained using the angle of attack and the polar data matches the lift obtained with the Kutta-Joukowski theorem.
This method is based on the work of van Garrel (2003) and is further detailed in the OLAF User's Guide and Theory Manual ).

110
Discretization and regularization studies were performed using the Big Adaptive Rotor (BAR) Downwind Rail Carbon (DRC) turbine (Bortolotti et al. (2021)), described in Table 1, to determine recommended values for various OLAF input parameters.
The studies were performed with uniform wind speeds ranging from 4 to 12 m/s, with the ranges for each study shown in Table 2. The studies were done sequentially as presented here. Meaning, the wake discretization study was the first to be completed, from which an optimal value was determined. This value was used as the wake discretization value for all remaining 115 studies. The near wake extent study was performed next, and this optimal value was used for the remaining studies, and so on. Output quantities of interest (QoIs) are summarized in Table 3. These quantities are used to evaluate the convergence and accuracy for a set of parameters. An overview of the results is included here; for detailed results, see the online OLAF documentation ).  The wake discretization and convergence studies observed how wind turbine QoIs varied with changing wake discretization, near wake extent, and far wake extent. These are all user-defined parameters for OLAF and are known to have a large effect on both wind turbine response and simulation time. The wake discretization is specified by the AeroDyn15 time step, using the equation dψ = dt × Ω. Because rotor speed is a turbine characteristic, dψ is defined using dt, which remains constant throughout the simulation. A smaller time step results in a finer wake discretization. The near wake extent is specified by the 125 number of revolution degrees to be used in the lattice near-wake computations. Similarly, the far wake extent is determined by specifying the remaining length of the turbine wake length used in the straight-line filaments for wake calculations. By reducing the time step or increasing the wake length, the simulation accuracy is expected to improved, along with the computational expense. Accuracy is assessed by computing the percent difference in the mean QoI value at each test point relative to the finest or largest value for the discretization and wake extent studies, respectively. To balance model accuracy with computational cost, 130 a convergence criteria of~2% was selected.
An optimal value was found across all wind speeds for each study, as reported in Table 4. Generator power convergence for each parameter is shown in Figure 2. The convergence is hard to observe for the wake discretization study. Results for fine wake discretization are likely to be a function of the wake regularization factor, which is a topic that we believe needs further research (see Section 2.2.2). The power converges for increasing near wake extend, since it is expected that a longer 135 wake captures more induction at the rotor (see the discussion related to Figure 4). Variability in relative error is observed for increased far wake extends, but the error generally decreases for longer far wakes. Computational time for variations in each parameter is shown in Figure 3. Computational time increased linearly with far wake extent (not shown). All simulations are run using a tree-code algorithm with a time complexity O(n log n). Note that the values reported in Table 4 were deemed optimal for the specific wind turbine used in this paper under uniform inflow wind, and are not necessarily optimal for other 140 configurations. However, they can serve as a starting point for other systems and inflows.
The wake length required for a given accuracy can be assessed by comparing the induced velocity at the rotor plane as obtained from a vortex cylinder of finite length (u L ) and from a vortex cylinder of infinite length (u ∞ ), the ratio of which is shown in Figure 4. Both velocities can be obtained using closed form formulae (Branlard and Gaunaa (2014)). It is seen that 99% of the axial induction is provided by the first 4 diameters of the wake of a vortex cylinder. This observation supports the 145 recommended value of 5D selected in Table 4.

Regularization
Wake regularization factor, blade regularization factor, and core spreading eddy viscosity were examined for the regularization study. The method was the same as the wake discretization and extent studies, except that results were compared against 150 those obtained using the actuator-line formulation of SOWFA (Churchfield et al. (2012)) to determine the optimal parameter values. This was done because it is not possible to know a priori which regularization parameter can be used as a reference for the convergence study. The questions remains as to whether the actuator-line simulations can be used as a reference. The different parameter values found to minimize the difference with SOWFA for different wind speeds are reported in Table 5.
Computational time was generally insensitive to parameter values, so is not reported. An example plot is shown in Figure 5    for the study of the wake regularization factor. Based on the results given in Table Table 5, no clear trend is found between the regularization factors and the wind speed. It appears that a reasonable approximation would consist of setting the regularization parameters to 0.5 dr, where dr is the blade spanwise discretization. More research is nevertheless needed, both for actuator line CFD and lifting-line based codes, to find the optimal values of these parameters.

160
In this work, the OpenFAST-FVW code was compared to traditional OpenFAST-BEM and large-eddy simulation results.
Comparison of the FVW method to the traditional BEM method in OpenFAST is performed using the BAR-DRC turbine 8 https://doi.org/10.5194/wes-2021-130 Preprint. Discussion started: 12 January 2022 c Author(s) 2022. CC BY 4.0 License. (Bortolotti et al. (2021)). Some characteristics of the turbines are given in Table 1. Large-eddy simulations were performed using the Simulator fOr Wind Farm Applications (SOWFA) (Churchfield et al. (2012)). SOWFA solves the filtered Navier-Stokes equations in primitive variables using a second-order control volume formulation. The wind turbines blades are modeled 165 as actuator lines (Sørensen and Shen (2002); Martínez-Tossas et al. (2015)). For all simulations, structural modeling was done using the structural module ElastoDyn from OpenFAST. For BAR-DRC comparisons, the turbine controller was not used and all simulations were performed at a set rotor speed and blade pitch.
Simulations were performed at a range of freestream velocities, turbulence intensities (TI), and yaw misalignment angles, as summarized in Table 6. For steady OpenFAST simulations, InflowWind was used to specify the hub-height inflow velocity and Table 6. Variation of inflow conditions for code comparisons.

Freestream Velocity TI Yaw Misalignment Shear
shear exponent. For the turbulent OpenFAST simulations, ambient wind inflow is generated synthetically by TurbSim (Jonkman (2014)), which creates two-dimensional (2D) turbulent flow fields. Turbulence is simulated using the Kaimal spectrum with exponential coherence model and the standard IEC turbulence model was used. The time-dependent 2D wind field is propagated along the wind direction at the mean wind speed of the midpoint of the field. Each case from the measurements is simulated in TurbSim and OpenFAST six times, using six different random turbulence seeds to generate the wind data at the spatially 175 coherent points. This was done to capture variability in the numerical simulation results associated with the uncertainty imposed by the unmeasured wind inflow.

Results
For the BAR-DRC turbine, the results focus on verification of statistical results between OLAF, BEM, and SOWFA. The QoIs were turbine performance, blade-root bending moments, and tower-base bending moments. Blade-span and azimuthal results 180 were also considered.

Code-to-Code Comparisons
For this code-to-code comparison, many QoIs were considered, as summarized in Table 7. The blade forces and induction quantities are only considered as blade-span quantities, whereas the others are also considered as time-and azimuthallyaveraged quantities. Though all these QoIs were considered, only a select few are presented for much of the analysis. All QoIs 185 are included in box-and-whisker plots. All inflow conditions are for steady inflow unless stated otherwise. Yaw misalignment cases were simulated in constant inflow (shear exponent = 0) and sheared and TI comparison cases were simulated with no yaw misalignment.
For all QoIs, the results from each computational method follow the same trend with changing yaw misalignment. The trend depends on the QoI, with most showing reduced values with increasing absolute yaw misalignment. Tower-base side/side and 195 yaw moments showed decreasing values with non-absolute increasing yaw misalignment, and blade-root pitching moment showed the opposite trend. Though the same trends were captured by each computational method, the actual results differ.
In particular, OLAF and SOWFA results show comparable values with an average percent difference of <5% for all yaw misalignment angles, including the most extreme. However, BEM results deviated significantly more from SOWFA results with increasing absolute yaw misalignment, with percent differences at yaw misalignment of ±30 • reaching up to 48.6%.

200
Alternatively, the percent difference between OLAF and SOWFA results does not appear to be dependent on yaw misalignment. Figure 7 are box and whisker plots of percent difference values between time-averaged mean and standard deviations, computed using Equation 4, of SOWFA results and OLAF or BEM results for all QoIs listed in Table 7 except the blade distributed quantities and those with means near zero.        percent difference for SOWFA is shown as zero because this is the highest fidelity code. This chart clearly shows the near 12× error reduction when using OLAF instead of BEM. OLAF does, however, require significantly more computational power than BEM, with OLAF simulations taking O(40 CPU hours to complete for the simulation specifications used here, whereas BEM requires minutes. When compared to SOWFA, however, OLAF is significantly faster, with SOWFA simulations for this 245 work requiring approximately 260 CPU hours to complete. As always, the choice of which computational method to use is a balancing act between computational cost and accuracy. These results indicate that, given the problem being considered, OLAF could be a reasonable middle ground between engineering models and CFD.

Shear Exponent
Shown in Figure 11 are the time-averaged quantities for several QoIs from OLAF and BEM simulations for a range of shear here at the zero-shear point. This is because of the complexities within SOWFA of running steady inflow with a specified shear exponent as well as the additional computational expense. For all results but tower-base yaw moment, the relative trends of OLAF and BEM are the same with changing shear exponent, though the percent difference between the results does increase slightly with increasing shear exponent. For all shear exponents, OLAF predicts higher loads for the considered components.

255
Tower-base yaw moment, however, increases at a sharper rate for BEM result compared to OLAF results. In fact, at low shear exponents OLAF predicts a higher tower-base yaw moment, while BEM predicts a higher yaw moment at a shear exponent of 0.2. This results in a reduced percent difference between the methods. These results are summarized in Figure 12  upper range of these values does fluctuate slightly. These results indicate that though shear exponent does have a slightly higher impact on time-averaged BEM results compared to OLAF results, it is to a much lesser extent than yaw misalignment. Figure 13 are azimuthal-averaged results for all computational methods as well as percent differences between 265 BEM and OLAF results. The overall shapes of these results are comparable to those shown for the varying yaw misalignment cases in Figure 9. For rotor torque and and tower-base fore/aft bending moment, there is a clear separation between OLAF and BEM results. For rotor torque, the results compare best at the dips corresponding to the blades passing behind the tower.

Shown in
For each method, there is minimal change to the results at this point in the blade path with changing shear exponent. This is expected, since the flow behind the tower is dominated by the tower wake and likely unaffected by the inflow shear exponent.

270
When the blades are not behind the tower, shear exponent has more of an effect on the results, though to a greater extent for BEM results, with a higher shear exponent resulting in reduced rotor torque. Minimal differences in tower-base fore/aft bending moment are seen for each method. Overall the moment increases with increasing shear exponent, though OLAF shows appreciably higher moment at all azimuthal angles at a shear exponent of 0.2. This results in higher percent difference at this inflow condition. When comparing results for out-of-plane blade-root bending moment, the largest difference are seen 275 when the blade passes behind the tower. At this location, OLAF predicts a higher bending moment compared to BEM. For both methods, this bending moment decreases with increasing shear exponent. The bending moment increases until the blade is above the turbine, with this maximum value increasing with increasing shear exponent. For all shear exponents, OLAF predicts a higher bending moment than BEM. Though these overall trends are comparable between OLAF and BEM, the average percent difference between these results does increase with increasing shear exponent, though to a greater extent when 280 the blade is in the tower wake. Differences in tower-base yaw moment are considerable between OLAF and BEM, with different azimuthal trends seen at all shear exponents. For all OLAF results, the yaw moment spikes when a blade is behind the tower, followed by a sharp decrease and then another smaller spike before another blade passes behind the turbine. This trend is seem for all shear exponents, though the sharp drop becomes shallower with increasing shear exponent. BEM results are also characterized by a spike when 285 a blade is behind the tower. However, at low shear exponents, the value of this spike is considerably less than that predicted by OLAF, and is preceded by a sharp drop to a value that remains constant until another blade passes behind the tower. This drop value is comparable to that predicted by OLAF, but the secondary bump that OLAF predicts is not present for the BEM results.
As shear exponent increases, the large spike predicted by BEM surpasses that predicted by OLAF, and the sharp drop is instead a gradual reduction until the minimum value is reached. Thus, while percent difference between the methods remains large for 290 all shear exponents, the primary takeaway is that changing shear exponent drastically changes the azimuthal trend predicted by BEM, whereas OLAF trend changes are more subtle. While these differences are not seen for most other QoIs, it is important to note that changing shear exponent can have a significantly different effect on OLAF compared to BEM results for certain important QoIs. Figure 14 are the time-averaged quantities for several QoIs from OLAF and BEM simulations for a range of TI values, along with percent differences of BEM results relative to OLAF results. In terms of mean results, varying turbulence intensity effects OLAF and BEM comparably. There is some change in percent difference for varying turbulence intensity, but it remains within a few percentage points. This is further supported in Figure 12, which shows box-and-whisker plots for each   soon as turbulence is introduced into the inflow, this relationship flips and BEM predicts significantly lower rotor torque and OoP blade-root bending moment for all azimuthal locations.

Conclusions
The purpose of this work was to determine the accuracy of an aeroelastically coupled FVW method under conditions known to be challenging to lower-fidelity aerodynamic methods. This was done by comparing OLAF results to high-fidelity SOWFA results and low-fidelity BEM results for a large range of TI, shear, and yaw misalignment conditions. Through these comparisons, no yaw misalignment. For OLAF results, this strong agreement was consistent for all yaw misalignment values, with percent difference values of time-averaged results remaining within 1.2% across all QoIs and varying little with yaw misalignment.
The BEM results, however, deviated significantly more from SOWFA results with increasing absolute yaw misalignment, with percent differences at yaw misalignment of ±30 • reaching up to 48.6%. These trends were true for standard deviations of the results as well, though the BEM results showed less of a change with increasing absolute yaw misalignment. Similar trends 320 were seen when comparing blade-span loads, with OLAF results remaining within 7.9% of SOWFA results across the blade span and for all yaw misalignment angles. Maximum BEM percent difference results increased from 22.1% for no yaw misalignment up to 56.6% for ±30 • yaw misalignment, with the differences beginning at 20% blade span and continuing outboard of the blade.
Consideration must also be given to modeling computational cost. When looking at percent differences of time-averaged 325 results for all QoIs, OLAF showed a near 12× error reduction when using OLAF instead of BEM. As with all methods, increased accuracy comes with increased computational cost, with BEM, OLAF, and SOWFA simulations taking on the order of minutes, 10 CPU hours, and 100 CPU hours to complete for the simulations performed in this work, respectively.
Author contributions. KS led the investigation of OLAF verification work and SOWFA and BEM comparison. BA completed the discretization and convergence studies. LM ran all SOWFA simulations. EB and NJ provided support for this project. KS prepared the article, with of Energy Office of Energy Efficiency and Renewable Energy, Wind Energy Technologies Office. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.
The research was performed using computational resources sponsored by the Department of Energy's Office of Energy Efficiency and