While Loop not converging??

7 次查看（过去 30 天）
Kevin
Kevin 2014-7-22
评论： Geoff Hayes 2014-7-22
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
  1 个评论
Geoff Hayes
Geoff Hayes 2014-7-22
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Michael Haderlein
Michael Haderlein 2014-7-22
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

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by