Investigation of a new model accounting for rotors of finite tip-speed ratio in yaw or tilt

The main results from a recently developed vortex model are implemented into a Blade Element Momentum(BEM) code. This implementation accounts for the effect of finite tip-speed ratio, an effect which was not considered in standard BEM yaw-models. The model and its implementation are presented. Data from the MEXICO experiment are used as a basis for validation. Three tools using the same 2D airfoil coefficient data are compared: a BEM code, an Actuator-Line and a vortex code. The vortex code is further used to validate the results from the newly implemented BEM yaw-model. Significant improvements are obtained for the prediction of loads and induced velocities. Further relaxation of the main assumptions of the model are briefly presented and discussed.


Introduction
The analytical studies from Glauert in 1926 [1] and Coleman in 1945 [2] form the basis of most yaw-models implemented in Blade Element Momentum(BEM) codes. Yet, these models strictly apply to rotors of infinite tip-speed ratios. Coleman used a skewed vortex cylinder to represent the wake behind a rotor. The underlying assumptions of his study are: potential flow, infinite number of blades, infinite tip-speed ratio and constant rotor circulation. Coleman focused on the axial velocity induced by the tangential vorticity component which is indeed the main source of induction for large tip-speed ratios. The ratio between the right cylinder induction and the skewed cylinder induction provides a correction factor that is applied in BEM codes. The lifting-line analogy of Glauert provides a similar correction factor.
In a recent work from the authors [3], Coleman's vortex system was extended to assess the effect of finite tip-speed ratio. Relaxing this sole assumption results in an extended vortex system made of a bound vortex disk, a root vortex and a vortex cylinder with tangential and longitudinal vorticity. Similar analysis were previously performed for subsets of this system by Snel and Schepers [4] and Burton et al. [5]. The system was studied by the current authors both for right [6] and yawed [3] configuration in order to provide correction factors to be applied in the BEM algorithm. The analysis was not restricted to the axial induction but considered also the tangential and radial velocity components. Analytical and semi-analytical formulations were provided for each velocity and vorticity components together with fitted models for engineering applications. The model may be applied both to yaw or tilt configurations. The current paper focuses on the application of this newly developed model. Results from BEM, Computational Fluid Dynamics (CFD) and vortex methods simulations are compared with measurements.
The structure of the paper is as follows. The vortex system used for the study is briefly presented. The induced velocity formulae for each vortex element are provided. The methodology for the BEM code implementation of the new yaw-model is then developed. Results obtained with the new yaw-model are compared to results obtained with Coleman's yaw-model. The model is validated using experimental data from the MEXICO experiment [7] and numerical results from actuator-line(AL) and vortex-method simulations. The limitations of the model are studied to reveal the effects of finite-number of blades, wake roll-up, shed-vorticity and radially varying circulation.

