function air_equilibrium()
% Air separation - phase equilibrium of air and the pressures of a double column
% bdrworks.org - Handbook - Air separation, worked model 1
%   Lee-Kesler vapour pressures (omega fitted to the normal boiling point), Raoult's law,
%   dew and bubble points of N2-Ar-O2 air, Clapeyron latent heats with Peng-Robinson volumes
% Save as air_equilibrium.m and run:  air_equilibrium      (MATLAB or GNU Octave)

air = [0.78118; 0.00934; 0.20957];  air = air/sum(air);         % N2, Ar, O2
fprintf('Lee-Kesler omega: N2 %.4f, Ar %.4f, O2 %.4f\n', wlk(1), wlk(2), wlk(3));
Tm = [63.15 80 100 120];  Pm = [0.1246 1.367 7.777 25.13];
for k = 1:4
    fprintf('  N2 %6.2f K: %7.4f bar vs %.4f bar (%+.2f %%)\n', Tm(k), psat(1, Tm(k)), Pm(k), 100*(psat(1, Tm(k))/Pm(k) - 1));
end
Hm = [2218+3248, 2466+2033, 2103+504.6];
for k = 1:3
    fprintf('  latent heat %d K: %.0f J/mol vs %.0f J/mol\n', Tm(k+1), latent(1, Tm(k+1)), Hm(k));
end
for P = [1 6 8]
    fprintf('air at %d bar: dew point %.2f K, bubble point %.2f K\n', P, dewT(air, P), bubT(air, P));
end
PL = 1.3;  TO2 = tsat(3, PL);
for dT = [1 1.5 2 2.5]
    fprintf('LP %.1f bar (O2 boils at %.2f K), dT = %.1f K -> HP column at %.2f bar\n', PL, TO2, dT, psat(1, TO2 + dT));
end
fprintf('latent heats in the main condenser: N2 at HP %.0f J/mol, O2 at LP %.0f J/mol\n', latent(1, TO2 + 2), latent(3, TO2));

Ps = linspace(1, 10, 40);
figure('Position', [60 100 1300 380]);
subplot(1,3,1);
plot(Ps, arrayfun(@(p) dewT(air, p), Ps), Ps, arrayfun(@(p) bubT(air, p), Ps), 'LineWidth', 2); hold on;
plot([1 6], [81.5 101], 'ks');
xlabel('pressure (bar)'); ylabel('temperature (K)'); title('Liquefying air'); legend('dew point', 'bubble point', 'quoted values'); grid on;
subplot(1,3,2);
aNO = arrayfun(@(p) psat(1, bubT(air, p))/psat(3, bubT(air, p)), Ps);
aAO = arrayfun(@(p) psat(2, bubT(air, p))/psat(3, bubT(air, p)), Ps);
plot(Ps, aNO, Ps, aAO, 'LineWidth', 2);
xlabel('pressure (bar)'); ylabel('relative volatility'); title('Volatilities fall with pressure'); legend('N_2/O_2', 'Ar/O_2'); grid on;
subplot(1,3,3);
Ts = linspace(70, 120, 100);
semilogy(Ts, arrayfun(@(t) psat(1, t), Ts), Ts, arrayfun(@(t) psat(3, t), Ts), 'LineWidth', 2); hold on;
semilogy([TO2 TO2+2], [PL psat(1, TO2+2)], 'ko-');
xlabel('temperature (K)'); ylabel('vapour pressure (bar)'); title('Coupling of the two columns'); legend('N_2', 'O_2'); grid on;
end

% ---------------------------------------------------------------- property model
function c = crit()
% Tc (K), Pc (bar), Tn (K), omega (Peng-Robinson): N2, Ar, O2
c = [126.2 34.00 77.35 0.038; 150.9 48.98 87.3 0.000; 154.6 50.43 90.2 0.022];
end

function v = lnPr0(Tr),  v = 5.92714 - 6.09648./Tr - 1.28862*log(Tr) + 0.169347*Tr.^6;  end
function v = lnPr1(Tr),  v = 15.2518 - 15.6875./Tr - 13.4721*log(Tr) + 0.43577*Tr.^6;  end

function w = wlk(i)
c = crit();  Trn = c(i,3)/c(i,1);
w = (log(1.01325/c(i,2)) - lnPr0(Trn))/lnPr1(Trn);
end

function p = psat(i, T)
c = crit();  Tr = T/c(i,1);
p = c(i,2)*exp(lnPr0(Tr) + wlk(i)*lnPr1(Tr));
end

function T = tsat(i, P)
c = crit();  T = fzero(@(t) psat(i, t) - P, [50 c(i,1) - 1e-6]);
end

function k = Kv(T, P),  k = [psat(1, T); psat(2, T); psat(3, T)]/P;  end
function T = dewT(y, P),  T = fzero(@(t) sum(y./Kv(t, P)) - 1, [55 150]);  end
function T = bubT(x, P),  T = fzero(@(t) sum(Kv(t, P).*x) - 1, [55 150]);  end

function L = latent(i, T)
% Clapeyron with Peng-Robinson vapour and liquid compressibility factors, J/mol
R = 8.314462618;  c = crit();  P = psat(i, T);
kap = 0.37464 + 1.54226*c(i,4) - 0.26992*c(i,4)^2;
a = 0.45724*(R*c(i,1))^2/(c(i,2)*1e5)*(1 + kap*(1 - sqrt(T/c(i,1))))^2;  b = 0.07780*R*c(i,1)/(c(i,2)*1e5);
A = a*P*1e5/(R*T)^2;  B = b*P*1e5/(R*T);
z = roots([1, -(1 - B), A - 3*B^2 - 2*B, -(A*B - B^2 - B^3)]);  z = sort(real(z(abs(imag(z)) < 1e-10)));
h = 1e-3;  dlnP = (log(psat(i, T + h)) - log(psat(i, T - h)))/(2*h);
L = R*T^2*(z(end) - z(1))*dlnP;
end