Comparison of a Coupled Near and Far Wake Model With a Free Wake Vortex Code

. This paper presents the integration of a near wake model for trailing vorticity, which is based on a prescribed wake lifting line model proposed by Beddoes, with a BEM-based far wake model and a 2D shed vorticity model. The resulting coupled aerodynamics model is validated against lifting surface computations performed using a free wake panel code. The focus of the description of the aerodynamics model is on the numerical stability, the computation speed and the accuracy of unsteady simulations. To stabilize the near wake model, it has to be iterated to convergence, using a relaxation factor that has 5 to be updated during the computation. Further, the effect of simplifying the exponential function approximation of the near wake model to increase the computation speed is investigated in this work. A modiﬁcation of the dynamic inﬂow weighting factors of the far wake model is presented that ensures good induction modeling at slow time scales. Finally, the unsteady airfoil aerodynamics model is extended to provide the unsteady bound circulation for the near wake model and to improve the modeling of the unsteady behavior of cambered airfoils. The model comparison with results from a free wake panel code 10 and a BEM model is centered around the NREL 5 MW reference turbine. The response to pitch steps at different pitching speeds is compared. By means of prescribed vibration cases, the effect of the aerodynamic model on the predictions of the aerodynamic work is investigated. The validation shows that a BEM model can be improved by adding near wake trailed vorticity computation. For all prescribed vibration cases with high aerodynamic damping, results similar to those obtained by the free wake model can be achieved in a small fraction of computation time with the proposed model. In the cases with low 15 aerodynamic damping, (2 v r ) , where ω is the angular velocity, c is the chord length, and v r is the relative ﬂow speed. For an airfoil pitching harmonically about the three-quarter chord point, the error has been estimated by Madsen and Gaunaa (2004) to be 10% at k = 0 . 1 and 100 % at k = 0 . 8 , which for the NREL 5 MW reference turbine at rated wind and rotor speed corresponds to frequencies of about 1.2 and 9.8 Hz at 60 m rotor radius with a chord of 2 m. Except for the ﬁrst ﬂapwise and 25 edgewise bending frequencies, most relevant modal frequencies for modern blades are between these values, which shows that it is important to include a modeling of the unsteady circulation.

is introduced to model these unsteady induction characteristics due to load changes by pitch, eigenmotion of the blades, turbulent inflow and shear. The accuracy of the computations is improved due to the added aerodynamic coupling between airfoil sections through the trailed vorticity, alleviating the limitations of the BEM strip theory. Especially in cases with large radial load gradients, for example close to trailing edge flaps or other aerodynamic devices or close to the blade root and tip, the cross sectional coupling will lead to an improved prediction of the steady and dynamic induction. The addition of the near 5 wake model in an aeroservoelastic code has an acceptable effect on the total computation speed. An aeroelastic simulation with the wind turbine code HAWC2, (Larsen and Hansen, 2007;Larsen et al., 2013;Kim et al., 2013), of the DTU 10 MW turbine (Bak et al., 2012) in normal operation with turbulent inflow takes roughly 10% (30 aerodynamic sections) to 40% (55 aerodynamic sections) longer if the near wake model is enabled than if a pure BEM model is used.
The coupled model using the modified BEM approach for the far wake has been proposed by Madsen and Rasmussen 10 (2004) and extended by Andersen (2010). Further improvement has been presented by Pirrung et al. (2012), where an iterative procedure was used to ensure convergence and avoid numerical instabilities of the NWM. An application of the coupled model to estimate the critical flutter speeds of the NREL 5MW turbine (Jonkman et al., 2009) also including blades with modified stiffness, has been described by Pirrung et al. (2014), where the coupled aerodynamics model has predicted 4-10 % higher critical flutter speeds than the unsteady BEM model in the aeroservoelastic wind turbine code HAWC2. 15 In the present paper, the iteration procedure of the NWM used by Pirrung et al. (2012) is presented in more detail, as well as a method to compute the necessary relaxation factor during a simulation, removing the need for additional input or very conservative relaxation factors that are independent of spatial and temporal discretization and increase the computation time.
Further, the NWM is simplified to accelerate the computations with small loss of accuracy of the unsteady results.
The dynamic responses to pitch steps and prescribed blade vibrations are validated by comparing them to results from 20 the more complex free wake code GENUVP (Voutsinas, 2006). The focus in the pitch step cases is the dynamic induction response, while the prescribed vibration cases are evaluated based on aerodynamic work during a period of oscillation. It is found that the coupled aerodynamic model is capable of producing results that agree much better with results obtained from the free wake code than the unsteady BEM model in most cases, without a dramatic increase in computation time. The more accurate computation of aerodynamic work can have a considerable impact on the aeroelastic response in the case where the 25 total damping is close to zero, such as for edgewise vibrations. This paper is structured as follows: In the next section a short description of the NWM and a previous implementation of the coupling to a far wake model and shed vorticity model are presented. In Sect. 4, modifications to far wake and shed vorticity model are proposed to improve the interaction of these models with the near wake model and to increase the accuracy of the dynamic lift computation for cambered airfoils. This is followed by a description of the iterative procedure to stabilize the near 30 wake model in Sect. 5. A way of simplifying the NWM to accelerate the computation is presented in Sect. 6. In Sect. 7 the free wake panel code used for validation of the coupled near and far wake model is briefly described. The effects of the model modifications and results from the code comparison are shown and discussed in Sect. 8.

