Help about BEM and NREL's Phase VI
Moderator: Bonnie.Jonkman

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Dear Jason,
I am currently using FASTV8.6, AeroDynV15.03. Now, I have downloaded OpenFASTV2.1.0 and going through the source code of AeroDyn in it. I see some differences compared to AeroDynV15.03. I could follow most of the details. If I understand it right, the starting point for the algorithm to predict induction factor requires a value of phi, Vx, Vy and many other variable, based on which a, a' are evaluated. Again based on the test region, sub_brent is called to find the phi_new. Again, the process is repeated until the residual is with in tolerance are met. My question is at the very beginning of the process, what value of "phi" is used. Is it atan2(Vx, VY); Where Vx = free strem wind velocity, Vy = rad_chord *Omega. are velocities without induction factors? Thanks.
Regards,
Kumara
I am currently using FASTV8.6, AeroDynV15.03. Now, I have downloaded OpenFASTV2.1.0 and going through the source code of AeroDyn in it. I see some differences compared to AeroDynV15.03. I could follow most of the details. If I understand it right, the starting point for the algorithm to predict induction factor requires a value of phi, Vx, Vy and many other variable, based on which a, a' are evaluated. Again based on the test region, sub_brent is called to find the phi_new. Again, the process is repeated until the residual is with in tolerance are met. My question is at the very beginning of the process, what value of "phi" is used. Is it atan2(Vx, VY); Where Vx = free strem wind velocity, Vy = rad_chord *Omega. are velocities without induction factors? Thanks.
Regards,
Kumara

 Posts: 5738
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Help about BEM and NREL's Phase VI
Dear Kumara,
The phi used at the beginning of the updatestates algorithm is the phi solved at the previous step. This phi is used as the starting guess for calculating phi at the next time step. Phi is calculated using free stream wind, rotor rotation, structural velocities (from aeroelastic deflection), and induction (a,a').
Best regards,
The phi used at the beginning of the updatestates algorithm is the phi solved at the previous step. This phi is used as the starting guess for calculating phi at the next time step. Phi is calculated using free stream wind, rotor rotation, structural velocities (from aeroelastic deflection), and induction (a,a').
Best regards,
Jason Jonkman, Ph.D.
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Dear Jason,
Sorry, it is still not clear to me. Assume that blades and tower are rigid (no structural velocities), pitch angle=0, wind velocity =6 m/s (uniform), initial angular velocity of rotor = 2 rpm. Now, while calculating induction factors at time zero, the Phi_guess (for the very first iteration) cannot contain induction factors, right? In my understanding it should be calculated as atan2(Vx, Vy). Right? (or some default starting value. )?
Regards,
Kumara
Sorry, it is still not clear to me. Assume that blades and tower are rigid (no structural velocities), pitch angle=0, wind velocity =6 m/s (uniform), initial angular velocity of rotor = 2 rpm. Now, while calculating induction factors at time zero, the Phi_guess (for the very first iteration) cannot contain induction factors, right? In my understanding it should be calculated as atan2(Vx, Vy). Right? (or some default starting value. )?
Regards,
Kumara

 Posts: 5738
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Help about BEM and NREL's Phase VI
Dear Kumara,
I'm not sure I understand what you want to know. Are you interested in the value of phi used/set in AD_Init, AD_UpdateStates, or AD_CalcOutput?
Best regards,
I'm not sure I understand what you want to know. Are you interested in the value of phi used/set in AD_Init, AD_UpdateStates, or AD_CalcOutput?
Best regards,
Jason Jonkman, Ph.D.
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Dear Jason.
My apologies for not making the question clear enough. I am looking for the initialization value for the phi. I think I found the answer in the source code (in the colored text below). Thank you!
!...............................................................................................................................
! if we haven't initialized z%phi, we want to get a better guess as to what the actual values of phi at t are:
!...............................................................................................................................
if (.not. OtherState%nodesInitialized) then
if (p%useInduction) then
do j = 1,p%numBlades ! Loop through all blades
do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements
NodeTxt = '(node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')'
call BEMT_UnCoupledSolve(z%phi(i,j), p%numBlades, p%kinVisc, AFInfo(p%AFIndx(i,j)), u1%rlocal(i,j), p%chord(i,j), u1%theta(i,j), &
u1%Vx(i,j), u1%Vy(i,j), p%useTanInd, p%useAIDrag, p%useTIDrag, p%useHubLoss, p%useTipLoss, p%hubLossConst(i,j), p%tipLossConst(i,j), &
p%maxIndIterations, p%aTol, OtherState%ValidPhi(i,j), m%FirstWarn_Phi, errStat2, errMsg2)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName//trim(NodeTxt))
if (errStat >= AbortErrLev) return
end do
end do
else
! We'll simply compute a geometrical phi based on both induction factors being 0.0
do j = 1,p%numBlades ! Loop through all blades
do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements
z%phi(i,j) = ComputePhiWithInduction(u1%Vx(i,j), u1%Vy(i,j), 0.0_ReKi, 0.0_ReKi)
end do
end do
end if
OtherState%nodesInitialized = .true. ! state at t+1
end if
My apologies for not making the question clear enough. I am looking for the initialization value for the phi. I think I found the answer in the source code (in the colored text below). Thank you!
!...............................................................................................................................
! if we haven't initialized z%phi, we want to get a better guess as to what the actual values of phi at t are:
!...............................................................................................................................
if (.not. OtherState%nodesInitialized) then
if (p%useInduction) then
do j = 1,p%numBlades ! Loop through all blades
do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements
NodeTxt = '(node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')'
call BEMT_UnCoupledSolve(z%phi(i,j), p%numBlades, p%kinVisc, AFInfo(p%AFIndx(i,j)), u1%rlocal(i,j), p%chord(i,j), u1%theta(i,j), &
u1%Vx(i,j), u1%Vy(i,j), p%useTanInd, p%useAIDrag, p%useTIDrag, p%useHubLoss, p%useTipLoss, p%hubLossConst(i,j), p%tipLossConst(i,j), &
p%maxIndIterations, p%aTol, OtherState%ValidPhi(i,j), m%FirstWarn_Phi, errStat2, errMsg2)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName//trim(NodeTxt))
if (errStat >= AbortErrLev) return
end do
end do
else
! We'll simply compute a geometrical phi based on both induction factors being 0.0
do j = 1,p%numBlades ! Loop through all blades
do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements
z%phi(i,j) = ComputePhiWithInduction(u1%Vx(i,j), u1%Vy(i,j), 0.0_ReKi, 0.0_ReKi)
end do
end do
end if
OtherState%nodesInitialized = .true. ! state at t+1
end if

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Dear Jason,
I have one more question. In the source code corresponding to the propeller brake region, (Copied below) after setting IsValidSolution =.false, why to assign value to a? Because at another portion in the source copied below (with blue color text) a and a' are assigned as zeros. ? Please let me know incase the question is not clear.
I have one more question. In the source code corresponding to the propeller brake region, (Copied below) after setting IsValidSolution =.false, why to assign value to a? Because at another portion in the source copied below (with blue color text) a and a' are assigned as zeros. ? Please let me know incase the question is not clear.
Code: Select all
else ! propeller brake
if ( EqualRealNos(k,1.0_ReKi) ) then
IsValidSolution = .false.
a = InductionLimit
else
a = k/(k1.0_ReKi)
end if
Code: Select all
[color=#0000BF] if (.not. IsValidSolution) then
a = 0.0_ReKi
ap = 0.0_ReKi
end if[/color]

 Posts: 5738
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Help about BEM and NREL's Phase VI
Dear Kumara,
I didn't write this routine, so, I'm not sure. Regardless, it looks like a and ap (a') are always set to zero if the solution is not valid.
Best regards,
I didn't write this routine, so, I'm not sure. Regardless, it looks like a and ap (a') are always set to zero if the solution is not valid.
Best regards,
Jason Jonkman, Ph.D.
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Dear Jason,
I went through the source code (inductionFactors subroutine, BEMTUncoupled module) again. Now I understand the reason for defining 'a' twice (before and after the evaluation of residual). It is because Residual calculations (which are performed in later part of the code) require the values of 'a' and 'kp' . Thus, 'a' value is assigned a large number (InductionLimit) even for the case of 'nosolution'. After the residual calculations, a and a' for the 'no solution' region are assigned 'zero' values for the next iterations. I don't know the logic for choosing these particular values though. Thanks.
Regards,
Kumara
I went through the source code (inductionFactors subroutine, BEMTUncoupled module) again. Now I understand the reason for defining 'a' twice (before and after the evaluation of residual). It is because Residual calculations (which are performed in later part of the code) require the values of 'a' and 'kp' . Thus, 'a' value is assigned a large number (InductionLimit) even for the case of 'nosolution'. After the residual calculations, a and a' for the 'no solution' region are assigned 'zero' values for the next iterations. I don't know the logic for choosing these particular values though. Thanks.
Regards,
Kumara

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Dear all,
Should the induction factors (a, a') and inflow angle (phi) obtained from FAST, satisfy the below equation always? I am using OpenFastV2.2.0.
[sin(phi) / (1  a) ]  [ (V_wind/V_r)*cos(x) / (1 + a') ] = 0
where V_r = r*Omega_r;
Regards,
Kumara
Should the induction factors (a, a') and inflow angle (phi) obtained from FAST, satisfy the below equation always? I am using OpenFastV2.2.0.
[sin(phi) / (1  a) ]  [ (V_wind/V_r)*cos(x) / (1 + a') ] = 0
where V_r = r*Omega_r;
Regards,
Kumara

 Posts: 5738
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Help about BEM and NREL's Phase VI
Dear Kumara,
I assume you meant to type "cos(phi)" instead of "cos(x)", in which case you seem to be restating equation (22) from our AIAA SciTech 2015 paper on the theory and validation of the BEM implementation within AeroDyn v15: https://www.nrel.gov/docs/fy15osti/63217.pdf. Because this equation is undefined at a=1 and/or a'=1, the current version of AeroDyn v15 uses a different form of this equation, but the equation is generally valid otherwise. Please note that instead of V_wind and V_r, our AIAA SciTech 2015 paper and the the AeroDyn v15 implementation use V_x and V_y, respectively, which can include blade motions other than r*omega.
I hope that helps.
Best regards,
I assume you meant to type "cos(phi)" instead of "cos(x)", in which case you seem to be restating equation (22) from our AIAA SciTech 2015 paper on the theory and validation of the BEM implementation within AeroDyn v15: https://www.nrel.gov/docs/fy15osti/63217.pdf. Because this equation is undefined at a=1 and/or a'=1, the current version of AeroDyn v15 uses a different form of this equation, but the equation is generally valid otherwise. Please note that instead of V_wind and V_r, our AIAA SciTech 2015 paper and the the AeroDyn v15 implementation use V_x and V_y, respectively, which can include blade motions other than r*omega.
I hope that helps.
Best regards,
Jason Jonkman, Ph.D.
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

 Posts: 51
 Joined: Mon Sep 12, 2016 8:34 pm
 Organization: IIT Bombay
 Location: Mumbai
Re: Help about BEM and NREL's Phase VI
Yes Jason, i meant "phi" not 'x'.
My main concern is how to deal with the case of a=1, and a' = 1; I understand your comments on Vx, Vy. Since I'm looking at the simple case of rigid blades, I'm using V_wind and r*Omega instead of Vx, Vy.
Could you please tell me how are these exceptional cases resolved in AeroDyn or FAST?
My understanding is that for a=1, the residual equation used is:  cos(phi) * (Vx/Vy) * (1  k_prime) not equation22.
But I can't understand how is the a' = 1 dealt. Thanks.
Regards,
Kumara
My main concern is how to deal with the case of a=1, and a' = 1; I understand your comments on Vx, Vy. Since I'm looking at the simple case of rigid blades, I'm using V_wind and r*Omega instead of Vx, Vy.
Could you please tell me how are these exceptional cases resolved in AeroDyn or FAST?
My understanding is that for a=1, the residual equation used is:  cos(phi) * (Vx/Vy) * (1  k_prime) not equation22.
But I can't understand how is the a' = 1 dealt. Thanks.
Regards,
Kumara

 Posts: 5738
 Joined: Thu Nov 03, 2005 4:38 pm
 Location: Boulder, CO
 Contact:
Re: Help about BEM and NREL's Phase VI
Dear Kumara,
AeroDyn uses 1/(1+a') = (1k'), which is equivalent based on Equation (21) in our AIAA SciTech 2015 paper linked above. Likewise, AeroDyn uses 1/(1a) = (1k), which is equivalent based on Equation (20).
Please note that if a'=1, then phi = +/90 degrees, based on the sign of Vx (see Figure 2). And if a=1, then phi = 0 or 180 based on the sign of Vy. If both a=1 and a'=1, then phi is undefined, but irrelevant because the dynamic pressure and aerodynamic loads are zero anyway (Equations (23)(27)).
Best regards,
AeroDyn uses 1/(1+a') = (1k'), which is equivalent based on Equation (21) in our AIAA SciTech 2015 paper linked above. Likewise, AeroDyn uses 1/(1a) = (1k), which is equivalent based on Equation (20).
Please note that if a'=1, then phi = +/90 degrees, based on the sign of Vx (see Figure 2). And if a=1, then phi = 0 or 180 based on the sign of Vy. If both a=1 and a'=1, then phi is undefined, but irrelevant because the dynamic pressure and aerodynamic loads are zero anyway (Equations (23)(27)).
Best regards,
Jason Jonkman, Ph.D.
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov
Senior Engineer  National Wind Technology Center (NWTC)
National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway  Golden, CO 80401
+1 (303) 384 – 7026  Fax: +1 (303) 384 – 6901
nwtc.nrel.gov
Return to “Rotor Aerodynamics”
Who is online
Users browsing this forum: No registered users and 1 guest