Wake-Skewed Angle Model

Albert.Meseguer · August 23, 2016, 10:02am

Hi everyone,

I am running some simulations where I compare the effect of the two available skewed-wake correction model in Aerodyn v15.03 (Uncoupled and Pitt and Peters model) . My case presents a wind turbine with no tilt neither blade pre-cone. I use an uniform wind speed without shear effect.

Equation 30 in the report " Development and Validation of a New Blade Element Momentum Skewed-Wake Model within Aerodyn" ( [url]http://www.nrel.gov/docs/fy15osti/63217.pdf[/url] ) , models the Pitt and Peters correction as function of the induction factor for each section:

Being X[size=50]0[/size] the angle between the vector normal to the rotor plane, in my case the yaw angle. I expect to obtain same results when using Uncoupled and Pitt and Peters models for yaw angle 0 ( the tangent term becomes 0). However, the models differ and I really do not understand why. I.e. I plot thrust of the turbine for both models:

Any idea why? Am I missing something?

Thank you,

Jason.Jonkman · August 23, 2016, 12:45pm

Dear Albert,

I have not tried to reproduce your results, but I quickly skimmed the AeroDyn BEMT algorithm and it is not obvious to me how changing input SkewMod between uncoupled (SkewMod=1) and Pitt and Peters (SkewMod=2) would have any influence on the results if the inflow-skew angle were in fact zero. Have you output the inflow-skew angle (AeroDyn output RtSkew) to verify that it is in fact zero? Have you also noticed the inflow angle (output BαNβPhi), angle of attack (output BαNβAlpha), or the axial and/or tangential induction factors (outputs BαNβAxInd and/or BαNβTnInd) changing for certain blade nodes between these cases?

Best regards,

Albert.Meseguer · August 24, 2016, 8:40am

Dear Jason,

Thank you for your fast reply.

In both simulations the inflow-skew angle is 0 in all the simulation (as expected) . I took a look into the aerodynamic parameters you mentioned in your post , for nodes close to the root both induction factors (B1N1AxInd , B1N1TnInd) are the same. Therefore , the inflow and attack angle (B1N1Alpha , B1N1Phi) present a really good agreement.

However, If I output other blade nodes , i.e. node 9 , the tangential induction factor value (B1N9TnInd) for the uncoupled option remains -1 all the simulation. On the other hand, the Pitt and Peters model is 0 the whole simulation. The axial induction factor (B1N9AxInd) has also discrepancies, the Uncoupled model predicts a value of +1 the complete simulation while the Pitt and Peters is 0 . This has an impact in inflow and attack angle since I’m considering the tangential induction factor in the calculation (TanInd = TRUE). This results in this discrepancies over time :

I do not understand why the induction factors are different if the module is only supposed to correct the yaw misalignment. I’m missing something probably.

Any idea?

P.S. : The first node “B1N1…” corresponds to the first node of the blade and “B1N9…” is the last node of the blade, defined in BlOutNd as 1,3,5,9,11,13,15,17,19.

Jason.Jonkman · August 24, 2016, 6:45pm

Dear Albert,

OK, thanks for looking into this more.

I looked in the source code (specifically SUBROUTINE ApplySkewedWakeCorrection() in source file BEMTUncoupled.f90) and noticed that if the axial and tangential induction factors (a and ap) are set to 1 or -1, respectively, both are reset to zero. This was added to avoid numerical problems in the calculation of the inflow-angle (phi) later in the code. However, earlier in the code, a=1 and ap=-1 is set at the root and tip when the root and tip loss models are enabled to ensure that no loads are computed at this nodes. So, in essence, the skewed-wake correction is over-writing the tip and hub loss effects at the root and tip in your case. Obviously, there is a problem in the logic here that we must fix. I’ve noted this as an issue that will be fixed in a future version.

Thanks for your help in debugging.

Best regards,

Saptarshi.Sarkar · June 20, 2018, 2:25pm

Dear Dr. Jonkman,

For the Pitt and Peter’s skewed wake correction is the angle between the vector normal to the rotor plane and the wind vector calculated at every blade section? What I mean is, under turbulent wind every blade section will have a different wind vector while the normal vector to the rotor plane remains same. In this situation should I calculate the angle for every section separately and use that angle for the skewed wake correction?

I would really appreciate if you could help clear this confusion.

Thanks a million,
Saptarshi

Jason.Jonkman · June 20, 2018, 7:53pm

Dear Saptarshi,

The skew angle in the Pitt and Peters skewed-wake correction is defined as the angle between the vector normal to the rotor plane and the vector of the wind averaged across the rotor. This angle is the same for all blades.

Best regards,

Saptarshi.Sarkar · June 21, 2018, 9:58am

Dear Dr. Jonkman,

Thanks a lot for clarifying this.

Saptarshi

Lorenzo.Cottura · September 3, 2020, 9:46am

Dear Jason,
I introduce myself to this topic because my question is about skew correction.
In particular, I was wondering the formulation of Xo (angle between the vector normal to the rotor plane) in case of non-constant wind and all 6 active platform DOFs.
Could you bring me the equation or a document in which it is reported?
Thanks for your help,
Lorenzo.

Jason.Jonkman · September 3, 2020, 2:54pm

Dear Lorenzo,

The equation for the skew angle is calculated as described, i.e. the angle between the vector normal to the rotor plane and the vector of the wind averaged across the rotor. In the source code, it is calculated as:

skew = ACOS( V_DiskAvg dot xhat_Disk / 2NORM ( V_DiskAvg ) )

where V_DiskAvg is a the relative velocity (wind minus structure) averaged over the rotor disk and xhat_Disk is a unit vector normal to the disk.

Best regards,

Lorenzo.Cottura · November 9, 2020, 6:02pm

Dear Jason,
I’m not sure what you mean by “wind averaged across the rotor”. Tell me if I understood correctly:
By defining V_DiskAvg = [V0_x, V0_y, V0_z] the vector of the wind speeds relative to the structure (wind minus structure) in the node under analysis, to be clear the same used in the BEM theory to obtain Vx and Vy according to the formula reported in “Development and Validation of a New Blade Element Momentum Skewed-Wake Model within AeroDyn” article where Vinf=V0_x.

The formula becomes:

skew = acos (V0_x / (sqrt (V0_x ^ 2 + V0_y ^ 2 + V0_z ^ 2)))

Confirm?
Thanks for the help,
Lorenzo.

Jason.Jonkman · November 9, 2020, 6:31pm

Dear Lorenzo,

Equation (28) from Ning et al’s paper is not implemented in FAST / AeroDyn, as this formula is only valid for a rigid wind turbine and uniform flow.

In AeroDyn, the V_DiskAvg is computed by spatially averaging the relative wind velocity at each aerodynamic analysis node across the rotor, i.e.:

V_DiskAvg = SUM( SUM( V_wind(j,k) - V_structural(j,k), j = 1, 2, …, NumBlNds ), k = 1, 2, …, NumBlades )/( NumBlades*NumBlNds )

I hope that helps.

Best regards,

Lorenzo.Cottura · November 9, 2020, 7:03pm

Dear Jason,
So V_DiskAvg becomes a constant vector for each time step like: V_DiskAvg = [V_DiskAvg_x, V_DiskAvg_y, V_DiskAvg_z]
and it will be
skew = acos (V_DiskAvg / (sqrt (V_DiskAvg_x^2 + V_DiskAvg_y^2 + V_DiskAvg_z^2)))
It is right?
Best regards,
Lorenzo

Jason.Jonkman · November 9, 2020, 7:22pm

Dear Lorenzo,

If xhat_Disk is directed along the x axis, then:

skew = acos (V_DiskAvg_x / (sqrt (V_DiskAvg_x^2 + V_DiskAvg_y^2 + V_DiskAvg_z^2)))

Best regards,

Lorenzo.Cottura · November 9, 2020, 7:50pm

Dear Jason,
Yes, I meant that, sorry for the typo.
Best regards,
Lorenzo.

Keerthana.Mohan · August 5, 2022, 2:59pm

Dear Jason,

This is follow-up to the questions on skewed wake model of OpenFAST. The constant used in the source code AeroDyn_IO.f90 of OpenFAST v3.1 for Pitt-Peters model is 15pi/32 (as also mentioned in the AeroDyn.dat file).

However, Ning’s paper mentions that a value of 15pi/64 in Eq. (30). May I know why this value is not implemented since it has been reported that the value of 15pi/64 gives a better comparison with literature?
When I changed the value to 15pi/64, I got abnormally high value of axial induction factor at skew angle of 30deg for 10 MW RWT in comparison with similar simulation with constant value of 15pi/32. The inflow conditions are steady wind with speed of 11 m/s, with all structural DOFs turned off. I am attaching a plot of the comparison of azimuthal variation of axial induction factors at a particular radial station. What could probably be the reason for this difference?

Jason.Jonkman · August 5, 2022, 10:15pm

Dear @Keerthana.Mohan,

We had originally implement Ning’s recommended value for SkewModFactor, but we found that the value of 15pi/32 gave results that better matched older versions of AeroDyn. In the end, we made SkewModFactor an input that that user can adjust as needed to match known solutions.

Different values of this parameter are reported in the literature. I suspect different values are needed based on how the BEM method is implemented. In AeroDyn, local values of induction at each blade node are always computed, resulting in sinusoidal variation of induction in skewed flow even with the Pitt and Peters model disabled, so, the Pitt and Peters model is really adding an additional correction.

Best regards,

Keerthana.Mohan · August 8, 2022, 8:55am

Dear Jason,

It makes more sense. The problem however seems to be getting the reference data for calibration of the coefficient.
Thanks a lot for taking the time to answer!

Miguel.Pereira · July 10, 2023, 11:09am

Dear Jason,

I am using a BEM code (CCBlade) with the Pitt/Peters correction, with a constant of 15pi/64. I am using the same framework where the equations are parameterized for phi, and solving the BEM consists of finding the roots of the function f:

However, for certain conditions (downwind positions, higher yaw angles and high local tip speed ratios) there are no solutions to the equations (with phi between pi/2 and 0).

Considering the NREL turbine, the following figure shows f0 as a function of phi (the solution is found where the function crosses the x axis) for an azimuth of 90 degrees (0º = 12 o’clock) and a radial position of r = 58.9m (NACA6_A17).

Similar conclusions can be drawned from increasing the local tip speed ratio, or moving from upwind positions to downwind positions.

How did you overcome this problem in FAST?

Best regards,
Miguel Pereira

Emmanuel.Branlard · July 18, 2023, 2:21am

Dear Miguel,

I would like to try to come up with a full answer for you, but right now I haven’t had the time to dive into the problem.
It is possible that the answer to your problem is by looking for solutions in different test regions. OpenFAST solves for phi in the following routine:

github.com

OpenFAST/openfast/blob/main/modules/aerodyn/src/BEMT.f90#L2123


      
                   end if
                end if
                
                
             else
                TestRegion = 2  ! No zero
             end if
             
          end function FindTestRegion
                
          subroutine BEMT_UnCoupledSolve(p, u, iBladeNode, jBlade, phi, AFInfo, ValidPhi, FirstWarn, ErrStat, ErrMsg)
          
             use mod_root1dim
             !use fminMod
             type(BEMT_ParameterType),intent(in  ) :: p
             type(BEMT_InputType),    intent(in  ) :: u
             integer(IntKi),          intent(in  ) :: iBladeNode         !< index for blade node
             integer(IntKi),          intent(in  ) :: jBlade             !< index for blade
             real(ReKi),             intent(inout) :: phi
             TYPE(AFI_ParameterType),INTENT(IN   ) :: AFInfo
             logical,                intent(inout) :: ValidPhi

And finds the test region using this function:

github.com

OpenFAST/openfast/blob/main/modules/aerodyn/src/BEMT.f90#L2038


      
          
                   test_lower(2) = -PiBy2 - BEMT_epsilon2
                   test_upper(2) = -pi + BEMT_epsilon2
                end if
          
             end if
             
             
          end subroutine GetSolveRegionOrdering
             
          integer function FindTestRegion(p, u, iBladeNode, jBlade, phiLower, phiUpper, AFInfo, &
                                  phiIn_IsValidSolution, phiIn, f_phiIn, f1, f2, errStat, errMsg) result(TestRegion)
          
             type(BEMT_ParameterType),intent(in  ) :: p
             type(BEMT_InputType),   intent(in   ) :: u
             integer(IntKi),         intent(in   ) :: iBladeNode         !< index for blade node
             integer(IntKi),         intent(in   ) :: jBlade             !< index for blade
             real(ReKi),             intent(inout) :: phiLower !intent "out" in case the previous solution can alter the test region bounds
             real(ReKi),             intent(inout) :: phiUpper !intent "out" in case the previous solution can alter the test region bounds
             !integer,                intent(in   ) :: numBladeNodes
             type(AFI_ParameterType),intent(in   ) :: AFInfo

Let me know what you find and if that’s indeed pointing you in the right direction.

Thanks

Emmanuel

Yangyang.Li · July 27, 2024, 12:08pm

Dear Jason,
I reviewed the formula for calculating V_DiskAvg provided above by you, and tested it with my case. However, I noticed a discrepancy in the Z direction’s RtVAvgzh. My calculated value is -0.52, while the AeroDyn 15 output is -1.05, resulting in a difference of 0.5. The detailed calculation results are as follows:

I attempted to perform a coordinate system transformation, converting the local blade coordinate system to the hub coordinate system, but the results were still incorrect. Additionally, I searched relevant information in this forum, but couldn’t find anything helpful.

My case is as follows:
a steady wind of 12 m/s without shear effect, NREL 5MW with zero yaw, PreCone(1-3)=-2.5, ShftTilt=5, RotSpeed=0,TMax=10, The above results are at t=0.

I hope you can provide the specific calculation formula for RtVAvgzh in the Z direction. If possible, it would be great to also include the formulas for BαNβUndy and BαNβUndz.

Any help would be greatly appreciated, Best Regards,
Yangyang. Li

Topic		Replies	Views
Yaw loads at skewed inflow Rotor Aerodynamics	1	6177	June 1, 2011
axial induction factor from AeroDyn13 Rotor Aerodynamics	8	6988	January 17, 2019
BEMT method in AeroDyn module Computer-Aided Engineering Software Tools	1	189	April 17, 2023
Generalized Dynamic Wake Theory and Implementation AeroDyn Development	13	97900	February 5, 2015
Direction of the incoming wind Rotor Aerodynamics	3	8655	March 26, 2012

Wake-Skewed Angle Model

Related topics