Original Model description
The structure of the previous implementation (Madsen and Rasmussen, 2004;Andersen, 2010) of the model is shown in Fig.   1. From the velocity triangle, denoted as V T , follows a geometric angle of attack (AOA) α QS and a relative velocity v r . An effective AOA α ef f is obtained through a 2D modeling of the shed vorticity effects, which is briefly described in Sect. 2.3.
This effective AOA is used to determine the aerodynamic forces and the thrust coefficient C T . The thrust coefficient leads to 5 a far wake induction factor a F W , requiring a coupling factor k F W as input. Section 2.2 contains the dynamic inflow model, using the weighting factors A 1 and A 2 , which is used to determine the unsteady far wake induction u F W,dyn .
Using this far wake induction, and the near wake induction from the previous time step, a new intermediate velocity triangle V T i is determined, with a new quasi steady AOA and relative velocity. These lead to the bound circulation Γ QS . The difference in Γ QS between adjacent blade sections, denoted as ∆Γ in the following, determines the trailed vorticity. In the next section it 10 is shown how the induced velocity W due to the near wake, which is added to u F W to obtain the total induced velocity u tot at each blade section, follows from the trailed vortices. The total induced velocity will then, in addition to the relative velocity due to blade motion and turbulence in the incoming wind, determine the velocity triangle after the time step ∆t.

Near wake model
The NWM enables a fast computation of the induction due to the trailed vorticity behind a rotor blade. The trailed wake can be 15 discretized into trailed vortex arcs from several positions on the blade, where each arc consists of a number of vortex elements.
The induction at a blade section due to each vortex element can be computed using the Biot-Savart law, but this computation is numerically expensive as the influence of each vortex element on the induction at each blade section has to be determined. Beddoes (1987) proposed to avoid these expensive computations by assuming that the trailed vorticity follows circular vortex arcs in the rotor plane and limiting the computation to a quarter rotation. In this quarter rotation, the axial induction dw from 20 a vortex element at a blade position is decreasing as the vortex element moves away from the blade, starting with a value dw 0 . This decreasing induction, following from the Biot-Savart law, is approximated by exponential functions: Figure 1. The previous implementation of the coupled near and far wake model, as described by Madsen and Rasmussen (2004) and Andersen (2010). The numbers in parenthesis refer to the equations in the following sections, [B] to the original model by Beddoes (1987). where Φ is a geometric factor depending on the radius from which the vortex is trailed and the distance between the vortex trailing point and the blade section where the induction is computed. The angle β determines how much the blade has rotated away from the vortex element. The numerically efficient trailing wake algorithm gives the induction W due to the trailed vorticity at time step i at a blade section s as: where N v is the number of vortex arcs trailed from the blade and W s,v is the induction due to a single vortex arc v at the blade section. It consists of components X i s,v and Y i s,v corresponding to both of the exponential terms in Eq. (1): where ∆Γ v is the trailed vortex strength, which depends on the radial difference in bound circulation between the blade sections adjacent to the vortex trailing point. The relative movement of the blade in the rotor plane during the time step at the vortex trailing point is denoted as ∆β v = (v r,in−plane /r)∆t. The in-plane velocity component perpendicular to the lifting line is denoted as v r,in−plane . Equations (3a) and (3b) show that the induction consists of a decreasing part of the induction at the previous time step, due to the previously trailed wake moving away from the blade, and the contributions from the newest 15 element, given by the D X,s,v and D Y,s,v terms. These equations are less time step sensitive and computationally faster than the original equations by Beddoes. They have been proposed by Pirrung et al. (2016), as well as a correction for the helix angle of the trailed vorticity. Pirrung et al. (2016) also describe how the tangential induction is computed based on the same approach.

