Rotor-integrated modeling of wind turbine aerodynamics

A novel approach for the modeling of rotor-integrated aerodynamic loads is suggested to answer the need for a comprehensive, insightful, and analytical actuator disc model. All the six degrees of freedom (including tangential components) are considered. It is shown that loads may be written as a quadratic form of a reduced six-component velocity vector at the hub. The individual contributions of lift and drag, azimuthal variations, as well as blade pitching and tip losses are isolated. Errors introduced by the necessary approximations are discussed, and parametric corrections are considered. Parameter identification methods are then suggested, and the performance of the resulting calibrated analytical models is assessed. Results show that the new modeling approach is able to accurately model both the mean values of the thrust and power coefficients and their derivatives with respect to tip-speed ratio and pitch angle across the full range of operating wind speeds. Furthermore, it is able to reconstruct the general rotor behavior with a minimal amount of information available. Tangential components are also well modeled, although they require the knowledge of airfoil properties. The model's architecture leaves room for extensions to dynamic flow, skewed flow, and azimuthal load variations. presents an alternative approach for the modeling of rotor-integrated aerodynamics, providing the “ best of both worlds ” : the simplicity of empirical actuator disc models with the fidelity of frequency-domain-oriented AHSE models. It is insightful, compact, nonlinear, and analytical, obtained through a comprehensive re-integration of the underlying physics and a minimal set of information about the rotor as input. may valuable for multidisciplinary or multiscale couplings when rotor-integrated loads and unnecessary complexity is undesir-able (e.g., aerodynamics and physical control aerodynamics not core competence). Example applications may be Including aerodynamics in hydrodynamic studies of floating wind turbines (FWTs): impact of tangential loads, 8 higher order coupling between aerodynamic and hydrodynamic loads design optimization of floating foundation and mooring system, interpreting wave tank testing of farm aerodynamics for control design purposes, (2) with mid-fidelity models (e.g., dynamic wake meandering) with possibility for stochastic (e.g., Monte Carlo-based) simulations in real time, or (3) high-fidelity large eddy simulations (LES) models for multiscale coupling with the atmosphere. 14 It is expected that a physics-based model order reduction approach will enable further physics-based simplifications in these couplings. and improvements of parametric

Look-up tables based on pre-computed thrust and power coefficients 1 have been commonly used as an alternative. [2][3][4] While modeling is indeed simplified at simulation time, insight and short development time are still lacking as this empirical modeling method amounts to re-using pre-run simulations from AHSE analytical models. Analytical frequency-domain-oriented modeling is another alternative, through tools such as HAWCStab2, 5 TURBU Offshore, 6 or more recently STAS 7 offering as-accurate-as-possible frequency-domain models linking structural blade dynamics with aerodynamics, hydrodynamics, and beyond. 7 Being an analytical variant of AHSE models, they feature the same level of complexity, but rely on linear system theory and its powerful set of tools to handle it. This paper presents an alternative approach for the modeling of rotor-integrated aerodynamics, providing the "best of both worlds": the simplicity of empirical actuator disc models with the fidelity of frequency-domain-oriented AHSE models. It is insightful, compact, nonlinear, and analytical, obtained through a comprehensive re-integration of the underlying physics and a minimal set of information about the rotor as input.
This may be valuable for multidisciplinary or multiscale couplings when rotor-integrated loads are sufficient and unnecessary complexity is undesirable (e.g., when turbine-level aerodynamics is not the central topic) and when physical understanding is desirable (for instance, for control design or when aerodynamics is not the core competence). Example applications may be • Including aerodynamics in hydrodynamic studies of floating wind turbines (FWTs): impact of tangential loads, 8 higher order coupling between aerodynamic and hydrodynamic loads 9,10 design optimization of floating foundation and mooring system, interpreting wave tank testing observations.
• Control design of FWTs: understanding the effect of controls (coupling between blade pitch, thrust, torque, and rotor speed) and the circulatory relationship between velocity and load components for instance to address pitch control-induced instability, 11 using individual pitch control (IPC) 12 and nonlinear control schemes.
• Farm simulations: wind turbine wakes being driven by rotor-integrated loads, the current approach is sufficient and may be coupled with (1) analytical models of wake dynamics 13 for a complete analytical model of farm aerodynamics for control design purposes, (2) with mid-fidelity models (e.g., dynamic wake meandering) with possibility for stochastic (e.g., Monte Carlo-based) simulations in real time, or (3) high-fidelity large eddy simulations (LES) models for multiscale coupling with the atmosphere. 14 It is expected that a physics-based model order reduction approach will enable further physics-based simplifications in these couplings.
This paper is an in-detail revision of the approach suggested by Pedersen,11,15 and complementary information may be found there. A number of novelties are brought here, from physical/mathematical explanations (load matrix derivation process and properties), to application to other rotors, through modeling rectifications and improvements (effect of blade twist and tip losses), and more reproducible and tangible (i.e., physics-based) parametric modeling and parameter identification methods.
The paper is structured as follows: Section 2 presents generalities about the current scope. Section 3 starts from 6-degree of freedom (DOF) rotor-integrated flow, translates it to airfoil kinematics, and combines this with a parametric, integrable form of blade-element momentum theory (BEMT) to derive an actuator disc load model, giving insight about which velocity component affects which load component through which phenomenon. Section 4 looks into detail at uncertainties lying in assumptions used along the way through comparison with engineering model results. Based on this, parametric corrections are suggested in Section 5. Section 6 presents various methods for parameter identification based on rotor specifications and engineering model simulations. Results are then presented and discussed in Section 7.

