While Loop not converging??

Hi,

I am using a while loop in order to perform a tolerance test. When tha calculated check value becomes greater than an assigned tolerance value, the loop should stop. However my check value appears to contain imaginary numbers which stops it from converging. Does anybody know how I could fix this? Or can anybody see any simple mistakes I have made as I am still relatively mnew to the world of Matlab? I have attached my code for clarity:

% Inputs
R=0.4;                                          % Radius of Rotor
B=3;                                            % Number of blades
U=1;                                            % Fluid velocity
Rho=998;                                        % Fluid Density
N=9;                                            % Number of Blade Elements
Cp_estimate=0.5;                                % Estimate power coefficient
Alpha_design=4;                                 % Design alpha
Cl_design=1;                                    % Design lift coefficient

% Variables

i=1;                                            % Counter
alpha_new=0;                                    % Initial value for alpha new
tolerance=0.00001;                                 % Tolerance Value
axial_induction=[0;0;0;0;0;0;0;0;0];
Check=1;                                        % Initial check value
axial_induction_old=0;                          % Initial value for old axial induction factor

Cl=[1.3; 1.1; 1; 0.9; 0.86; 0.83; 0.8; 0.75; 0.5];                                          % Lift Coefficients
Cd=[0.027; 0.024; 0.02; 0.019; 0.018; 0.016; 0.013; 0.012; 0.01];                           % Drag Coefficients

r_local=R/9*(1:9)';
r_over_R=r_local / R;

ii=0;                                               % Counting Interval for integral values

for TSR=0:0.1:10                                     % TSR from 1 to 10
    disp(TSR)                                                                    
    Check=1;
ii=ii+1;

TSR_local=r_over_R .* TSR;    
Phi=(2/3)*atan(1./TSR_local);
C=((8*pi.*r_local) ./ (B.*Cl_design)).*(1-cos(Phi));
sigma=(B*C) ./ (pi.*r_local.*2);

axial_induction= 1 ./ (((4.*(sin(Phi).^2)) ./ (sigma.*Cl_design.*cos(Phi)))+1);
angular_induction= (1-(3*axial_induction)) ./ ((4.*axial_induction)-1);
angular_velocity=TSR.*U./R;

     while abs(Check)>=tolerance

relative_wind = atan(U.*(1-axial_induction))./(angular_velocity*R).*(1+angular_induction);

F=(2/pi) .* acos(exp(-(((B/2) .* (1-(r_over_R))) ./ ((r_over_R) .* sin(relative_wind)))));      % Tip Loss Factor

C_T=(sigma .* ((1-axial_induction).^2) .* ((Cl.*cos(relative_wind))+(Cd.*sin(relative_wind)))) ./ ((sin(relative_wind)).^2);     

axial_induction_old = axial_induction;        
TSR_local = TSR .* (r_local./R);                                                                 % Local Tip Speed Ratio
Phi = (2/3) .* atan(1./TSR_local);                                                               % Angle of Relative Fluid

for i=1:length(C_T)
   if  C_T(i) > 0.96;
      axial_induction(i) = 1 / (((4*F(i)*cos(relative_wind(i))) / (sigma(i)*Cl(i)))-1);

     else
         axial_induction(i) = 1 / (1+(4*F(i)*(sin(relative_wind(i))^2)) / (sigma(i)*Cl(i)*cos(relative_wind(i)))); 
     end;
  end;

D=(8./(TSR.*9)).*(F.*(sin(Phi).^2).*(cos(Phi)-((TSR_local).*(sin(Phi)))).*(sin(Phi)+((TSR_local).*(cos(Phi)))).*(1-(Cd./Cl).*cot(Phi)).*(TSR_local.^2));
Cp=sum(D);

Diff=axial_induction-axial_induction_old;
Check=max(Diff(:));
Check

     end

store_sigma(:,ii)=sigma;
store_Phi(:,ii)=Phi;
store_TSR_local(:,ii)=TSR_local;
store_axial_induction(:,ii)=axial_induction;
store_angular_induction(:,ii)=angular_induction;
store_axial_induction_old(:,ii)=axial_induction_old;
store_relative_wind(:,ii)=relative_wind;
store_Check(:,ii)=Check;
store_Diff(:,ii)=Diff;
store_Cp(:,ii)=Cp;
store_TSR(:,ii)=TSR;
store_F(:,ii)=F;

end

[Geoff Hayes]

Kevin - you state that When tha calculated check value becomes greater than an assigned tolerance value, the loop should stop whereas your code is doing the opposite

 while abs(Check)>=tolerance

continuing so long as the check value is greater than the assigned tolerance.

Note that when TSR is zero (on the first iteration of the outer for loop) there is a division by zero error which leads to a couple of vectors being set to all Inf and (later) all NaN values. You may want to skip this first case.

Your code is getting imaginary numbers at the calculation of F

 F=(2/pi) .* acos(exp(-(((B/2) .* (1-(r_over_R))) ./ ...
      ((r_over_R) .* sin(relative_wind)))));

Has this equation been written correctly?

Answer 1

[Michael Haderlein]

If you have imaginary numbers, maybe the problem comes from the acos in the calculation of F? acos is rad based (not degree). Could that be the error? If you want to use degree, use the acosd function.

If that's not the problem, you should check stepwise where the imaginary numbers appear first (use breakpoints and F10).

Best,

Michael