Coupling to far wake model
The NWM, which only computes a fraction of the total rotor induction, is complemented by a modified BEM model for the far 20 wake. The total induced velocity at a blade section is computed as where u F W is the far wake component of the induced velocity and W is the near wake component, cf. Eq. (2).
The far wake component u F W is computed based on the BEM model implementation in HAWC2 that uses a polynomial to relate the local thrust coefficient with the axial induction factor at each annular element. To account for the near wake induction, 25 the far wake induction factor a F W is computed without a tip loss correction and based on a thrust coefficient C T that is reduced by the coupling factor k F W . The coupling factor is automatically adjusted during the computation to closely match the thrust of a reference BEM model, where the induction factor a ref is computed including tip loss effects (Pirrung et al., 2016). To account for the far wake dynamics, this work uses the dynamic inflow model implemented in HAWC2. Two parallel first order low pass filters are applied on the quasi steady induced velocities u F W,QS = a F W u ∞ from the BEM model: In a pure BEM computation and the previous far wake model implementation, the factors A i are A 1 = 0.6 and A 2 = 0.4. They are used to divide the induction into a faster and slower reacting part, corresponding to a faster time constant τ 1 and the slower time constant τ 2 . Both time constants are a function of radius and mean loading. The constants A i and τ i have been tuned to the actuator disc simulations of step changes in uniform loading (Sørensen and Madsen, 2006).

20
where the superscript i denotes the time step and c the chord length. Further, α QS is the quasi steady angle of attack resulting from the velocity triangle at the blade section and v r denotes the corresponding relative velocity.

Model overview
The structure of the current implementation of the coupled near and far wake model is shown in Fig. 3. The changes to the previous implementation, cf. Fig. 1 are: -The weighting factors A i of the far wake dynamic inflow are adjusted during the computation to account for the induction computed by the near wake model, which is explained in Sect. 4.1.

-
The trailed vorticity is no longer based on the quasi steady bound circulation Γ QS , but instead on a dynamic bound circulation Γ dyn . The computation of the dynamic bound circulation is shown in Sect. 4.2.1.
-The near wake induction is computed in an iteration loop, which is detailed in Sect. 5.
-The coupling factor is no longer needed as input, but instead continually updated during the computation, as described by Pirrung et al. (2016).

10
-The trailed vorticity is assumed to follow helix arcs to account for the downwind convection of the trailed vorticity. To achieve this, Φ is multiplied with a correction function f , depending on the blade section and vortex trailing point, as well as the helix angle at which the vortex is trailed (Pirrung et al., 2016).
-The computation of α ef f according to shed vorticity effects is improved for cambered airfoils, which is explained in Sect. 4.2.2.
The factors are continuously updated during the computations. A first order low pass filter with the far wake time constant τ 2 of the dynamic inflow model is applied on A 1,F W to make sure this model does not introduce unphysical rapid induction variations due to instantaneous changes of the weighting factors.

Unsteady circulation computation
The influence of shed vorticity on the bound circulation buildup has to be be considered when determining the strength of the trailed vortices of the NWM. Joukowski's relation between quasi steady lift L QS and circulation Γ QS , which has been used by Madsen and Rasmussen (2004) and Andersen (2010) to determine the bound vorticity, is not valid 20 for unsteady conditions. The error of calculating the circulation based on the unsteady lift at an airfoil section depends on the reduced frequency k = ωc/(2v r ), where ω is the angular velocity, c is the chord length, and v r is the relative flow speed.
For an airfoil pitching harmonically about the three-quarter chord point, the error has been estimated by Madsen and Gaunaa (2004) to be 10% at k = 0.1 and 100 % at k = 0.8, which for the NREL 5 MW reference turbine at rated wind and rotor speed corresponds to frequencies of about 1.2 and 9.8 Hz at 60 m rotor radius with a chord of 2 m. Except for the first flapwise and 25 edgewise bending frequencies, most relevant modal frequencies for modern blades are between these values, which shows that it is important to include a modeling of the unsteady circulation. In this paper, the step response of the circulation is approximated by the three term indicial function used by Madsen and Gaunaa (2004).
The algorithm is implemented analogue to the computation for the effective angle of attack in Equations (8)- (11): where the quasi steady circulation is computed from the quasi steady lift coefficient using Eq. (14).

