Works out phase splits, bubble and dew points, enthalpy, entropy and volume of a mixture for the flowsheet and the calculators.
How the precision is kept
Everything is 64-bit. The flashes are iterations with tight stopping rules (1e-13 on ln K, 1e-13 relative on a saturation bracket), so the arithmetic is not what limits them. The results rest on a gas constant typed to ten digits, textbook component data of three to six digits, and correlations (Antoine, Watson, Rackett, Wilson's K estimate) for what the data does not give.
Open item. In the bubble and dew point scan, throw away an incipient phase that has collapsed onto the feed and confirm the bracket with the stability test, so a flip of the single-fluid flag is not taken for a saturation point.
Calculations under this policy
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.
The ideal-gas enthalpy of a component, with its formation enthalpy, at a temperature.
H = h_f + R (A (T − T₀) + B⁄2 (T² − T₀²) + C⁄3 (T³ − T₀³) − D (1⁄T − 1⁄T₀))
The integral of the heat-capacity polynomial is written as differences of powers of T and 298.15, so close to 298.15 K the sensible part carries fewer digits. The formation enthalpy is added to it.
Watson's scaling with the exponent 0.38. When the value at the boiling point is missing it is estimated from the critical constants with typed coefficients (3.978, 3.958, 1.555, 1.07), and a missing boiling point is taken as 0.65 T_c.
process/thermo.js
lines 53 to 59
CadCraft 0.8.2
/** Heat of vaporisation J/mol at T: Watson scaling from the value at Tb (Chen's estimate when Tb or the value is missing). */
dHvap(i, T) {
const c = this.c[i]; if (T >= c.Tc) return 0;
const Tb = c.Tb || 0.65 * c.Tc, Tbr = Tb / c.Tc;
const h = c.dHvapTb || R * c.Tc * Tbr * (3.978 * Tbr - 3.958 + 1.555 * Math.log(c.Pc / 1e5)) / (1.07 - Tbr);
return h * Math.pow((1 - T / c.Tc) / (1 - Tbr), 0.38);
}
Liquid molar volume
64-bit
The saturated liquid volume of a component for the ideal and activity models.
V = R T_c ÷ P_c × Z_ra^(1 + (1 − T_r)^(2⁄7)), Z_ra = 0.29056 − 0.08775 ω
The Rackett correlation; its two typed coefficients set the accuracy. The reduced temperature is capped at 1.
The molar enthalpy, entropy and volume of a liquid or vapor of given composition.
EOS: H = Σ xᵢ H_ig,ᵢ + departure; V = ZRT ÷ P. Activity liquid: H = Σ xᵢ (H_ig,ᵢ − ΔH_vap,ᵢ) + H^E, H^E = −RT² Σ xᵢ ∂ln γᵢ⁄∂T
Sums over the components in index order. For the activity models the excess enthalpy uses a central difference of ln γ with a step of 0.001 × T, so it carries the truncation and rounding error of that difference.
process/thermo.js
lines 67 to 89
CadCraft 0.8.2
/** Enthalpy, entropy, volume of one phase. type 'V' | 'L'. EOS phases take the root of that type (or the only root). */
phase(T, P, x, type) {
const n = this.n; let H = 0, S = 0, V, Z = null;
if (this.eos) {
const r = this.eos.roots(T, P, x), zs = r.z; if (!zs.length) throw new Error('no EOS root');
Z = type === 'L' ? zs[0] : zs[zs.length - 1];
const d = this.eos.departure(T, P, x, Z, r);
for (let i = 0; i < n; i++) { H += x[i] * this.hIg(i, T); S += x[i] * this.sIg(i, T); if (x[i] > 0) S -= R * x[i] * Math.log(x[i]); }
H += d.H; S += d.S - R * Math.log(P / P0); V = Z * R * T / P;
return { H, S, V, Z, liquidLike: this.eos.liquidLike(Z, r.beta) };
}
for (let i = 0; i < n; i++) { H += x[i] * this.hIg(i, T); S += x[i] * this.sIg(i, T); if (x[i] > 0) S -= R * x[i] * Math.log(x[i]); }
if (type === 'V') { S -= R * Math.log(P / P0); return { H, S, V: R * T / P, Z: 1, liquidLike: false }; }
// activity / ideal liquid: ideal-gas enthalpy minus the heat of vaporisation, plus the excess terms
let vol = 0, lg = null;
for (let i = 0; i < n; i++) { H -= x[i] * this.dHvap(i, T); S -= x[i] * (this.dHvap(i, T) / T + R * Math.log(this.psatIdeal(i, T) / P0)); vol += x[i] * this.vLiq(i, T); }
if (this.model !== 'ideal') {
lg = this.lnGamma(x, T); const dT = 1e-3 * T, a = this.lnGamma(x, T + dT), b = this.lnGamma(x, T - dT);
let HE = 0, GE = 0; for (let i = 0; i < n; i++) { HE += -R * T * T * x[i] * (a[i] - b[i]) / (2 * dT); GE += R * T * x[i] * lg[i]; }
H += HE; S += (HE - GE) / T;
}
return { H, S, V: vol, Z: P * vol / (R * T), liquidLike: true };
}
K-values
64-bit
ln K of every component between a liquid and a vapor composition.
EOS: ln Kᵢ = ln φᵢ(liquid) − ln φᵢ(vapor); activity: ln Kᵢ = ln γᵢ + ln(P_sat,ᵢ ÷ P), or ln(Hᵢ ÷ P) for a Henry solute
A difference of two logarithms of fugacity coefficients. Near a critical point the two are nearly equal and ln K loses relative digits.
process/thermo.js
lines 99 to 119
CadCraft 0.8.2
/** ln K = ln(x-phase fugacity ratio): liquid x against vapour y. Returns { lnK, trivial } where trivial = ± 1 when the roots say
* "this is a single dense/dilute fluid" (see saturation()). */
lnK(T, P, x, y) {
const n = this.n, out = new Array(n);
if (this.eos) {
const rl = this.eos.roots(T, P, x), rv = this.eos.roots(T, P, y);
const zl = rl.z[0], zv = rv.z[rv.z.length - 1];
const ll = this.eos.lnphi(T, P, x, zl, rl), lv = this.eos.lnphi(T, P, y, zv, rv);
for (let i = 0; i < n; i++) out[i] = ll[i] - lv[i];
let trivial = 0;
if (rv.z.length === 1 && this.eos.liquidLike(zv, rv.beta)) trivial = -1;
else if (rl.z.length === 1 && !this.eos.liquidLike(zl, rl.beta)) trivial = 1;
return { lnK: out, trivial };
}
const lg = this.lnGamma(x, T);
for (let i = 0; i < n; i++) {
const H = this.henryAt(i, T);
out[i] = H != null ? Math.log(H / P) : lg[i] + Math.log(this.psatIdeal(i, T) / P);
}
return { lnK: out, trivial: 0 };
}
Vapor fraction (Rachford–Rice)
64-bit
The vapor fraction that balances a feed for given K-values.
Σ zᵢ (Kᵢ − 1) ÷ (1 + V (Kᵢ − 1)) = 0, 0 ≤ V ≤ 1
Newton steps kept inside a bracket that shrinks every pass, with bisection when a step leaves it. Stops when |f| < 1e-15 or the bracket is below 1e-16, at most 300 passes.
process/thermo.js
lines 122 to 138
CadCraft 0.8.2
/** Vapour fraction in [0,1] for feed z and K; returns { V, state: 'two'|'liquid'|'vapor' }. */
rachfordRice(z, K) {
const n = z.length; let s1 = 0, s2 = 0;
for (let i = 0; i < n; i++) { s1 += z[i] * K[i]; s2 += z[i] / K[i]; }
if (s1 <= 1) return { V: 0, state: 'liquid' };
if (s2 <= 1) return { V: 1, state: 'vapor' };
const f = (V) => { let s = 0; for (let i = 0; i < n; i++) s += z[i] * (K[i] - 1) / (1 + V * (K[i] - 1)); return s; };
let lo = 0, hi = 1, V = 0.5; // f(0) > 0 and f(1) < 0 here, f decreasing
for (let it = 0; it < 300; it++) {
const fv = f(V);
if (fv > 0) lo = V; else hi = V;
if (Math.abs(fv) < 1e-15 || hi - lo < 1e-16) break;
let d = 0; for (let i = 0; i < n; i++) { const a = K[i] - 1; d -= z[i] * a * a / ((1 + V * a) ** 2); }
const Vn = V - fv / d; V = (Vn > lo && Vn < hi) ? Vn : 0.5 * (lo + hi);
}
return { V, state: 'two' };
}
Phase stability test
64-bit
Whether a mixture at a temperature and pressure stays one phase or splits.
Successive substitution from a vapor-like and a liquid-like start, up to 500 passes each, to 1e-10 on ln Y. A trial within 1e-4 of the feed counts as the trivial solution only if the passes converged.
process/thermo.js
lines 141 to 162
CadCraft 0.8.2
/** → { stable, trials[] }; if unstable, `trial` is the composition of the incipient phase and `kind` 'V' or 'L'. */
stability(T, P, z) {
const e = this.eos, n = this.n, r = e.roots(T, P, z);
let best = r.z[0], bg = Infinity;
for (const Z of r.z) { const g = e.gibbs(T, P, z, Z, r); if (g < bg) { bg = g; best = Z; } }
const lz = e.lnphi(T, P, z, best, r), d = z.map((v, i) => (v > 0 ? Math.log(v) + lz[i] : -Infinity));
const K = this.wilsonK(T, P);
for (const kind of ['V', 'L']) {
let Y = z.map((v, i) => (kind === 'V' ? v * K[i] : v / K[i]));
let y = normed(Y), conv = false;
for (let it = 0; it < 500; it++) {
y = normed(Y); const ry = e.roots(T, P, y), Zy = kind === 'V' ? ry.z[ry.z.length - 1] : ry.z[0], ly = e.lnphi(T, P, y, Zy, ry);
const Yn = z.map((v, i) => (v > 0 ? Math.exp(d[i] - ly[i]) : 0)); let err = 0;
for (let i = 0; i < n; i++) err = Math.max(err, Math.abs(Math.log((Yn[i] + 1e-300) / (Y[i] + 1e-300))));
Y = Yn; if (err < 1e-10) { conv = true; break; }
}
const S = sum(Y), yy = normed(Y);
if (maxAbs(yy, z) < 1e-4 && conv) continue; // trivial solution (the feed itself)
if (S > 1 + 1e-8) return { stable: false, trial: yy, kind };
}
return { stable: true, root: best };
}
Two-phase flash iteration
64-bit
The vapor fraction and the two compositions at a temperature and pressure.
V from Rachford–Rice; xᵢ = zᵢ ÷ (1 + V (Kᵢ − 1)); yᵢ = Kᵢ xᵢ; K ← exp(ln K(x, y)) until max |Δ ln K| < 1e-13
Successive substitution, at most 2000 passes. The pass count is returned, but nothing marks a result that used all 2000 without meeting the tolerance.
process/thermo.js
lines 165 to 180
CadCraft 0.8.2
/** Successive substitution VLE from initial K; returns { V, x, y, K, iter, trivial }. */
_ss(T, P, z, K) {
const n = this.n; let V = 0.5, x = z.slice(), y = z.slice(), it = 0;
for (it = 0; it < 2000; it++) {
const rr = this.rachfordRice(z, K); V = rr.V;
x = normed(z.map((v, i) => v / (1 + V * (K[i] - 1)))); y = normed(x.map((v, i) => K[i] * v));
if (rr.state === 'liquid') { x = z.slice(); y = normed(z.map((v, i) => K[i] * v)); }
if (rr.state === 'vapor') { y = z.slice(); x = normed(z.map((v, i) => v / K[i])); }
const { lnK } = this.lnK(T, P, x, y); let err = 0; const Kn = new Array(n);
for (let i = 0; i < n; i++) { Kn[i] = Math.exp(lnK[i]); err = Math.max(err, Math.abs(lnK[i] - Math.log(K[i]))); }
K = Kn;
if (err < this.tolSS) break;
}
let tr = 0; for (let i = 0; i < n; i++) tr = Math.max(tr, Math.abs(Math.log(K[i])));
return { V, x, y, K, iter: it, trivial: tr < 1e-6 };
}
Liquid–liquid split
64-bit
Two liquid phases and their fractions for an activity model.
trial phases from each component at 0.9; split when the tangent-plane distance < −1e-9; K ← exp(ln γ(x₂) − ln γ(x₁)) until max |Δ ln K| < 1e-12
The phase fraction is found by 200 bisections on [0, 1], which reaches 64-bit resolution. A split is reported only when the fraction is inside (1e-12, 1 − 1e-12) and the two compositions differ by more than 1e-6.
process/thermo.js
lines 222 to 253
CadCraft 0.8.2
/** Liquid–liquid split (activity models). Phase 'liquid' = the one richer in the first seed key, 'liquid2' the other. */
flashLL(T, P, z) {
if (this.eos || this.model === 'ideal') throw new Error('liquid–liquid flash needs an activity model');
z = normed(z); const n = this.n, lz = this.lnGamma(z, T);
let best = null; const trials = [];
for (let k = 0; k < n; k++) {
let w = z.map((_, i) => (i === k ? 0.9 : 0.1 / (n - 1)));
for (let it = 0; it < 500; it++) {
const lw = this.lnGamma(w, T), Y = z.map((v, i) => Math.exp(Math.log(v) + lz[i] - lw[i])), wn = normed(Y);
const e = maxAbs(wn, w); w = wn; if (e < 1e-12) break;
}
const lw = this.lnGamma(w, T); let tpd = 0; for (let i = 0; i < n; i++) if (w[i] > 0) tpd += w[i] * (Math.log(w[i]) + lw[i] - Math.log(z[i]) - lz[i]);
if (maxAbs(w, z) > 1e-4) { trials.push({ w, tpd }); if (!best || tpd < best.tpd) best = { w, tpd }; }
}
if (!best || best.tpd > -1e-9) return this._single(T, P, z, 'L');
// successive substitution on K_i = x1_i / x2_i with x1 the seed-rich phase
let other = z.slice(), dmax = 0; for (const t of trials) { const d = maxAbs(t.w, best.w); if (d > dmax && d > 1e-3) { dmax = d; other = t.w; } }
let K = best.w.map((v, i) => v / other[i]), V = 0.5, x1 = best.w, x2 = z.slice();
for (let it = 0; it < 3000; it++) {
const f = (v) => { let s = 0; for (let i = 0; i < n; i++) s += z[i] * (K[i] - 1) / (1 + v * (K[i] - 1)); return s; };
let lo = 0, hi = 1; if (f(0) <= 0) { V = 0; break; } if (f(1) >= 0) { V = 1; break; }
for (let j = 0; j < 200; j++) { V = 0.5 * (lo + hi); if (f(V) > 0) lo = V; else hi = V; }
x2 = normed(z.map((v, i) => v / (1 + V * (K[i] - 1)))); x1 = normed(x2.map((v, i) => K[i] * v));
const g1 = this.lnGamma(x1, T), g2 = this.lnGamma(x2, T); let err = 0; const Kn = K.slice();
for (let i = 0; i < n; i++) { const kn = Math.exp(g2[i] - g1[i]); err = Math.max(err, Math.abs(Math.log(kn / K[i]))); Kn[i] = kn; }
K = Kn; if (err < 1e-12) break;
}
if (V <= 1e-12 || V >= 1 - 1e-12 || maxAbs(x1, x2) < 1e-6) return this._single(T, P, z, 'L');
const p1 = this.phase(T, P, x1, 'L'), p2 = this.phase(T, P, x2, 'L');
return { T, P, z, V: 0, phases: [{ type: 'liquid', frac: V, x: x1, H: p1.H, S: p1.S, v: p1.V }, { type: 'liquid2', frac: 1 - V, x: x2, H: p2.H, S: p2.S, v: p2.V }],
H: V * p1.H + (1 - V) * p2.H, S: V * p1.S + (1 - V) * p2.S, v: V * p1.V + (1 - V) * p2.V, state: 'liquid-liquid' };
}
Saturation residual
64-bit
How far a mixture is from its bubble or dew point at a trial temperature or pressure.
bubble: Σ Kᵢ zᵢ − 1; dew: 1 − Σ zᵢ ÷ Kᵢ
The incipient composition is iterated to 1e-13, at most 300 passes, starting from the composition kept from the previous trial point. When the equation-of-state roots report a single fluid, the function returns that flag (±1) in place of the residual.
process/thermo.js
lines 256 to 273
CadCraft 0.8.2
/** r > 0 means the mixture is "more vapour" than saturated (T above / P below the saturation point). */
_satResidual(kind, T, P, z, st) {
const n = this.n;
let x, y;
if (kind === 'bubble') { x = z; y = st.y || normed(this.wilsonK(T, P).map((k, i) => k * z[i])); }
else { y = z; x = st.x || normed(this.wilsonK(T, P).map((k, i) => z[i] / k)); }
let s = 0, tr = 0;
for (let it = 0; it < 300; it++) {
const { lnK, trivial } = this.lnK(T, P, x, y); tr = trivial;
const K = lnK.map(Math.exp);
if (kind === 'bubble') { const t = z.map((v, i) => K[i] * v); s = sum(t); const yn = normed(t); const e = maxAbs(yn, y); y = yn; if (e < 1e-13) break; }
else { const t = z.map((v, i) => v / K[i]); s = sum(t); const xn = normed(t); const e = maxAbs(xn, x); x = xn; if (e < 1e-13) break; }
}
st.x = x; st.y = y;
if (tr !== 0 && (kind === 'bubble' || kind === 'dew')) { st.trivial = true; return tr; }
st.trivial = false;
return kind === 'bubble' ? s - 1 : 1 - s;
}
Bubble and dew points
64-bit
The bubble or dew temperature at a pressure, or pressure at a temperature.
scan a geometric grid for a sign change of the residual, then bisect on √(a b) until b ÷ a − 1 < 1e-13
The bracket is refined in the logarithm of the unknown to a relative 1e-13. The scan takes the first sign change it meets; because the residual can be the single-fluid flag and the incipient composition is carried from point to point, that sign change can be a flip of the flag and not a zero of the sum.
process/thermo.js
lines 274 to 300
CadCraft 0.8.2
/** Bubble or dew point. kind 'bubble'|'dew'; fixed 'T' (value K → returns P) or 'P' (value Pa → returns T). → { T, P, x, y } */
saturation(kind, fixed, value, z, guess = null) {
z = normed(z); const st = {}, isT = fixed === 'P'; // solving for T when P is fixed
const res = (v) => { const T = isT ? v : value, P = isT ? value : v; return { r: this._satResidual(kind, T, P, z, st), T, P }; };
const sgn = isT ? 1 : -1; // r increases with T, decreases with P
let lo = null, hi = null; // lo: r*sgn<0 side, hi: r*sgn>0
const grid = [];
if (guess) { for (let k = 0; k < 40; k++) { grid.push(guess * Math.pow(1.08, k)); grid.push(guess / Math.pow(1.08, k + 1)); } }
else { const a = isT ? 30 : 1, b = isT ? 2500 : 1e10, f = isT ? 1.04 : 1.12; for (let v = a; v < b; v *= f) grid.push(v); }
if (guess) {
for (const v of grid) { const q = res(v).r * sgn; if (q < 0) { if (lo == null || Math.abs(Math.log(v / guess)) < Math.abs(Math.log(lo / guess))) lo = v; } else if (hi == null || Math.abs(Math.log(v / guess)) < Math.abs(Math.log(hi / guess))) hi = v; if (lo && hi) break; }
} else {
let prev = null;
for (const v of grid) { const q = res(v).r * sgn; if (prev !== null && prev < 0 && q >= 0) { lo = prev_v(grid, v); hi = v; break; } prev = q; }
}
if (lo == null || hi == null) return null;
if (lo > hi) [lo, hi] = [hi, lo];
// bisect in the log of the variable, keeping the side signs
const sLo = res(lo).r * sgn; let a = lo, b = hi, sa = sLo;
for (let it = 0; it < 200; it++) {
const m = Math.sqrt(a * b), sm = res(m).r * sgn;
if (sm < 0) { a = m; sa = sm; } else b = m;
if (b / a - 1 < 1e-13) break;
}
const v = Math.sqrt(a * b), q = res(v);
return { T: q.T, P: q.P, x: st.x, y: st.y, kind };
}
Temperature for a given enthalpy or entropy
64-bit
The temperature at which a PT flash has the wanted enthalpy or entropy.
false position inside a bracket; stop when b − a < 1e-10 K or |f| < 1e-9 (|target| + 1)
A bracket is searched around the guess in widening steps, then the whole range 20 to 3500 K. Each trial runs a full PT flash. The residual tolerance is 1e-9 times the size of the target plus 1.
process/thermo.js
lines 306 to 328
CadCraft 0.8.2
// ---- PH, PS, and a general enthalpy helper -----------------------------------------------------------------------------------
_solveT(P, z, target, prop, T0g, opts = {}) {
const f = (T) => this.flashPT(T, P, z)[prop] - target;
let a = opts.Tmin || 20, b = opts.Tmax || 3500, g = T0g && T0g > a && T0g < b ? T0g : 300;
// bracket around the guess first (cheap), then fall back to the full range
let fa = null, fb = null;
for (const step of [5, 25, 100, 400, 1500]) {
const lo = Math.max(opts.Tmin || 20, g - step), hi = Math.min(opts.Tmax || 3500, g + step), flo = f(lo), fhi = f(hi);
if (flo <= 0 && fhi >= 0) { a = lo; b = hi; fa = flo; fb = fhi; break; }
if (step === 1500) { a = opts.Tmin || 20; b = opts.Tmax || 3500; fa = f(a); fb = f(b); }
}
if (fa > 0 || fb < 0) return null;
for (let it = 0; it < 200; it++) {
const m = 0.5 * (a + b), fm = f(m);
if (fm === 0) { a = b = m; break; }
// regula falsi proposal kept inside the bracket
let t = a - fa * (b - a) / (fb - fa); if (!(t > a && t < b)) t = m;
const ft = f(t);
if (ft < 0) { a = t; fa = ft; } else { b = t; fb = ft; }
if (b - a < 1e-10 || Math.abs(ft) < 1e-9 * (Math.abs(target) + 1)) { a = b = t; break; }
}
return 0.5 * (a + b);
}
Flash at a vapor fraction
64-bit
The state at a pressure and a wanted vapor fraction.
bisect T between the bubble and dew temperatures until the bracket is below 1e-9 K
The stop is on the temperature bracket, not on the vapor fraction, so the fraction returned is the one the PT flash gives at that temperature. A pure component uses the lever rule at its saturation temperature.
process/thermo.js
lines 345 to 353
CadCraft 0.8.2
/** Flash with a vapour fraction specified (0 = bubble, 1 = dew) at P. */
flashPVF(P, vf, z) {
z = normed(z); if (this.n === 1) { const pd = this._pureDomeVF(P, vf); if (pd) return pd; }
const b = this.bubbleT(z, P), d = this.dewT(z, P); if (!b || !d) return null;
if (vf <= 0) return this.flashPT(b.T, P, z); if (vf >= 1) return this.flashPT(d.T, P, z);
let lo = b.T, hi = d.T;
for (let i = 0; i < 100; i++) { const m = 0.5 * (lo + hi), f = this.flashPT(m, P, z); if (f.V < vf) lo = m; else hi = m; if (hi - lo < 1e-9) break; }
return this.flashPT(0.5 * (lo + hi), P, z);
}
Saturation temperature of a pure component
64-bit
The temperature at which a pure component's vapor pressure equals the pressure.
bisect T on [20 K, T_c − 1e-6] until the bracket is below 1e-11 K
Plain bisection, at most 200 halvings; each trial computes a vapor pressure, which for the equation-of-state models is itself an iteration to 1e-12.
process/thermo.js
line 344
CadCraft 0.8.2
_pureTsat(P) { const f = this.eos ? (T) => this.eos.psat(0, T) : (T) => this.psatIdeal(0, T); let lo = 20, hi = this.c[0].Tc - 1e-6; if (f(hi) < P || f(lo) > P) return null; for (let i = 0; i < 200; i++) { const m = 0.5 * (lo + hi); if (f(m) < P) lo = m; else hi = m; if (hi - lo < 1e-11) break; } return 0.5 * (lo + hi); }
R, used by the property package, the equations of state, the reactors and the dynamic simulation.
R = 8.314462618 J/(mol K)
Typed to ten significant digits. The value is fixed by definition as 8.31446261815324, so every result that multiplies by R carries ten digits from it, not fifteen. mechanism.js in the same engine uses the full value.
process/data.js
line 17
CadCraft 0.8.2
export const R = 8.314462618; // J/(mol K), CODATA 2018
Component data to SI
64-bit
A component's critical pressure, molar mass and Antoine constants in the units the engine works in.
P_c (Pa) = P_c (bar) × 1e5; M (kg/mol) = M ÷ 1000; A′ = A + ln 1000; C′ = C − 273.15
Exact powers of ten, one logarithm and one subtraction. The data it converts is the textbook table above it, of three to six significant digits.
process/data.js
lines 58 to 64
CadCraft 0.8.2
/** A component in SI units: Pc Pa, MW kg/mol, Antoine as ln P(Pa) = A' - B/(T[K] + C'). Anything missing stays null. */
export function normalise(c) {
const o = { ...c, Pc: c.Pc * 1e5, MWkg: c.MW / 1000 };
if (c.antoine) o.antoineSI = { A: c.antoine.A + Math.log(1000), B: c.antoine.B, C: c.antoine.C - 273.15 };
o.cp = { A: 0, B: 0, C: 0, D: 0, ...c.cp };
return o;
}
In detail
Method
Rachford–Rice by Newton inside a bracket; successive substitution on K; Michelsen's stability test for the equation-of-state models; liquid–liquid split for activity models; bubble and dew points by a grid scan and bisection in the logarithm of the unknown; PH and PS flashes by false position inside a bracket.
Number types
64-bit numbers in plain arrays. R = 8.314462618 J/(mol K), T₀ = 298.15 K, P₀ = 1e5 Pa.
Tolerances
K iteration: max |Δ ln K| < 1e-13, 2000 passes. Rachford–Rice: |f| < 1e-15 or bracket < 1e-16, 300 passes. Stability: 1e-10 on ln Y, 500 passes; a trial closer than 1e-4 to the feed is trivial; unstable when ΣY > 1 + 1e-8. Saturation: inner 1e-13, 300 passes; bracket to 1e-13 relative, 200 halvings. PH/PS temperature: bracket below 1e-10 K or residual below 1e-9 × (|target| + 1). Vapor-fraction flash: 1e-9 K. Liquid–liquid: 1e-12.
Limits
The header of data.js says the component table was typed without the book open and is to be checked before design use. The K iteration and the stability test return their last iterate without a flag when their pass limits are reached. In the saturation scan a point where the equation of state reports a single fluid returns a flag of ±1 in place of a residual.
Approach: 64-bit, stated limit
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.