Course AE4-135 Rotor and Wake Aerodynamics
Press right arrow to continue.
On the bottom right there are four directional arrows.
At the lower right corner, you can find the slide number (horizontal number. vertical number).
Press bottom arrow to test the vertical slide motion.
Use the side arrows to move horizontally.
To be able to program a BEM code for application to an horizontal axis wind turbine...
BEM segments the blades into blade elements, to calculate 2D loads in each element.
The circular region swept by the blade elements at a given radial position defines an annulus.
Actuator disk momentum theory is applied at each annulus to determine induction as a function of the loading on the blade segments.
BEM uses an iterative approach:
The figure on our right shows the equivalent actuator disk that represents the area swept by the blades; in detail, the annulus swept by the blade elements.
Velocity perceived at blade element: axial, tangential and total.
Forces at blade element (two dimensional): lift, drag, axial and azimuthal.
Note that all forces are in 2D, and have dimensions $N/m$.
$C_l$ and $C_d$ can be retrieved from airfoil polar data, as a function of angle of attack $\alpha$.
Press bottom arrow to see an example of an airfoil polar.
The axial induction $a$ is a function of the axial loading on the annulus. The thrust coefficient is given by
$C_T = \frac{F_{axial}N_{blades} \delta r}{\frac{1}{2}\rho U_\infty^2 A_a} $, with $A_a = 2 \pi r \delta r$
with $C_T=4 a \left( 1 -a \right)$ and $a = \frac{1}{2}\left( 1 - \sqrt{1-C_T} \right)$
The azimuthal induction $a'$ is a function of the tangential loading on the annulus (torque), given by $ \delta Q = F_{azim} r \delta r N_{blades} $
From momentum theory, we obtain $\delta Q=\rho \left(2\pi r\delta r \right)U_{\infty}\left(1-a\right) \left(2\Omega a'r^{2}\right)$
resulting in $a' \Omega r = \frac{F_{azim} N_{blades}}{2 \rho \left(2\pi r \right)U_{\infty}\left(1-a\right) }$ and $a' = \frac{F_{azim} N_{blades}}{2 \rho \left(2\pi r \right)U_{\infty}^2 \left(1-a\right) \lambda \frac{r}{R} }$
Our original assumption of the relation between loading on the actuator and induction is given by the equations below.
$C_T=4 a \left( 1 -a \right)$
and
$a = \frac{1}{2}\left( 1 - \sqrt{1-C_T} \right)$
Above a certain loading on the streamtube, the wake will enter "turbulent wake state" or even "vortex ring state". The relations above are no longer valid, and should be corrected.
Press bottom arrow to see Glauert's correction for heavily loaded rotors.
One proposed correction was defined by Glauert, and it is given by the equations below.
$ C_T= \left\{ {\begin{array}{*{20}{c}} {4 a \left( 1 -a \right) \rm{,~for~} a< 1- \frac{\sqrt{C_{T_1}}}{2} }\\ {C_{T_1}-4 \left(\sqrt{C_{T_1}} -1 \right) \left( 1 -a \right) \rm{,~for~} a \ge 1- \frac{\sqrt{C_{T_1}}}{2} } \end{array}} \right. $
$ a = \left\{ {\begin{array}{*{20}{c}} { \frac{1}{2} - \frac{\sqrt{1-C_{T}}}{2} \rm{,~for~} C_T < C_{T_2} }\\ {1+\frac{C_T-C_{T_1} }{4\sqrt{C_{T_1}} - 4} \rm{,~for~} C_T \ge C_{T_2} } \end{array}} \right. $
and
$C_{T_1}=1.816$ and $C_{T_2}= 2\sqrt{C_{T_1}} - C_{T_1}$
Press bottom arrow to see the implementation of the function.
Press bottom arrow to see the source code of the function.
function thrust_coefficient_from_induction_Gluert_correction(a) {
// calculates thrust coefficient as a function of induction factor a
var CT=[];
var CT1=1.816;
a1=1-Math.sqrt(CT1)/2;
for (var i = 0; i < a.length; i++) {
if (a[i] < a1) {
temp = 4*a[i]*(1-a[i]);
} else {
temp = CT1-4*(Math.sqrt(CT1)-1)*(1-a[i]);
};
CT.push(temp);
};
return CT;
};
The actuator disk model assumes an uniform distribution of loading over the actuator. A finite number of blades leads to a concentration of loading on the blades, and the shedding of concentrated tip and root vortices. The concentrated tip vortex will cause an increased induction at the tip and the root, which is also experienced over the remaining of the blade.
Press bottom arrow to see Prandtl's correction for finite number of blades.
Prandtl proposed a tip and root correction of induction, which increases induction as a function of the distance to the tip and root, and the proximity of the vortices to the blade, a function of the number of blades and the local tip speed ratio. The tip and root corrections are given by the expressions:
$ f_{tip \left( \mu \right)} =\frac{2}{\pi} \arccos{\left( e^{ \left(-\frac{B}{2} \left(\frac{1-\mu}{\mu} \right) \sqrt{\left(1+\frac{\lambda^2 \mu^2}{\left(1 -a \right)^2} \right)} \right)} \right)} $
$ f_{root \left( \mu \right)} =\frac{2}{\pi} \arccos{\left( e^{ \left(-\frac{B}{2} \left(\frac{\mu-\mu_{root}}{\mu} \right) \sqrt{\left(1+\frac{\lambda^2 \mu^2}{\left(1 -a \right)^2} \right)} \right)} \right)} $
where $\mu = r/R$ (radial position r over radius R), and $\mu_{root}$ is the location of the root vortex
The tip and root effect are combined in a total correction $ f_{total \left( \mu \right)} = f_{tip \left( \mu \right)} * f_{root \left( \mu \right)}$
The corrected values of induction are given by $a=\frac{a}{f_{total}}$ and $a'=\frac{a'}{f_{total}}$
Press bottom arrow to see the implementation of the function.
Press bottom arrow to see the source code of the function.
function calculatePrandtlTipRootCorrection(r_R, rootradius_R, tipradius_R, TSR, NBlades, axial_induction) {
// applies Prandtl's tip and root correction to the induction vector aind
var mu; // temp variabl
var Ftip; // tip correction
var Froot; // tip correction
var result;
var temp1;
temp1 = -NBlades/2*(tipradius_R-r_R)/r_R*Math.sqrt( 1+ (Math.pow(TSR*r_R,2))/(Math.pow(1-axial_induction,2)) )
Ftip = 2/Math.PI*Math.acos(Math.exp(temp1));
if (isNaN(Ftip)) {
Ftip = 0;
};
temp1 = NBlades/2*(rootradius_R-r_R)/r_R*Math.sqrt( 1+ (Math.pow(TSR*r_R,2))/(Math.pow(1-axial_induction,2)) )
Froot = 2/Math.PI*Math.acos(Math.exp(temp1));
if (isNaN(Froot)) {
Froot = 0;
};
result= {Ftotal: (Froot*Ftip), Ftip: Ftip ,Froot: Froot};
return result;
}
BEM assumes that all annuli are independent; we will solve each annulus independently. For each annulus, the iterative algorithm on the right is applied.
Suggestion 1: initialize the iteration of by seeding the values of $a=0.3$ and $a'=0$.
Suggestion 2: to help convergence: at each iteration $i$, correct the new values of induction by $a_i= 0.25a_i+0.75a_{i-1}$ and ${a'}_i= 0.25{a'}_i+0.75{a'}_{i-1}$.
Suggestion 3: if required, bound axial induction to $a<0.95$.
Suggestion 4: bound cases of "not a number" and "divide by zero".
Press bottom arrow to see an example of code that represents this algorithm.
function solveStreamtube(Uinf, r1_R, r2_R, rootradius_R, tipradius_R , Omega, Radius, NBlades ){
// solve balance of momentum between blade element load and loading in the streamtube
// input variables:
// Uinf - wind speed at infinity
// r1_R,r2_R - edges of blade element, in fraction of Radius ;
// rootradius_R, tipradius_R - location of blade root and tip, in fraction of Radius ;
// Radius is the rotor radius
// Omega -rotational velocity
// NBlades - number of blades in rotor
// initialize properties of the blade element, variables for output and induction factors
var r_R = (r1_R+r2_R)/2; //centroide
var Area = Math.PI*(Math.pow(r2_R*Radius,2)-Math.pow(r1_R*Radius,2)); // area streamtube
var a = 0.3; // axial induction factor
var anew; // temp new axial induction factor
var aline = 0.; // tangential induction factor
var Urotor; // axial velocity at rotor
var Utan; // tangential velocity at rotor
var loads; // normal and tangential loads 2D
var load3D = [0 , 0]; // normal and tangential loads 3D
var CT; //thrust coefficient at streamtube
var Prandtl; // Prandtl tip correction
// iteration cycle
var Niterations =100; // maximum number of iterations
var Erroriterations =0.00001; // error limit for iteration rpocess, in absolute value of induction
for (var i = 0; i < Niterations; i++) {
///////////////////////////////////////////////////////////////////////
// this is the block "Calculate velocity and loads at blade element"
///////////////////////////////////////////////////////////////////////
Urotor = Uinf*(1-a); // axial velocity at rotor
Utan = (1+aline)*Omega*r_R*Radius; // tangential velocity at rotor
// calculate loads in blade segment in 2D (N/m)
loads = loadBladeElement(Urotor, Utan, r_R);
load3D[0] =loads[0]*Radius*(r2_R-r1_R)*NBlades; // 3D force in axial direction
load3D[1] =loads[1]*Radius*(r2_R-r1_R)*NBlades; // 3D force in azimuthal/tangential direction (not used here)
///////////////////////////////////////////////////////////////////////
//the block "Calculate velocity and loads at blade element" is done
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// this is the block "Calculate new estimate of axial and azimuthal induction"
///////////////////////////////////////////////////////////////////////
// calculate thrust coefficient at the streamtube
CT = load3D[0]/(0.5*Area*Math.pow(Uinf,2));
// calculate new axial induction, accounting for Glauert's correction
anew = induction_from_thrust_coefficient_Gluert_correction([CT]);
// correct new axial induction with Prandtl's correction
Prandtl=calculatePrandtlTipRootCorrection(r_R, rootradius_R, tipradius_R, Omega*Radius/Uinf, NBlades, anew);
if (Prandtl.Ftotal < 0.0001) { Prandtl.Ftotal = 0.0001; } // avoid divide by zero
anew = anew/Prandtl.Ftotal; // correct estimate of axial induction
a = 0.75*a+0.25*anew; // for improving convergence, weigh current and previous iteration of axial induction
// calculate aximuthal induction
aline = loads[1]*NBlades/(2*Math.PI*Uinf*(1-a)*Omega*2*Math.pow(r_R*Radius,2));
aline =aline/Prandtl.Ftotal; // correct estimate of azimuthal induction with Prandtl's correction
///////////////////////////////////////////////////////////////////////////
// end of the block "Calculate new estimate of axial and azimuthal induction"
///////////////////////////////////////////////////////////////////////
// test convergence of solution, by checking convergence of axial induction
if (Math.abs(a-anew) < Erroriterations) {
i=Niterations; // converged solution, this is the last iteration
}
};
// we have reached a solution or the maximum number of iterations
// returns axial induction factor a, azimuthal induction factor a',
// and radial position of evaluations and loads
return [a , aline, r_R, loads[0] , loads[1]];
};
Scroll the button below to change tip speed ratio. Tip speed ratio = 6.
Press bottom arrow to see other examples.
Scroll the button below to change tip speed ratio. Tip speed ratio = 6.
$~~~~C_T$= 0
$~~~~C_{Protor}$= 0
Suggestion: compare $F_{axial}$ and $F_{azim}$. Derive the relation.
Note: $F_{axial}$ and $F_{azim}$ are non-dimensioned by $\frac{1}{2} \rho U_\infty^2 R$
Scroll the button below to change tip speed ratio. Tip speed ratio = 6.
$~~~~C_T$= 0
$~~~~C_{Protor}$= 0
Suggestion: compare $\Gamma$ and $C_T$. Derive the relation.
Note: The circulation $\Gamma$ is non-dimensioned by $\frac{ \pi U_\infty^2 }{N_{Blades} \Omega}$
To be able to program a BEM code for application to an horizontal axis wind turbine...
Use the simulation tools to compare the relations between the two terms of induction, loading, circulation, and tip speed ratio.
Derive the relations between the two terms of induction, loading, circulation, and tip speed ratio.