While Loop not converging??
1 回表示 (過去 30 日間)
古いコメントを表示
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
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
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
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Creating and Concatenating Matrices についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!