Dear Wesley,
I’m sorry, but I don’t understand your questions. Please clarify.
Best regards,
Sir,
I’m doing BEM using Excel for initial testing.
I notice for the blade element at the tip (i.e. r=R), the axial and tangential induction factor becomes 1 and -1 after the first run. This becomes a problem when the second run starts as the Relative Velocity at that element is zero. Any idea about this one sir?
Kind regards,
Wesley
Dear Wesley,
I’m assuming this Excel file you are referring to has your own implementation of BEM?
When looking at the limit of the BEM solution with the Prandtl tip-loss model as r → R at the tip, I would expect that a → 1 and a’ → 0, which is what we have implemented in AeroDyn v15. (AeroDyn v14 and earlier never allowed nodes to lie right at the tip.)
In our AeroDyn v15 implementation, we simply set a = 1 and a’ = 0 for any nodes lying right at the tip when Prandtl tip loss is enabled. You can do the same in your Excel implementation of BEM, or you can probably get around any problems right at the tip by ensuring that nodes are not placed within some tolerance of the tip location.
Best regards,
Sir, thanks for this!
Yes sir, the Excel file has my own implementation (i.e. following the usual procedure of BEM as explained in several texts).
The value of a, a’ you provided sir are irrespective of the empirical correction used concerning turbulent wake state (e.g. Glauert, Buhl)?
As for the matlab implementation, the path to follow sir is using newton raphson?
Kind regards,
Wesley
Dear Wesley,
Regarding the limiting values at the tip, “yes,” that is correct.
There are different methods to solve the BEM equations. Newton iterations are one method. In AeroDyn v15, we use Brent’s method.
Best regards,
I see,
thanks sir!
Kind regards,
Wesley
Sir,
Good day!
Sir, I’m using trapezoidal rule in the computation of Tangential Force/Torque/Power.
The formula requires indicating segments (i.e. distance between two radial positions, r(n+1)-r(n)).
My computation started with r=1.257m (since this is the section where s809 is utilized).
Sir, if I will get the Tangential force for the first blade element (n=1), any idea what would the segment be?
Will the segment of the first blade element covers be r=1.257 and r=1.343?
Kind regards,
Wesley
Dear Wesley,
I’m sorry, but I don’t understand your question. Are you using AeroDyn v15, an older version of AeroDyn, or something else? What does your nodal distribution look like?
Best regards,
Dear Jason,
I have implemented axial induction factor calculation in MATLAB as given in the paper “Ning, Andrew, et al. “Development and validation of a new blade element momentum skewed-wake model within AeroDyn.” 33rd Wind Energy Symposium. 2015”.
I get good match for AxInd and TanInd between MATLAB code and FAST at the middle section of the blade, but at the tip , FAST reports, AxInd =1 and TanInd =-1; (I chose Wakemod=1, included tip, hub loss, drag is taken in the induction factor calculation)
But my matlab implementation gives values different from 1 and -1. From reading this post, I realized that AxInd and TanInd are hardcoded to be 1 and -1 respectively near tip. Is it incorrect to proceed with the values given by the algorithm given in the paper? I understand the problem with the implementation in paper right at the tip. Could you please elaborate on the implications of this on the further calculations.
Also, I don’t see hard coded values for AxInd and TanInd in the source code (given below). Where is this implemented in the code?
[code] ! compute axial induction factor
if (phi > 0.0_ReKi) then ! momentum/empirical
! update axial induction factor
if (k <= 2.0_ReKi/3.0_ReKi) then ! momentum state
if ( EqualRealNos( k, -1.0_ReKi) ) then
k = k - 0.1_ReKi ! Need to bump k to avoid singularities
end if
a = k/(1.0_ReKi+k)
else ! Glauert(Buhl) correction
g1 = 2.0_ReKi*F*k - (10.0_ReKi/9-F)
g2 = 2.0_ReKi*F*k - (4.0_ReKi/3-F)*F
g3 = 2.0_ReKi*F*k - (25.0_ReKi/9-2*F)
if (abs(g3) < 1e-6_ReKi) then ! avoid singularity
a = 1.0_ReKi - 1.0_ReKi/2.0/sqrt(g2)
else
a = (g1 - sqrt(g2)) / g3
end if
end if
else ! propeller brake region (a and ap not directly used but update anyway) !bjj: huh? when k is slightly larger than 1, a is definitely getting used (and causing issues)…
if (k > 1.0_ReKi .and. .not. EqualRealNos(k, 1.0_ReKi) ) then
!if (sigma_pcn > Fsphi) then
a = k/(k-1.0_ReKi) !sigma_pcn / (sigma_pcn - Fsphi ) !
! axial induction is blowing up, so I'm putting a band-aid here. BJJ 25-Feb-2016
a = min(a, 10.0_ReKi )
else
a = 0.0_ReKi ! dummy value
end if
end if
! compute tangential induction factor
if ( cphi==0.0_ReKi ) then ! We don’t want NaN here
kp = HUGE(kp)
else
kp = sigma_p*ct/4.0_ReKi/F/sphi/cphi
end if
! Per conversation with Rick, we should only trigger this if phi = 0 , so we will return predefined values as if phi=0.0
if (EqualRealNos(kp, 1.0_ReKi)) then
fzero = 0.0_ReKi
a = 0.0_ReKi
ap = 0.0_ReKi
return
end if
ap = kp/(1.0_ReKi-kp)
! tangential induction is blowing up, so we’re putting a band-aid here. GJH, JMJ, BJJ 1-Sep-2015
if ( abs(ap) > 10.0_ReKi ) then
ap = sign( 10.0_ReKi, ap )
end if
!bjj: 3-jun-2015: TODO: was able to trigger divide-by-zero here using ccBlade_UAE.dvr without tiploss or hubloss
if (.not. wakerotation) then
ap = 0.0_ReKi
kp = 0.0_ReKi
end if[/code]
Regards,
Kumara
Dear Kumara,
The BEM solution of AeroDyn generally follows Ning et al’s paper from AIAA SciTech 2015, but some of the details have been changed to ensure AeroDyn is robust for all cases, e.g., aero-elastic interaction under turbulent wind inflow.
Which version of FAST / AeroDyn are you using? In the current version of OpenFAST / AeroDyn v15, I would expect the axial induction to be 1.0 and the tangential induction to be 0.0 at the root and tip, when hub and tip losses are enabled. This is calculated in routine BEMTUncoupled/BEMTU_InductionWithResidual(). I would recommend upgrading if you are not at least using OpenFAST 0.1 or newer because AeroDyn v15 was upgraded and made more robust between FAST v8.16 and OpenFAST v0.1.
We need to develop documentation with the most up-to-date implementation of BEM, based on updates to Ning et al’s paper, but we do not currently have the resources to do so. Until then, the source code is probably your best resource as to the details of the implementation.
Best regards,
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
Dear Kumara,
The phi used at the beginning of the update-states 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 aero-elastic deflection), and induction (a,a’).
Best regards,
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
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,
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
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.
[code]else ! propeller brake
if ( EqualRealNos(k,1.0_ReKi) ) then
IsValidSolution = .false.
a = InductionLimit
else
a = k/(k-1.0_ReKi)
end if
[/code]
if (.not. IsValidSolution) then
a = 0.0_ReKi
ap = 0.0_ReKi
end if
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,
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
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 = rOmega_r;
Regards,
Kumara
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: 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,