Runs the units of a flowsheet in order, converges recycle loops and design specifications, and works out each unit's outlet streams, heat and work.
How the precision is kept
Streams carry 64-bit flows, compositions, temperature, pressure, enthalpy and entropy. A unit operation is a material and energy balance in direct arithmetic followed by a flash from the property package. Recycles are iterated on tear streams until the scaled change is below 1e-9 by default, and design specifications are solved by a bracketed root or Newton; an unconverged loop or specification is reported as a problem.
Calculations under this policy
15 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.
A stream rebuilt from the component flows, temperature and pressure the recycle iteration is working on.
nᵢ = max(xᵢ, 0); F = Σ nᵢ; z = n ÷ F; state from a PT flash
Negative component flows are set to zero before they are summed. A total at or below 1e-300 gives an empty stream at the vector's temperature and pressure.
process/flowsheet.js
lines 88 to 91
CadCraft 0.8.2
const fromVec = (ids, x) => ids.forEach((id, k) => {
const base = k * (nc + 2), n = x.slice(base, base + nc).map((v) => Math.max(v, 0)), F = sum(n);
S[id] = F > 1e-300 ? fromFlash(F, th.flashPT(x[base + nc], x[base + nc + 1], n.map((v) => v / F)), id) : { ...emptyStream(th, id), T: x[base + nc], P: x[base + nc + 1] };
});
Recycle convergence test
64-bit
Whether the tear streams have stopped changing from one pass to the next.
The change of every component flow, temperature and pressure is scaled by its own size plus a floor, so trace components are judged against 1e-6 of the largest feed. The default tolerance is 1e-9 and the limit 300 passes.
process/flowsheet.js
lines 92 to 99
CadCraft 0.8.2
const scale = tears.flatMap(() => [...new Array(nc).fill(1e-6 * Fscale), 1, 1e3]);
// initial values: the feeds' T and P when known, empty flow
let x = toVec(tears); fromVec(tears, x);
const G = (xv) => { fromVec(tears, xv); g.forEach(exec); return toVec(tears); };
let xPrev = null, gPrev = null, err = Infinity, it = 0, H = null;
for (it = 1; it <= maxIter; it++) {
const gx = G(x); err = 0; for (let i = 0; i < x.length; i++) err = Math.max(err, Math.abs(gx[i] - x[i]) / (Math.abs(x[i]) + scale[i]));
if (err < tol) { x = gx; break; }
Wegstein acceleration
64-bit
The next guess of each tear variable from its last two passes.
s = (g − g_prev) ÷ (x − x_prev); q = clamp(s ÷ (s − 1), −5, 0); x_next = q x + (1 − q) g
Applied to each variable on its own. A change in x below 1e-300 falls back to direct substitution. s close to 1 would divide by nearly zero, but the clamp on q bounds the step.
process/flowsheet.js
lines 102 to 104
CadCraft 0.8.2
else if (method === 'wegstein') {
xn = x.map((xi, i) => { const dx = xi - xPrev[i]; if (Math.abs(dx) < 1e-300) return gx[i]; const s = (gx[i] - gPrev[i]) / dx, q = Math.max(-5, Math.min(0, s / (s - 1))); return q * xi + (1 - q) * gx[i]; });
} else { // Broyden's good method on F = x - G(x)
Broyden acceleration
64-bit
The next guess of all tear variables together, from an updated inverse Jacobian.
H ← H + (Δx − H ΔF) (Δxᵀ H) ÷ (Δxᵀ H ΔF); x_next = x − H F, F = x − g(x)
The rank-one update is skipped when its denominator is below 1e-300 in absolute value. The variables are not scaled, so flows, temperatures and pressures enter the update with their own magnitudes.
process/flowsheet.js
lines 105 to 113
CadCraft 0.8.2
const Fx = x.map((v, i) => v - gx[i]);
if (H == null) { H = x.map((_, i) => x.map((__, j) => (i === j ? 1 : 0))); }
else {
const dx = x.map((v, i) => v - xPrev[i]), dF = Fx.map((v, i) => v - (xPrev[i] - gPrev[i]));
const Hd = H.map((row) => row.reduce((s, h, j) => s + h * dF[j], 0)), dxH = x.map((_, j) => dx.reduce((s, d, i) => s + d * H[i][j], 0)), den = dx.reduce((s, d, i) => s + d * Hd[i], 0);
if (Math.abs(den) > 1e-300) for (let i = 0; i < x.length; i++) for (let j = 0; j < x.length; j++) H[i][j] += (dx[i] - Hd[i]) * dxH[j] / den;
}
const step = H.map((row) => row.reduce((s, h, j) => s + h * Fx[j], 0));
xn = x.map((v, i) => v - step[i]);
Design specifications
64-bit
The value of a unit parameter that makes a stream property hit its target.
residual = (value − target) ÷ (|target| + 1e-6); one bounded spec: scan 48 intervals, then a bracketed root to 1e-12; otherwise Newton to 1e-10
Every residual evaluation solves the whole flowsheet, so the specification can be no closer than the recycle tolerance inside it. An unmet specification marks the result as not converged.
process/flowsheet.js
lines 137 to 161
CadCraft 0.8.2
/** Solve with design specs; each spec varies one unit parameter. Returns the usual result plus `specs: [{ vary, value, residual }]`. */
export function solve(def) {
if (!def.specs || !def.specs.length) return solveFlowsheet(def);
const th = thermoFor(def), sp = def.specs;
const ov = (xv) => { const o = {}; sp.forEach((s, k) => { (o[s.vary.unit] ||= {})[s.vary.param] = xv[k]; }); return o; };
const x0 = sp.map((s) => s.vary.x0 ?? (def.units.find((u) => u.id === s.vary.unit)?.p?.[s.vary.param]) ?? 0.5 * ((s.vary.min ?? 0) + (s.vary.max ?? 1)));
let last = null;
const F = (xv) => { last = solveFlowsheet(def, ov(xv)); return sp.map((s) => { const v = propOf(th, last, s.target); return (v - s.target.value) / (Math.abs(s.target.value) + 1e-6); }); };
const feas = (xv) => sp.every((s, k) => (s.vary.min == null || xv[k] >= s.vary.min) && (s.vary.max == null || xv[k] <= s.vary.max));
let sol = null;
if (sp.length === 1 && sp[0].vary.min != null && sp[0].vary.max != null) { // one spec with limits: bracket by a scan (flat regions defeat Newton), then Illinois
const { min, max } = sp[0].vary, N = 48; let prev = null;
for (let k = 0; k <= N && !sol; k++) {
const v = min + (max - min) * k / N, f = F([v])[0];
if (f === 0) sol = { x: [v], ok: true };
else if (prev && prev.f * f < 0) { const r = illinois((t) => F([t])[0], prev.v, v, { tol: 1e-12 }); sol = { x: [r], ok: r != null }; }
prev = { v, f };
}
}
if (!sol) sol = newtonFD(F, x0, { tol: 1e-10, maxIter: 60, feasible: feas, rel: 1e-6 });
const final = solveFlowsheet(def, ov(sol.x));
final.specs = sp.map((s, k) => ({ unit: s.vary.unit, param: s.vary.param, value: sol.x[k], target: s.target, achieved: propOf(th, final, s.target) }));
if (!sol.ok) { final.converged = false; final.problems.push(problem('spec.unconverged', 'error', 'A design specification could not be met within the limits on the varied parameter.', sp.map((s) => ({ kind: 'unit', id: s.vary.unit })))); }
return final;
}
Property a specification reads
64-bit
The stream or unit value a design specification is compared with.
flow of a component = F × zᵢ
Values are read from the solved result without rounding; a component flow is one multiplication.
process/flowsheet.js
lines 127 to 135
CadCraft 0.8.2
export function propOf(th, res, t) {
const s = res.streams[t.stream], p = t.prop;
if (p === 'T' || p === 'P' || p === 'F' || p === 'vf' || p === 'H') return s ? s[p] : NaN;
if (p.startsWith('x:')) return s.z[th.index(p.slice(2))];
if (p.startsWith('flow:')) return s.F * s.z[th.index(p.slice(5))];
if (p.startsWith('duty:')) return (res.units[p.slice(5)] || {}).Q;
if (p.startsWith('work:')) return (res.units[p.slice(5)] || {}).W;
throw new Error('unknown spec property ' + p);
}
A feed's total flow and mole fractions from component flows, mass flows or a total with fractions.
F = Σ nᵢ; zᵢ = nᵢ ÷ F; from mass: nᵢ = mᵢ ÷ Mᵢ
The flows are added in component order and each is divided by the sum. Given fractions are divided by their own sum, so they need not add to exactly 1.
process/stream.js
lines 20 to 28
CadCraft 0.8.2
/** A feed from flows { id: mol/s } (or F and z) at T and P. */
export function feedStream(th, spec, id = '') {
let F, z;
if (spec.flows) { z = th.ids.map((c) => spec.flows[c] || 0); F = sum(z); z = z.map((v) => v / F); }
else if (spec.massFlows) { const m = th.ids.map((c, i) => (spec.massFlows[c] || 0) / th.c[i].MWkg); F = sum(m); z = m.map((v) => v / F); }
else { F = spec.F; z = spec.z.slice(); const s = sum(z); z = z.map((v) => v / s); }
const fl = spec.vf != null ? th.flashPVF(spec.P, spec.vf, z) : th.flashPT(spec.T, spec.P, z);
return fromFlash(F, fl, id);
}
Mass flow of a stream
64-bit
The mass flow in kg/s from the molar flow and composition.
ṁ = F × Σ zᵢ Mᵢ
One sum of products and one multiplication. The molar masses are the table's values of five significant digits.
Each phase's molar volume is evaluated again from the property package at the stream's temperature and pressure and the terms are added.
process/stream.js
line 32
CadCraft 0.8.2
export const volFlow = (th, s) => { let v = 0; for (const p of s.phases) { const ph = th.phase(s.T, s.P, p.x, p.type === 'vapor' ? 'V' : 'L'); v += s.F * p.frac * ph.V; } return v; }; // m3/s
Mixer
64-bit
The outlet of several streams mixed at the lowest inlet pressure.
nᵢ = Σ F zᵢ; H_flow = Σ F H; outlet from a flash at that enthalpy
Component flows and enthalpy flows are added stream by stream in inlet order. The outlet temperature comes from the enthalpy flash, started from the flow-weighted mean temperature.
process/stream.js
lines 42 to 49
CadCraft 0.8.2
mixer(th, ins, p = {}) {
const live = ins.filter((s) => !isZero(s)), n = th.n;
if (!live.length) return { outs: [emptyStream(th)], Q: 0, W: 0, notes: [] };
const nf = new Array(n).fill(0); let Hf = 0, Tg = 0, Fsum = 0;
for (const s of live) { for (let i = 0; i < n; i++) nf[i] += s.F * s.z[i]; Hf += s.F * s.H; Tg += s.F * s.T; Fsum += s.F; }
const F = sum(nf), z = nf.map((v) => v / F), P = p.P != null ? p.P : Math.min(...live.map((s) => s.P));
return { outs: [atEnthalpy(th, F, z, P, Hf, Tg / Fsum)], Q: 0, W: 0, notes: [] };
},
Splitter
64-bit
Outlets that share one stream by given fractions.
F_k = F × f_k ÷ Σ f
Every fraction is divided by the sum of the fractions, so the outlets add up to the inlet to within rounding whatever was typed. A sum more than 1e-9 from 1 adds a note.
process/stream.js
lines 50 to 54
CadCraft 0.8.2
splitter(th, [s], p) {
const fr = p.fractions; if (!fr) throw new Error('splitter needs fractions');
const tot = sum(fr), outs = fr.map((f) => (isZero(s) ? emptyStream(th) : { ...s, F: s.F * f / tot, z: s.z.slice(), phases: s.phases.map((q) => ({ ...q, x: q.x.slice() })) }));
return { outs, Q: 0, W: 0, notes: Math.abs(tot - 1) > 1e-9 ? ['split fractions summed to ' + tot + ' and were normalised'] : [] };
},
Pump
64-bit
The shaft work and outlet state of a pump.
w = V_liquid (P_out − P_in) ÷ η; W = F w; outlet at H_in + w
The liquid molar volume is taken at the inlet and held constant over the pressure rise. The efficiency is 0.75 unless given.
process/stream.js
lines 65 to 73
CadCraft 0.8.2
pump(th, [s], p) {
if (isZero(s)) return { outs: [emptyStream(th)], Q: 0, W: 0, notes: [] };
const eff = p.eff || 0.75, notes = [];
if (s.vf > 0) notes.push('pump inlet has vapour (vapour fraction ' + s.vf.toFixed(3) + ')');
let v = 0; for (const q of s.phases) v += q.frac * th.phase(s.T, s.P, q.x, 'L').V;
const w = v * (p.Pout - s.P) / eff; // J/mol
const out = atEnthalpy(th, s.F, s.z, p.Pout, s.F * (s.H + w), s.T);
return { outs: [out], Q: 0, W: s.F * w, notes, eff, idealW: s.F * v * (p.Pout - s.P) };
},
Compressor and expander
64-bit
The work and outlet state from an isentropic step and an efficiency.
h_s from a flash at P_out and the inlet entropy; compressor: h = H_in + (h_s − H_in) ÷ η; expander: h = H_in + η (h_s − H_in)
The isentropic enthalpy comes from the entropy flash, whose temperature is found to its own tolerance. The work is the enthalpy difference times the flow.
Two outlets that each take a set fraction of every component.
aᵢ = F zᵢ fᵢ; bᵢ = F zᵢ (1 − fᵢ); Q = F_a H_a + F_b H_b − F H
The two outlet flows of a component are formed separately from f and 1 − f, so their sum equals the inlet to within one rounding. The duty is the difference of enthalpy flows.
process/stream.js
lines 100 to 108
CadCraft 0.8.2
/** Component splitter: split { id: fraction to outlet 1 }, the rest to outlet 2; both leave at T, P (or Tout1/Tout2, P). */
compsplit(th, [s], p) {
if (isZero(s)) return { outs: [emptyStream(th), emptyStream(th)], Q: 0, W: 0, notes: [] };
const n = th.n, a = new Array(n).fill(0), b = new Array(n).fill(0);
for (let i = 0; i < n; i++) { const f = p.split[th.ids[i]] ?? 0, m = s.F * s.z[i]; a[i] = m * f; b[i] = m * (1 - f); }
const mk = (m, T) => { const F = sum(m); if (F <= 0) return emptyStream(th); const z = m.map((v) => v / F); return fromFlash(F, th.flashPT(T, p.P != null ? p.P : s.P, z)); };
const o1 = mk(a, p.T1 != null ? p.T1 : s.T), o2 = mk(b, p.T2 != null ? p.T2 : s.T);
return { outs: [o1, o2], Q: o1.F * o1.H + o2.F * o2.H - s.F * s.H, W: 0, notes: [] };
},
Two-stream heat exchanger
64-bit
The outlet states, duty, log-mean temperature difference, UA and smallest approach of a counter-current exchanger.
duty = F_hot (H_in − H_out); LMTD = (ΔT₁ − ΔT₂) ÷ ln(ΔT₁ ÷ ΔT₂), or ΔT₁ when |ΔT₁ − ΔT₂| < 1e-9; UA = duty ÷ LMTD
Just outside the 1e-9 K window the log-mean formula divides a small difference by a small logarithm and loses digits. The smallest approach is the minimum over 41 sampled points along the duty, so a pinch between samples can be missed. When the two end differences have opposite signs the logarithm has a negative argument and the LMTD is not a number; UA is then left empty.
process/stream.js
lines 109 to 128
CadCraft 0.8.2
/** Two-stream heat exchanger. ins [hot, cold]; spec one of { hotOutT, coldOutT, duty, approach (hot out - cold in... min) } plus dPhot, dPcold. */
hx(th, [hot, cold], p) {
if (isZero(hot) || isZero(cold)) return { outs: [hot, cold].map((s) => (isZero(s) ? emptyStream(th) : s)), Q: 0, W: 0, notes: [] };
const Ph = hot.P - (p.dPhot || 0), Pc = cold.P - (p.dPcold || 0);
let ho, co, duty;
if (p.hotOutT != null) { ho = fromFlash(hot.F, th.flashPT(p.hotOutT, Ph, hot.z)); duty = hot.F * (hot.H - ho.H); co = atEnthalpy(th, cold.F, cold.z, Pc, cold.F * cold.H + duty, cold.T); }
else if (p.coldOutT != null) { co = fromFlash(cold.F, th.flashPT(p.coldOutT, Pc, cold.z)); duty = cold.F * (co.H - cold.H); ho = atEnthalpy(th, hot.F, hot.z, Ph, hot.F * hot.H - duty, hot.T); }
else if (p.duty != null) { duty = p.duty; ho = atEnthalpy(th, hot.F, hot.z, Ph, hot.F * hot.H - duty, hot.T); co = atEnthalpy(th, cold.F, cold.z, Pc, cold.F * cold.H + duty, cold.T); }
else throw new Error('heat exchanger needs hotOutT, coldOutT or duty');
// profile check: temperatures along the duty, counter-current
const N = 40, prof = []; let minDT = Infinity;
for (let k = 0; k <= N; k++) {
const q = duty * k / N, h = atEnthalpy(th, hot.F, hot.z, Ph, hot.F * hot.H - duty + q, ho.T), c = atEnthalpy(th, cold.F, cold.z, Pc, cold.F * cold.H + q, cold.T);
prof.push([q, h.T, c.T]); minDT = Math.min(minDT, h.T - c.T);
}
let lm = 0, cnt = 0; const dT1 = ho.T - cold.T, dT2 = hot.T - co.T;
lm = Math.abs(dT1 - dT2) < 1e-9 ? dT1 : (dT1 - dT2) / Math.log(dT1 / dT2);
const notes = []; if (minDT < 0) notes.push('temperature cross: the hot side is colder than the cold side (minimum difference ' + minDT.toFixed(2) + ' K)');
return { outs: [ho, co], Q: 0, W: 0, notes, duty, lmtd: lm, minApproach: minDT, UA: lm > 0 ? duty / lm : null, profile: prof, internalQ: duty };
},
In detail
Method
Sequential-modular solution with Tarjan's strongly connected components; tear streams accelerated by bounded Wegstein or Broyden's method; design specifications by a 48-point scan and the Illinois root finder, or a finite-difference Newton.
Number types
64-bit numbers in plain arrays and objects.
Tolerances
Recycle: max |g(x) − x| ÷ (|x| + scale) < 1e-9, with scale 1e-6 × the largest feed flow for component flows, 1 K for temperature and 1000 Pa for pressure; 300 passes. Wegstein factor q kept in [−5, 0]. Specification: residual (value − target) ÷ (|target| + 1e-6); Illinois to 1e-12 or Newton to 1e-10 with a step of 1e-6. Split fractions that do not sum to 1 within 1e-9 are normalized with a note.
Limits
The heat exchanger's temperature profile is sampled at 41 points. Pump work uses the inlet liquid volume as constant over the pressure rise.
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.