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

1 view (last 30 days)
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)')

Accepted Answer

Alan Stevens
Alan Stevens on 8 Aug 2020
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)')
  1 Comment
asd ad
asd ad on 8 Aug 2020
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

Sign in to comment.

More Answers (1)

Alan Stevens
Alan Stevens on 8 Aug 2020
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.
  1 Comment
asd ad
asd ad on 8 Aug 2020
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

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!