How to remove the zeros from the array upon breaking out of the for loop?

That's the code I've written, supposedly i want the for loop to break once the dp(i+1) starts becoming negative, but apparently what's happening is that it's filling the rest of the array with zeros due to the zeros array, my objective is to have the array stopped with last value once that condition is met without the zeros filling the rest of the array.

clear all
close all

RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]


Tb = 293.15;

pinf =RHo*2.31; % it's constant because Tb is assumed to be constant
cinf = (pinf*1000)*(M) / ((Ru)*(Tb)); %it's constant because pinf and Tb are constants




dt=0.001; % time step
tend = 3; %total time
num_steps = tend/dt; %number of timesteps
Tp = zeros(1,num_steps); %array to store Temperature of droplet at every timestep
dp = zeros(1,num_steps); %array to store diameter at every timestep
mp = zeros(1,num_steps); %array to store mass at every timestep
tp =zeros (1,num_steps);
Tp(1) = 293;
dp(1) = 40000e-9; %initial diameter value of the droplet
mp(1) = (pi/6) *den*(dp(1))^3; %initial value of the masss of the droplet
tp(1)= 0;

for i =1:num_steps-1
    

    
    Kr = exp ((4*st*M)/(Ru*den*dp(i)*Tp(i)));
    psat= RHp*(exp(16.7-(4060/(Tp(i) -37)))); %it's constant because Tp is constant (Relative humidity* Psat)
    pd = Kr*psat;
    cp = ((pd*1000)*(M))/((Ru)*(Tp(i)));
    Kn = (2*64*10^-9)/(dp(i));
    Cm = (1+Kn)/(1+(((4/3)+0.377)*Kn)+((4*Kn^2)/3));
    mpdot=(2*pi)*D*Cm*dp(i)*(cp-cinf);
    ddot = (-4*D*Cm*(cp - cinf))/(den*dp(i));
    
   
    
    Tp(i+1) = (((-L*Cm*D*(cp-cinf))-(kair*(Tp(i)-Tb)))/(den*c*(dp(i)^(2))*(1/12)*(1/dt))) +Tp(i);
    dp(i+1) = dp(i) + dt*ddot;
    
    if (dp(i+1) < 0)
        
        break; 
    end
    
    
    mp(i+1) = ((den*pi/6)*(dp(i+1))^(3));
    tp(i+1) = tp(i) + dt;
     
   
   
    
end
figure(1)
plot (tp,dp)
xlabel('time [s]'), ylabel('diameter [nm]')

figure(2)
plot(tp,mp)
xlabel('time [s]'), ylabel('mass[Kg]')

figure(3)
plot(tp,Tp)
xlabel('time [s]'), ylabel('Temp[C]')

upon plotting the arrays, i don't want that converting to zero line, it's misleading and i do not want to edit the plotting command for this specific case, i'm looking for a generic solution

Answer 1

[Davide Masiello]

Would this suit your purpose?

RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]
Tb = 293.15;
pinf =RHo*2.31; % it's constant because Tb is assumed to be constant
cinf = (pinf*1000)*(M) / ((Ru)*(Tb)); %it's constant because pinf and Tb are constants
dt=0.001; % time step
tend = 3; %total time
num_steps = tend/dt; %number of timesteps
Tp = zeros(1,num_steps); %array to store Temperature of droplet at every timestep
dp = zeros(1,num_steps); %array to store diameter at every timestep
mp = zeros(1,num_steps); %array to store mass at every timestep
tp =zeros (1,num_steps);
Tp(1) = 293;
dp(1) = 40000e-9; %initial diameter value of the droplet
mp(1) = (pi/6) *den*(dp(1))^3; %initial value of the masss of the droplet
tp(1)= 0;
for i =1:num_steps-1
    
    
    Kr = exp ((4*st*M)/(Ru*den*dp(i)*Tp(i)));
    psat= RHp*(exp(16.7-(4060/(Tp(i) -37)))); %it's constant because Tp is constant (Relative humidity* Psat)
    pd = Kr*psat;
    cp = ((pd*1000)*(M))/((Ru)*(Tp(i)));
    Kn = (2*64*10^-9)/(dp(i));
    Cm = (1+Kn)/(1+(((4/3)+0.377)*Kn)+((4*Kn^2)/3));
    mpdot=(2*pi)*D*Cm*dp(i)*(cp-cinf);
    ddot = (-4*D*Cm*(cp - cinf))/(den*dp(i));
    
   
    
    Tp(i+1) = (((-L*Cm*D*(cp-cinf))-(kair*(Tp(i)-Tb)))/(den*c*(dp(i)^(2))*(1/12)*(1/dt))) +Tp(i);
    dp(i+1) = dp(i) + dt*ddot;
    
    if (dp(i+1) < 0)
        
        break; 
    end
    
    
    mp(i+1) = ((den*pi/6)*(dp(i+1))^(3));
    tp(i+1) = tp(i) + dt;
     
   
   
    
