Appendices
Appendix A: Complete BEM code
%%%%%% BEM code, with Empirical model for the axial induced velocities in yaw,
including the 3D corrections %%
addpath BL
close all
clear all
%input quantities on geometric conditions
Rt=2.25; %[m]
Rh=0.21; %[m]
omega=424.4*2*pi/60; %[rad/s]
U_inf=24; %[m/s]
density=1.225; %[kg/m3]
B=3;
N=8;
angular_increment=10;
pitch=-2.3; % in degrees
yaw=30; % in degrees
tolerance = 0.0005;
k_visc=15.1e-6;
Nrev=3 ; % number of revolutions considered
DS=1; % flag to activate the Dynamic Stall model
C3D=1; % flag to activate the 3D correction of the force coefficients
CTL=1; % flag to activate the Tip loss correction
% derived quantities
lbd=Rt*omega/(U_inf*cos(yaw*pi/180));
dt=angular_increment/(omega*180/pi);
%Derived quantities
element_length=(Rt-Rh)/N;
r=Rh+(Rt-Rh)/(2*N):element_length:Rt; % the actual radius of each element
section
r_real_mat=[r r r]; % for the 3 blades
r_frac=(r-Rh)/(Rt-Rh); % fraction of the tip radius
r_mat=[r_frac r_frac r_frac]; % for 3 blades, i.e. relative radius from
the hub
r_frac_real=(r)/(Rt); % actual fraction of the tip radius
r_mat_real=[r_frac_real r_frac_real r_frac_real]; % for 3 blades
span_MEX_a=[0.25 0.35 0.6 0.82 0.92]; % spanwise positions of the MEXICO,
divided by Rt
x_positions_for_interp=[0.21 0.23 0.235 0.3 0.45 0.675 0.9 1.025 1.125 1.225
1.35 1.475 1.575 1.675 1.8 2.025 2.165 2.193 2.222 2.25];
c_points=[0.195 0.195 0.09 0.09 0.24 0.207 0.178 0.166 0.158 0.15 0.142 0.134
0.129 0.123 0.116 0.102 0.092 0.082 0.056 0.011];
c=interp1(x_positions_for_interp,c_points,r,'cubic');
c_mat=[c c c];
c_MEX=interp1(x_positions_for_interp,c_points,span_MEX_a*Rt,'cubic');
twist_points=[16.4 16.4 16.4 16.4 16.4 12.1 8.3 7.1 6.1 5.5 4.8 4 3.7 3.2 2.6
1.5 0.7 0.469 0.231 0];
twist=interp1(x_positions_for_interp,twist_points,r,'cubic');
twist_mat=[twist twist twist];
ReMEX=(sqrt((omega*(span_MEX_a*(Rt))).*(omega*(span_MEX_a*(Rt)))+
(U_inf*cos(yaw*pi/180))^2)).*c_MEX/k_visc;
%%% setting up the azimuth angle, because it is a constant matrix
for j=1:1:1+Nrev*(360/angular_increment)
for i=1:1:B*N
az(i,j)=(j-1)*angular_increment;
if i>N
az(i,j)=120+(j-1)*angular_increment;
end
if i>2*N
az(i,j)=240+(j-1)*angular_increment;
end
while az(i,j)>360
az(i,j)=az(i,j)-360;
end
end
end
% initializing variables for the BL model
% these are the initial values that describe the state
variables=zeros(1,B*N,14);
variables(1,:,6)=1.2; %Cn
variables(1,:,5)=1; % f parameter
variables(1,:,9)=1; % f' parameter
%%% we iterate for the average axial induced velocity
for j=1:1:1+Nrev*(360/angular_increment)
if j==1 % in the first iteration the value of the axial induction factor is
assumed
a0_av=0.6-0.02*U_inf*cos(yaw*pi/180); % from normal values
a_av=a0_av*0.99;
count_jota=0;
else % if it is not the first time instant, we simply take the induction
factor from the previous iteration
a0_av=a(j-1);
a_av=0.99*a(j-1);
end
if rem(j,16)==0 % to keep track of the iteration
count_jota=j
end
while abs(a0_av-a_av)>tolerance
for i=1:1:B*N
Vax(i)=emp_model2(U_inf,yaw,a_av,r_mat(i),az(i,j)); % axial
velocity, from the empirical model
Vtg(i)=r_real_mat(i)*omega-
U_inf*sin(yaw*pi/180)*cos(az(i,j)*pi/180); % tangential velocity, including
skewed inflow
phi(i)=(180/pi)*atan(Vax(i)/Vtg(i));
if Vax(i)<0 % these lines are included to prevent the iteration from
diverging at higher yaw errors
phi(i)=0;
end
if Vtg(i)<0
phi(i)=90;
end
alpha(i,j)=phi(i)-twist_mat(i)-pitch;
Vtot(i)=sqrt(Vax(i)^2+Vtg(i)^2);
%%%%% to compute the Cl,Cd we must choose the airfoil section %%%%
airfoil_type_v(i)=2; % the transition between airfoils is measured
with tip radius fraction, from
if r_mat(i)<0.4485 % the hub
airfoil_type_v(i)=1;
end
if r_mat(i)>0.6691
airfoil_type_v(i)=3;
end
airfoil_type=airfoil_type_v(i);
%% getting the critical normal force coefficient
Cn1(i)=getCnCrit(airfoil_type,r_mat(i),c_mat(i),twist_mat(i)+pitch);
%%% using 2D aerodynamic coefficients
if airfoil_type==1
twodimadd=DUW;
elseif airfoil_type==2
twodimadd=RISO;
else
twodimadd=NACA;
end
Cl_mat(i)=interp1(twodimadd(:,1),twodimadd(:,2),alpha(i,j),'cubic','
extrap');
Cd_mat(i)=interp1(twodimadd(:,1),twodimadd(:,3),alpha(i,j),'cubic','
extrap');
%%%
%%% to compute the 3 dimensional aerodynamic coefficients
if C3D==1
Cl_mat(i)=get3DCl(airfoil_type,r_mat(i),c_mat(i),alpha(i,j));
Cd_mat(i)=get3DCd(airfoil_type,r_mat(i),c_mat(i),twist_mat(i)
+pitch,alpha(i,j));
end
Cn_mat(i)=Cl_mat(i)*cos(alpha(i,j)*pi/180)+Cd_mat(i)*sin(alpha(i,j)*
pi/180); % computation of force normal to the airfoil
%%% the dynamic stall is introduced
if DS==1;
if j>1
intermediate=BLnew(alpha(i,j),alpha(i,j-
1),Vtot(i),dt,c_mat(i),Cn1(i),airfoil_type,variables(j-
1,i,:),r_mat(i),twist_mat(i)+pitch);
Cn_mat(i)=intermediate(14);
for index=1:1:length(intermediate)
variables(j,i,index)=intermediate(index);
end
end
end
Ctg_mat(i)=Cl_mat(i)*sin(alpha(i,j)*pi/180)-
Cd_mat(i)*cos(alpha(i,j)*pi/180); % computation of the forcel tang to the
airfoil
%%% including the tip loss, from Shen et al., applied in the
%%% normal and tangential coefficients
if CTL==1;
Cn_mat(i)=TipCorrection(Cn_mat(i),lbd,B,phi(i),Rt,r_mat_real(i))
;
Ctg_mat(i)=TipCorrection(Ctg_mat(i),lbd,B,phi(i),Rt,r_mat_real(i
));
end
%the momentum balance is done with the corrected load
%coefficients
%Cx_mat(i)=Cl_mat(i)*cos(phi(i)*pi/180)+Cd_mat(i)*sin(phi(i)*pi/180)
; % computation of force perpend to rotor plane
Cx_mat(i)=Cn_mat(i)*cos((phi(i)-
alpha(i,j))*pi/180)+Cd_mat(i)*sin((phi(i)-alpha(i,j))*pi/180); % computation of
force perpend to rotor plane
T_mat(i)=Cx_mat(i)*0.5*density*(Vax(i)^2+Vtg(i)^2)*c_mat(i)*element_
length;
Nf_mat(i)=Cn_mat(i)*0.5*density*(Vax(i)^2+Vtg(i)^2)*c_mat(i)*element
_length;
Tf_mat(i)=Ctg_mat(i)*0.5*density*(Vax(i)^2+Vtg(i)^2)*c_mat(i)*elemen
t_length;
end
T(j)=sum(T_mat); % here we sum the contribution of all
the blades to the axial force
Ct(j)=T(j)/(pi*(Rt^2)*0.5*density*(U_inf*cos(yaw*pi/180))^2);
if Ct(j)>0.8889 % we obtain 'a' from Ct, using the
Glauert correction for
roots_a=roots([3 -5 4 -Ct(j)]); % turbulent wake state
a(j)=1;
for r=1:1:3
if abs(roots_a(r))<a(j)
a(j)=abs(roots_a(r));
end
end
else
a(j)=(1-(1-Ct(j))^0.5)/2; % otherwise
end
a0_av=a_av; % update the average induction factor
a_av=0.9999*a_av+0.0001*a(j);
end
for k=1:1:length(span_MEX_a) % interpolation of the forces at the
MEXICO span positions
NfBEM(j,k)=(1/element_length)*interp1(r_mat_real(1:N),Nf_mat(1:N),span_M
EX_a(k),'cubic');
CnBEM(j,k)=interp1(r_mat_real(1:N),Cn_mat(1:N),span_MEX_a(k),'cubic');
if j>1
%CnBEM2(j,k)=interp1(r_mat(1:N),Cn_mat2(1:N),span_MEX_a(k),'cubic');
end
TfBEM(j,k)=(1/element_length)*interp1(r_mat_real(1:N),Tf_mat(1:N),span_M
EX_a(k),'cubic');
CtgBEM(j,k)=interp1(r_mat_real(1:N),Ctg_mat(1:N),span_MEX_a(k),'cubic');
alphaMEX(j,k)=interp1(r_mat_real(1:N),alpha(1:N,j),span_MEX_a(k),'cubic'
);
end
end
%
span=1; % we choose the span to plot
grid on
hold on
axis ([0 360 -0.1 (max(CnBEM(:,span))+0.5)])
plot(az(1,1:1+360/angular_increment),CnBEM(((Nrev-
1)/Nrev)*length(CnBEM(:,1)):end,span),'r');
surf(CnBEM(((Nrev-1)/(Nrev))*length(CnBEM(:,1)):end,:));
Appendix B - DS model main code
%%% this function computes the normal force coefficient using the BL dynamic
stall
%%% model for a given aoa history
function [variables]=
BLnew(alfa,alfa_old,V,dt,chord,Cn1,airfoil,variables,radius,t_pitch)
%derived quantity
delta_s=V*dt/(chord/2); % semichords travelled in each time step
% steady data
[cn_alfa, alfa_t, f_t, cn0] = steady_t_e_sep(airfoil,chord,radius,t_pitch);
save f_coeff cn_alfa alfa_t f_t cn0
pachacha=[alfa_t f_t];
% variables needed from previous condition
X_old=variables(1);
Y_old=variables(2);
Dl_old=variables(3);
Df_old=variables(4);
f_prime_old=variables(5);
cn_old=variables(6);
cn_prime_old=variables(7);
tv=variables(8);
f_double_prime_old=variables(9);
cnv_old=variables(10);
cv_old=variables(11);
delta_alfa_old=variables(12);
D_imp_old=variables(13);
% if flag==1 % on the first time instant initial values are assigned
% [X_old Y_old Dl_old Df_old f_prime_old cn_old cn_prime_old tv
f_double_prime_old cnv_old cnv cv_old delta_alfa_old
D_imp_old]=inicio1(cn_alfa,alpha_mean);
% else
% load cenas1 % loading the values calculated in the previous iteration
% end
% geometry
c = chord; % chord length [m]
% Attached flow constants
%
A1 = 0.3;
A2 = 0.7;
b1 = 0.14;
b2 = 0.53;
etha = 0.95; % recovery factor
cn_1 = Cn1; % critical value for the leading edge pressure
% time-constants
Tp = 1.5; % peak pressure - cn lag
Tf = 5.0; % boundary layer - peak pressure lag
Tv =6.0; % vortex decay constant
Tv1 = 5.0; % trailing edge position, in semichords
delta_alfa=alfa-alfa_old;
[alfa_e, cn, cc, X_old, Y_old, D_imp_old,cn_i,cnc]=unsteady_attached(V, c, alfa,
dt, X_old, Y_old, alfa_old, A1, A2, b1, b2,
cn_alfa,cn0,delta_alfa,delta_alfa_old,D_imp_old);
[alfa_f, cn_prime, cn_f, cc_f, f_double_prime, Dl_old, Df_old,
f_prime_old,cn_old]=unst_t_e_sep(Dl_old, Df_old, f_prime_old, cn, cn_old, dt,
Tp, Tf, V, c, cn_alfa, alfa_e, etha, tv, Tv1,delta_s,cn0,cn_i);
[tv, f_double_prime_old]=leading_edge_sep(cn_prime, f_double_prime,
f_double_prime_old, cn_1, dt, tv, V, c);
[cnv,cnv_old,cv_old]=vortex_lift(cnc, dt, Tv, V, c, f_double_prime, cnv_old,
cv_old, tv, Tv1,delta_s);
delta_alfa_old=delta_alfa; % updating of the variable
cn_tot = cn_f + cnv;
save cenas1 X_old Y_old Dl_old Df_old f_prime_old cn_old cn_prime_old tv
f_double_prime_old cnv_old cv_old delta_alfa_old D_imp_old
variables=[X_old Y_old Dl_old Df_old f_prime_old cn_old cn_prime_old tv
f_double_prime_old cnv_old cv_old delta_alfa_old D_imp_old cn_tot cn cn_i cnc
cn_f];
end
Appendix C - Results obtained at other data points and
spanwise stations
Cn vs Azimuth for U=18m/s and Beta=15deg for 25% span Cn vs Azimuth for U=18 m/s and Beta=15 deg for 35% span
1.8
2.5
1.6
2 1.4
1.2
1.5
1
Model Model
Cn
Cn
Model_noDS Model_noDS
0.8
MEXICO MEXICO
1
0.6
0.4
0.5
0.2
0
0
0 45 90 135 180 225 270 315 360
0 45 90 135 180 225 270 315 360
Azimuth angle Azimuth angle
Cn vs Azimuth for U=18m/s and Beta=30 deg at 60% span
Cn vs Azimuth for U=18m/s and Beta=30 deg at 82% span
1.6 1.4
1.4
1.2
1.2
1
0.8
0.8 Model
Cn
Model
Cn
Model_noDS Model_noDS
MEXICO 0.6
MEXICO
0.6
0.4
0.4
0.2 0.2
0 0
0 45 90 135 180 225 270 315 360 0 45 90 135 180 225 270 315 360
Azimuth angle Azimuth angle
Appendix D - Results of validation cases of the 2D BL DS
model
Appendix E - Convergence of BEM code including the DS
model
Convergence of Model w ith Nr Elements at 180 deg azimuth Cn Convergence w ith Angular Increment at 240 deg azimuth
1.7
1.4
1.5
1.2
1.3
1
0.8 25%
1.1 25%
35%
Cn
35% 0.9
Cn
0.6
60% 60%
0.4 0.7 82%
82%
92% 0.5 92%
0.2
0 0.3
6 8 10 12 14 16 18 20 22 24 26 22 20 18 16 14 12 10 8 6 4 2
Nr Elements Angular Increment (deg)
Cn Convergence w ith Angular Increment at 0 deg azimuth
3.8
3.3
2.8
2.3 25%
35%
Cn
1.8
60%
1.3 82%
0.8 92%
0.3
22 20 18 16 14 12 10 8 6 4 2
Angular Increment (deg)
Appendix F - MEXICO Cp distributions at the 35% spanwise
station
Cp variation at data point 152
Cp variation at data point 153
Cp variation at data point 160 Cp variation at data point 167
Appendix G - DS model results at very high reduced
frequencies