Resources Numerical accuracy Calculation policy
Sizes a column by shortcut methods and solves stage-by-stage columns, absorbers and strippers for flows, temperatures, compositions and duties.
The shortcut formulas are single 64-bit expressions on K-values from the property package; two of them are empirical fits. The rigorous columns are iterations: tridiagonal component balances by the Thomas algorithm, bubble points and energy balances, repeated until temperatures and vapor flows stop changing. Each result reports its component and energy closure so the reader can see how far the iteration got.
Open item. Evaluate Kremser's equation near A = 1 from its series or with expm1 and log1p so it does not divide two cancelled differences.
18 calculations in 2 source files. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests tools/process-test/test.mjs
The least number of stages at total reflux for a split of two key components.
N_min = ln((d_LK ÷ d_HK)(b_HK ÷ b_LK)) ÷ ln α
Two logarithms and a division. When the relative volatility is close to 1 the denominator is close to zero and the result is very sensitive to α.
export function fenske(alpha, dLK, dHK, bLK, bHK) { return Math.log((dLK / dHK) * (bHK / bLK)) / Math.log(alpha); }The stage ratio Y for a reflux ratio X above the minimum.
Y = 0.75 (1 − X^0.5668)
An empirical fit with two typed constants; the fit sets the accuracy.
export function eduljee(X) { return 0.75 * (1 - Math.pow(X, 0.5668)); }The same stage ratio by the fit of Molokanov and others.
Y = 1 − exp((1 + 54.4X) ÷ (11 + 117.2X) × (X − 1) ÷ √X)
An empirical fit. X at or below 0 gives 1 and X at or above 1 gives 0, so the division by √X is never by zero.
export function molokanov(X) { return X <= 0 ? 1 : X >= 1 ? 0 : 1 - Math.exp((1 + 54.4 * X) / (11 + 117.2 * X) * (X - 1) / Math.sqrt(X)); }The stages needed at a reflux ratio, from the minimum stages and minimum reflux.
X = (R − R_min) ÷ (R + 1); N = (N_min + Y) ÷ (1 − Y)
One expression on the correlation's Y. As Y approaches 1 (reflux at its minimum) the denominator approaches zero.
/** Stages N from Nmin, R, Rmin by the chosen correlation. */
export function gilliland(Nmin, R, Rmin, kind = 'eduljee') {
const X = (R - Rmin) / (R + 1), Y = kind === 'molokanov' ? molokanov(X) : eduljee(X); return (Nmin + Y) / (1 - Y);
}The Underwood root between the key volatilities and the minimum reflux ratio from it.
Σ αᵢ z_F,ᵢ ÷ (αᵢ − θ) = 1 − q; R_min = Σ αᵢ x_D,ᵢ ÷ (αᵢ − θ) − 1
Bisection between the two key volatilities, started a relative 1e-9 inside each pole, until the bracket is below an absolute 1e-15 or 300 halvings. The sum has poles at the bracket ends, so terms are large and of opposite sign near them. No sign change in the bracket returns null.
/** Underwood: alphas (relative to the heavy key = 1), feed mole fractions zF, thermal condition q, distillate fractions xD → { theta, Rmin }.
* The root searched is the one between alpha_HK and alpha_LK (two adjacent keys; light components only in D, heavy only in B). */
export function underwood(alpha, zF, q, xD, iLK, iHK) {
const f = (th) => { let s = 0; for (let i = 0; i < alpha.length; i++) s += alpha[i] * zF[i] / (alpha[i] - th); return s - (1 - q); };
let lo = alpha[iHK] + 1e-9 * Math.abs(alpha[iHK] || 1), hi = alpha[iLK] - 1e-9 * Math.abs(alpha[iLK]);
if (!(f(lo) * f(hi) < 0)) return null;
for (let i = 0; i < 300; i++) { const m = 0.5 * (lo + hi); if (f(lo) * f(m) <= 0) hi = m; else lo = m; if (hi - lo < 1e-15) break; }
const theta = 0.5 * (lo + hi); let s = 0; for (let i = 0; i < alpha.length; i++) s += alpha[i] * xD[i] / (alpha[i] - theta);
return { theta, Rmin: s - 1 };
}The fraction of a component absorbed or stripped in N stages.
(A^(N+1) − A) ÷ (A^(N+1) − 1); N ÷ (N + 1) when |A − 1| < 1e-12
Numerator and denominator each subtract nearly equal numbers when A is close to 1. The exact limit is used only within 1e-12 of 1, so just outside that window the quotient is formed from two differences that have lost most of their digits.
export function kremser(A, N, mode = 'absorb') { // absorption factor A = L/(K V) (stripping factor S = K V/L)
if (Math.abs(A - 1) < 1e-12) return N / (N + 1);
return (Math.pow(A, N + 1) - A) / (Math.pow(A, N + 1) - 1);
}q of a feed: the fraction that joins the liquid below the feed stage.
q = (H_dew − H_feed) ÷ (H_dew − H_bubble)
A ratio of enthalpy differences from two saturation flashes. If either flash fails q is taken as 1.
/** Thermal condition of the feed q = (H_sat vapour - H_feed)/(H_sat vapour - H_sat liquid) at the feed pressure. */
export function feedQ(th, s) {
const v = th.flashPVF(s.P, 1, s.z), l = th.flashPVF(s.P, 0, s.z); if (!v || !l) return 1;
return (v.H - s.H) / (v.H - l.H);
}How every component divides between distillate and bottoms, and the volatilities used for it.
dᵢ ÷ bᵢ = (d_HK ÷ b_HK) αᵢ^N_min; αᵢ = √((Kᵢ ÷ K_HK)_top (Kᵢ ÷ K_HK)_bottom)
Three fixed passes: the volatilities are re-evaluated at the dew point of the distillate and the bubble point of the bottoms and averaged geometrically. There is no convergence test on the passes.
for (let pass = 0; pass < 3; pass++) {
Nmin = 0; const dHK = f[iHK] * (1 - p.recHK), bHK = f[iHK] * p.recHK, dLK = f[iLK] * p.recLK, bLK = f[iLK] * (1 - p.recLK);
Nmin = fenske(alpha[iLK], dLK, dHK, bLK, bHK);
d = new Array(n); b = new Array(n);
for (let i = 0; i < n; i++) {
if (i === iLK) { d[i] = dLK; b[i] = bLK; } else if (i === iHK) { d[i] = dHK; b[i] = bHK; }
else { const r = (dHK / bHK) * Math.pow(alpha[i], Nmin); d[i] = f[i] * r / (1 + r); b[i] = f[i] - d[i]; }
}
const D = sum(d), B = sum(b), xD = d.map((v) => v / D), xB = b.map((v) => v / B);
const dew = th.dewT(xD, P, bub.T), bot = th.bubbleT(xB, P, bub.T);
if (dew && bot) {
const Kt = Kat(dew.T, dew.x, xD), Kb = Kat(bot.T, xB, bot.y);
alpha = Kt.map((k, i) => Math.sqrt((k / Kt[iHK]) * (Kb[i] / Kb[iHK])));The split of stages above and below the feed (Kirkbride) and the condenser and reboiler duties.
N_R ÷ N_S = ((z_HK ÷ z_LK)(B ÷ D)(x_B,LK ÷ x_D,HK)²)^0.206; Q_c = (R + 1) D (H_liquid − H_vapor); Q_r from the overall enthalpy balance
Kirkbride's exponent 0.206 is an empirical constant. The reboiler duty is what closes the overall balance, a sum and difference of enthalpy flows that can be much larger than the duty itself.
const NR = (zF[iHK] / zF[iLK]) * (B / D) * Math.pow(xB[iLK] / xD[iHK], 2); const ratio = Math.pow(NR, 0.206);
const Nfeed = N != null ? { rectifying: N * ratio / (1 + ratio), stripping: N / (1 + ratio) } : null;
const dist = th.flashPVF(P, 0, xD), vapD = th.flashPVF(P, 1, xD), bot = th.flashPVF(P, 0, xB);
const outD = fromFlash(D, dist), outB = fromFlash(B, bot);
const Qc = R != null ? (R + 1) * D * (dist.H - vapD.H) : null, Qr = Qc != null ? D * dist.H + B * bot.H - s.F * s.H - Qc : null;The liquid flow of each component on each stage of the rigorous column.
l_{j−1} − (1 + S_j) l_j + S_{j+1} l_{j+1} = −f_j, S_j = K_j V_j ÷ L_j; solved by the Thomas algorithm
Forward elimination and back substitution without pivoting. Stripping factors can span many orders of magnitude between light and heavy components. A negative flow from the solve is set to zero.
for (let i = 0; i < n; i++) { // tridiagonal for component i (Thomas)
const a = new Array(N).fill(0), b = new Array(N), c = new Array(N).fill(0), r = new Array(N);
for (let j = 0; j < N; j++) {
const S = K[j][i] * V[j] / (j === N - 1 ? Math.max(Ftot - D, 1e-30) : L[j]);
b[j] = j === 0 ? -(1 + S / (R + 1)) : -(1 + S); r[j] = -f[j][i];
if (j > 0) a[j] = 1;
if (j < N - 1) c[j] = K[j + 1][i] * V[j + 1] / (j + 1 === N - 1 ? Math.max(Ftot - D, 1e-30) : L[j + 1]);
}
for (let j = 1; j < N; j++) { const m = a[j] / b[j - 1]; b[j] -= m * c[j - 1]; r[j] -= m * r[j - 1]; }
const sol = new Array(N); sol[N - 1] = r[N - 1] / b[N - 1]; for (let j = N - 2; j >= 0; j--) sol[j] = (r[j] - c[j] * sol[j + 1]) / b[j];
for (let j = 0; j < N; j++) lij[i][j] = Math.max(sol[j], 0);
}The vapor flow leaving each stage and the liquid flow below it, and the test that ends the iteration.
V_{j+1} = (V_j (H_j − h_j) − L_{j−1} (h_{j−1} − h_j) − (H_F,j − F_j h_j)) ÷ (H_{j+1} − h_j); stop when max |ΔT| < tol and max |ΔV| ÷ |V| < 1e-9
Each flow is computed from the stage above, so a rounding error is carried down the column. A flow that comes out at or below 1e-12 is replaced by 1e-12 of the feed. The default tolerance on temperature is 1e-8 K.
const Vn = new Array(N), Ln = new Array(N); Vn[0] = (R + 1) * D; let Lprev = R * D, hprev = h0;
for (let j = 0; j < N - 1; j++) {
const rhs = Vn[j] * (HV[j] - hL[j]) - Lprev * (hprev - hL[j]) - (HF[j] - Fj[j] * hL[j]);
let vnext = rhs / (HV[j + 1] - hL[j]); if (!(vnext > 1e-12)) vnext = 1e-12 * Ftot; Vn[j + 1] = vnext;
Ln[j] = Lprev + vnext - Vn[j] + Fj[j]; if (!(Ln[j] > 1e-12)) Ln[j] = 1e-12 * Ftot; Lprev = Ln[j]; hprev = hL[j];
}
Ln[N - 1] = Ftot - D;
dV = 0; for (let j = 0; j < N; j++) dV = Math.max(dV, Math.abs(Vn[j] - V[j]) / (Math.abs(V[j]) + 1e-12)); for (let j = 0; j < N; j++) { V[j] = Vn[j]; L[j] = Ln[j]; }
if (dT < tol && dV < 1e-9 && iter > 1) break;The condenser and reboiler duties and how well the component and energy balances close.
Q_c = V₁ (H₁ − h₀); component closure = max |Σf − d − b| ÷ F; energy closure = |in − out| ÷ (|in| + 1)
The closures are reported with the result, so a column that stopped short of its tolerance shows it in these two numbers. The energy closure is relative to the energy in, with 1 W added to the denominator.
const QcDuty = V[0] * (HV[0] - h0); // heat removed (W, positive)
const Hbot = hL[N - 1], Qr = (Ftot - D) * Hbot + V[N - 1] * HV[N - 1] - (N > 1 ? L[N - 2] * hL[N - 2] : 0) - HF[N - 1]; // heat added to the reboiler stage
const dist = fromFlash(D, th.flashPT(Tc, Pj[0], y[0])), bot = fromFlash(Ftot - D, th.flashPT(T[N - 1], Pj[N - 1], x[N - 1]));
const dFlow = lij ? y[0].map((v) => v * D) : null, bFlow = x[N - 1].map((v) => v * (Ftot - D));
let comp = 0; for (let i = 0; i < n; i++) comp = Math.max(comp, Math.abs(f.reduce((s, r) => s + r[i], 0) - dFlow[i] - bFlow[i]) / Ftot);
// energy closure: feeds + Qr = D h0 + B hB + Qc (condenser removes QcDuty from the vapour V1 to saturated liquid)
const eIn = sum(HF) + Qr, eOut = D * h0 + (Ftot - D) * Hbot + QcDuty - (R + 1) * D * 0 ;Tests tools/process-test/test-chem4.mjs
The solution of a tridiagonal system, used for the component balances.
forward: m = a_j ÷ b_{j−1}; b_j −= m c_{j−1}; r_j −= m r_{j−1}; back: x_j = (r_j − c_j x_{j+1}) ÷ b_j
The Thomas algorithm on copies of the diagonal and right side. There is no pivoting and no test for a zero pivot.
function thomas(a, b, c, r) {
const N = b.length, bb = b.slice(), rr = r.slice();
for (let j = 1; j < N; j++) { const m = a[j] / bb[j - 1]; bb[j] -= m * c[j - 1]; rr[j] -= m * rr[j - 1]; }
const x = new Array(N); x[N - 1] = rr[N - 1] / bb[N - 1]; for (let j = N - 2; j >= 0; j--) x[j] = (rr[j] - c[j] * x[j + 1]) / bb[j];
return x;
}The local K-value and enthalpy models the inner loop uses in place of the property package.
ln K_b = A + B ÷ T; B = (ln K_b(T + 0.5) − ln K_b(T)) ÷ (1 ÷ (T + 0.5) − 1 ÷ T); h = h_a + h_b T with h_b = (h(T + 0.5) − h(T)) ÷ 0.5
Slopes are one-sided differences over 0.5 K: each subtracts two nearly equal property values. The models are refitted at every outer pass.
const k0 = cth.lnK(T[j], Pj[j], x[j], y[j]).lnK, k1 = cth.lnK(T[j] + dTfd, Pj[j], x[j], y[j]).lnK;
const w = normed(y[j].map((v) => v + 1e-300));
let lkb = 0, lkb1 = 0; for (let i = 0; i < n; i++) { lkb += w[i] * k0[i]; lkb1 += w[i] * k1[i]; }
B[j] = (lkb1 - lkb) / (1 / (T[j] + dTfd) - 1 / T[j]); A[j] = lkb - B[j] / T[j];
lnK0.push(k0); alpha.push(k0.map((v) => Math.exp(v - lkb)));
const l0 = cth.phase(T[j], Pj[j], x[j], 'L').H, l1 = cth.phase(T[j] + dTfd, Pj[j], x[j], 'L').H, v0 = cth.phase(T[j], Pj[j], y[j], 'V').H, v1 = cth.phase(T[j] + dTfd, Pj[j], y[j], 'V').H;
hb[j] = (l1 - l0) / dTfd; ha[j] = l0 - hb[j] * T[j]; Hb[j] = (v1 - v0) / dTfd; Ha[j] = v0 - Hb[j] * T[j]; hL[j] = l0; HV[j] = v0;The temperature and vapor composition of each stage from the simple K model.
K_b = 1 ÷ Σ αᵢ xᵢ; T = B ÷ (ln K_b − A); yᵢ = αᵢ xᵢ, normalized
A closed-form bubble point on the fitted model, so no inner iteration is needed for the temperature. 1e-300 is added to each component flow before normalizing so an absent component does not divide by zero.
let dT = 0;
for (let j = 0; j < N; j++) {
const xj = normed(lij.map((row) => row[j] + 1e-300)); x[j] = xj;
let s = 0; for (let i = 0; i < n; i++) s += alpha[j][i] * xj[i]; // bubble point on the simple model: Kb = 1 / sum(alpha x)
const lkb = -Math.log(s), tn = B[j] / (lkb - A[j]); dT = Math.max(dT, Math.abs(tn - T[j])); T[j] = tn;
y[j] = normed(alpha[j].map((al, i) => al * xj[i]));
}When the inner and outer loops of the inside-out column end.
inner: max |ΔT| < 1e-10 K and max |ΔV| ÷ |V| < 1e-11; outer: max |ΔT| < tol and max |ΔV| ÷ |V| < 1e-8
The inner loop is successive substitution, at most 200 passes per outer pass. The outer tolerance on temperature is 1e-8 K by default, and the column is reported converged when the last outer change was below ten times that and 1e-6 on the flows.
let dV = 0; for (let j = 0; j < N; j++) { dV = Math.max(dV, Math.abs(Vn[j] - V[j]) / (Math.abs(V[j]) + 1e-12)); V[j] = Vn[j]; L[j] = Ln[j]; }
if (dT < 1e-10 && dV < 1e-11) break;
}
dTouter = 0; dVouter = 0;
for (let j = 0; j < N; j++) { dTouter = Math.max(dTouter, Math.abs(T[j] - Told[j])); dVouter = Math.max(dVouter, Math.abs(V[j] - Vold[j]) / (Math.abs(Vold[j]) + 1e-12)); }
if (dTouter < tol && dVouter < 1e-8 && outer > 1) break;The vapor and liquid flow leaving each stage of an absorber or stripper.
V_j = Σᵢ l_ij K_ij V_j ÷ L_j; L_j = F(1..j) − V₁ + V_{j+1}
The vapor flow is a sum of component flows in index order; the liquid flow is a running total balance. Both are floored at 1e-12 of the feed, and the update is damped by 0.7 after the third pass.
const vj = lj.map((l, i) => l * K[j][i] * V[j] / L[j]); Vn[j] = Math.max(sum(vj), 1e-12 * Ftot); yn.push(normed(vj.map((v) => v + 1e-300)));
}
// total balance gives L; V is the sum of the vapour component flows (sum-rates)
const Ln = new Array(N); let acc = 0; for (let j = 0; j < N; j++) { acc += Fj[j]; Ln[j] = Math.max(acc - Vn[0] - (j > 0 ? 0 : 0) + (j > 0 ? 0 : 0), 1e-12 * Ftot); }
for (let j = 0; j < N - 1; j++) Ln[j] = Math.max(sum(Fj.slice(0, j + 1)) - Vn[0] + Vn[j + 1] + 0, 1e-12 * Ftot); // L_j = F(1..j) - V_1 + V_{j+1}
Ln[N - 1] = Math.max(Ftot - Vn[0], 1e-12 * Ftot);
const dampV = iter < 4 ? 1 : 0.7; dV = 0;
for (let j = 0; j < N; j++) { const v = V[j] + dampV * (Vn[j] - V[j]); dV = Math.max(dV, Math.abs(v - V[j]) / (Math.abs(V[j]) + 1e-12)); V[j] = v; L[j] = L[j] + dampV * (Ln[j] - L[j]); }The temperature of every stage from the stage energy balances.
E_j = L_{j−1} h_{j−1} + V_{j+1} H_{j+1} + H_F,j + Q_j − L_j h_j − V_j H_j = 0; J ΔT = −E with heat capacities by a 0.5 K difference
One Newton step per pass on a tridiagonal Jacobian, solved by the dense linear solver. Steps are capped at 30 K and halved for the first five passes. The loop ends when the step is below 1e-6 K and the flow change below 1e-8.
const E = new Array(N), Jm = Array.from({ length: N }, () => new Array(N).fill(0)), dd = 0.5;
for (let j = 0; j < N; j++) {
hL[j] = hh(j, T[j]); HV[j] = HH(j, T[j]);
const cpl = (hh(j, T[j] + dd) - hL[j]) / dd, cpv = (HH(j, T[j] + dd) - HV[j]) / dd;
const inL = j > 0 ? L[j - 1] * (hh(j - 1, T[j - 1])) : 0, inV = j < N - 1 ? V[j + 1] * HH(j + 1, T[j + 1]) : 0;
E[j] = inL + inV + HF[j] + Qs[j] - L[j] * hL[j] - V[j] * HV[j];
Jm[j][j] = -(L[j] * cpl + V[j] * cpv);
if (j > 0) Jm[j][j - 1] = L[j - 1] * ((hh(j - 1, T[j - 1] + dd) - hh(j - 1, T[j - 1])) / dd);
if (j < N - 1) Jm[j][j + 1] = V[j + 1] * ((HH(j + 1, T[j + 1] + dd) - HH(j + 1, T[j + 1])) / dd);
}
const dTv = solveLinear(Jm, E.map((v) => -v)); dT = 0;
if (dTv) for (let j = 0; j < N; j++) { const s = Math.max(-30, Math.min(30, dTv[j])) * (iter < 6 ? 0.5 : 1); T[j] += s; dT = Math.max(dT, Math.abs(s)); }
else dT = Infinity;
if (dT < tol * 100 && dV < 1e-8 && iter > 2) break;The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.