Unsteady aerodynamics of cambered airfoils
10 Any change in bound circulation Γ, which is a function of v r C L , cf. Eq. (14) should lead to the corresponding shed vorticity.
The implementation of the shed vorticity model according to Hansen et al. (2004), cf. Equations (9-11) is based on the term α QS v r . The camber of the airfoil is neglected in this computation of the shed vorticity effects. We propose in this work to replace α QS in Equations (9 to 11) by α QS,camber , with    In the unsteady circulation computation described in the previous section, the camber is accounted for through the quasi steady circulation Γ QS , which is based on the lift coefficient, cf. Eq. (14).

Iteration scheme
The NWM can become numerically unstable depending on the time step, operating point of the turbine, blade geometry and 5 radial calculation point distribution (Pirrung et al., 2012). Fig. 7 shows the maximum time step where a stable computation is possible for a fine and coarse geometry definition, shown in Fig. 6, of the NREL 5 MW blade. The coarse geometry definition is a blade geometry typically distributed for BEM computations and the fine distribution is more suitable for computations with higher fidelity codes. The aerodynamic calculation points and vortex trailing points follow a cosine distribution, which means they are placed at equi-angle increments. The time steps have been determined in a numerical experiment, where the time step 10 has been decreased until large oscillations of the induction disappear. The results are accurate to the first significant digit. It can be seen that the finer blade geometry leads to a more stable computation. This can be explained by the smoother blade tip, where the blade chord is approaching zero. Thus the radial circulation gradient at the very blade tip is smaller and the vortex strength of the tip vortex is distributed to several weaker trailed vortices in the tip region that are less likely to cause numerical instabilities. In a coupled aeroelastic simulation, the small stable time steps for resolutions of 30 to 60 points would lead to a 15 very slow computation especially in case of the coarser blade geometry.
The numerical instability which occurs at larger time steps can be explained as follows: The axial induction due to trailed vortices typically reduces the angle of attack at a blade section, which in attached flow leads to a reduced lift. In the original  implementation of the NWM the constant circulation trailed during a time step is only depending on the flow conditions at the blade at the beginning of a time step. Thus a longer time step will lead to a bigger induction and thus a further reduction in lift in the next time step. If the time step is too large, the induction can become big enough to create a negative lift in the next time step, that is bigger in absolute value than the previous positive lift. This in turn leads to stronger trailed vortices of opposite sign, which will cause even bigger induced velocities in the opposite direction, which again leads to stronger vortices.

5
To stabilize the NWM the balance between trailed vortex strength based on the sectional circulation and the induced velocities are iterated to equilibrium in each time step, which removes the need for small time steps to stabilize the aerodynamics model. The iteration is structured as follows: 1 The quasi-steady circulation is computed according to Joukowski's law using the velocity triangle at the airfoil section based on the induction from the last iteration. The BEM model for the far wake is excluded from this iteration procedure. The AOA and relative velocity used to compute the far wake induction are the values from the converged iteration in the previous time step. This is accelerating the computation and is feasible because the near wake effects are on a much faster time scale than the dynamic inflow effects in the BEM model.

Estimation of the necessary relaxation factor
In the following, an estimation of the relaxation factor for a blade section is described. A conservative estimation is based on 5 the least stable case which is characterized by the following properties: -One single blade section with one vortex trailing from each side. Adjacent sections would tend to have similar circulations and therefore reduce the vortex strengths and the corresponding induction at the blade section. The trailed vortices on both sides of the section depend only on the bound circulation Γ of that section.
-The lift coefficient is linearly dependent on the angle of attack, C L = 2πα. A reduced but still positive gradient due to 10 stall would stabilize the model.
-No prior trailed vorticity is present. It would stabilize the model, because the induction would not only be determined by the momentary circulation at the section, but also by the decaying influence of the wake trailed before. If the model converges in the very first time step, with a given induction at the section from the previous iteration then the iterations will also converge with prior trailed vorticity.

15
-The helix angle at which the vortices are trailed is assumed to be small. Thus all the induction due to trailed vorticity is assumed to be axial induction.
With these assumptions, the downwash after a time step ∆t can be determined by summing up the contributions of the newest element, cf. the right terms in Equations (3a) and (3b), for both adjacent vortices:

