How to get iterations to converge in while loop

6 views (last 30 days)
Hi,
I am trying to write a code to calculate the induction factors along a marine turbine blade, using the Blade Element Momentum Theory. The blade is divided into 9 elements and I am trying to calculate the induction factors at each element. The iteration process is to have the difference between (a_old and a_new) and (aa_old and aa_new) to be less than the tolerance, before it should move onto the next blade element and the process repeated until all the induction factors are calculated for each element. The code seems to calculate along the blade, but as soon as the tolerance is met for one element it jumps out of the loop. What do i need to do so that the while loop condition is met for each blade elemnt. Delta_aa is very close to the tolerance but delta_a is nowhere near the tolerance criteria.
clc
clear
R=0.4; %Radius Rotor
B=3; %Number of blades
U=1.73; %Fluid velocity
rho=998; %Fluid density
Ne=9; %Number of blade elements
tolerance=0.001; %Tolerance for induction factors difference
dr=0.04; %Element width
BladeElementRadii=([0.08 0.12 0.16 0.2 0.24 0.28 0.32 0.36 0.4]); %Blade element radii
BladeElementPitch=([20 14.5 11.1 8.9 7.4 6.5 5.9 0.4 0]);
%Blade element pitch
BladeElementPitch_rad=BladeElementPitch.*(pi()/180);
Chord_ratio=([0.125 0.1156 0.1063 0.0969 0.0875 0.0781 0.0688 0.0594 0.05]);
%Blade element chord lengths
Chord=Chord_ratio.*0.4;
%Angle of attack
AoA=[-10;-9;-8;-7;-6;-5;-4;-3;-2;-1;0;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20];
%Drag coefficients for corresponding angle of attack
Drag=
[0.021;0.019;0.018;0.017;0.015;0.015;0.015;0.015;0.015;0.015;0.015;0.022;0.032;0.037;0.041;0.045;0.05;0.056;0.059;0.068;0.082;0.085;0.089;0.105;0.125;0.143;0.161;0.176;0.191;0.195;0.2];
%Lift coefficients for corresponding angle of attack
Lift=[-0.5;-0.42;-0.35;-0.22;-0.08;0;0.11;0.21;0.31;0.42;0.5;0.58;0.67;0.76;0.86;0.95;1.07;1.14;1.22;1.26;1.32;1.35;1.37;1.36;1.35;1.33;1.32;1.31;1.29;1.28;1.27];
for Angular_Velocity=1:0.1:45
TSR=(Angular_Velocity*R)/U;
Local_Tip_Speed_Ratio=(Angular_Velocity*BladeElementRadii)./U;
a_old=0;
aa_old=0;
delta_a=1;
delta_aa=1;
itr=0;
%Step a: Calculate the solidity
solidity=zeros(1,9);
solidity=(B*Chord)./(2*pi().*BladeElementRadii);
for Ne=1:9
while delta_a & delta_aa > tolerance
itr=itr+1;
% Step 1: Calculate angle of relative fluid
phi=atan((1-a_old)./((1+aa_old).*(Local_Tip_Speed_Ratio)));
phi_degrees=phi./(pi()/180);
% Step 2: Calculate Prandtls Tip loss correction
F=(2/pi)*acos(exp(-(((B/2).*(R-(BladeElementRadii)))./((BladeElementRadii).*sin(phi)))));
% Step 3: Calculate the angle of attack
alpha=phi-BladeElementPitch_rad;
alpha_degrees=alpha.*(180/pi());
% Calculate Lift and Drag coefficients for angle of attack >20 degrees
if alpha_degrees >20
C_d=2.*sin(phi_degrees).^2;
C_l=2.*sin(phi_degrees).*cos(phi);
else C_d=interp1(AoA,Drag,(alpha));
C_l=interp1(AoA,Lift,(alpha));
end
% Step 4: Calculate the normal and tangential load coefficients
C_n=(C_l.*cos(phi))+(C_d.*(sin(phi)));
C_t=(C_l.*sin(phi))-(C_d.*(cos(phi)));
%Step 5: Calculate the solidity
solidity=(B.*Chord)/(2*pi().*BladeElementRadii);
% Step 6: Calculate new values for induction factors
a_new=1./((1+((4.*sin(phi).^2)./(solidity*C_n))));
aa_new=1./(((4.*F.*sin(phi).*cos(phi))./(solidity.*C_t))-1);
% Step 7: Calculate the difference between itterations
delta_a=abs(a_new-a_old);
delta_aa=abs(aa_new-aa_old);
a_old=a_new; %Set the old axial induction factor to the new value
aa_old=aa_new; %Set the old angular induction factor to the new value
end
end
dT=4*pi().*BladeElementRadii.*rho*U.^2.*a_new.*(1-a_new).*dr.*F; %Thrust
dM=4*pi().*BladeElementRadii.^3.*rho.*Angular_Velocity.*aa_new.*(1-a_new).*dr.*F; %Torque
dP=dM*Angular_Velocity; %Power
M=sum(dM);
P=sum(dP);
P_a=0.5*rho*pi()*R^2*U^3;
C_p=P/P_a;
end

Answers (1)

Enrique Montero
Enrique Montero on 23 Aug 2018
Did you solve the problem? I have to develop the same code and i am having the same problem, i would really appreciate it if you share it! thank you beforehand.

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!