Why is my loop not ending?

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 ?

Answer 1

[Image Analyst]

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')

[smith]

Thank you!

Answer 2

[Alan Stevens]

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?

[smith]

No, this is a diameter, it should'nt be in negative, this line should get me to the new diameter of the particle at each iteration, until this diameter reaches zero theoritically. but it is reaching a place where it goes imaginary so the loop isn't really stopping.

[Alan Stevens]

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]

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 ?