How to get iterations to converge in while loop

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

Answer 1

[Enrique Montero]

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.