end

%% trimming the vectors %%
tp = tp(1:i);
dp = dp(1:i);
mp = mp(1:i);
Tp = Tp(1:i);
%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
plot (tp,dp)
xlabel('time [s]'), ylabel('diameter [nm]')

figure(2)
plot(tp,mp)
xlabel('time [s]'), ylabel('mass[Kg]')

figure(3)
plot(tp,Tp)
xlabel('time [s]'), ylabel('Temp[C]')

[smith]

Hello, I made a little edit on my if statement and applied your approach, look what happened to the temperature curve, thoughts?

clear all
close all

RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]


Tb = 293.15;

pinf =RHo*2.31; % it's constant because Tb is assumed to be constant
cinf = (pinf*1000)*(M) / ((Ru)*(Tb)); %it's constant because pinf and Tb are constants




dt=0.00000000001; % time step
tend = 4.9e-5; %total time
num_steps = tend/dt; %number of timesteps
Tp = zeros(1,num_steps); %array to store Temperature of droplet at every timestep
dp = zeros(1,num_steps); %array to store diameter at every timestep
mp = zeros(1,num_steps); %array to store mass at every timestep
tp =zeros (1,num_steps);
Tp(1) = 293;
dp(1) = 100e-9; %initial diameter value of the droplet
mp(1) = (pi/6) *den*(dp(1))^3; %initial value of the masss of the droplet
tp(1)= 0;

for i =1:num_steps-1
  
    Kr = exp ((4*st*M)/(Ru*den*dp(i)*Tp(i)));
    psat= RHp*(exp(16.7-(4060/(Tp(i) -37)))); %it's constant because Tp is constant (Relative humidity* Psat)
    pd = Kr*psat;
    cp = ((pd*1000)*(M))/((Ru)*(Tp(i)));
    Kn = (2*64*10^-9)/(dp(i));
    Cm = (1+Kn)/(1+(((4/3)+0.377)*Kn)+((4*Kn^2)/3));
    mpdot=(2*pi)*D*Cm*dp(i)*(cp-cinf);
    ddot = (-4*D*Cm*(cp - cinf))/(den*dp(i));
   
    Tp(i+1) = (((-L*Cm*D*(cp-cinf))-(kair*(Tp(i)-Tb)))/(den*c*(dp(i)^(2))*(1/12)*(1/dt))) +Tp(i);
    dp(i+1) = dp(i) + dt*ddot;
    
   if (ddot > 0)
        
        break; 
    end
  
    mp(i+1) = ((den*pi/6)*(dp(i+1))^(3));
    tp(i+1) = tp(i) + dt;
   
end

tp = tp(1:i);
dp = dp(1:i);
mp = mp(1:i);
Tp = Tp(1:i);

figure(1)
plot (tp,dp)
xlabel('time [s]'), ylabel('diameter [nm]')

figure(2)
plot(tp,mp)
xlabel('time [s]'), ylabel('mass[Kg]')

figure(3)
plot(tp,Tp)
xlabel('time [s]'), ylabel('Temp[C]')

[Davide Masiello]

Well, different results are to be expected if you change the condition for breaking the loop. The time span is also much smaller.

[smith]

The new condition serves my objective more accurately, my question is about the shape of the plot of the temperature at the end, do you have any clue why's it behaving this way?

