VanDerWaal.m

%Van der Waals
clear;
clc;

Tvdw = linspace(0,3.5,150);
V = linspace(0.25,5,202);

[TT,VV] = meshgrid(Tvdw,V);

P = 8/3.*(TT./(VV-1/3))-3./(VV.^2);
Pc = 8/3.*(1./(V-1/3))-3./(V.^2);


[Cplot,h1] = contourf(VV,P,TT,50);
##clabel(Cplot,h1);
colormap("white");
axis([0,5,0,1.5]);
grid minor;
xlabel("Reduced Volume V/V_c");
ylabel("Reduced Pressure P/P_c");

hold on

%the following equations are from %[1]

dd = linspace(0.001,25,200);
T = [length(dd)];
vg = [length(dd)];
vl = [length(dd)];
p = [length(dd)];

for ii = 1:length(dd)

d = dd(ii);
x = dd(ii);

vg(ii) = -1/6*(4*x*exp(2*x) - exp(4*x) + 1)*exp(x)/(x*exp(3*x) +...
x*exp(x) - exp(3*x) + exp(x)) + 1/3;

vl(ii) =-1/6*(4*x*exp(2*x) - exp(4*x) + 1)*exp(-x)/(x*exp(3*x) + x*exp(x) -
exp(3*x) + exp(x)) + 1/3;

T(ii) = -27/4*((4*d*exp(2*d) - exp(4*d) + 1)*((4*d*exp(2*d) - exp(4*d) + 1)*exp(-d)...
/(d*exp(3*d) + d*exp(d) - exp(3*d) + exp(d)) - 2)^2*exp(d)/(d*exp(3*d) + d*exp(d) -...
exp(3*d) + exp(d)) + (((4*d*exp(2*d) - exp(4*d) + 1)*exp(d)/(d*exp(3*d) + d*exp(d) -...
exp(3*d) + exp(d)) - 2)^2 + 4*(4*d*exp(2*d) - exp(4*d) + 1)*exp(d)/(d*exp(3*d) +...
d*exp(d) - exp(3*d) + exp(d)) - 4)*((4*d*exp(2*d) - exp(4*d) + 1)*exp(-d)/...
(d*exp(3*d) + d*exp(d) - exp(3*d) +...
exp(d)) - 2) + 2*((4*d*exp(2*d) - exp(4*d) + 1)*exp(d)/...
(d*exp(3*d) +...
d*exp(d) - exp(3*d) + exp(d)) - 2)^2 +...
4*(4*d*exp(2*d) - exp(4*d) + 1)*exp(d)/(d*exp(3*d) + d*exp(d) - exp(3*d) +...
exp(d)) - 8)/(((4*d*exp(2*d) - exp(4*d) +...
1)*exp(-d)/(d*exp(3*d) + d*exp(d) - exp(3*d) + exp(d)) - 2)^2*((4*d*exp(2*d) -...
exp(4*d) + 1)*exp(d)/(d*exp(3*d) + d*exp(d) - exp(3*d) + exp(d)) - 2)^2);

p(ii) = 8*T(ii)/(3*vg(ii)-1)-3/vg(ii)^2;
endfor

#you will have to adjust the "33" if you change the resolution of V
area(V(33:end),Pc(33:end),'facecolor','red','facealpha',0.2);
area(V(1:33),Pc(1:33),'facecolor','blue','facealpha',0.2);
area(vg,p,'facecolor','magenta','facealpha',0.75);
area(vl,p,'facecolor','magenta','facealpha',0.75);

h2 = plot(vg,p,'r-+');
h3 = plot(vl,p,'b-x');
h4 = plot(V,Pc,'k-*');

vg50 = (vg+vl)/2;
vg75 = (3*vg+vl)/4;
vg25 = (vg+3*vl)/4;
vg10 = (vg+9*vl)/10;
vg90 = (vg*9+vl)/10;

plot(vg25,p,'k:.');
plot(vg50,p,'k--');
plot(vg75,p,'k:.');
plot(vg10,p,'k:');
plot(vg90,p,'k:');

hold off
legend([h2,h3,h4],'Vapor Locus','Liquid Locus','Supercritical Locus');

%test properties
T_test = 0.5
V_test = 0.5

Tcrit = 154.6; %K
Pcrit = 5043; %kPa

vg_test = interp1(T,vg,T_test,"cubic")
vl_test = interp1(T,vl,T_test,"cubic")

p_test = 8*T_test/(3*vg_test-1)-3/vg_test^2

gas = (V_test-vl_test)/(vg_test-vl_test);
liquid = 1 - gas;

disp(['Gas = ',num2str(gas*100),'%']);
disp(['Liquid = ',num2str(liquid*100),'%']);
disp(['Temperture = ',num2str(T_test*Tcrit),'K']);
disp(['Pressure = ',num2str(p_test*Pcrit),'kPa']);


%References

%[1] https://opencommons.uconn.edu/cgi/viewcontent.cgi?article=1107&context=chem_educ