20
where the subscript v denotes the vortex further inboard (v = 1) and outboard (v = 2) of the section with the bound circulation Γ. The subscript i denotes the iteration. Because the tangential induction is neglected, ∆β is only a function of the rotation speed of the turbine and the time step. Thus Γ is the only variable in Eq. (21) that depends on the induction at the section: where v ∞ is the free wind speed and the step response function from Eq. (15) evaluated at half a time step gives Γ dyn /Γ QS because we consider the first time step, thus a buildup of the circulation from zero.

11
Wind the iterative process converges. As seen in Fig. 8, the gradient of the curves is almost independent of W i−1 . The gradients are 10 negative because induction reduces the angle of attack. Therefore an approximation of condition (25) can be used: This gradient can be derived from Equations (21) and (24) as:

15
The gradient is mainly depending on the time step and point density (through B 1 and B 2 ) and the rotational speed. Instead of reducing time step and point density until a simulation is stable, which can lead to time steps orders of magnitude smaller than commonly used in aeroelastic codes and low spatial resolution, a relaxation factor f r can be introduced, so that: The derivative of this downwash with regard to the old downwash is: 5 For the minimum relaxation factor r, that allows for a stable computation (dW i,r /dW i−1 = −1), follows: which can be determined depending on the time step ∆t, the point distribution, and the number of points on the blade.
In the initial phase of the simulation, the maximum relaxation factor for all blade sections can be quickly determined by setting W i−1 = 0 in Eq. (27) and looping through the sections. The highest necessary relaxation factor for one section that 10 has been found is then used for the whole blade. As the simulation continues, the relaxation factor can be updated whenever there are big changes in rotational speed, induction, or blade pitch. If the relaxation factor is updated every several time steps, determining the relaxation factor takes negligible computation time. Choosing a slightly more conservative relaxation factor than what has been estimated will ensure stability also in different conditions than the ones the factor was based on.

15
In this section, an approach to accelerate the model is presented. The number of exponential terms used to approximate the decreasing induction with increasing distance from the blade in Eq. (1) is reduced to one. Using only one exponential term removes the Y w component in the near wake algorithm, Eq. (3b) and thus halves the computation time.
The reduced approximation function is defined as: 20 The values of A * and Φ * are found by solving the following equations: Equation (33)  into their new positions. The layout of the modelling is shown in Fig. 10. Details of the model can be found in Voutsinas (2006).
Since GENUVP is defined as a potential flow solver, the loads need correction in order to account for viscous effects. This is done by means of the generalized ONERA unsteady aerodynamics and dynamic stall model (Petot, 1989). The potential load  provided by the free wake model (Riziotis and Voutsinas, 1997).
In the case of aeroelastic coupling, the aerodynamic part will receive the deformed geometry and the deformation velocity and feedback the loading. The deformed geometry as well as the deformation velocities are introduced into the boundary 10 conditions and therefore the flow is accordingly adjusted.
The GENUVP free wake code has been thoroughly validated over the past years against measured data both on wind turbines and helicopter rotors in the framework of numerous EU funded projects. Blade loads and wake velocities comparisons against measurements have been performed on the MEXICO rotor in the context of Innwind.eu project (Madsen et al., 2015).
Moreover, detailed blade load calculations have been performed for the NREL test rotor and results have been compared to 15 experimental data (NREL experiment) and CFD computations (Chassapoyiannis and Voutsinas). Extensive validation of the code has been also performed in the framework of the HeliNovi project where aerodynamic and structural loads, wake velocities and elastic deflections have been compared to tunnel measurements on a BO105 helicopter model (Dieterich et al., 2005).

Results
In the following section, the effectiveness of the iteration procedure and the estimation of the relaxation factor are demonstrated for a horseshoe vortex. Then in Sect. 8.2 the unsteady induction predicted by the coupled near and far wake model is compared with results from an unsteady BEM model and the free wake code described in Sect. 7. Pitch steps and prescribed vibrations of the blades of the NREL 5 MW reference turbine are investigated. The comparison in Fig. 12 shows that the estimated relaxation factor is conservative, but the safety margin towards unstable computation is smaller in the 25 m/s case.

