Works out compressibility, fugacity coefficients, enthalpy and entropy departures and vapor pressure from the Peng–Robinson and Soave–Redlich–Kwong equations.
How the precision is kept
The cubic is solved in closed form in 64-bit and each root gets Newton steps on the original polynomial. What limits the digits is the model's own constants, typed to four or five significant digits (0.45724, 0.07780, 0.42748, 0.08664), the gas constant taken to ten digits, and component data of three to six digits.
Open item. Write the Peng–Robinson and Soave–Redlich–Kwong constants to full 64-bit precision instead of four or five digits.
Calculations under this policy
8 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.
Ω_a and Ω_b are typed to four or five significant digits, so every a and b, and everything computed from them, carries no more than that from the model. σ and ε for Peng–Robinson come from Math.SQRT2 at full precision.
process/eos.js
lines 11 to 15
CadCraft 0.8.2
export const MODELS = {
pr: { name: 'Peng–Robinson', sigma: 1 + Math.SQRT2, eps: 1 - Math.SQRT2, oa: 0.45724, ob: 0.07780, zc: 0.3074,
kappa: (w) => (w <= 0.49 ? 0.37464 + 1.54226 * w - 0.26992 * w * w : 0.379642 + 1.48503 * w - 0.164423 * w * w + 0.016666 * w * w * w) },
srk: { name: 'Soave–Redlich–Kwong', sigma: 1, eps: 0, oa: 0.42748, ob: 0.08664, zc: 1 / 3, kappa: (w) => 0.480 + 1.574 * w - 0.176 * w * w },
};
Real roots of a cubic
64-bit
The real roots of z³ + c₂z² + c₁z + c₀ = 0 in ascending order.
In the one-root branch one of the two cube roots subtracts nearly equal numbers when √disc is close to |q| ÷ 2. The acos argument is clamped to [−1, 1]. Each root then gets up to three Newton steps on the original cubic. The branch thresholds of 1e-30 are absolute.
process/eos.js
lines 17 to 34
CadCraft 0.8.2
/** Real roots of z^3 + c2 z^2 + c1 z + c0 = 0, ascending (trigonometric/Cardano, no iteration). */
export function cubicRoots(c2, c1, c0) {
const p = c1 - c2 * c2 / 3, q = 2 * c2 * c2 * c2 / 27 - c2 * c1 / 3 + c0, sh = -c2 / 3;
const disc = q * q / 4 + p * p * p / 27;
let roots;
if (disc > 1e-30) { // one real root
const s = Math.sqrt(disc), u = Math.cbrt(-q / 2 + s), v = Math.cbrt(-q / 2 - s);
roots = [u + v + sh];
} else if (p === 0 && q === 0) roots = [sh];
else {
const m = 2 * Math.sqrt(-p / 3), arg = Math.max(-1, Math.min(1, (3 * q) / (p * m))), th = Math.acos(arg) / 3;
roots = [m * Math.cos(th) + sh, m * Math.cos(th - 2 * Math.PI / 3) + sh, m * Math.cos(th - 4 * Math.PI / 3) + sh].sort((a, b) => a - b);
if (disc > -1e-30 && Math.abs(roots[0] - roots[1]) < 1e-9) roots = [roots[0], roots[2]];
}
// one Newton polish each, so a result is the same on every engine
const f = (z) => ((z + c2) * z + c1) * z + c0, d = (z) => (3 * z + 2 * c2) * z + c1;
return roots.map((z) => { for (let i = 0; i < 3; i++) { const dd = d(z); if (dd === 0) break; const zn = z - f(z) / dd; if (!Number.isFinite(zn)) break; z = zn; } return z; }).sort((a, b) => a - b);
}
Attraction parameter of each component
64-bit
a(T) and its temperature derivative for every pure component.
The derivative is the analytic one, not a difference. Close to the critical temperature 1 − √(T ÷ T_c) subtracts nearly equal numbers, but it is then a small correction to 1.
process/eos.js
lines 46 to 55
CadCraft 0.8.2
/** pure-component a_i(T), da_i/dT */
pure(T) {
const a = [], da = [];
for (let i = 0; i < this.n; i++) {
const c = this.comps[i], sq = 1 + this.kap[i] * (1 - Math.sqrt(T / c.Tc)), al = sq * sq;
a.push(this.ac[i] * al);
da.push(-this.ac[i] * this.kap[i] * sq / Math.sqrt(T * c.Tc));
}
return { a, da };
}
Mixture parameters
64-bit
a, da/dT and b of a mixture of given composition.
a = Σᵢ Σⱼ yᵢ yⱼ √(aᵢ aⱼ) (1 − kᵢⱼ); b = Σ yᵢ bᵢ
A double sum over the components in index order; all terms have the same sign unless a kᵢⱼ exceeds 1.
process/eos.js
lines 56 to 67
CadCraft 0.8.2
mix(T, y) {
const { a, da } = this.pure(T), n = this.n; let am = 0, dam = 0, bm = 0; const aij = Array.from({ length: n }, () => new Array(n));
for (let i = 0; i < n; i++) {
bm += y[i] * this.b[i];
for (let j = 0; j < n; j++) {
const s = Math.sqrt(a[i] * a[j]), k = 1 - this.kij[i][j]; aij[i][j] = s * k;
am += y[i] * y[j] * s * k;
if (s > 0) dam += y[i] * y[j] * k * 0.5 * s * (da[i] / a[i] + da[j] / a[j]);
}
}
return { a: am, da: dam, b: bm, aij };
}
Compressibility roots
64-bit
The physical compressibility factors of a mixture at a temperature and pressure.
β = bP ÷ RT; A = aP ÷ (RT)²; Z³ + c₂Z² + c₁Z + c₀ = 0 with c from β, A, ε, σ; keep Z > β (1 + 1e-12)
The coefficients are formed in one expression each and passed to the cubic solver. A root is kept only when it exceeds β by a relative 1e-12.
process/eos.js
lines 68 to 75
CadCraft 0.8.2
/** All real compressibility roots above beta (physical), ascending, for the composition y. */
roots(T, P, y) {
const mx = this.mix(T, y), beta = mx.b * P / (R * T);
const A = mx.a * P / (R * T * R * T), u = this.m.eps + this.m.sigma, w = this.m.eps * this.m.sigma;
const c2 = u * beta - 1 - beta, c1 = w * beta * beta - (1 + beta) * u * beta + A, c0 = -(1 + beta) * w * beta * beta - A * beta;
const zs = cubicRoots(c2, c1, c0).filter((z) => z > beta * (1 + 1e-12));
return { z: zs, beta, A, mx };
}
Fugacity coefficients
64-bit
ln φ of every component for one compressibility root.
ln φᵢ = bᵢ⁄b (Z − 1) − ln(Z − β) − a ÷ (bRT) × (2 Σⱼ yⱼ aᵢⱼ ÷ a − bᵢ⁄b) × I, I = ln((Z + σβ) ÷ (Z + εβ)) ÷ (σ − ε)
ln(Z − β) loses digits for a dense liquid root where Z is only slightly above β. For Soave–Redlich–Kwong ε is 0 and the logarithm is of 1 + β ÷ Z.
process/eos.js
lines 76 to 85
CadCraft 0.8.2
/** ln of the fugacity coefficients of every component for root Z. */
lnphi(T, P, y, Z, r = null) {
r = r || this.roots(T, P, y); const { beta, mx } = r, s = this.m.sigma, e = this.m.eps;
const I = Math.log((Z + s * beta) / (Z + e * beta)) / (s - e), qv = mx.a / (mx.b * R * T), out = [];
for (let i = 0; i < this.n; i++) {
let sum = 0; for (let j = 0; j < this.n; j++) sum += y[j] * mx.aij[i][j];
out.push(this.b[i] / mx.b * (Z - 1) - Math.log(Z - beta) - qv * (2 * sum / mx.a - this.b[i] / mx.b) * I);
}
return out;
}
Enthalpy and entropy departures
64-bit
How far a real phase's enthalpy and entropy are from the ideal gas at the same temperature and pressure.
H − H_ig = RT (Z − 1) + (T da/dT − a) ÷ b × I; S − S_ig = R ln(Z − β) + (da/dT) ÷ b × I
Closed forms on the same I as the fugacity coefficients. R is the ten-digit constant from data.js.
process/eos.js
lines 86 to 91
CadCraft 0.8.2
/** departures at (T,P) for root Z: { H, S } in J/mol and J/(mol K) (H - Hig, S - Sig at the same T and P). */
departure(T, P, y, Z, r = null) {
r = r || this.roots(T, P, y); const { beta, mx } = r, s = this.m.sigma, e = this.m.eps;
const I = Math.log((Z + s * beta) / (Z + e * beta)) / (s - e);
return { H: R * T * (Z - 1) + (T * mx.da - mx.a) / mx.b * I, S: R * Math.log(Z - beta) + mx.da / mx.b * I };
}
Vapor pressure of a pure component
64-bit
The pressure at which the liquid and vapor roots of one component have equal fugacity.
P ← P × exp(ln φ_L − ln φ_V) until |ln φ_L − ln φ_V| < 1e-12
Successive substitution from Wilson's estimate, at most 100 passes. While only one root exists the pressure is moved by 10 % per pass. If the passes run out the last pressure is returned with no flag.
process/eos.js
lines 96 to 109
CadCraft 0.8.2
/** Pure-component vapour pressure at T (Pa) from equal fugacity of the two roots (i = component index), or null above Tc. */
psat(i, T, guess = null) {
const c = this.comps[i]; if (T >= c.Tc) return null;
const sub = new Cubic([c], this.model), y = [1];
let P = guess || c.Pc * Math.exp(5.373 * (1 + c.omega) * (1 - c.Tc / T));
for (let it = 0; it < 100; it++) {
const r = sub.roots(T, P, y);
if (r.z.length < 2) { P *= r.z.length && sub.liquidLike(r.z[0], r.beta) ? 0.9 : 1.1; continue; }
const zl = r.z[0], zv = r.z[r.z.length - 1], lf = sub.lnphi(T, P, y, zl, r)[0], vf = sub.lnphi(T, P, y, zv, r)[0];
const d = lf - vf; if (Math.abs(d) < 1e-12) return P;
P *= Math.exp(d * (zv > 0 ? 1 : 1)); // phi_L/phi_V > 1 means P is too low
}
return P;
}
In detail
Method
Generic cubic form of Smith, Van Ness and Abbott (eq. 3.41 and Table 3.1, cited in the file) with van der Waals one-fluid mixing and binary interaction coefficients; Cardano and trigonometric roots; vapor pressure by successive substitution on equal fugacity.
Number types
64-bit numbers in plain arrays.
Tolerances
Cubic discriminant: above 1e-30 one real root, otherwise three (absolute). Two lowest roots closer than 1e-9 are merged when the discriminant is within 1e-30 of zero. Roots are kept only above β × (1 + 1e-12). Vapor pressure stops at |ln φ_L − ln φ_V| < 1e-12 or after 100 passes.
Limits
The equation-of-state constants are rounded forms of the model's exact values. psat returns its last pressure without a flag when 100 passes do not reach the tolerance.
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.