function air_column()
% Air separation - single-column nitrogen generator (cryogenic), Wang-Henke bubble-point method
% bdrworks.org - Handbook - Air separation, worked model 2
%   N equilibrium stages, saturated air vapour F = 1 below stage N, gaseous N2 product D,
%   constant molar overflow (V = 1, L = 1 - D), Raoult K-values from Lee-Kesler vapour pressures.
%   For fixed temperatures each component balance is a tridiagonal system (Thomas algorithm).
% Save as air_column.m and run:  air_column      (MATLAB or GNU Octave)

P = 8;  z = [0.78118; 0.00934; 0.20957];  z = z/sum(z);
[Dmax, xs] = pinch(P, z);
fprintf('pinch: at most D = %.3f mol N2 per mol air (%.1f %% recovery); bottom liquid %.1f %% O2\n', Dmax, 100*Dmax/z(1), 100*xs(3));
Ns = [20 30 40 50 60];  rec = zeros(size(Ns));  ar = rec;
for k = 1:numel(Ns)
    D = D_for_ppm(Ns(k), P, z, 1.0);
    [T, X, Y] = column(Ns(k), D, P, z);
    rec(k) = 100*D*Y(1,1)/z(1);  ar(k) = 1e6*Y(1,2);
    fprintf('N = %d: 1 ppm O2 at D = %.4f, N2 recovery %.1f %%, Ar in product %.0f ppm, bottom liquid %.1f %% O2\n', ...
            Ns(k), D, rec(k), ar(k), 100*X(end,3));
end
N = 40;  D = D_for_ppm(N, P, z, 1.0);  [T, X, Y] = column(N, D, P, z);
Tc = dewT(Y(1,:)', P);
Pev = fzero(@(p) dewT(X(end,:)', p) - (Tc - 2), [1 P]);
fprintf('design N = 40: N2 condenses at %.2f K; the crude liquid must evaporate at or below %.2f bar\n', Tc, Pev);

figure('Position', [60 100 1300 380]);
subplot(1,3,1);
Ds = linspace(0.30, Dmax - 0.002, 20);
for Nn = [20 30 40 60]
    ppm = zeros(size(Ds));
    for k = 1:numel(Ds), [~, ~, Yk] = column(Nn, Ds(k), P, z); ppm(k) = 1e6*Yk(1,3); end
    semilogy(100*Ds/z(1), ppm, 'LineWidth', 2); hold on;
end
xlabel('N_2 recovery (%)'); ylabel('O_2 in product (ppm)'); title('Purity against recovery, 8 bar');
legend('20 stages', '30 stages', '40 stages', '60 stages'); grid on;
subplot(1,3,2);
semilogx(X(:,3), 1:N, X(:,2), 1:N, 'LineWidth', 2); set(gca, 'YDir', 'reverse');
xlabel('liquid mole fraction'); ylabel('stage (1 = top)'); title('Profiles, 40 stages'); legend('O_2', 'Ar'); grid on;
subplot(1,3,3);
plot(Ns, rec, 'o-', 'LineWidth', 2); xlabel('equilibrium stages'); ylabel('N_2 recovery at 1 ppm O_2 (%)'); grid on;
end

% ---------------------------------------------------------------- column model
function [T, X, Y] = column(N, D, P, z)
L = 1 - D;  T = dewT(z, P)*ones(N, 1);  w = 1;  dTold = Inf;
for it = 1:3000
    K = Kmat(T, P);  X = zeros(N, 3);
    for i = 1:3
        a = L*ones(N, 1);  a(1) = 0;
        b = -(L + K(:,i));  b(1) = -(L + D*K(1,i));      % top: reflux = condensed top vapour
        c = [K(2:end,i); 0];
        d = zeros(N, 1);  d(end) = -z(i);
        X(:,i) = thomas(a, b, c, d);
    end
    X = max(X, 1e-300);  X = X./sum(X, 2);
    Tn = bubT_all(X, P, T);  dT = max(abs(Tn - T));
    if dT < 1e-9, T = Tn; break; end
    if it > 1 && dT > dTold, w = max(0.5*w, 0.05); end   % damp oscillations near the pinch
    T = T + w*(Tn - T);  dTold = dT;
end
Y = Kmat(T, P).*X;
end

function x = thomas(a, b, c, d)
n = numel(d);  cp = zeros(n,1);  dp = zeros(n,1);
cp(1) = c(1)/b(1);  dp(1) = d(1)/b(1);
for i = 2:n
    m = b(i) - a(i)*cp(i-1);  cp(i) = c(i)/m;  dp(i) = (d(i) - a(i)*dp(i-1))/m;
end
x = zeros(n,1);  x(n) = dp(n);
for i = n-1:-1:1, x(i) = dp(i) - cp(i)*x(i+1); end
end

function [Dmax, xs] = pinch(P, z)
Td = dewT(z, P);  xs = z./Kv(Td, P);  xs = xs/sum(xs);
Dmax = 1 - z(3)/xs(3);
end

function D = D_for_ppm(N, P, z, ppm)
f = @(D) log(1e6*top_O2(N, D, P, z)/ppm);
Dhi = pinch(P, z) - 0.002;
if f(Dhi) < 0, D = Dhi; else, D = fzero(f, [0.20 Dhi]); end
end

function y = top_O2(N, D, P, z)
[~, ~, Y] = column(N, D, P, z);  y = Y(1,3);
end

% ---------------------------------------------------------------- properties (Lee-Kesler, Raoult)
function c = crit(),  c = [126.2 34.00 77.35; 150.9 48.98 87.3; 154.6 50.43 90.2];  end   % N2, Ar, O2
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 g = dlnpsat(i, T)
c = crit();  Tr = T/c(i,1);
g = ((6.09648./Tr.^2 - 1.28862./Tr + 6*0.169347*Tr.^5) + wlk(i)*(15.6875./Tr.^2 - 13.4721./Tr + 6*0.43577*Tr.^5))/c(i,1);
end
function k = Kv(T, P),  k = [psat(1, T); psat(2, T); psat(3, T)]/P;  end
function K = Kmat(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_all(X, P, T)
% bubble points of all stages at once, Newton's method on ln(sum K x) = 0
for it = 1:60
    K = Kmat(T, P);  s = sum(K.*X, 2);
    g = sum(K.*X.*[dlnpsat(1, T), dlnpsat(2, T), dlnpsat(3, T)], 2)./s;
    dT = max(min(-log(s)./g, 10), -10);  T = T + dT;
    if max(abs(dT)) < 1e-11, break; end
end
end