{'fraction_non_alphanumeric': 0.02217741935483871, 'fraction_numerical': 0.0065609622744669215, 'mean_word_length': 5.004308576118178, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic model for aerodynamic force dynamics on wind turbine blades in unsteady wind inflow

The paper presents a stochastic approach to estimate the aerodynamic forces with local dynamics on wind turbine blades in unsteady wind inflow. This is done by integrating a stochastic model of lift and drag dynamics for an airfoil into the aerodynamic simulation software AeroDyn. The model is added as an alternative to the static table lookup approach in blade element momentum (BEM) wake model used by AeroDyn. The stochastic forces are obtained for a rotor blade element using full field turbulence simulated wind data input and compared with the classical BEM and dynamic stall models for identical conditions. The comparison shows that the stochastic model generates additional extended dynamic response in terms of local force fluctuations. Further, the comparison of statistics between the classical BEM, dynamic stall and stochastic models' results in terms of their increment probability density functions gives consistent results.


I. INTRODUCTION
Wind energy being safe, significant and fundamental for economical as well as social development is getting great interest and successfully penetrating the energy market today. According to Global Wind Energy Council 2012 projections, the power generation from wind is expected to reach twice the present global installed capacity by the end of 2017 [1].
In the success of modern wind energy, aerodynamic research has a significant role [2]. The aerodynamic models for wind turbines are used to study the acting loads on rotor blades in response to wind inflow. At present, various engineering as well as computational fluid dynamics (CFD) models exist to foresee the performance of wind turbines. In computational terms, the investigations suggest wide choice of engineering methods in particular the well-known blade element momentum (BEM) method [3,4]. The high fidelity CFD is yet extremely costly [3] and needs faster and bigger memory computers to achieve acceptable computational efficiency [5]. For distributed loads on wind turbine rotor blades, the aerodynamic models are combined with dynamic analysis codes such as FAST [6], YawDyn [7], ADAMS/WT [8], SIMPACK [9], DHAT [10], FLEX5 [11] etc. The aerodynamic models used by these dynamic analysis codes or other similar codes are fundamentally based on simple lookup tables, which contain mean static characteristics for an airfoil at constant angles of attack (AOAs) [12,13] and disregard the information on system local dynamics. Even the dynamic models such as Beddoes-Leishman dynamic stall model use static airfoil coefficients which are modified according to AOA and its rate of variation and mostly disobey the local force dynamics.
In this contribution, a new concept based on a stochastic approach has been integrated into the aerodynamic model AeroDyn [12] as an alternative to traditional table lookup method used by the classical BEM model. The concept represents a stochastic lift and drag model, which provides the lift and drag forces with local dynamics under unsteady wind inflow conditions. The model estimates the lift and drag coefficients numerically based on the local AOA [14]. The proposed approach thus is a stochastic alternative to the classical BEM and Beddoes-Leishman dynamic stall models.
The scope of this paper is to prove the concept by integrating the stochastic model into AeroDyn and showing that the newly developed concept extracts more load information on rotor blades compared to traditional approaches. A future aim is to achieve a complete stochastic rotor model, which could provide the full local loading information on blades in a stochastic sense, leading to an optimum rotor design. Such aerodynamic model could be combined with a wind energy converter (WEC) model to obtain a stochastic rotor model. It is necessary to mention here that this contribution is a primary step towards the final goal and for simplicity reasons disregards the tower shadow and other related effects at this stage.
The paper is structured such that Section II provides a short introduction to the AeroDyn, the input files and the elemental forces. Section III introduces the stochastic lift and drag model. Section IV presents the stochastic model integration into AeroDyn and the results achieved by classical BEM, dynamic stall and stochastic model. The final Section V summarizes the outcome and outlook of the work.

II. AERODYN
AeroDyn is an add-in software consisting a set of routines to execute aerodynamic computations for horizontal axis wind turbines (HAWTs). It can be interfaced with number of dynamic analysis codes in particular with FAST, SymDyn, YawDyn and ADAMS/WT to carry out the aeroelastic simulations. These aeroelastic simulation codes differ only in structural dynamics, the aerodynamic calculations are identical for them. AeroDyn has no stand alone executable functionality and is invoked by a dynamic analysis code. However, its separate aerodynamic forces output file for any element can be obtained by setting an option from "NOPRINT" to "PRINT" in AeroDyn primary input data file [15]. It creates one output file in one simulation time for one selected blade element only.
AeroDyn requires information about wind turbine geometry, blade element velocity, element location, airfoil aerodynamic data, operating conditions and wind input [12]. Based on given information, it computes the corresponding elemental aerodynamic forces and delivers to the aeroelastic simulation program to estimate the distributed forces on wind turbine blades. AeroDyn uses different models to perform aerodynamic calculations for aeroelastic simulations of HAWTs; however, for current computations, the BEM and dynamic stall models are used. The BEM model is a well-known classical approach used by different wind turbine designers with various corrections, whereas the dynamic stall model is based on the semi-empirical Beddoes-Leishman model which is especially important for yawed wind turbines. A detailed description of the classical BEM and dynamic stall models used for present computations can be found in Moriarty and Hansen [12]. The BEM model is used for estimation of the steady forces (with static table lookup approach) and the stochastic forces (with stochastic model addition).

A. AeroDyn input files
AeroDyn requires three input files to perform the aerodynamic calculations. These are, the primary data file, the airfoil data file and the wind file. Latter two are called through paths provided in primary data file. The primary data file consists of different options of models for calculation of flow influence. Additionally, it carries the information on wind reference height (hub height), convergence tolerance for induction factors, air density, kinematic air viscosity, time interval for aerodynamic calculations, and the number of blade elements per blade. For each blade element it is provided with the location, twist, span and chord length. For more clarity of the element related parameters; see blade segment terminology given in Figure 1.
The airfoil data file contains two dimensional static airfoil characteristics. It carries lift, drag and pitching moment (optional) coefficients for range of AOAs. Besides, it comprises some additional parameters pertaining to dynamic stall model. The airfoil characteristics for intermediate AOAs are obtained by linear interpolation.
The wind files in AeroDyn are used in two formats, one the hub-height wind and other the full field turbulence, which are created by either measurements or simulations. The hub-height wind files are the simple ones containing either steady or time varying wind data. The full field turbulence is generated by TurbSim program [16,17], which creates two files, one the binary wind data file and other the summary file. The created wind data corresponds to all three wind components changing in time and space. The wind is sampled at frequency of 20 Hz and the turbulence is generated by a square grid spread over the whole rotor area. The velocity components at each point of the square grid are provided as function of time. The components of velocity at each blade element are obtained with linear interpolation (in terms of averaging or smoothing) in time and space [12]. For more details; see [16,17].

B. Elemental forces
To determine the elemental forces on a rotating blade, the BEM method is applied on blade section as shown in Figure 2. The illustration shows the local velocities with flow angles and the acting forces on the element.
In Figure 2(a), φ represents the flow angle, V Total the relative speed, α the AOA and β the combination of pitch and twist angles. The flow angle φ is the angle between the relative speed and the plane of rotation, whereas the AOA α is the angle between the relative speed and the chord of the blade element. The parameter U ∞ denotes the free stream wind velocity, Ω the blade rotational speed and r the local radius of the blade element. The variables v e−op and v e−ip are the out-of-plane and in-plane element velocities, respectively, originating from blade structural deflections under pronounced rotation. When the blade rotational speed is very small, the latter velocities are ignored.
The terms U ∞ (1 − a) and Ωr(1 +á) are the effective axial wind and tangential blade speeds, respectively. The parameters a andá are the axial and tangential induction factors, where a represents the amount of reduction in axial wind speed when approaching the blade andá the amount of rotational acceleration to the blade caused by induced wake rotation opposite to the rotor rotation [18].
The induction factors are estimated using an iterative process described in the flow chart Figure 3. After initialization, the algorithm iteratively finds values fulfilling the condition expressed for Tol. The parameter Tol is an acceptable tolerance allowed around the true values of axial and tangential induction factors. The process repeats for each element. Once the induction factors, inflow angles and AOAs converged to their final values, the acting forces are estimated for the element.
In flow chart Figure 3, the function F represents the tip and root loss correction factors in combination, C T the thrust loading on the element, σ the local solidity of the rotor, C n the normal force coefficient and C t the tangential force coefficient. The parameter a skew is the local element induction factor for skewed wake, R the blade radius, χ the wake angle and ψ the blade azimuth angle [19]. For more details and derivations of the relations given in flow chart Figure 3; see [12,18,[20][21][22][23][24][25][26][27].

III. STOCHASTIC LIFT AND DRAG MODEL
A stochastic model of the lift and drag dynamics has been developed for an advanced characterization of forces on wind turbine airfoils under unsteady wind inflow. The model parameters are derived from dynamic measurements of lift and drag forces performed for FX 79-W-151A airfoil in a wind tunnel. The turbulent inflow with intensity of 4.6% was generated using a fractal square grid, which produces the flow with typical intermittent velocity fluctuations commonly known for free field wind situations; see [14]. The model extracts most of the information available in the system dynamics in terms of a dynamic response. The model reads [14] X model (k) = where the first part of the equation represents the basic model and the second the extension which accounts for the oscillation effects with amplitude modulation (breathing) contained in the lift and drag time series. The parameter X represents the lift and drag coefficients, A the constant to fix the oscillation amplitude for lift and drag coefficients, k the discrete time variable and T the most dominant oscillation period of lift and drag coefficients. The exponential function in equation (1) controls the breathing of oscillation along the lift and drag time series, where k o is half the average breathing length, k = (k mod k o ) and S is described as where n = 0, 1, 2, . . ..
The basic model X langevin (k) in equation (1) corresponds to a first order stochastic differential equation termed as Langevin equation, cf. [28], expressed in discrete form as where D (1) (X, α) and D (2) (X, α) are the drift and diffusion functions, also known as first and second Kramers-Moyal coefficients. Here τ is the integration time step and α the mean AOA that varies slowly compared to the fluctuations caused by turbulent inflow conditions. The parameterF is the correction factor for diffusion function to incorporate the model extension and Γ(k) the Gaussian white noise termed as Langevin force [28] with mean value Γ(k) = 0 and variance Γ 2 (k) = 2.
The D (1) (X, α) reflects the deterministic part of the system, whereas D (2) (X, α) quantifies the amplitude of the stochastic fluctuations. These functions for lift and drag coefficients are parameterized as [14] D (1) Here m is the slope of drift function, X the lift or drag coefficient, X o the stable fix point in X where drift function is zero and β the constant diffusion function. The optimized values for these parameters as function of mean AOAs are given in Tables I and II. The intermediate values can be obtained by linear interpolation. For more details and validation of the model; see [14].

IV. MODEL INTEGRATION INTO AERODYN
The stochastic lift and drag model described in Section III is integrated into AeroDyn to obtain aerodynamic forces with local dynamics on the rotating blade. The model is added in form of a routine based on model equation (1) as an alternate to static airfoil data files used by the classical BEM model in AeroDyn. The characteristics of model related parameters given in Tables I  and II are imported as input files to the model routine. The axial and tangential induction factors are estimated (using mean lift and drag coefficients) through an iterative procedure described in Figure 3.

A. Numerical setup
Before starting the aerodynamic computations, the AeroDyn input files described in Section II A are set up Assume initial values for a andá Calculate flow angle by φ = Calculate new a andá by   The aerodynamic calculations are performed for a threebladed rotor having radius of 13.76 m. The blade is divided into 10 equal segments along the radius. The local de-sign specifications for each element are given in Table  III containing the element nodal radius (from blade hub centre to the centre of element), twist angle, span and chord length. All three blades have the same distribution yielding identical aerodynamics for same pitch angle -1 • . The elemental pitch can be obtained by adding the blade pitch angle to the element's local twist angle. The com- II. CD model parameters for AOAs 0 • to 30 • . The parameter m is the slope of drift function, CDo the stable fix point, β the optimized diffusion function,F the correction factor for diffusion function to incorporate the model extension, A the constant to fix the oscillation amplitude, T the most dominant oscillation period and ko half the average breathing length.  putations are carried out for all 10 elements; however, the results here are presented for element number 5 only to avoid the repetition of similar statistics.
Since the model is based on aerodynamic characteristics of an FX 79-W-151A, for comparison purpose the airfoil data file of AeroDyn is replaced with measured aerodynamic characteristics of this airfoil (for measurement details see [14]).
For wind input, the full field turbulence simulated binary wind data file is used. The file represents all three components of the wind vector created with TurbSim program. The components are variable in time and space obtained at each element each moment by subroutine interpolations. A summary of the meteorological parameters of the wind data file is given in Table IV.   It has to be stressed that the model parameters have been derived from wind tunnel measurements in fractal square grid generated stationary turbulence, while the simulation uses synthetic turbulent inflow according to Table IV. To compensate at least partially for the different flow situations, the diffusion function of the basic model equation (3) is multiplied with a correction factor obtained by dividing synthetic wind input interpolated hub-height turbulence intensity with fractal square grid generated turbulence intensity. Due to the high complexity of turbulence, this can nevertheless not ensure completely comparable flow situations.
As described in Section II, AeroDyn can not be operated independently, but has to be initiated by a dynamic analysis code. Here it is invoked by FAST, which uses a combined modal and multi-body dynamics representation. In FAST, the wind turbine blades and tower are modeled by applying the linear representation considering small deflections with mode shapes (degrees of freedom (DOF)) listed in Table V [6,29]. FAST collects the basic information from its primary input file containing the details of wind turbine operating conditions and the basic geometry. For additional information such as blade properties, tower properties, furling properties, wind time histories and the aerodynamic characteristics, FAST reads some supplementary files [6].
Once FAST calls AeroDyn and exchanges the information on the model status including elemental pitch and velocities, AeroDyn starts to compute the elemental aero-dynamic forces. The velocity components are expressed normal and tangential to the plane of rotation.

B. Results
Results are achieved for a blade element following the numerical setup described in Section IV A. The force calculations are performed using the classical BEM, dynamic stall and stochastic models according to procedure given in flow chart Figure 4. The TurbSim generated synthetic wind is used as an input to AeroDyn. Simulations are performed for 10 realizations of 10 minutes of wind input. The synthetic wind, containing irregular speed and irregular turbulence level, led to strong fluctuations in local wind components and AOA as expected; see Figure 5. The figure portrays the complex fluctuations of the axial velocity component experienced by the blade element and its effect on AOA dynamics; compare [30]. The variations in velocity component are proportional to the AOA, i.e., the higher the wind fluctuations, the higher the AOA variation. A pronounced oscillation at T = 1.13 s is visible in Figure 5, which possibly stems from boundary layer shear effects, as the tower-blade interaction was not included in the simulations [31].
The blade aerodynamics is function of the AOA, therefore, a change in AOA means a change in aerodynamic forces; compare Figure 5(b) with Figure 6. Figure 6 shows the resulting aerodynamic forces behavior for the selected element, where the force coefficient signals obtained with classical BEM represent the mean force dynamics as expected, being dependent on mean aerodynamic characteristics of the airfoil. The force coefficient signals achieved with the dynamic stall model represent the force dynamics with small fluctuations around the mean, while the force coefficient signals contributed by the stochastic model represent the force dynamics with extended local dynamics around the mean compared to the dynamic stall model. The means of the dynamic stall and stochastic models' local force dynamics along the signals match almost perfectly with the classical BEM model force signals.
In order to investigate the quality of the stochastic model results, their statistical properties are compared with the classical BEM and dynamic stall models' results in terms of increment probability density functions (PDFs). The increments in our case are the differences of an aerodynamic force over a specific time lag described as [32] The increment statistics are basically two-point statistics, which determine the nature of a parameter variation against a selected time lag. Using equation (6), the increment signals of the aerodynamic forces are derived for different time lags. Increment PDFs of the stochastic, dynamic stall and classical BEM models' results for different time lags are shown in Figures 7, 8 and 9, respectively. For quantitative comparison, further kurtosis and standard deviations of the increment signals have been estimated. The kurtosis is calculated using the expression  where γ 2 is the excess kurtosis, (x −x) 4 the fourth moment around the mean and σ 4 x the square of the variance of the probability distribution. The γ 2 = 0 resemble a normal distribution, γ 2 > 0 a distribution with sharp peak and long heavy tails, and γ 2 < 0 a distribution with round peak and short light tails. The excess kurtosis is a measure of the deviation of PDF shape from the normal distribution. Figure 7 presents the increment PDFs for stochastic model results, where all four aerodynamic force coefficients look similar at all three time lags except some differences at the tails. All PDFs almost resemble the normal distribution shape up to ±3σ; compare with added Gaussian fit. At the tails, the PDFs deviate from the Gaussian distribution and show intermittent behavior. To quantify for these shapes, the kurtosis and standard deviations of the increment signals have been estimated as given in Table VI. The estimations show slight to pronounced higher positive kurtosis for all four aerodynamic force coefficient signals at the selected time lags, meaning that the increment PDFs have slight to pronounced heavier tails. The main reason for this effect seem to be the non-linear lift and drag characteristics of the airfoil over AOA. The standard deviations of the increment signals in each force coefficient case at all three time lags are close to each other in magnitude; see Table VI. Figure 8 presents the increment PDFs for dynamic stall model results. In this case, the increment PDFs except those of the drag coefficient, resemble similar shape with minor differences at the tails. The drag increment PDFs present heavier tails compared to other three force coefficients. However, the increment PDFs of all forces correspond fairly to the normal distribution up to ±3σ (compare with added Gaussian fit) at all three time lags. The estimated kurtosis values given in Table VII indicate very slight to pronounced heavier tails for all force coefficient signals except almost no deviations in case of the lift, normal force and tangential force coefficients at higher time lags. The contribution of pronounced intermittency in the force signals is believed to stem from the same phenomenon as described in stochastic model case. The standard deviations of the all three increment signals in each force case depict increase with an increase in time lag; see Table VII.
Similarly, Figure 9 shows the increment PDFs for classical BEM model results. Here, except drag coefficient, the increment PDFs of other three force coefficients portray similar shape for same time lags. The increment PDFs of the lift, normal force and tangential force coefficients, resemble fairly the normal distribution shape up to ±3σ (compare with added Gaussian fit), whereas the drag coefficient correspond to the intermittent shape. The evaluated kurtosis given in Table VIII indicate strong heavier tails for drag coefficient increment PDFs at all three time lags. The other three force coefficients show slight to pronounced heavier tails for all three time lags except negligible lighter tails for higher time lag in case of the tangential force coefficient. The contribution of strong intermittency in the drag coefficient as well as slight to pronounced intermittency in other force coefficients can probably be addressed to the same phenomenon as in the stochastic model case. The standard deviations of the force coefficient signals follow the same behavior as in the dynamic stall model case above; compare Tables VII and VIII.

V. DISCUSSION AND OUTLOOK
The stochastic model of the lift and drag dynamics for an airfoil FX 79-W-151A described in Section III has been integrated to AeroDyn in the context of BEM wake model. The forces are obtained through the classical BEM, dynamic stall and stochastic models used by AeroDyn. For wind input, the full field turbulence simulated binary wind data is used. The forces are estimated for 10 blade elements; however, the results are analyzed here for one element only to avoid the repetition of similar statistics. With this ansatz it could be shown that local, short-time fluctuations of aerodynamic forces can be integrated as a stochastic model into a current WEC simulation tool.
The comparison of force time series given in Figure 6 show that the classical BEM model force coefficient signals represent the mean dynamics as expected, being dependent on mean aerodynamic characteristics of the airfoil. The dynamic stall model forces reflect small fluctuations around the mean, whereas the stochastic model introduces additional extended local force dynamics around the mean. The means of the local force dynamics achieved with dynamic stall and stochastic models seem to coincide with the classical BEM model force coefficient signals. Nevertheless, the stochastic model contributes additional extended dynamic response in terms of local force fluctuations, which is neglected by the classical BEM model fully and partially by the dynamic stall model.
The statistics of the stochastic and dynamic stall models presented in Figures 7 and 8   also resemble fairly a normal distribution up to ±3σ. Both the stochastic and dynamic stall models' force increment PDFs yield slight to pronounced intermittency due to the non-linear lift and drag characteristics of the airfoil over AOA. The magnitude of the standard deviations of the stochastic model force increment signals is significantly higher that of the dynamic stall model; compare Tables VI and VII. That is because the dynamic stall model forces possess smaller magnitudes of local fluctuations compared to the stochastic model. Moreover, the behavior of the standard deviations of the increment signals is different in both the stochastic and dynamic stall model cases. The standard deviations in the stochastic model case seem nearly stable for all time lags, whereas for the dynamic stall model case they are proportional to time lags. Similarly, the comparison of statistics between stochastic and classical BEM models given in Figures 7 and 9 also reflect good consistency except for the drag coefficient. The increment PDFs of the stochastic and classical BEM models' force coefficients resemble normal distributions up to ±3σ besides the BEM model drag coefficient. The intermittency in this range in the stochastic model drag coefficient is mostly covered by Gaussian fluctuations. The strong intermittency in the BEM model drag coefficient as well as slight to pronounced intermittency in other force coefficients in both the stochastic and classical BEM model cases can probably be addressed to the same reason described above. The magnitude of the stochastic model force standard deviations is larger than the classical BEM model forces; compare Tables VI and VIII. The reason is the extended local force dynamics in case of stochastic model. The behavior of the standard deviations of the increment signals between stochastic and classical BEM models is very similar as observed between stochastic and dynamic stall models.
As described in Section IV A the stochastic model parameters have necessarily been obtained under different conditions than the simulations carried out here. As a first compensation, the diffusion function has been scaled according to the interpolated hub-height turbulence intensity of wind input taken for the present computations. This can; however, not guarantee completely comparable flow situations. Therefore, the comparability of results between the stochastic model and the dynamic stall and BEM models is limited at the current state. For better quantitative comparison of stochastic versus dynamic stall and BEM models, in the future the stochastic model parameters have to be derived from measurements in more realistic flow situations at wind turbine airfoils. Also a comparison to load measurements at a real WEC should be performed.
Additionally, at the present state, the stochastic model could not generate the expected intermittency in the resulting forces, which is probably washed-out by the Gaussian noise term of the model. To improve this point, multiplicative noise may be introduced to achieve intermittency in the forces corresponding to intermittent properties of atmospheric flows, which is currently work in progress.
The stochastic model is being developed to extract and provide more complete local loading information on wind turbine blades, which could lead to an optimized rotor design under turbulent wind inflow. Here we have shown that the local force dynamics can be provided by a stochastic model integrated into existing BEM codes used by aerodynamic models such as AeroDyn. Further work will have to include more realistic model parameterization and corresponding experiments to allow for quantitative evaluation of the results.