| Actuator disc modeling
Let C τ be a hub-based 6-DOF dimensionless load vector including the thrust and torque coefficients C T and C Q , plus the shear force C Fy , C Fz and out-of-plane moments C My ,C Mz coefficients. Let ν be a hub-based 6-DOF dimensionless velocity vector (which includes the tip-speed ratio λ) and β the blade pitch angle. The overall goal is to derive a model of the form where f is an analytical continuous function. The approach may be divided in four steps as illustrated in Figure 1: (1) averaging the flow to a hub-based 6-DOF velocity vector, (2) re-characterize the flow at airfoil by means of airfoil kinematics, (3) adapt common airfoil aerodynamics (i.e., blade-element theory) to (4) integrate loads over the rotor. The loop is then closed by a rotor-integrated inflow model derived in Step (1) and a controller model feeding back the rotor speed and blade pitch angle.
Empirical models also follow this approach, but Steps [2] to [4] as well as inflow modeling are hidden inside the engineering model used to compute the load coefficients. The function f in (1) is then a set of pre-computed empirical values. To reduce the set dimension and hence the computational burden in practice, empirical models often use a linearization around the operating point n: providing a quasi-static model of aerodynamic loads.

| Scope
The problem may be divided into three parts: 1. Derivation of an analytical model of the form (1) and quasi-static analysis through comparison of its linearized form (2) with empirical models 2. Adaptation of flow averaging methods and assessment of their ability to capture wind turbulence 16 and dynamic and skewed inflow 15,17 3. Application to rigid-body response modeling of FWTs under wind-wave loading. Benchmarking the new model with empirical models, engineering models, and possibly other simplified models.
To allow for an in-detail study, the scope of this paper has been restricted to Point 1; the references given in Point 2 may be used to efficiently complement the load model. The following limitations apply: 1. The focus is set on low-frequency rotor-integrated aerodynamic loads. Higher frequency local processes (turbulence-induced radial load variations, elastic blade dynamics, blade sweeping-3p-frequency loads) are not considered.
2. Only operating conditions are considered. Although all velocity and load components are looked at, the mean flow is assumed purely axial. The vertical component arising from the rotor tilt is overlooked for simplicity; its impact on the validity of the various models is assumed negligible.
The effect of tangential velocities on axial loads is negligible and will not be treated, for the sake of conciseness rather than due to modeling limitations.
3. Only modern horizontal-axis wind turbines (HAWTs) are considered. This implies a control system of the variable-speed variable-pitch (VSVP) kind, with power production optimization in the below-rated wind regime and constant power and rotor speed in the above-rated wind speed regime. Nevertheless, no assumption is made on the upwind or downwind configuration, nor on the number of blades.

| Case study
Three reference turbines are used as examples throughout this paper. Working with dimensionless quantities enables comparing rotors of various sizes. In addition to the well-known NREL 5MW 18 and DTU 10MW, 19 the common research 20MW 20 (CR 20MW) turbine is added to the list of VSVP HAWT rotors under study. The three rotors have been designed independently, ensuring that no direct upscaling which would have skewed the comparison was used. Relevant rotor information is shown in Table 1.

| Coordinate systems
This paper makes use of four different coordinate systems illustrated in Figure 2: • Non-rotating rotor-plane coordinates (superscript h ): body-fixed and located at the hub with shaft-aligned South-West-Up orientation • Rotating rotor-plane coordinates (superscript r ): cylindrical system r, θ,z ð Þrotating with Blade 1.
• Airfoil coordinates (superscript a ): local system following the blade cross-section with z axis aligned with blade-pitch axis pointing outward, y axis aligned with chord line toward the trailing edge.
• Local flow coordinates (superscript f ): airfoil coordinates rotated about the pitch axis by the angle of attack (AoA).
Coordinate systems are linked through rotation matrices multiplied in cascade. Starting from non-rotating rotor-plane coordinates, azimuthal rotation by the angle θ leads to rotating (cylindrical) rotor-plane coordinates, as in (3). Airfoil coordinates are then obtained via the blade coning angle ϕ c , then the sum of the blade pitch β and twist angles ψ, as in (4). Finally, local flow coordinates are obtained via the AoA α, as in (5).
Blade coning and pre-bend have only a minimal effect on aerodynamics. 19 They are neglected in the following, greatly simplifying modeling.

| Rotor-integrated flow
Let v be linear velocities and ω angular velocities. The incoming flow is expressed in rotor-plane coordinates; the rotations by the tilt angle ϕ t in Figure 2 and possible floating platform motions are implicit and removed for clarity. The dimensionless flow vector ν is then defined as with U ∞ the undisturbed relative wind velocity, λ = ΩR U∞ the tip-speed ratio, R the rotor radius. Fluctuating wind and structural velocities are represented by σ x for the linear axial velocity component, by σ y,z , ζ y,z for tangential components, and neglected for the angular axial component. a a 0 ½ T isthe vector of induction factors (rotor-integrated inflow model presented in Figure 1), that may be modeled by the standard model of Pitt and Peters. 21 Steady-state momentum balance links the thrust coefficient with the axial induction factor a x =a and the angular moment coefficients with their respective tangential angular inductions factors: where B is a coefficient representing tip losses (standard, simplistic tip-loss model from helicopter theory 22 ). As in BEMT, an empirical correction for highly loaded rotors is applied. The in-plane angular inflow, commonly represented by the factor a 0 in BEMT, would read 1 and is neglected in Pitt and Peters model.

