Why is my loop not ending?

3 ビュー (過去 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 ?

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

製品


リリース

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by