Resources Numerical accuracy Calculation policy
Splits a mixture into a vapor and two liquids at a temperature and pressure with an activity-coefficient model.
64-bit throughout. The phase fractions come from a damped Newton iteration on the multiphase Rachford–Rice equations and the compositions from successive substitution on the K-values, stopped at 1e-11 on ln K. The result says whether that tolerance was met and gives the iso-fugacity and balance residuals.
4 calculations in 1 source file. 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-dec.mjs
The Newton step for the two phase fractions.
det = a d − b c; x = (r₁ d − b r₂) ÷ det; y = (a r₂ − c r₁) ÷ det
Cramer's rule. Only a determinant below 1e-300 (absolute) is refused, so a nearly singular Hessian is solved as it stands; the step is then checked by the line search.
const solve2 = (a, b, c, d, r1, r2) => { const det = a * d - b * c; return Math.abs(det) < 1e-300 ? null : [(r1 * d - b * r2) / det, (a * r2 - c * r1) / det]; };The vapor and second-liquid fractions for the current K-values.
gⱼ = Σ zᵢ (Kⱼᵢ − 1) ÷ Dᵢ = 0, Dᵢ = 1 + β_v (K_v,ᵢ − 1) + β_l (K_l,ᵢ − 1)
Newton with the analytic Hessian. A step is halved until both fractions are non-negative, their sum is at most 1 and every Dᵢ is above 1e-12. A phase at zero whose gradient is not positive is held at zero.
for (let k = 0; k < 200; k++) {
const D = den(bv, bl);
let g1 = 0, g2 = 0, h11 = 0, h12 = 0, h22 = 0;
for (let i = 0; i < n; i++) {
const a = Kv[i] - 1, b = Kl[i] - 1, q = z[i] / (D[i] * D[i]);
g1 += z[i] * a / D[i]; g2 += z[i] * b / D[i]; h11 -= q * a * a; h12 -= q * a * b; h22 -= q * b * b;
}
const stuckV = bv <= 1e-14 && g1 <= 0, stuckL = bl <= 1e-14 && g2 <= 0;
let dv = 0, dl = 0;
if (stuckV && stuckL) break;
if (stuckV) dl = h22 ? -g2 / h22 : 0; else if (stuckL) dv = h11 ? -g1 / h11 : 0; else { const s = solve2(h11, h12, h12, h22, -g1, -g2); if (!s) break; [dv, dl] = s; }
let t = 1, moved = false;
for (let j = 0; j < 60; j++) {
const nv = Math.max(bv + t * dv, 0), nl = Math.max(bl + t * dl, 0), Dn = den(nv, nl);
if (nv + nl <= 1 && Dn.every((v) => v > 1e-12)) { bv = nv; bl = nl; moved = true; break; }
t *= 0.5;
}
if (!moved || (Math.abs(g1) < 1e-14 && Math.abs(g2) < 1e-14) || Math.abs(t * dv) + Math.abs(t * dl) < 1e-16) break;
}The liquid compositions and the new K-values after each pass.
x₁,ᵢ = zᵢ ÷ Dᵢ; x₂,ᵢ = x₁,ᵢ K_l,ᵢ; error = max |ln(K_new ÷ K)|
Compositions are normalized to sum to 1 after each update. The loop ends when the largest change in ln K is below the tolerance (1e-11 by default).
const D = den(bv, bl), nx1 = norm(z.map((v, i) => v / D[i])), nx2 = norm(nx1.map((v, i) => v * Kl[i]));
const nKv = kv(nx1), nKl = bl > 1e-12 ? kl(nx1, nx2) : Kl;
let err = 0; for (let i = 0; i < n; i++) err = Math.max(err, Math.abs(Math.log(nKv[i] / Kv[i])), bl > 1e-12 ? Math.abs(Math.log(nKl[i] / Kl[i])) : 0);
x1 = nx1; x2 = nx2; Kv = nKv; Kl = nKl;
if (err < tol) { converged = true; break; }The largest iso-fugacity difference between the phases and the largest component-balance error, as the calculator reports them.
iso = max |ln(xᵢγᵢ)₁ − ln(xᵢγᵢ)₂|, |ln(yᵢP) − ln(xᵢγᵢP_sat,ᵢ)|; balance = max |Σ_phases frac × xᵢ − zᵢ|
Computed from the returned phases, independently of the iteration's own stopping test. A component with zero mole fraction in a phase would give the logarithm of zero.
/** Iso-fugacity across the phases that exist (max |ln| difference) and the component balances: what tests and the calculator report. */
export function vlleResiduals(th, r) {
const g = (x) => th.lnGamma(x, r.T), liq = r.phases.filter((p) => p.type.startsWith('liquid')), v = r.phases.find((p) => p.type === 'vapor');
let iso = 0, bal = 0;
if (liq.length === 2) { const a = g(liq[0].x), b = g(liq[1].x); for (let i = 0; i < th.n; i++) iso = Math.max(iso, Math.abs(Math.log(liq[0].x[i]) + a[i] - Math.log(liq[1].x[i]) - b[i])); }
if (v && liq.length) { const a = g(liq[0].x); for (let i = 0; i < th.n; i++) iso = Math.max(iso, Math.abs(Math.log(v.x[i] * r.P) - Math.log(liq[0].x[i]) - a[i] - Math.log(th.psatIdeal(i, r.T)))); }
for (let i = 0; i < th.n; i++) bal = Math.max(bal, Math.abs(r.phases.reduce((s, p) => s + p.frac * p.x[i], 0) - r.z[i]));
return { iso, bal };
}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.