Presentation of the model
The vortex model used in this study comprises: a (skewed) semi-infinite vortex cylinder with tangential and longitudinal vorticity, a bound vortex disk, and a root vortex. The coordinate and vortex systems are illustrated in Figure 1. The model is described in details in a previous study by the authors [3]. The normal to the rotor plane is directed along e z . The wake axis is directed along e ζ which is obtained by rotation of e z by an angle χ around the y-axis. The angle χ is referred to as the skew angle. The skew angle is slightly larger than the yaw angle, as sketched in Figure 2. Relations between the two angles can be found in the literature [8,5,9] or can be computed numerically (see section 4). The total circulation of the rotor is noted Γ tot . The circulation is positive for the convention of a wind turbine rotating along e ψ . The root vortex is directed along the wake axis such as its intensity is Γ r = −Γ tot e ζ . The vortex cylinder consists of two continuous and constant vorticity components in the tangential and longitudinal direction, respectively noted γ t and γ l and given by: where h is the pitch of the helical wake. The bound vorticity is assumed uniformly distributed on the rotor disk. Its influence is not considered in the current study since it is zero in the rotor plane. The notation u •, will refer to the •-component of the velocity induced by the vortex element , where • ∈ {x, y, z, r, ψ} and ∈ {r, t, l}. The velocity field induced by the full vortex system in the rotor plane is: Coleman et al. [2] and Castles and Durham [10] studied the axial velocity induced by γ t , i.e. u z,t . The different velocity components from the full vortex system were studied in a previous study from the authors [3]. Analytical and semi-analytical solutions to the Biot-Savart integrals were obtained and engineering models provided for r < R. These models are summarized below. The inductions from the tangential vorticity are modelled as: where u z,0 = γ t /2. The notation ψ m = ψ − ψ 0 is introduced to conveniently change between different conventions for the origin of the variable ψ. The case ψ 0 = 0 corresponds to the convention of this document ( Figure 1). The case ψ 0 = −π/2 should be used if the azimuth angle is zero for a blade that is pointing upwards, i.e. along y. Different expressions for the envelope factor F t exist [3]. The induced velocities u z,t and their exact envelope are plotted in Figure 2. In the current paper, the function F t will be approximated by its tangent at the origin: F t (r) =r/2, withr = r/R. Under this assumption, the left expression of Equation 3 is equivalent to Glauert's result [1, p. 572]. The velocities induced by the root vortex can be Wake axis Figure 2: Top-view of the wake cross-section laying in the plane y = 0, data for χ = −30 • . The velocities induced by the tangential vorticity along the fore-aft diameter of the rotor (z = 0) are shown. Their relative fluctuations around the mean induction, or envelope, has a slope of tan χ 2 . The distribution of velocity induces normal forces that can possibly (e.g. below stall) generate a restoring moment. determined exactly using the induction formula from a semi-infinite vortex filament [6]: u z,r (r, ψ, χ) = Γ r sin ψ m sin χ 4πr(1 − cos ψ m sin χ) , u ψ,r (r, ψ, χ) = Γ r cos χ 4πr(1 − cos ψ m sin χ) .
The following engineering model was suggested for the influence of the longitudinal tip-vorticity: where the dependencies in (r, ψ, χ) were dropped to shorten notations and with: 3. Implementation Analyses of the right vortex cylinder [6] reveal the following link between the induced velocities from actuator-disk(AD) momentum theory and the vortex model: with Ω the rotational speed of the AD, U 0 the incoming free stream, a the axial induction factor positive along e z and a the tangential induction factor negative along e ψ .
To implement BEM yaw-models based on vortex theory results it will be assumed that the differences in induced velocities between the right and skewed vortex systems reflects the change of induced velocities that should apply to correct momentum theory from the right case to the skewed case. Two approaches may be chosen for the BEM implementation of this yaw-model [3]. In this study, it is first chosen to use a single vortex cylinder representative of the rotor. This is the option chosen when using Coleman or Glauert yaw-models. The alternative approach, consisting in using a superposition of concentric cylinders [11], is discussed in section 5.
The values of u z,t and u ψ,r for χ = 0 were seen to be related to the AD momentum theory inductions. Using Equation 3 and Equation 4, the ratios between the induced velocities for a given χ and for χ = 0 provides the following factors: The above factors can be applied to the momentum theory inductions. They will imply a redistribution of the inductions over the AD without changing their azimuthal averages. The induced velocities u z,r , u ψ,t and u ψ,l are added to aU 0 R z,t and a ΩrR ψ,r to form the total induced velocities according to Equation 2.
To further close the system, the vorticity distributions γ l , γ t and Γ r need to be computed, which in turn require the knowledge of Γ tot , h and χ. The skew angle can be determined using different methods (see e.g. [12, p. 99]), it is here assumed given. The value Γ tot is taken as the sum of the average blade circulations: where V rel,⊥ , c, C l are respectively the norm of the relative velocity in the airfoil plane, the chord and the lift along each of the B blades. Alternatively, under the assumption of constant circulation and high tip-speed ratio, λ, the circulation of the rotor may be obtained from the total thrust coefficient over the rotor, C T , as: Γ tot = πRU 0 λ C T . Equation 9 is yet more precise. If the axial and tangential induction in the wake are neglected, the pitch of the helix is simply the distance run by a tip-vortex in one rotor rotation: h = 2πR λ cos χ. In general, the pitch of the helix is related to the inflow angle distribution on the rotor [13]. Assuming high tip speed ratio (a = 0) and 1D momentum theory other measures of the pitch can be obtained [14].

