Why is my loop not ending?

2 ビュー (過去 30 日間)
smith
smith 2022 年 9 月 23 日
コメント済み: smith 2022 年 9 月 23 日
So, I wrote a code to run iterations for a water particle as it undergoes mass transfer, to monitor the diameter size.
%%% DEFINITIONS %%%
clear
T = 20; %% degrees
t =0.01; %% delta t
RH_out = 0.5; %% 50 %
RH_surf = 1; %% 100%
dp_nano = 200; %% first diameter of the particle in nanometers
cs = 0.0172; %% kg/m3
cinf = (RH_out/RH_surf)*cs;
row = 1000; %% kg/m3 density of water at room temp
D = 2.5*10^-5 %% diffusivity coefficient of water m2/s at room temp
Ru = 8.314; %% kg m2/ mol.s2.K universal gas constant
st = 0.073; %% N/m surface tension at room temp for water
M= 0.018 ;%% molecular mass for water droplet kg/ mol
%%%% INITIAL CONDITIONS %%%%
dp(1) = dp_nano*10^-7; %% diameter in meters
m_particle(1) = (pi/6) * row * (dp(1))^3; %% in kgs
mdot(1) = (pi/6)*D*dp(1) *(cs-cinf) %% Kg/sec
m(1) = (pi/6)*row*(dp(1))^3; %% kg
K = exp ((4*st*M)/(Ru*row*dp(1)* T+273));
i = 0;
newdp = dp(1);
t2(1) = 0
while newdp > 0
i = i+1
newdp = newdp
m(i+1) = m(i) - mdot(i)*t
dp(i+1)= ((6*m(i+1))/(pi*row))^(1/3)
mdot(i+1) = (pi/6)*D*dp(i+1) *(cs-cinf)
t2(i+1) =t2(i) +0.01
newdp = newdp + dp(i+1)
end
When i run the code, the arrays i get from dp, become complex and infinite, my aim is for it just to reach zero. why is it not stopping at zero ?

採用された回答

Image Analyst
Image Analyst 2022 年 9 月 23 日
You could quit the loop if it becomes complex. Not sure why that happens though, but it does. For a while loop you must always set up a failsafe which will limit the number of iterations so you don't get an infinite loop. I show you how to do that.
%%% DEFINITIONS %%%
clear
T = 20; %% degrees
t =0.01; %% delta t
RH_out = 0.5; %% 50 %
RH_surf = 1; %% 100%
dp_nano = 200; %% first diameter of the particle in nanometers
cs = 0.0172; %% kg/m3
cinf = (RH_out/RH_surf)*cs;
row = 1000; %% kg/m3 density of water at room temp
D = 2.5*10^-5 %% diffusivity coefficient of water m2/s at room temp
Ru = 8.314; %% kg m2/ mol.s2.K universal gas constant
st = 0.073; %% N/m surface tension at room temp for water
M= 0.018 ;%% molecular mass for water droplet kg/ mol
%%%% INITIAL CONDITIONS %%%%
dp(1) = dp_nano*10^-7; %% diameter in meters
m_particle(1) = (pi/6) * row * (dp(1))^3; %% in kgs
mdot(1) = (pi/6)*D*dp(1) *(cs-cinf) %% Kg/sec
m(1) = (pi/6)*row*(dp(1))^3; %% kg
K = exp ((4*st*M)/(Ru*row*dp(1)* T+273));
i = 0;
newdp = dp(1);
t2(1) = 0
% Set up a failesafe to bail out if we hit too many iterations.
maxIterations = 100; % Way more than you think it would ever need.
loopCounter = 0;
% Now loop until we obtain the required condition: a random number equals exactly 0.5.
% If that never happens, the failsafe will kick us out of the loop so we do not get an infinite loop.
while newdp > 0 && loopCounter < maxIterations
i = i+1
newdp = newdp
m(i+1) = m(i) - mdot(i)*t
dp(i+1)= ((6*m(i+1))/(pi*row))^(1/3)
mdot(i+1) = (pi/6)*D*dp(i+1) *(cs-cinf)
t2(i+1) =t2(i) +0.01
newdp = newdp + dp(i+1)
% Quit loop if it becomes complex.
if imag(newdp) ~= 0
break;
end
% Increment our failsafe.
loopCounter = loopCounter + 1;
end
plot(dp, 'b-');
xlabel('index')
ylabel('newdp')
  1 件のコメント
smith
smith 2022 年 9 月 23 日
Thank you!

サインインしてコメントする。

その他の回答 (1 件)

Alan Stevens
Alan Stevens 2022 年 9 月 23 日
Shouldn't
dp(i+1)= ((6*m(i+1))/(pi*row))^(1/3);
be
dp(i+1)= -((6*m(i+1))/(pi*row))^(1/3);
i.e. have a negative sign in front?
  3 件のコメント
Alan Stevens
Alan Stevens 2022 年 9 月 23 日
A little more like this perhaps:
%%% DEFINITIONS %%%
clear
T = 20; %% degrees
t =10^-7; %% delta t %%%%%%%%%%%%%%%%%%%%%%%%%
RH_out = 0.5; %% 50 %
RH_surf = 1; %% 100%
dp_nano = 200; %% first diameter of the particle in nanometers
cs = 0.0172; %% kg/m3
cinf = (RH_out/RH_surf)*cs;
row = 1000; %% kg/m3 density of water at room temp
D = 2.5*10^-5; %% diffusivity coefficient of water m2/s at room temp
Ru = 8.314; %% kg m2/ mol.s2.K universal gas constant
st = 0.073; %% N/m surface tension at room temp for water
M= 0.018 ;%% molecular mass for water droplet kg/ mol
%%%% INITIAL CONDITIONS %%%%
dp(1) = dp_nano*10^-9; %%%%%%%%%%%%%% diameter in meters ????? 1 nanometer = 10^-9m
m_particle(1) = (pi/6) * row * (dp(1))^3; %% in kgs
mdot(1) = (pi/6)*D*dp(1) *(cs-cinf); %% Kg/sec
m(1) = (pi/6)*row*(dp(1))^3; %% kg
K = exp ((4*st*M)/(Ru*row*dp(1)* (T+273))); %%%%%%%%%%%%% This doesn't seem to be used anywhere
i = 0;
newdp = dp(1);
t2(1) = 0;
while dp(i+1) > 0
i = i+1;
% newdp = newdp
m(i+1) = m(i) - mdot(i)*t;
dm = mdot(i)*t; %%%%%%%%%%%%%%%%%%%%% mass loss
dp(i+1)= max(dp(i)-((6*dm)/(pi*row))^(1/3),0); %%%%%%%%%%%%%% diameter loss
mdot(i+1) = (pi/6)*D*dp(i+1) *(cs-cinf);
t2(i+1) =t2(i) + t;
%newdp = newdp + dp(i+1);
end
subplot(2,1,1)
plot(t2,dp*10^9, 'o'), grid
xlabel('time [s]'), ylabel('diameter [nm]')
subplot(2,1,2)
plot(t2,m*10^18, '+'), grid
xlabel('time [s]'), ylabel('mass [kg*10^{18}]')
smith
smith 2022 年 9 月 23 日
Hello, Thank you for your answer, but I've got a question;
Why did you choose this specific delta t, or time step ? and if I chose to change it, shouldn't the shape profile stays the same ?

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

製品


リリース

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by