function air_membrane()
% Air separation - membrane nitrogen generator (hollow fibres) and the energy of the three routes
% bdrworks.org - Handbook - Air separation, worked model 3
%   binary O2/N2, local permeate from phi(1-a) y^2 + [1 + (a-1)(x+phi)] y - a x = 0
%   complete mixing and cross flow (dn/dx = n/(y - x)); energy w = R T0 ln(pf/p0)/eta per mol air
% Save as air_membrane.m and run:  air_membrane      (MATLAB or GNU Octave)

xf = 0.2095;  PF = 12.35;  PP = 1.013;  phi = PP/PF;  R = 8.314462618;  T0 = 298.15;
tab = [0.005 21.7 0.73; 0.01 27.3 0.58; 0.02 35.0 0.45; 0.04 46.1 0.34; 0.06 55.2 0.29; 0.08 63.6 0.26];
err = @(a) sum((arrayfun(@(k) 100*recov(crossflow(a, phi, tab(k,1), xf), tab(k,1), xf), 1:6)' - tab(:,2)).^2);
afit = fminbnd(err, 2.5, 8);
fprintf('cross flow fitted to Table 7: alpha = %.2f (rms %.1f points)\n', afit, sqrt(err(afit)/6));
for k = 1:6
    x = tab(k,1);
    fprintf('%4.1f %% O2: Table 7 %5.1f, cross(fit) %5.1f, cross(alpha 6) %5.1f, complete(alpha 6) %6.1f\n', 100*x, tab(k,2), ...
            100*recov(crossflow(afit, phi, x, xf), x, xf), 100*recov(crossflow(6, phi, x, xf), x, xf), 100*recov(complete(6, phi, x, xf), x, xf));
end
fprintf('complete mixing, alpha 6, phi -> 0: the retentate cannot fall below %.2f %% O2\n', 100*xf/(6 - 5*xf));
wair = R*T0*log(PF/PP);  eta = zeros(6,1);
for k = 1:6
    eta(k) = kwh(wair*(1 - tab(k,1))/(tab(k,2)/100*(1 - xf)))/tab(k,3);
end
fprintf('Table 7 energies imply an isothermal compression efficiency of %.2f-%.2f\n', min(eta), max(eta));
wcol = R*T0*log(8.5/1.013)/(0.640*(1 - xf));
fprintf('single column, 64 %% recovery: isothermal work %.3f kWh/Nm3 -> Table 5 (0.25) implies eta = %.2f\n', kwh(wcol), kwh(wcol)/0.25);
fprintf('minimum work, 1 ppm N2 at 8 bar, 64 %% recovery: %.4f kWh/Nm3\n', kwh(wmin(1e-6, 8, 0.64, xf)));
fprintf('minimum work, 1 %% O2 at 11.35 bar, 27 %% recovery: %.4f kWh/Nm3\n', kwh(wmin(0.01, 11.35, 0.273, xf)));
prod = 1000/3600/0.022414;
A = area(6, 40, PF, PP, 0.01, xf)*prod/crossflow(6, phi, 0.01, xf);
fprintf('alpha 6, 40 GPU fibres: %.0f m2 for 1000 Nm3/h at 1 %% O2\n', A);

xs = linspace(0.004, 0.10, 40);
figure('Position', [60 100 1300 380]);
subplot(1,3,1);
plot(100*xs, arrayfun(@(x) 100*recov(crossflow(6, phi, x, xf), x, xf), xs), ...
     100*xs, arrayfun(@(x) 100*recov(crossflow(afit, phi, x, xf), x, xf), xs), 'LineWidth', 2); hold on;
plot(100*tab(:,1), tab(:,2), 'ks');
xlabel('O_2 in product (%)'); ylabel('N_2 recovery (%)'); title('Membrane: purity against recovery');
legend('cross flow, \alpha = 6', 'cross flow, fitted \alpha', 'Table 7'); grid on;
subplot(1,3,2);
loglog(100*tab(:,1), tab(:,3), 'ks', [3 0.1], [0.2 0.5]*1.2505, 'd', [1e-4 1e-4], [0.25 0.22], 'o'); hold on;
yps = logspace(-6, -0.9, 30);
loglog(100*yps, arrayfun(@(y) kwh(wmin(y, 8, 1e-6, xf)), yps), 'k--');
xlabel('O_2 in product (%)'); ylabel('kWh per Nm^3 N_2'); title('Energy of the three routes');
legend('membrane (Table 7)', 'PSA', 'cryogenic (Table 5)', 'minimum work'); grid on;
subplot(1,3,3);
xa = linspace(0.005, 0.08, 20);
semilogy(100*xa, arrayfun(@(x) area(6, 40, PF, PP, x, xf)*prod/crossflow(6, phi, x, xf), xa), 'LineWidth', 2);
xlabel('O_2 in product (%)'); ylabel('membrane area (m^2)'); title('Area for 1000 Nm^3/h'); grid on;
end

function y = ylocal(x, a, phi)
A = phi*(1 - a);  B = 1 + (a - 1)*(x + phi);  C = -a*x;
y = -2*C/(B + sqrt(B^2 - 4*A*C));
end

function nr = crossflow(a, phi, xr, xf)
opts = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);
[~, n] = ode45(@(x, n) n/(ylocal(x, a, phi) - x), [xf xr], 1, opts);
nr = n(end);
end

function nr = complete(a, phi, xr, xf)
y = ylocal(xr, a, phi);  th = (xf - xr)/(y - xr);
if th >= 0 && th <= 1, nr = 1 - th; else, nr = NaN; end
end

function r = recov(nr, xr, xf),  r = nr*(1 - xr)/(1 - xf);  end
function e = kwh(J),  e = J/0.022414/3.6e6;  end

function A = area(a, QO2gpu, pf, pp, xr, xf)
% membrane area per mol/s of feed air, cross flow
GPU = 3.3464e-10;  QO2 = QO2gpu*GPU;  QN2 = QO2/a;  phi = pp/pf;
f = @(x, u) [u(1)/(ylocal(x, a, phi) - x); ...
             -(u(1)/(ylocal(x, a, phi) - x))/(QO2*(pf*1e5*x - pp*1e5*ylocal(x, a, phi)) + QN2*(pf*1e5*(1 - x) - pp*1e5*(1 - ylocal(x, a, phi))))];
opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-12);
[~, u] = ode45(f, [xf xr], [1; 0], opts);
A = u(end, 2);
end

function w = wmin(yp, Pp, Rn, xf)
% ideal work, J per mol of product (O2 fraction yp at Pp bar), N2 recovery Rn, reject at 1 bar, T0 = 298.15 K
RT = 8.314462618*298.15;  np = Rn*(1 - xf)/(1 - yp);  nw = 1 - np;  yw = (xf - np*yp)/nw;
s = @(y, P) sum(y(y > 0).*log(y(y > 0)*P));
w = RT*(np*s([yp 1-yp], Pp) + nw*s([yw 1-yw], 1) - s([xf 1-xf], 1))/np;
end