Application and results
Measurement data and numerical tools The MEXICO rotor [7] will be used for application of the model. Data from the test case 2.1 of the IEA Task 29, first post-processed by Pascal [15] and finalized by Schepers et al. [16], will be used in the current study. The coordinate convention for the experimental setup is found in [17]. The wind speed is U 0 = 15 m/s, the yaw angle is θ yaw = −30 • using the sign convention of Figure 2, the blade pitch is −2.3 • [17], the rotational speed is Ω = 424.5 RPM and the air density ρ = 1.237 kg/m 3 . The tip-speed ratio is λ = 6.67 and the resulting thrust coefficient is C T = 0.75.
Results from the BEM yaw-model are compared with numerical simulations from geometryresolving CFD, AL-CFD and a free-wake vortex-code. Results from the in-house flow solver EllipSys3D [18,19] were presented in a previous study by Sørensen [20] using geometry-resolving CFD. An AL simulation of the rotor using EllipSys3D [21] was performed for the current study using a Cartesian structured mesh of 14 millions cells. The extent of the domain in each direction is (16R, 16R, 16R) and the finest cell size at the rotor is equal to R/55. The computational time step is chosen as dt = 6.7 × 10 −4 R/U 0 . The vortex code simulations are performed using the lifting-line formulation of the Omnivor library [22]. The BEM-, AL-and vortex code use the same airfoil data [23]. Vortex code simulations were run with 40 sections per blade, 50 time steps per revolution (dt = 0.02R/U 0 ) and 35 rotor revolutions (t = 33R/U 0 ), resulting in 4 × 10 5 vortex elements. The BEM code used is part of the aerodynamic module of hawc2 [24]. The angle χ is determined from the free-wake code by evaluating the mean slope of the wake axis between z = 0 and z = 6R. The value 6R is chosen since previous studies showed reliable performance of the vortex code compared to AD simulations up to this distance [22]. The skew angle found, χ = −36 • , is temporarily prescribed in the BEM code implementation. This value is seen to be consistent in this case with the model χ = 1.2θ yaw [25]. The slope was determined using the vortex-code and not the CFD results since the former will be later used as a reference.
The performance of the different codes are evaluated by comparison with the measurement data in Figure 3. The BEM code labelled "Coleman/Glauert", only accounts for the  redistribution of axial induced velocities using the factor R z,t . It corresponds to a standard BEM code yaw-model implementation. Offsets with the measurements were observed in previous studies [17,20,26]. Priority should therefore be given to the azimuthal trends when comparing with the measurement data. Yet, the three codes that rely on airfoil coefficients can be compared directly. The mid-part of the blade is well described by all codes. In particular, the CFD and vortex codes capture a linear trend in F nc between 90 • and 200 • which is present in the measurements. The vortex and AL code, which both use the same airfoil coefficient data, are in reasonable agreement at the blade root and in excellent agreement elsewhere. Complex aerodynamic effects take place at the root and affects the airfoil performances so that this specific location will be treated in a dedicated paragraph. Beyond these effects, the causes for the differences between the AL and vortex code results at this location are currently under active investigation. The vortex code will be used as a reference to evaluate the performance of the new yaw-model since the latter is also vortex-based. Results using the new yaw-model were included in Figure 3 and labelled "BEM-new". Further investigation of the new model will follow and its overall performance will be evaluated. Correlation and least-squared measures will show that the performance of the baseline BEM code can be further improved using the new yaw-model.
Application of the new yaw-model Velocities and loads obtained using the yaw-model described in this paper are compared with the vortex code results in Figure 4. The new model performs best when the corrections to the inductions are applied outside of the BEM loop. The two implementations are referred to as "BEM new(out.)" and "BEM new(in.)". In Figure 4, the concentrated root vortex, u z,r , is seen to have a strong influence on the axial velocity. It is expected that an implementation using a superposition of vortex cylinders will moderate the influence of the root vortex (see section 5). The terms acting in the tangential direction were all observed to give a significant contribution to u ψ . They can be listed from the most influential to the least as: u ψ,t , R ψ,r and u ψ l . The amplitude of u ψ with the new model is in good agreement with the one obtained from the vortex code due to the contribution of u ψ,t . The induced velocity in the tangential direction is nevertheless small compared to the wind contribution U 0 sin θ yaw cos ψ.  Inboard section The effect of dynamic stall and stall delay were not included for the results of Figure 4. Results from the vortex code and the new yaw-model with dynamic stall [27] are compared to the measurements in Figure 5 for the inboard part of the blade. It is possible to assess the part of the dynamic stall which is due to the vorticity shedding by comparing with results from the vortex code without shed vorticity 1 . The absence of shed vorticity modelling in the BEM code was likely to be the cause for the increased loading seen around the azimuth 270 • . It was originally thought that the increased loading towards the azimuth 90 • could be explained by the root vortex. From Equation 4 and Figure 4, it is seen that this effect cannot be explained simply by a concentrated root vortex. Since the effect is captured by the vortex code, it is possible that a superposition of vortex cylinders will also capture the extra loading occurring around 90 • . This will indeed be confirmed in the following section. At the inboard part of the blade, the results for the new yaw-model are clearly in better agreement with the measurements when the effect of dynamic stall is modelled. Apart from this area, the effect of the dynamic stall is not as pronounced.  Figure 5: Effect of dynamic stall on the normal load at r/R = 0.25. Dynamic stall modelling for both BEM and Vortex code is necessary to capture the aerodynamic performances in the in-board part.

Relaxation of assumptions and extension of the model
The model presents four main assumptions which are: an infinite number of blades, a rigid wake, a constant circulation over the radius, and a constant circulation over the azimuth. The relaxation of the three first assumptions are presented here. The vortex code Omnivor with its prescribed-circulation and no-shedding formulation is used to study the relaxation of the two first assumptions.

