Implementation of Buhl's Correction

Is this correct way of implementing Buhl’s correction?

F1=(2./pi).(acos(exp(-(((n./2)(R-r(k)))./(r(k).sind(phae))))))
F2= ((2./pi).
(acos(exp(-((n./2).*(r(k)-Rb))/((Rb).sind(phae))))))
F=((F1).
(F2));

CT=((Cx*((1-a)^2)/(sind(phae))^2))((nc(k))/(2pir(k))) eq (1)

if CT<=(.96*F)    
                 
    a=1./((((4.*F.*(sind(phae)).^2))./(sigma.*Cl.*cosd(phae)))+1)          eq(2)    (From BEMT)

else

    a=((18.*F)-20-3.*((CT.*(50-36.*F)+(12.*F.*(3.*F-4))).^0.5))./(36.*F-(50))     eq(3)

I am confused if i should use the value of CT from eq(1) for substituting in eq (3) or should i use following equation.

CT= (8/9)+ (((4.*F)-(40/9)).*a)+(((50/9)-(4.F)).(a.^2)) eq(4) (Buhl’s equation).

If if use eq(1) for calculating CT my solution converges for very small range of tip speed ratio.

Dear Zeeshan,

You should use Eq. (1) to calculate CT for use in Eq. (3). Eq. (1) is the thrust coefficient according to blade-element theory, with no assumption about momentum theory. Eq. (4) is equivalent to Eq. (3), but solved for CT instead of “a”.

Of course, there are other BEM solution procedures that include the Glauert correct. The approach you’ve described is equivalent to the one described in the 2005 AeroDyn theory manual: nrel.gov/docs/fy15osti/63217.pdf.

Best regards,

Thank you Jason. I have been applying Glauret’s correction as you described in your comment but this method converges only for very small variation in tip speed ratio and pitch angles. Values become complex most probably due to negative value of CT(from eq 1) in calculating axial induction factor through Buhl’s expression(eq 3) during iterations. Any idea why this problem arises?
However i have been able to solve this problem by using the approach shown in the attachment for value of ac=.2(also advised in the book). But i had to use the absolute value of “K”. This approach gives convergence for a large range and graphs of Coefficient of performance vs tip speed ratio obtained match with the graphs given in several books. But i am not sure if this approach is mathematically and aerodynamically correct.

Dear Zeeshan,

I’m not sure I understand your comment about “values become complex most probably due to negative value of CT(from eq 1) in calculating axial induction factor through Buhl’s expression(eq 3) during iterations.” Eq. (3) is only used when CT > 0.96F, so, CT cannot be negative in Eq. (3).

I’m not familiar with the solution approach in “Aerodynamics of Wind Turbines” by Hansen, so, perhaps someone else on the forum can comment about that. I’m most familiar with the approach summarized in our AIAA SciTech 2015 paper linked above, which is the approach we’ve implemented in AeroDyn v15 (although AeroDyn v15 includes a few modifications to the procedure described in the paper). I do see that “K” in your Eq. (6.45) is the inverse of “K” in Eq. (8) in our AIAA SciTech 2015 paper. In the later, “K” is allowed to go negative, although numerical problems arise at K = -1.

Best regards,

Sorry for not stating it properly. What i meant to say was that if i use the algorithm decribed by you in your reply, the value of Induction factor becomes complex number( imaginary number) during iterative process.
Thank you again for your valuable time. I am looking forward to reply from other members.