| Airfoil flow
Equations (6) and (7) linearize flow variations across the rotor disc and shorten its effective radius by a factor B to model tip losses. Under these approximations, the relative air velocity at the airfoil derives from the relative air velocity at the hub through airfoil kinematics. Kinematics are modeled by Jacobian matrices J linking linear and angular velocities using cross-product operations and coordinate transforms presented in Section 3.1. The translation from ν to blade-element relative velocities in rotating rotor coordinates v r a is given by (9): with I n the n identity matrix, Sð Þ the cross-product-equivalent matrix giving SðxÞy = x × y for any x, y ð Þ in ℝ 3 × 3 , and r r a ðrÞ the position vector from hub to airfoil in rotating rotor-plane coordinates (i.e., following the blade's azimuthal rotation).

| Airfoil loads
Once the airfoil flow is modeled, standard blade-element theory may be used to derive airfoil loads. The AoA α being given by the elementary aerodynamic load vector in local flow coordinates (airfoil coordinates rotated by α) reads: with ρ the air density, c(r) the radially varying chord length, and C D (r,α),C L (r,α) the drag and lift coefficients dependent on the radially varying airfoil profile and the AoA. See Figure 2 for an illustration. Note that, consistently with blade-element theory, v f a does not include the radial component of the velocity. The local aerodynamic moments dm f (only non-zero in z) are negligible for rotor-integrated loads.
Equation (11) is not convenient for analytical modeling as it is based on a look-up table of pre-computed airfoil coefficients. Lift may be instead expressed in terms of the circulation Γ through the Kutta-Jukowsky theorem. An equivalent quantity called dissipation Δ is created on purpose for drag (without physical interpretation). The relation reads Further, assuming that the rotor has been designed to optimize power production enables assuming that the lift coefficient varies linearly with the AoA. 1 Linear lift is equivalent to assuming the airfoil behaves as an ideal flat plate of unknown, equivalent chord length c e , 23 yielding where α 0 is the zero-lift AoA and ∂CL ∂α is the slope of the linear relationship. By means of trigonometric relations, we then get Regarding the drag coefficient, it may be approximated as independent of the AoA. Dissipation reads then with C D the average drag coefficient.
Note that these approximations are reasonable for small angles and hence also in above-rated wind speeds as blade pitching makes the AoA even smaller. At inner blade sections, however, demanding structural constraints and smaller contributions to the total torque (due to lower relative velocities) leads in practice to sub-optimal aerodynamic design, and the approximations become questionable. This will be treated in Section 4.2.

| Rotor-integrated loads
The rotor-integrated load vector derives from kinematic relations where N b is the number of blades, r 0 the radius of the blade root (i.e., hub diameter), and A the rotor area. Similar trigonometric relations to those used to derive (14) simplify the result of the multiplication of R r f by v f a , replacing the square root relationship in v f a by trigonometric ones.
These simplifications are conveniently obtained by the use of symbolic computing, and mathematical details will not be given here.
To establish the link between load and velocity coefficients using (16), a parametric azimuthal expansion of Γ and Δ is defined as Where lift is linear, (14) holds, and the identification of Γ ij simply reduces to higher order terms being zero. Regarding drag, Taylor expansion may be used to model Δ ij from (15): showing nonlinear variation with ν in addition to non-zero higher order terms. The same would be observable at inner blade sections where lift is not linear.
The blade pitch angle β may also be expanded azimuthally: with β 00 characterizing collective pitch control and higher order terms IPC. Assuming small IPC terms, Taylor expansion may be used include cyclic pitch variations into equivalent cyclic circulatory and dissipative variations b After setting a truncation order, the IPC-modified version of (17) and its equivalent for Δ may be inserted into (16) in a symbolic solver. After summing over the blades, C τ yields the following azimuthal expansion: The sought rotor-integrated load coefficient vector of the form (1) is finally obtained by extracting and simplifying the constant term 1 in (20).
It takes then a matrix product form, namely, with the Kronecker product operator, and the Hadamard (element-wise) product and division operators, respectively, and 1 m×n the m×n uniform matrix of ones. If the truncation order is set to 1 and the combination of IPC with tangential flow and drag is neglected, G reads 2   6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  5 ð23Þ which gives insight about which velocity and IPC component contributes to which load component through which phenomenon. As an example, looking at axial components only (first and fourth rows and columns in 23), we have showing that the cross-product-based contribution of global circulation Γ 00 is proportional to rotor speed for thrust and to axial relative wind speed for torque, and inversely for the dot-product-based contribution of global dissipation Δ 00 .
Note that as expected, the tangential load components show a symmetry relating the horizontal and vertical directions (second and third rows in 23 for in-plane forces, fifth and sixth rows for out-of-plane moments). For the sake of conciseness, only the horizontal components will be treated in the following, the extension to vertical components being straight forward. The notations C Fy = C F , C My = C M , σ y = σ , and ζ y = ζ will be used.
IPC is put aside for the rest of this paper.

| RELATION TO ENGINEERING MODELS
Rough approximations were introduced in Section 3. First in Section 3.3, the inflow (induced) velocity at each airfoil arose from a simple translation of a rotor-integrated value, and not from local momentum balance (or more advanced theories); tip-losses were given a simplistic model. In Section 3.4, the drag coefficient has been approximated as constant, and the lift coefficient as linear also at inner blade sections. Finally, azimuthal expansions have been truncated. To bring more insight with parametric modeling in perspective, these approximations are brought to light in this section and compared with engineering models for the reference turbines introduced in Section 2.3.
The inflow is calculated through a verified replica of FAST's BEMT algorithm 24 or variants of it; derivatives are computed numerically by finite differences (perturbation method) with fixed tip-speed ratio and pitch angle; steady-state values are computed using documented characteristics of the turbine controller for each reference rotor.

| Inflow modeling
As its name states, BEMT is based on local momentum balance for each blade element, dividing the rotor disc into elementary annuli. Using an iterative procedure based on a blade-element version of (7), BEMT gives elementary linear (also called axial, confusingly) a r and angular ar 0 (also called tangential) axial induction factors.

| Axial components
The rotor-averaged model in Section 3.2 assumed a constant linear induction factor a r ≈a over the blade and a negligible angular induction factor a 0 r ≈ 0 , which is a questionable approximation. To be consistent with BEMT, G should have been weighted by elementary induction factors, that is, To see the effect of the approximation, let us compare quantitatively the terms which link axial loads to axial velocities. The terms linking global circulation to tangential loads G 25 may also be included in the analysis. Results for all wind speeds and rotors are presented in Figure 3, showing a strong similarity. Regarding G 14 and the tangential terms, similarity was expected as the angular induction factor ar 0 is only significant-and still small-at inner blade sections, so weighting with r would remove any discrepancy. However, this is not applicable to 1− ar 1− a which, though approximately constant in the optimized power region, shows significant radial variations in outer blade sections at above-rated wind speeds, as seen in Figure 6. The similarity between G 41 and the other global circulation terms may therefore seem surprising. This is actually not a coincidence as it stems from energy conservation, circulation being loss-less by nature: links the circulatory power, thrust, and torque coefficients-respectively C Pc , C Tc , and C Qc -and their elementary versions (denoted with the subscript r). Knowing we obtain The same reasoning may be applied to cyclic circulatory components using conservation of cyclic energy. This result is important. Beside enabling modeling several components of G simultaneously, it makes the circulatory part of G skew-symmetric, implying strong physical and mathematical properties. It also enables modeling circulatory rotor-integrated loads independently on the circulation itself. As radial variation of inflow does not interfere in the rotor integration process, the circulatory terms of G will directly inherit modeling properties from Γ.
This analysis did not consider the effect of inflow variations in the modeling of the circulation itself. In (18), ν 1 should have been weighted by 1− ar 1− a and ν 4 by 1+ar 0 . While the latter has negligible effect, the former induces a slight error shown on the right-hand side of Figure 3 (for G 41 only, the other global circulatory terms behaving similarly), where linearized inflow means no radial weighting. Errors of the order +/− 5% are depicted.
Note the larger discrepancies at low wind speeds, where rotors are highly loaded. There, BEMT breaks down as the wake becomes turbulent.
Corrections such as Glauert's are used by engineering codes, but their extension from rotor to annuli is questionable due to turbulent mixing. On Effect of linearizing radial inflow variations and linearizing lift on global circulation the other hand, the accuracy of the simplistic tip-loss model in (7) is expected to worsen for highly loaded rotors. Invoking energy conservation there is therefore not legitimate. Given that low wind speeds are statistically less important to the overall turbine response, choice has been made not to put emphasis on discrepancies in these environmental conditions in the following.
A similar analysis is carried out for dissipation and shown in Figure 4. As expected, since energy conservation inherently does not hold, radial inflow variations induces significant modeling errors in the modeling and integration of axial global dissipative terms.

| Tangential components
Regarding tangential components, things are more complicated. By nature, BEMT is not capable of modeling non-uniform inflow (from yawed/ tilted inflow or non-uniform relative wind field), nor the effect of a finite number of blades, and even less in a turbulent wake state. 1

| Tip losses
Tip losses represent the effect of the finite number on the inflow (hub losses are neglected). BEMT is inherently not capable of capturing them, and Prandtl's correction is used as a standard. 1 As for a and a r , an elementary equivalent to B may be defined as and used in the iterative calculation of a r . Using Prandtl's theory, B may be related to an effective blade length R eff = RB, which yields B = 0:975 for an ideal blade design in optimal conditions. 1 In practice, no such direct relation can be used to derive B and a directly from B r and a r , but their general behavior should be similar. Figure 6 plots B 2 r at a selection of radii, showing a roughly linear trend against the inverse of the tip-speed ratio questioning the use of a constant value for B as done in helicopter theory. 22

| Airfoil coefficients
The slope of the linear lift coefficient versus AoA might be obtained from curve fitting. At inner blade sections, the slope is small, resulting in large values of −α 0 which violate the small angle assumption which is essential in the derivation of (13). An equivalent slope might then be derived as where α 0min is a threshold value (in the following taken as −8 degrees). α is a representative value of the AoA, see Section 6.1 for a practical estimation. Figure 7 shows

| General effect of dissipation on load coefficients
Crudely modeling aerodynamic coefficients may however not compromise the integrity of the modeling approach: • As the rotor design is optimized, the contribution of lift is expected to be far greater than that of drag. • Equations (21) and (23) show that aerodynamic coefficients are weighted by a power of r. This power is non-zero except for the relation between linear velocity and force components. Therefore, despite a chord length decreasing with radius, inner blade sections contribute only little to most load components.
It is well known that drag affects the mean value of the torque coefficient 2 and to a lower extent the mean value of the thrust coefficient at high wind speeds. The effect on derivatives is shown in Figure 8. In most cases, dissipative contributions do not exceed a few percents of circulatory contributions whenever the total contribution is of significance. The three exceptions are dCQ dλ , dCF dζ , and dCF dσ . Higher order terms in the azimuthal expansion of dissipation showed insignificant.

| PARAMETRIC MODELING
The approximations enlightened in Section 4 inhibit the direct use of relations (18) and (19) in (23) to compute G and hence any practical use of the model given by (21). This section aims at re-using the essence of these relations combined with parametric corrections informed by the results of Section 4, this for each component of G.  Correction candidates for (18) and (19) are linear weightings by the tip-speed ratio λ = − ν 4 and the induction factor a = 1−ν 1 . Their mean values are plotted in Figure 9. Corrections on terms involving in-plane velocities use the tip-speed ratio, while corrections on terms involving outof-plane velocities use the induction factor. Searching for a trade-off between fidelity and complexity (number of parameters), variants may be thought according to the need, and the identification procedure adapted. The following realization avoids any blade-element inflow modeling.

| Global circulatory terms
As a first step, we decouple the pitch angle β and the radially varying effective twist angle ψ+α 0 in (14). As seen in Figure 7, ψ+α 0 is naturally small in outer blade sections. β takes also only moderate values (not exceeding 30 degrees). Expanding (14) by trigonometric relations and neglecting insignificant terms enables the following approximation: This structure is preserved through integration. It enables separating the modeling of the circulation between a cosine and a sine term of the pitch angle, denoted with the subscripts C and S in the following. Looking first at the global component, F I G U R E 9 BEMT and AADT results for axial components Section 4.1.1 showed that the effect of radial variation in linear inflow should be corrected when integrating (32), as well as the effect of linearizing lift. The following parametric model is therefore suggested: with k Γ , k λ , and ℓ parameters. ℓ has not been made dimensionless on purpose. It may be seen as the radial location where an equivalent airfoil representing the entire blade may be defined. 11

| Global axial dissipative terms
Section 4.2 showed that modeling dissipative effects directly from (16) is challenging. Regarding G 44 , rotor integration includes weighting by r 2 .
Given that −ν 4 is from about 3 to about 12 times larger than ν 1 , it is reasonable to assume r 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ν 2 1 + r R ν 4 À Á 2 q ≈ r 2 ν 4 . This is confirmed by comparing trends of − ν 4 in Figure 9 and G 44 in Figure 4. It leads to with k ΔQ a parameter.
Regarding G 11 , the approximation does not hold. First, there is no radial weighting. Second, drag has a non-negligible effect on thrust only at high wind speeds, where the ratio − ν4 ν1 is lowest. Third, as seen in Figure 4, radial inflow variations accentuate the trend and rather indicate a proportional relationship with a, plus an offset. This offset appears to dominate the overall effect at high wind speeds: Keeping in mind that drag has only a limited effect on thrust, limiting complexity is preferred and a model based purely on the undisturbed wind velocity squared is suggested as with k ΔT a parameter.

| Tip losses
In light of the results of Section 4.1.3, the following model is suggested for the tip-loss factor: with B ? and k B parameters and λ ? the optimal tip-speed ratio.

| Cyclic circulatory terms
Beside for dCF dσ , direct integration of cyclic circulation from (18) does not induce much error when using the Pitt and Peters inflow model (see Figure 5). Still, a parametric correction enabling calibration against other inflow models is suggested as (splitting between sine and cosine terms as in 33): with k ζ C i , k σC i , k ζ S i , k σS i parameters, and γ ζ C i , γ σC i , γ σS i coefficients reading: πc e r i + 1 dr πc e r i dr: Note that direct integration is obtained when parameter values are zero, with results shown in Figure 5.

| Tangential dissipative terms
G 42 may be simplified using 1 R r 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi yielding a constant dG42 dσ as a fair approximation given the overall effect of dissipation. Regarding G 22 , Figure 4 indicates a relationship with the tipspeed ratio consistent with its tangential nature. The combined effect on the dissipative load coefficient variations may then be approximated as with k Fσ and k Fζ parameters. However, seeking for a minimal number of parameters, the dissipative effects may be included in equivalent circulatory parameters by relating (40) and (37) through (23): In order for the merging between dissipative and circulatory effects to be consistent, these equivalent parameters should be constant. Keeping in mind that circulatory effects are dominating, case studies showed that the relative variations across the wind speed range of interest are unlikely to exceed 20%. Considering the already present uncertainties in dissipation modeling, this is deemed reasonable. This is even more true for equivalent angular coefficients g k σC 2 and g k ζ C 2 as drag effects have been shown negligible there.

| Summary
Wrapping up the results of this section, we obtain to be used in (21). Note that G may be written as a sum of a skew-symmetric (circulatory part) and a diagonal (dissipative part) matrix, consistently with the properties of lift and drag. Note also that G is-with the exception of G 11 -linearly varying with ν, consistently with the quadratic nature of aerodynamic loads.
Regarding inflow, the quasi-static induction factors are given through momentum balance by a a 0 with a standard C T -based empirical estimation of a x =a for highly-loaded rotors, using B 2 ? + k B

| Method
Multiple alternatives may be considered for parameter identification, depending on the data that are available and the desired level of complexity and fidelity. Previous studies were based either on general rotor specifications only 11 (for axial components only), or on extensive curve fitting. 15 In both cases, models featured fewer parameters than suggested in Section 5, and only one rotor to validate against. The method suggested in this paper aims at maximizing applicability while preserving the physics-based approach. General conservation principles (momentum and energy) are combined with rotor design specification (optimal tip-speed ratio, power coefficient, and pitch angle) and performance curves (thrust and power coefficient and pitch angle). Depending on availability, alternative closure methods may use engineering model simulations, airfoil data, semi-empirical approximations, or manual tuning.
A procedure making use of information at two particular wind speeds is suggested: an arbitrary wind speed in the upper below-rated region where power is optimized and an arbitrary wind speed in the far above-rated region, close to cut-out. Quantities at these points will be denoted by the subscripts ? and +, respectively.
• In the optimized power production region, a ? ≈a r? ≈ 1 3 over the blade consistently with the Jukowsky representation of optimal rotors. 25 The optimal tip-speed ratio λ ? is typically given in the rotor specification. The pitch angle β ? is assumed 0 (any offset may be included in the twist).
Equations (10) to (23) may then be used to compute α ? , C L? , C D? , Γ ? , Δ ? and up to G ? . The optimal power coefficient C P? is also typically given (or refined from simulations).
• In the far above-rated region, a + may be neglected (a + U ∞ = 25m=s ð Þ ≈0:019 to 0.026 for the considered rotors) or computed from the thrust coefficient using (7). The latter, as well as the pitch angle β + , may be obtained from the rotor's performance curves or from simulations. The tip-speed ratio trivially reads λ + = ΩrtdR U∞ + with Ω rtd the rated rotor speed. The rated mechanical torque Q rtd is derived from rated power and efficiency.
The angle α ? is taken as a reference in the airfoil coefficient approximation procedure of Section 4.2: C D = C D? , α = α ? , c e , and α 0 are fitted over the range 0 α ? ½ .
Using engineering models, equilibrium values are found from closed-loop simulations with wind turbine controller active. Derivatives are computed about fixed equilibrium values of λ and β using the perturbation method. Variations in λ, σ, and ζ are modeled by modifying the wind input file. The current identification procedure does not make use of zero-drag simulations to isolate circulatory contribution, nor of frozen wake modeling to isolate load and inflow variations.

| Axial components
The corresponding parameters are k Γ , k λ , ℓ, k ΔQ , k ΔT , B ? , and k B , requiring a system of seven equations. This number may be raised to nine to refine a ? and a + , thus considered as parameters. A first identity is given by momentum balance in the optimized power region, reading involving k Γ , k λ , B ? , and a ? . A second identity is given by energy conservation, relating power losses to the optimized power coefficient and its loss-free limit (the famous Betz limit 4a ? 1 −a ? ð Þ 2 = 16 27 ): involving k ΔQ , B ? , and a ? . Information about the variation of the power coefficient with tip-speed ratio in the optimized power region is a good complement to (46) to match derivatives 3 . It reads and relates k ΔQ , k B and a ? . dCP dλ ? and dCT dλ ? may be calculated using engineering model simulations. Alternatively, assuming a ? = 1 3 and dCP dλ ? = 0 from optimal design removes this need by reducing (47) to Performance curves also readily give C T? , which gives B ? through (45) given a ? . Two more identities are furnished by 3 The derivative with respect to pitch angle dCP dβ ?
≈0 could in theory be used, but appears to be ill-posed and lead to the nonphysical result ℓ = 0 refining a + and linking the angular and linear drag parameters k ΔQ and k ΔT with the tip-loss parameters B ? and k B given C T + . To close the identification of the loss parameters k ΔQ , k ΔT ,B ? , k B À Á , k ΔQ may be calculated by direct integration: Regarding loss-free parameters, one more identity is provided by the torque coefficient near cut-out given the pitch angle: involving the parameters k Γ , k λ , ℓ, and k ΔQ . It may be appealing to use one (or more through curve fitting) extra point(s) in the power curves. However, it appears that an apparently good match of the C T , C P curves does not guarantee a satisfactory modeling of the variations with respect to β and λ. In fact, several sets of parameters are in appearance solutions. One piece of information of different nature is required to close the system of equations linking k Γ , k λ , and ℓ without redundancy. This information may be the variation of C T or C P with λ or β at one of the particular wind speeds (except the already used dCP dλ ? and discussed dCP dβ ? ), obtained from engineering model simulations. A good choice appears to be dCT dλ ? , yielding by use of (49) involving k Γ , k λ , B ? , and a ? . If simulations are not a possibility, an alternative may be derived from airfoil data. As radial inflow variations do not affect Γ 00s , and assuming the error induced by linearizing lift is reasonable (see Figure 3), direct integration leads to linking k Γ and ℓ.
If airfoil data are not available, some estimation of the parameters must be used to close the system. The least important parameter k λ may be set to 0. Alternatively, ℓ may be estimated by assuming that blade pitching takes action at the aerodynamic center of the blade. 11 If rotor performance curves are not available, further approximations on the tip-loss and linear dissipative parameters may be used. 11 Table 2 lists the approximations that may be made. Manual tuning may then be thought provided that some fitting criterion is available.

| Tangential components
Tangential components being secondary in rotor design, no relevant information is expected to be found in published documentation. Recalling the results of Figure 5 presented in Section 4.2, direct integration leads to good results. Variations of tangential load coefficients with tangential velocity components at each particular wind speed are used, giving each parameter an equation, namely,  2   6  6  6  6  6  4   3   7  7  7  7  7  5 , , where the circulatory partial derivatives are derived from the total derivatives by removal of dissipative contributions and modification by the inflow model of Section 3.2: : ð57Þ It is here assumed that if engineering model simulations can be used, airfoil data is available. If this is not the case, the tunable parameters may be set to γ ζ C i , γ σC i , k ζ S i , γ σS i , leaving g k ζ C i , g k σC i , k σS i as zero.

| COMPARATIVE ANALYSIS
This section applies the parameter identification methods presented in Section 6 to the 6-DOF simplified load model suggested in Section 5.7, named here analytical actuator disc theory (AADT), and compares its performance with BEMT. Method 6 is similar as Method 1 except that angular drag is estimated directly from airfoil drag coefficients. The identities of Section 6.2 used for each method are summed up in Table A1T3, and the corresponding parameter values are given in Table A2T4. U ∞+ has been set to 20m/s-for a better average fit in above rated wind speeds-for Methods 1, 5, and 6 and to 25m/s-where errors in the estimation of a + have less effect-for Methods 2, 3, and 4. General results are plotted in Figure 9. Figure 10 shows the effect of the identification method on relative errors with respect to BEMT for axial quantities. The errors related to quantities used as direct constraints in parameter identification such as C T? or C P+ are by definition zero for all methods and therefore not shown in the figure. The errors are split between the optimized power production (subscript ?) and above-rated wind speed ranges (rated wind speed not included, subscript +). Average values over the range are used 4 , that is, δ X = P are by definition close to 0, absolute instead of relative errors are used and should be related to Figure 9 to be meaningful. δdC T dλ + should also be taken cautiously as dCT dλ crosses zero in the above-rated range.

| Axial components
The following observations come up: • Method 1 matches mean quantities almost perfectly. Most derivative quantities show an excellent match too. This shows that AADT's level of accuracy is comparable to that of look-up tables (regarding quasi-static quantities).
• Power coefficient derivatives tend to be slightly less accurately predicted than their thrust coefficient counterparts. The discrepancies observed seem to be rotor-dependent and hence not systematic. Although wind speeds beyond the optimized power region were set out of focus in Section 4.1.1, the match is still very good for most quantities.
• Method 3, although performing less accurately than Method 1 in general (with remarkable exceptions, such as for derivative quantities of the NREL 5-MW rotor), still provides an excellent way of transforming C T ,C P curves into an analytical model and hence obtaining their derivatives for free.
• The only slight discrepancy observed between the results of Methods 3 and 4 is also promising. With only basic information about any rotor, most of its rotor-integrated loads may be reconstructed. Some errors are observed in the equilibrium pitch angle at high wind speeds. It should be however noted that this error is at maximum 2 degrees (for the DTU 10-MW, which in general shows the poorest match for the notsimulation-based Methods 2 to 4), which is still comparable to possible pitch actuator errors.
• Surprisingly, using airfoil data does not add accuracy in general, when comparing Method 2 with Method 3. To the contrary, it constrains parameters to somewhat erroneous values, instead of letting the identification procedure adapt them to available error-free data. Still, using (52) and (55) should be kept in the toolbox as it may increase accuracy in situations not presented here.
• Assuming a constant tip-loss factor does not affect the modeling of mean values, but does decrease fidelity for derivative quantities when comparing Methods 1 and 5. Although not shown here, the detrimental effect of forcing k B to 0 is greater for not-simulation-based methods and would significantly affect the three latter remarks above.
• Refining drag modeling using airfoil data (Method 6), though improving the modeling of the contribution of drag itself (see below), does not seem to significantly increase overall fidelity (it may even decrease it in some cases).
• Variations of parameter values between rotors are slight, as expected from design optimality. Including more rotors would be beneficial to identify whether changes are more related to the design method or to the power rating.
The modeling of the contribution of drag to mean thrust and power coefficients is looked more into detail in Figure 11. Looking first at the thrust, it is seen that the model suggested by (35) is accurate for high wind speeds if calibrated against the "real" dissipative contribution output by engineering model simulations (i.e., deduced from zero-drag simulations). Empirical modeling as given in 2 when matching dissipative thrust only also gives fair results. To the contrary, calibrating against BEMT simulations greatly overestimates dissipative thrust. This is to compensate for inaccuracies in the modeling of the circulatory contribution. Choice has to be made between modeling accurately dissipative only or total thrust; hence, the suggestion of two empirical values for k ΔT in Table 2. However, the empirical value suggested to match total thrust appears to significantly diverge from that identified by BEMT simulations. This does not have a large effect on general results.
Regarding torque and power, it is seen that the estimation of C Q Δ? using airfoil data through (52) (Method 6) is accurate and provides the best identification method for k ΔQ (except for the DTU 10-MW rotor where using Method 1 is-by chance-spot-on), although it does not improve the general results as outlined above. The remaining uncertainty of order of magnitude 10% to 30% in angular drag (unveiled by zero-drag simulations) comes from the parametric modeling structure itself.

| Tangential components
Five identification methods are considered: • Using only input data to compute γ σ,ζ S,C i coefficients from (39) without correction, k σ,ζ S,C i parameters being zero (Method 1, same as Direct Integration in Figure 5).
• Not using input data and identifying directly the γ σ,ζ S,C i coefficients (and k ζ S i ) as parameters using BEMT (Method 4) or GDW simulations (Method 5).
U + ∞ has been set to 20m/s. Results are shown in Figure 12 and Table A3T5. Methods 2 and 3 efficiently correct and adapt AADT to their respective inflow models, which considers the use of higher fidelity ones. Unlike for axial components, larger errors are observed at far belowrated wind speeds, as expected from Section 4.1.1. Elsewhere, AADT is accurate enough to be used in lieu of empirical look-up tables. The poor performance of Methods 4 and 5, although limited by the parameter identification constraints, shows inconsistent trends witnessing an inability of capturing the underlying phenomena. The effect of in-plane velocity on out-of-plane moment does not seem to involve these intricate phenomena, and all methods agree well. Looking at parameter values, similarities between rotors are again evident. Correlation with turbine rating is not obvious, although trends might be seen for some parameters.
F I G U R E 1 0 Error with respect to BEMT, relative in percent or absolute where specified. The method number is printed on the top of each bin F I G U R E 1 1 Relative error of drag modeling with respect to BEMT simulations. The different curves correspond to different identification methods. The DTU 10-MW "dashed" curve on the right-hand plot is actually a dotted curve coinciding with a dashed-dotted one

| CONCLUSION
A novel approach for actuator disc modeling of wind turbine aerodynamics has been suggested in this paper. It enables finding the 6-DOF rotorintegrated load components (including thrust, torque, in-plane forces, and out-of-plane moments) from a reduced 6-DOF hub-based velocity vector. In this manner, it resembles the standard, empirical modeling approach based on pre-computed load coefficients from engineering model simulations. A blunt comparison between the results of the two approaches for three different rotors of different ratings shows that • AADT is able to reconstruct almost perfectly the thrust and power coefficient tables and their derivatives with respect to tip-speed ratio and pitch angle across the entire operating wind speed range. To achieve this, AADT requires only the mean coefficient values at two particular wind speeds and the derivatives with respect to tip-speed ratio at one of these speeds.
• If no information about derivatives are available, and even if only general rotor specification is given, AADT is expected to be able to capture most of the axial behavior of modern wind turbines.
• Regarding tangential loads, AADT may be calibrated against simulations using either BEMT or GDW to accurately reconstruct tangential load coefficient derivatives. This requires the knowledge of airfoil properties (chord length, twist angle, and lift coefficient).
The derivation process also led to a better understanding of wind turbine aerodynamics: • The dependency of the various load components on the various velocity components takes a generic, interpretable matrix form. Circulatory (lift) and dissipative (drag) as well as global (averaged azimuthally) and cyclic (varying azimuthally) contributions may be clearly separated. IPC may be included.
• The skew symmetry of the circulatory part of this matrix form is linked to the lossless nature of lift through energy conservation. 11 • Skew symmetry for axial components highlights the dependency of thrust on the rotor speed and on torque on the axial flow velocity. 11 • The matrix form enables to clearly grasp the effect of blade pitching separately to that of velocity. 11 • Tip losses may be given a rotor-integrated version based on helicopter theory and Prandtl's tip-loss correction for BEMT.   Using empirical, engineering, frequency-domain-oriented or high-fidelity models does not lead as clearly to these conclusions. In addition, AADT features other advantages: • Its mathematical nature makes it designated for control and state-estimation applications. The skew-symmetric nature of the dominating circulatory part implies strong physical properties that may be used to prove the stability of control laws. 11 • Only quasi-static results have been considered in this paper. Moving to properly speaking aerodynamics, beside inherently coupling thrust and rotor dynamics, AADT may readily be combined with dynamic wake models (Pitt and Peters or other actuator disc models 11,26 ) which is not the case for empirical models (frozen wake modeling needed). The ability to calibrate AADT against higher-fidelity inflow theories should also be highlighted.
• As physics are better represented, linearizing AADT is expected to give more accurate results than when using empirical models. Although it does not reach the level of detail of state-of-the-art frequency-domain-oriented models, it may be used in complement when insight is needed and when rotor-integrated loads are sufficient.
This paper does not treat the dynamic performance of AADT. After having added a dynamic inflow and a turbulent wind model, comparison with engineering models should be done through wind-wave simulations on FWTs.
AADT is an approach rather than a final model. Its promising results motivate the development of extensions to skewed flow and bladesweeping (3p) frequency loads. It might then be coupled with a model for the first tower-bending eigenmode -which is also only dependent on rotor-integrated loads-to provide an alternative to empirical-based elastic models for control applications. 2 Adding skewed flow (using for instance the analytical model suggested in, 17 linking AADT to Pedersen's Dynamic Vortex Theory 11 ) would enable AADT to be used for control and state estimation of large wind parks using nacelle yaw control.
Beside validation of dynamic performance and development of extensions, further work includes the application of AADT to more rotors to strengthen this paper's conclusions and possibly draw trends relating parameter values to rotor properties. Improvements may be made through a better understanding of the turbulent wake state, but this is a thorny subject. A better modeling of drag may also be achieved. Parametric corrections on tangential components involved a significant part of trial-and-error and would also benefit from more understanding. Parametric corrections for IPC should also be suggested, based on results for tangential loads.