Comparison of the coupled model with a BEM model and a free wake panel code
The predicted force responses to pitch steps (Sect. 8.2.1) and blade vibrations (Sect. 8.2.2) are investigated in the following.
All computations use the refined blade model shown in Fig. 6. The blade has been discretized using 40 radial aerodynamic stations in the BEM based codes, with corresponding 41 vortices in the coupled model trailed from root, tip and in between stations. For the lifting surface free wake simulations, the blade has been discretized using 35 span wise and 11 chord wise 5 grid lines. Compared to the the faster models, the resolution was mainly reduced close to the blade root.

Pitch steps
Pitch steps with stiff blades have been performed, where the NREL 5 MW reference turbine is operating at a wind speed of 8 m/s and a rotation speed of 9.2 rpm. The turbine starts with blades that are pitched by 5 degrees to feather. The time steps are 0.054 seconds (120 steps per revolution) for GENUVP and 0.05 seconds for the other models. After 60 seconds simulated 10 time, the blades are pitched to 0 degrees at constant pitch rate in 1 or 4 seconds. The forces are normalized to compare the dynamics of the pitch response, such that the force before the pitch step is 0 and the force 45 seconds after the pitch step is 1. Figure 13 show the axial force response at a position at mid-blade and close to the blade tip for the fast pitch step. The free wake code predicts a slower force response during the pitch step than the BEM model. The results of the coupled model during the pitch step lie in between the other codes. In the free wake code results, some oscillations due to the changing wake The results of the slower pitch step in Fig. 14 show less oscillations of the free wake code results. In this case, the predicted results from the coupled model clearly agree better with the free wake code than the BEM results, both on the slope during the pitching motion and on the predicted overshoot.
Axial force distributions for a partial pitch comparison at 8 m/s are shown in Fig. 15. In this case, only the outer half of the blade is pitched from five to zero degrees during 1 second. As shown in the left plot of Fig. 15, the load smoothing effect of 5 the cross sectional coupling at the mid blade due to the trailed vorticity is predicted by both the coupled aerodynamics model and the free wake code to a similar degree. In the right plot of Fig. 15, the time history of the axial force is compared at 34.3 % radius. The constant force predicted by the BEM model on this non-pitching part of the blade is not included in this comparison. increasing after 1.6 seconds, corresponding to a quarter revolution at 9.2 rpm, while the force predicted by GENUVP continues to increase. Figure 15 thus illustrates that the coupled model can predict a change in loading that can not be computed based on BEM theory, but the restriction to a quarter revolution is limiting in this case. The difference in the overshoot prediction at 21.6 meters radius amounts to roughly 70 N/m.

