Plotting a contour plot for a variable dependednt on the while loop

Hello everyone,

I'm currently running a code to see how long does it take for a range of volume to reach 0 for different relative humidities and I'm trying to plot a contour plot for the values with Volume on the y-axis and Relative humidity on the x axis and the time taken on the right side. For some reason, my volume doesn't seem to converge and the code won't stop running. Does anyone know why?

Thanks in advance

clear all 
close all
clc

%% Inputs

IR = 1e-3; %initial radius [m]
Rho = 1000; %density [kg/m^3]
T = 23.5; %temperature [celsius]

% Loop 1
aInitial = 1;
aStep = 1;
aMax = 25;
VMin = 1e-12; %minimum volume [m^3]
VMax = 2e-9; %maximum volume [m^3]

% Loop 2
bInitial = 1;
bStep = 1;
bMax = 25;
ReHuMin = 0.0; %minimum relative humidity [m]
ReHuMax = 0.9; %maximum relative humidity [m]

% Loop 3
c = 1;
tInitial = 0; %initial time [s]
tStep = 1; %final time [s]



D_T = 2.5e-4*exp(-684.15/(T+273.15)); %coefficient [m^2/s]
c_sat = (9.99e-7)*T^3 - (6.94e-5)*T^2 + (3.2e-3)*T - 2.87e-2; %concentration [kg/m^3]

%% Computing