[Davide Masiello]

That could have to do with the equations you have implemented into code, so we would have to dwell into the physical correctness of your mathematics and the correctness of the impekentation itself, but much more information is needed. I suggets you open a different question about it.

Answer 2

[prasanna balaji]

A perfect solution would be to check for negative values in the arrays and either set a minimum value or remove them before plotting.

%Check for negative values in dp and set a minimum value

dp(dp<0) = min(dp(dp>0));

%Check for negative values in mp and set a minimum value

mp(mp<0) = min(mp(mp>0));

%Check for negative values in Tp and set a minimum value

Tp(Tp<0) = min(Tp(Tp>0));

figure 1

plot (tp,dp)

xlabel('time [s]'), ylabel('diameter [nm]')

ylim([min(dp), max(dp)])

figure 2

plot(tp,mp)

xlabel('time [s]'), ylabel('mass[Kg]')

ylim([min(mp), max(mp)])

figure 3

plot(tp,Tp)

xlabel('time [s]'), ylabel('Temp[C]')

ylim([min(Tp), max(Tp)])

[smith]

I changed my if statement, it's no longer about the negative in these arrays, it's about a negative iterated value "ddot"

clear all
close all

RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]


Tb = 293.15;

pinf =RHo*2.31; % it's constant because Tb is assumed to be constant
cinf = (pinf*1000)*(M) / ((Ru)*(Tb)); %it's constant because pinf and Tb are constants




dt=0.00000000001; % time step
tend = 4.9e-5; %total time
num_steps = tend/dt; %number of timesteps
Tp = zeros(1,num_steps); %array to store Temperature of droplet at every timestep
dp = zeros(1,num_steps); %array to store diameter at every timestep
mp = zeros(1,num_steps); %array to store mass at every timestep
tp =zeros (1,num_steps);
Tp(1) = 293;
dp(1) = 100e-9; %initial diameter value of the droplet
mp(1) = (pi/6) *den*(dp(1))^3; %initial value of the masss of the droplet
tp(1)= 0;

for i =1:num_steps-1
  
    Kr = exp ((4*st*M)/(Ru*den*dp(i)*Tp(i)));
    psat= RHp*(exp(16.7-(4060/(Tp(i) -37)))); %it's constant because Tp is constant (Relative humidity* Psat)
    pd = Kr*psat;
    cp = ((pd*1000)*(M))/((Ru)*(Tp(i)));
    Kn = (2*64*10^-9)/(dp(i));
    Cm = (1+Kn)/(1+(((4/3)+0.377)*Kn)+((4*Kn^2)/3));
    mpdot=(2*pi)*D*Cm*dp(i)*(cp-cinf);
    ddot = (-4*D*Cm*(cp - cinf))/(den*dp(i));
   
    Tp(i+1) = (((-L*Cm*D*(cp-cinf))-(kair*(Tp(i)-Tb)))/(den*c*(dp(i)^(2))*(1/12)*(1/dt))) +Tp(i);
    dp(i+1) = dp(i) + dt*ddot;
    
   if (ddot > 0)
        
        break; 
    end
  
    mp(i+1) = ((den*pi/6)*(dp(i+1))^(3));
    tp(i+1) = tp(i) + dt;
   
end



figure(1)
plot (tp,dp)
xlabel('time [s]'), ylabel('diameter [nm]')
figure(2)
plot(tp,mp)
xlabel('time [s]'), ylabel('mass[Kg]')
figure(3)
plot(tp,Tp)
xlabel('time [s]'), ylabel('Temp[C]')

Answer 3

[fi]

You can remove the range of zeros at the end of an array like this:

A(800:end) = []; % Removes all elements at the end of A, starting at index 800

This of course means you'll need to know at which index you zeros start. Since your for loop stopped there, in your case, you could just use the loop index i for that:

for i =1:num_steps-1
    % ...
    Tp(i+1) = % ...
    if %...
        break;
    end
end

% Tp(i+1) was the last element written to, so we need to
% remove all elements starting with Tp(i+2)
Tp(i+2:end) = [];