5
The aerodynamic response to blade vibrations is investigated for normal operation at 8 and 25 m/s. The corresponding rotor speeds are 9.2 rpm and 12.1 rpm and the pitch angles 0 degrees and 23.2 degrees, respectively. The force response is compared in terms of radial distributions of aerodynamic work during one oscillation, where a positive aerodynamic work corresponds to a positive aerodynamic damping of the vibration. The mode shapes are chosen as the first and second structural mode shapes of the NREL 5 MW reference turbine blade at stand still, cf. Fig. 16. To simplify the comparison, the vibrations have 10 been prescribed as collective in-plane or out-of-plane vibrations. The frequencies, amplitudes and time steps used for the computations are shown in Table 1, as well as the modal masses that are used for damping estimations. For the first modes only results at the larger amplitudes are shown in the following. Investigating the results at the smaller amplitudes leads to the same conclusions.
In the BEM and coupled model, the blade section velocities due to the vibrations are added to the relative wind speed. The  been seen to that extend for the first edgewise cases, cf. Figures 17 and 19. A reason for this might be that the second edgewise case is computed with fewer time steps per period of oscillation to limit the computational cost.
For easier evaluation of the force response differences, the aerodynamic work can be expressed in terms of a damping ratio of a respective blade mode. Because the computations have been based on prescribed purely in-plane and out-of-plane structural mode shapes, these dampings do not correspond to any aeroelastic blade modes. For a single degree of freedom system with 5 the modal mass m and frequency f , given in Table 1, the damping ratio ξ and logarithmic decrement δ are: where A is the amplitude and W aero the aerodynamic work per oscillation period. The estimated logarithmic decrements according to Eq. (36) corresponding to the first flap motion at 8 m/s with an amplitude of 0.5 m, cf. Fig. 17, are 334 % for the BEM results, 300 % for the coupled model and 292 % for the free wake code results. Flapwise modes are highly damped and 10 thus these changes of the damping will not significantly alter the blade fatigue loads. On the other hand, lower aerodynamic damping of flapwise blade motion will correspond to a lower aerodynamic damping of tower fore-aft motion and might thus lead to increased tower fatigue loads. It is expected that this lower aerodynamic damping is balanced to some degree because the near wake effects reduce the aerodynamic excitation due to atmospheric turbulence.
Aerodynamic damping estimations for the first in-plane vibrations at 8 m/s at 1 m amplitude are shown in Table 2. The damping has been estimated in four cases, which differ in the airfoil polars and the modeling of the camber effect on the unsteady airfoil aerodynamics. Comparison of the first two cases of Table 2 shows that the induced drag caused by airfoil camber in the shed vorticity modeling results in a damping of roughly 0.7 % logarithmic decrement. According to the BEM and coupled model results in case (3) and (4) the airfoil drag increases the logarithmic decrement by about 0.3 %. The trailed vorticity decreases the absolute value of the damping by roughly 0.14 % log. dec. Further, comparing columns (2) and (4), the 20 combined influence of airfoil polars with lift coefficients other than 2π and drag is close to three times larger in the free wake code computations, which is caused by the different unsteady drag modeling. The small differences in estimated logarithmic decrement can have an impact on loads and stability computations for edgewise modes with a very low aeroelastic damping.
In the out-of-plane prescribed vibration cases investigated, the trailed vorticity reduces the aerodynamic work. Further a previous study by Pirrung et al. (2014) showed that the trailed vorticity effects will delay the onset of flutter towards higher 25 rotor speeds. This is in agreement with findings on the influence of shed vorticity, which leads to both a decrease of the flapwise damping and increased flutter speeds of a vibrating 2D blade section (Hansen, 2007).
(1) C ′ L = 2π, CD = 0, no camber (2) C ′ L = 2π, CD = 0 (3) NREL CL, CD = 0 (4) NREL CL and CD BEM -1. It has been shown that the acceleration of the model by reducing the number of exponential functions in the trailing wake 5 approximation from two to one is possible with negligible effect on the results. The approach presented here does not change the steady results predicted by the NWM.
An iteration scheme to stabilize the model has been presented. It applies a relaxation factor that is computed dynamically based on the blade discretization and the operating point of the turbine. To evaluate the computed relaxation factors, minimum necessary relaxation factors have been determined by trial-and-error and the estimated factors are found to be conservative.

10
The iterative process enables stable computations without the need for very small time steps and reduces oscillations of the near wake induction.
The 2D shed vorticity modeling, based on thin airfoil theory, has been extended by including the unsteady effects on the bound circulation. Further it has been found that it is necessary to include airfoil camber in the modeling of the influence of varying inflow velocity on the dynamic angle of attack to obtain good results if the direction of vibration is close to parallel to 15 the inflow direction.
A comparison of pitch step responses of the NREL 5 MW reference turbine using the coupled near and far wake model, a BEM model based on the aerodynamics model in HAWC2 and the free wake panel code GENUVP has been presented. The trailed vorticity modeling in the coupled model gives results closer to the free wake code than the BEM model during the pitching motion, and for a slow pitching rate a clear improvement is seen in the computation of the overshoot. Fast pitch rates 20 resulted in oscillations due to the motion of the wake in the free wake code, which could not be achieved in the coupled model due to the prescribed wake assumption. The response to a partial pitch of the outer half of the blade demonstrated the cross sectional aerodynamic coupling, which will have an influence on the load distribution in the presence of trailing edge flaps.
The coupled model agreed better than the BEM model with the free wake code in all prescribed vibration cases investigated.
The main improvement due to the trailed vorticity is found close to the tip of the blade, even in case of the higher modes 25 investigated. The work response to the edgewise vibrations has been found to be difficult to model if the direction of vibration is close to parallel to the inflow direction. The results in this case compare much better if no drag forces are computed. If drag is included, the coupled model still compares well with the free wake code close to the blade tip, but there are larger deviations in the results of all models further inboard. In general, the simulations agreed better for out-of-plane vibrations than in-plane vibrations.

30
The implementation of the coupled near and far wake model presented here delivers promising results and will be further investigated and validated against computational fluid dynamics results and measurements in future work. In particular the more accurate prediction of aerodynamic work for edgewise vibrations is considered to be important for stability analyses and load predictions due to the low aeroelastic damping typically associated with these vibrations.