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

Question

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]

clear all 
close all
clc

%% Inputs

IR = 1e-3; %initial wetted 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]

%% Initialising

D_T = 2.5e-4*exp(-684.15/(T+273.15)); %diffusion coefficient [m^2/s]
c_sat = (9.99e-7)*T^3 - (6.94e-5)*T^2 + (3.2e-3)*T - 2.87e-2; %saturation 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 contact angle [radians]
            
            CAD(c+1,a,b) = CAR(c+1,a,b)*180/pi; %new contact degrees [radians]
            
            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 (m^3)')
            

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)')

Alan Stevens · Accepted Answer

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

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

采纳的回答

1 个评论 显示 -1更早的评论 隐藏 -1更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论 隐藏 -1更早的评论

类别

标签

另请参阅

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论