Wake distortion
The theoretical results of the model were derived for an infinite number of blades. A number of blades of 10 was enough to reproduce these results within more than 99% accuracy using the vortex code without wake relaxation. The effect of wake distortion computed with 10 blades is shown in Figure 6 for the axial velocity. The wake distortion introduces phase and amplitude shifts along the blade which could potentially be modelled in a BEM extension. The rigid-wake assumption appears justified to a first order approximation since the main azimuthal and radial trends are captured. Figure 7 where the velocities are normalized by u z,0 and u ψ,0 = Γ tot /4πr. The yaw-model derived for an infinite number of blades can be applied with sufficient accuracy to three-bladed rotors. Figure 7 shows the induction that can be obtained if the tip-loss factors from Prandtl(Pr.) or 1 The simulation is not physical but gives insight on the shed vorticity contribution.  Glauert(Gl.) [14] are used to correct the yaw-model for a number of blade of 1. From the figure, the Glauert formula appears satisfactory to model the effect of finite number of blades on the axial induction but over-predicts the changes in tangential induction. This over-prediction of the tip-loss is consistent with previous results from the authors [14,13]. Radially varying circulation The effect of radially varying circulation is added by using a superposition of concentric vortex-cylinder systems. This approach was used by Heyson in 1956 [11] and is here extended to comprise all the vortex components. Each system consists of a root vortex, a bound vortex disk and a cylinder of radius r i . The upperscript is used to index the concentric cylinders and subscripts will be used to index the control points which are assumed to lay in-between two successive cylinders. To a first approximation all systems can be considered to have the same skew angle. Different pitch values should be used for each system to ensure the proper convection velocity of each vortex system [6], such that:

Finite number of blades Computations performed for a lower number of blades are shown in
A constant pitch value simplifies the derivations but violates the kinematic condition across the vortex cylinders. The intensity of each vortex system, notedΓ i , is determined from the bound circulation according to the decomposition presented in Figure 8. The total induction from all the vortex systems for the -component at a control point j is:  Similar to the implementation presented in section 3, the momentum theory inductions are corrected by comparing the inductions from the right and the skewed cylinder models. The velocities induced by a superposition of right cylinders with constant pitch values reduce to [6]: where the first equality of each term is non-trivial but results from properties of the right cylinder and a sum over theΓ i . The correction factors applied to the momentum theory inductions are then: The influence of the other components are added to the total induction (see section 3). For simplicity the influence of the root vortex may be computed at once. Indeed, the sum over all the root vortices of intensities −Γ i is −Γ 1 . Further, in the common case where the circulation is zero about the root of the blade due to the use of cylindrical profile, the few first inner cylinders will have no influence. In particular, Γ 1 = 0 and the root vortex can be neglected for all vortex systems. The formulae presented in section 2 can directly be applied to compute the induced velocities inside of each cylinder. Unlike the right vortex cylinder, the induced velocities are not zero outside of the cylinder and models should thus be derived for this condition. Preliminary results using only the inner induced velocities are shown in Figure 9 to illustrate the applicability of the method. The error introduced by using only the inner velocities is expected to be small closer to the root since few cylinders have a radius lower than the control point at this location. The comparison of the constant-circulation yaw-model and the varying-circulation yaw-model reveals a clear improvement when using the second model since this model shows closer agreement with the vortex code results. The increased normal loading around the 90 • -azimuth is well captured by the model. Yet, further investigation using both inner and outer velocities are required to fully assess the performance of the model.

Conclusions
The implementation of a new yaw-model based on recently developed vortex-theory results was presented. The model relaxes the assumption of infinite tip-speed ratios used in current BEM yaw-models. Three tools using the same 2D airfoil coefficient data were used in this study. Good agreement was found between the actuator-line and the vortex code simulations. The latter was used as a reference to validate the new vortex-based yaw-model. The induced velocities and loads obtained with the new BEM yaw-model showed better agreement with the vortex code  Figure 9: Normal and tangential loads with respect to the chord line for r/R = 0.25. The yaw-model using a superposition of vortex cylinders shows better agreement with the vortex code.
than the standard yaw-model implementation. The influence of a single root vortex was not seen to be the source of the shift of loading towards lower azimuth at the inboard part of the blade but this shift was captured using a superposition of cylinders. The new yaw-model combined with the effect of dynamic stall showed reasonable agreement for the inner part of the blade which is usually difficult to model. The limitations due to the assumptions of infinite number of blades and rigid-wake were briefly discussed. Accounting for these effects do not present major difficulties. The assumption of rigid wake was seen to be justified to a first order approximation while the effect of finite number of blades can be consistently modelled using Glauert tip-loss factor. The assumption of constant circulation over the span was successfully relaxed and preliminary results showed improved accuracy towards the root. Future work will further address this topic and investigate the relaxation of the assumption of constant circulation along the blade azimuth.