for a = aInitial:aStep:aMax
    for b = bInitial:bStep:bMax
        Vol(a,b) = VMin + (VMax - VMin)*(a-aInitial)/(aMax);
        h(a,b) = ((sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(2/3) - pi^(2/3)*IR^2)/(pi^(1/3)*(sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(1/3));
        CAR(a) =  2*atan(h(a,b)/IR);
        CAD(a) = CAR(a)*180/pi;
        
        RH(a,b) = ReHuMin + (ReHuMax - ReHuMin)*(b - bInitial)/(bMax);
        
        h(c,a,b) = h(a,b);
        CAR(c,a,b) = CAR(a);
        Vol(c,a,b) = Vol(a,b);
        timeSec(c,a,b) = tInitial;
        
        while Vol(c,a,b) > 0
            
            M_dot(c,a,b) = -pi*IR*D_T*(1 - RH(a,b))*c_sat*(0.27*CAR(c,a,b)^2+1.30); %mass flow
            
            Mkg(c,a,b) = M_dot(c,a,b)*tStep; %mass loss at each time step [kg]
            
            Vm3(c,a,b) = Mkg(c,a,b)/Rho; %volume loss [m^3]
            
            Vol(c+1,a,b) = Vol(c,a,b) + Vm3(c,a,b); %new volume [m^3]
            
            h(c+1,a,b) = ((sqrt(pi^2*IR^6 + 9*(Vol(c+1,a,b))^2) + 3*(Vol(c+1,a,b)))^(2/3) - pi^(2/3)*IR^2)/(pi^(1/3)*(sqrt(pi^2*IR^6 + 9*(Vol(c+1,a,b))^2) + 3*(Vol(c+1,a,b)))^(1/3)); %new height [m]
            
            CAR(c+1,a,b) =  2*atan(h(c+1,a,b)/IR); %new angle [radians]
            
            CAD(c+1,a,b) = CAR(c+1,a,b)*180/pi; %new angle [degrees]
            
            timeSec(c+1,a,b) = tInitial + tStep*(c - 1);
            
            c = c + 1;
            
        end
        
        EvaporationTime(a,b) = timeSec(c,a,b);
        
        c = 1;
        
    end
end

%% Plotting

contourf(RH*100,Vol*1e9,EvaporationTime)
c = colorbar;
c.Label.String = 'Evaporation Time (s)';
xlabel('Relative Humidity (%)')
ylabel('Volume (mm^3)')
            

Answer 1

[Alan Stevens]

I'm not sure if I've interpreted your requirements correctly, but what about the following:

%% Inputs

IR = 1e-3; %initial radius [m]
Rho = 1000; %density [kg/m^3]
T = 23.5; %temperature [celsius]

% Loop 1
aInitial = 1;
aStep = 1;
aMax = 25;
VMin = 1e-12; %minimum volume [m^3]
VMax = 2e-9; %maximum volume [m^3]

% Loop 2
bInitial = 1;
bStep = 1;
bMax = 25;
ReHuMin = 0.0; %minimum relative humidity [m]
ReHuMax = 0.9; %maximum relative humidity [m]

% Loop 3
c = 1;
tInitial = 0; %initial time [s]
tStep = 1; %final time [s]



D_T = 2.5e-4*exp(-684.15/(T+273.15)); %coefficient [m^2/s]
c_sat = (9.99e-7)*T^3 - (6.94e-5)*T^2 + (3.2e-3)*T - 2.87e-2; %concentration [kg/m^3]

%% Computing

for a = aInitial:aStep:aMax
    Va = VMin + (VMax - VMin)*(a-aInitial)/(aMax);
    for b = bInitial:bStep:bMax
        V(a,b) = Va;  %Initial volume
        Vol(a,b) = Va;
        h(a,b) = ((sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(2/3) - pi^(2/3)*IR^2)/(pi^(1/3)*(sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(1/3));
        CAR(a,b) =  2*atan(h(a,b)/IR);
        CAD(a,b) = CAR(a,b)*180/pi;
        
        RH(a,b) = ReHuMin + (ReHuMax - ReHuMin)*(b - bInitial)/(bMax);
        
        timeSec(a,b) = tInitial;
        
         Vold = 0; c = 0;
         
        while Vol(a,b)>VMin && c<10000
            Vold = Vol(a,b);
            M_dot(a,b) = -pi*IR*D_T*(1 - RH(a,b))*c_sat*(0.27*CAR(a,b)^2+1.30); %mass flow
            
            Mkg(a,b) = M_dot(a,b)*tStep; %mass loss at each time step [kg]
            
            Vm3(a,b) = Mkg(a,b)/Rho; %volume loss [m^3]
            
            Vol(a,b) = Vol(a,b) + Vm3(a,b); %new volume [m^3]
            Vol(a,b) = max(Vol(a,b), VMin);
            
            h(a,b) = ((sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(2/3) - pi^(2/3)*IR^2)/(pi^(1/3)*(sqrt(pi^2*IR^6 + 9*(Vol(a,b))^2) + 3*(Vol(a,b)))^(1/3)); %new height [m]
            
            CAR(a,b) =  2*atan(h(a,b)/IR); %new angle [radians]
            
            CAD(a,b) = CAR(a,b)*180/pi; %new angle [degrees]
            
            timeSec(a,b) = tStep*c;
            
          
            c = c + 1;
        end
        
        EvaporationTime(a,b) = timeSec(a,b);
        
       
        
    end
end

%% Plotting


contourf(RH*100,V*1e9,EvaporationTime)  % plot against Initial volume not final volume
k = colorbar;
k.Label.String = 'Evaporation Time (s)';
xlabel('Relative Humidity (%)')
ylabel('Initial Volume (mm^3)')
figure
surf(RH*100,V*1e9,EvaporationTime)
k = colorbar;
k.Label.String = 'Evaporation Time (s)';
xlabel('Relative Humidity (%)')
ylabel('Initial Volume (mm^3)')
zlabel('Evaporation Time (s)')

[asd ad]

Thanks a lot Alan Stevens. It works! You're a lifesaver. I've been at it for the past 2 days and didn't know where I'm going wrong

Answer 2

[Alan Stevens]

Your code doesn't stop because Vol never reaches zero.

What does the c index do for you?

Your contour function won't plot while you have RH and EvaporationTime as 2d matrices with Vol as a 3d matrix.

[asd ad]

The cindex is just for the time step of the while loop. I might have coded it wrong as well. Is there a better way to do it?

Thanks