Resources Numerical accuracy Calculation policy
Works out the orifice area a relief valve needs for gas, liquid, steam or two-phase flow, the heat input from a fire, the back pressure in the outlet line and flare header, flare radiation distance and the pressure of a vessel as it is vented.
64-bit arithmetic. The orifice equations are closed forms; several results come from iterations with their own stop tests (viscosity correction, critical pressure ratio, outlet-line pressure) and the depressuring curve from a fixed-step integration. US-unit equations are converted with constants of eight to ten digits or exact ones.
Open item. Interpolate the time at which the depressuring curve crosses the target pressure, and take the minimum temperature over every step instead of every tenth.
17 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-equip.mjs, tools/process-test/test-dec.mjs
Converts the pool-fire coefficient from Btu/h per ft^1.64 to W per m^1.64.
C_SI = C_US × 0.29307107 × (1 ÷ 0.3048²)^0.82
The Btu/h to watt factor is written to eight significant digits (0.29307107); the foot is the exact 0.3048. The conversion is computed once when the file loads.
const FIRE_C = (c) => c * BTU_H * Math.pow(1 / (FT * FT), 0.82); // W per m^1.64The orifice area a relief valve needs for a gas or vapor flow.
A = W ÷ (Kd·Kb·Kc × G); critical: G = P₁·√(k·M ÷ Z·R·T) × (2 ÷ (k + 1))^((k+1)/(2(k−1)))
64-bit closed forms. Flow is treated as critical unless the back-pressure ratio is above (2 ÷ (k + 1))^(k/(k−1)); with no back pressure given it is always critical. In the subcritical branch the factor 1 − r^((k−1)/k) loses digits as the ratio approaches 1.
/** Gas/vapour required area (m²). W kg/s, P1 relieving pressure Pa absolute (set + overpressure + atmosphere), P2 total back pressure Pa abs (default: critical flow assumed). */
export function gasArea({ W, k, M, T, Z = 1, P1, P2 = 0, Kd = 0.975, Kb = 1, Kc = 1 }) {
const rc = Math.pow(2 / (k + 1), k / (k - 1)), r = P2 / P1, critical = !(r > rc), notes = [];
let flux;
if (critical) flux = P1 * Math.sqrt(k * M / (Z * R * T)) * Math.pow(2 / (k + 1), (k + 1) / (2 * (k - 1)));
else { flux = P1 * Math.sqrt(M / (Z * R * T)) * Math.sqrt(2 * k / (k - 1) * Math.pow(r, 2 / k) * (1 - Math.pow(r, (k - 1) / k))); notes.push('subcritical flow: the sizing used the isentropic expansion; check the valve maker\'s subcritical method (conventional valves are often limited)'); }
const A = W / (Kd * Kb * Kc * flux); return { A, flux, critical, rc, notes, source: SOURCES.gas };
}The orifice area for a liquid flow, with the viscosity correction.
A = Q ÷ (Kd·Kw·Kc·Kv) × √(ρ ÷ 2ΔP); Kv = 1 ÷ (0.9935 + 2.878 ÷ √Re + 342.75 ÷ Re^1.5)
The viscosity correction is a fixed-point loop on the area: at most 40 passes, stopped when the area changes by less than 1e-12 relative. Nothing reports a loop that ends at the 40th pass.
/** Liquid relief area (m²). Q m³/s, dP = set pressure + overpressure − back pressure (Pa). Iterates the viscosity correction. */
export function liquidArea({ Q, rho, mu = null, dP, Kd = 0.65, Kw = 1, Kc = 1 }) {
let A = Q / (Kd * Kw * Kc) * Math.sqrt(rho / (2 * dP)), Kv = 1, Re = null;
if (mu != null) for (let i = 0; i < 40; i++) { const Do = Math.sqrt(4 * A / Math.PI); Re = 4 * rho * Q / (Math.PI * mu * Do); Kv = 1 / (0.9935 + 2.878 / Math.sqrt(Re) + 342.75 / Math.pow(Re, 1.5)); const An = Q / (Kd * Kw * Kc * Kv) * Math.sqrt(rho / (2 * dP)); if (Math.abs(An - A) < 1e-12 * A) { A = An; break; } A = An; }
return { A, Kv, Re, source: SOURCES.liquid };
}The pressure ratio at which a flashing two-phase flow chokes.
η² + (ω² − 2ω)(1 − η)² + 2ω²·ln η + 2ω²(1 − η) = 0
Bisection on η between 1e-9 and 1 − 1e-12, 200 halvings, which is past 64-bit resolution. It returns null when the function has the same sign at both ends.
/** Critical pressure ratio ηc of the (flashing) omega method: ηc² + (ω² − 2ω)(1 − ηc)² + 2ω² ln ηc + 2ω²(1 − ηc) = 0. */
export function omegaCritical(omega) {
const f = (e) => e * e + (omega * omega - 2 * omega) * (1 - e) ** 2 + 2 * omega * omega * Math.log(e) + 2 * omega * omega * (1 - e);
let lo = 1e-9, hi = 1 - 1e-12; if (f(lo) * f(hi) > 0) return null;
for (let i = 0; i < 200; i++) { const m = (lo + hi) / 2; if (f(lo) * f(m) <= 0) hi = m; else lo = m; } return (lo + hi) / 2;
}The mass flux of a flashing mixture through a nozzle.
G = √(P₀·ρ₀) × √(2[ω·ln(1 ÷ η) − (ω − 1)(1 − η)]) ÷ (ω(1 ÷ η − 1) + 1)
Closed form at η = max(Pa ÷ P₀, η_c). The bracket under the root is a difference of two terms that are both small when η is close to 1.
/** Mass flux (kg/m²/s) of a flashing two-phase mixture through a nozzle: P0 Pa, rho0 kg/m³ at P0, back pressure Pa. */
export function twoPhaseFlux({ P0, rho0, omega, Pa }) {
const ec = omegaCritical(omega), eta = Math.max(Pa / P0, ec ?? 0), choked = ec != null && Pa / P0 <= ec;
const G = Math.sqrt(P0 * rho0) * Math.sqrt(2 * (omega * Math.log(1 / eta) - (omega - 1) * (1 - eta))) / (omega * (1 / eta - 1) + 1);
return { G, etaC: ec, eta, choked, source: SOURCES.omega };
}The heat a wetted vessel absorbs in a pool fire.
Q = C × F × A^0.82
One Math.pow in 64-bit. A note is added above 260 m² wetted area, where other coefficients apply and are not used.
/** Heat absorbed by a wetted vessel in a pool fire: W. F the environment factor (1 bare vessel; <1 with credit-worthy insulation). */
export function fireHeat({ A, F = 1, drainage = 'adequate' }) {
const c = FIRE_COEFF[drainage]; if (!c) throw new Error("drainage: 'adequate' (drainage and prompt fire fighting) or 'poor'");
const notes = []; if (A > 260) notes.push('wetted area above about 260 m²: API 521 gives other coefficients for large areas (not applied)');
return { Q: c * F * Math.pow(A, 0.82), coeff: c, notes, source: SOURCES.fire };
}The vessel surface in contact with liquid up to the fire height.
vertical: π·D·level + k·D²; horizontal: r·θ·L + 2k·D²·(h ÷ D), θ = 2·acos(1 − h ÷ r)
The level is capped at the fire height (7.6 m by default) and at the diameter but is not floored at 0, so a negative level gives acos an argument above 1 and the result is NaN. The heads of a horizontal vessel are counted in proportion to the level, which the code marks as an approximation; the head factor 1.084 is a rounded constant.
/** Wetted surface (m²) of a vessel to a given height above the fire source (default 7.6 m = 25 ft), liquid level included. cyl: { type: 'vertical'|'horizontal', D, L (tan-tan), liquidLevel (m from the bottom tangent), headFactor (area/D², 2:1 ellipsoidal 1.084) } */
export function wettedArea({ type, D, L, liquidLevel, headFactor = 1.084, fireHeight = 7.6 }) {
const top = Math.min(liquidLevel, fireHeight);
if (type === 'vertical') { const cyl = Math.PI * D * Math.max(0, Math.min(top, L)); return cyl + headFactor * D * D; }
const r = D / 2, h = Math.min(top, D), th = 2 * Math.acos(1 - h / r); // horizontal: segment of the circumference wetted
return r * th * L + 2 * headFactor * D * D * (h / D); // heads wetted in proportion to the level (approximation)
}The flow through broken tube ends from the high-pressure side to the low-pressure side.
Q = Cd × A × √(2(P_h − P_l) ÷ ρ), A = n × π·d² ÷ 4
One 64-bit expression; the pressure difference is formed first.
/** Tube rupture: liquid or gas through two tube ends (double-ended), orifice flow with discharge coefficient Cd. High pressure side Ph, low Pl (Pa). → m³/s (liquid) */
export function tubeRupture({ di, nEnds = 2, Ph, Pl, rho, Cd = 0.61 }) { const A = nEnds * Math.PI * di * di / 4; return { Q: Cd * A * Math.sqrt(2 * (Ph - Pl) / rho), A }; }The smallest standard orifice whose area covers the required area.
in² = m² ÷ 0.0254²
The required area is converted with the exact inch and compared with a table of standard areas in in²; the first one at least as large is returned, converted back to m².
export function orificeLetter(areaM2) { const a = areaM2 / (0.0254 * 0.0254); for (const [k, v] of Object.entries(ORIFICES)) if (v >= a) return { letter: k, area: v * 0.0254 * 0.0254, source: SOURCES.orifices }; return { letter: null, area: null, notes: ['larger than orifice T: use more than one valve'] }; }The distance from a flare at which the heat flux falls to an allowable level.
D = √(τ·F·Q ÷ (4π·K))
One 64-bit expression. The radiant fraction defaults to 0.2 and the allowable flux to 4730 W/m².
/** Radiation distance: D = sqrt(tau F Q / (4 pi K)); Q heat release W, F radiant fraction (0.1–0.4 by gas; the caller's), tau transmissivity (≤ 1), K allowable W/m². */
export const flareRadiation = ({ Q, tau = 1, F = 0.2, K = 4730 }) => ({ D: Math.sqrt(tau * F * Q / (4 * Math.PI * K)), source: SOURCES.radiation });The back pressure at each relief valve as a fraction of its set pressure.
fraction = (P_back − 101325) ÷ P_set
The pressures come from the gas network solver. Atmospheric pressure is the constant 101325 Pa, not an input. The stored fraction is unrounded; toFixed appears only in the message.
/**
* Flare header: a gas network where the relief valves put mass flow in (demand negative) and the flare tip is a fixed-pressure node.
* net = same as solveNetwork with fluid.kind 'gas'. Returns the network solution plus, for each node named in `valves`, the back pressure and its fraction of the set pressure.
*/
export function flareHeader(net, valves = []) {
const r = solveNetwork(net), problems = r.problems.slice();
const bp = valves.map((v) => { const Pb = r.nodes[v.node].P, f = (Pb - 101325) / v.set; if (f > (v.limit ?? 0.1)) problems.push({ id: 'process:relief.backpressure:' + v.node, code: 'relief.backpressure', severity: 'error', text: `Valve at ${v.node}: back pressure ${(Pb / 1e3).toFixed(0)} kPa abs is ${(100 * f).toFixed(0)} % of the set pressure (limit ${(100 * (v.limit ?? 0.1)).toFixed(0)} % for a conventional valve).`, about: [{ id: v.node }], source: { kind: 'engine', name: 'process' } }); return { node: v.node, Pback: Pb, fraction: f }; });
return { ...r, backPressures: bp, problems };
}The orifice area for steam relief.
A[in²] = W[lb/h] ÷ (51.5 × P₁[psia] × Kd·Kb·Kc·Kn·Ksh)
kg/s to lb/h divides by the exact 0.45359237, Pa to psi divides by 6894.757293 (ten digits) and in² to m² multiplies by the exact 0.00064516. The Napier correction is applied above 1500 psia in the code, while the comment above the function says 1515 psia.
/**
* Steam relief area (m²), API 520 Part I steam equation in US units converted: A[in²] = W[lb/h] / (51.5 P1[psia] Kd Kb Kc Kn Ksh).
* Kn corrects the Napier equation above 1515 psia: Kn = (0.1906 P1 − 1000) / (0.2292 P1 − 1061) up to 3200 psia (typed without the standard open: check). Ksh, the superheat correction,
* is read from the standard's table by the person (1 for saturated steam). W kg/s, P1 Pa absolute relieving pressure.
*/
export function steamArea({ W, P1, Kd = 0.975, Kb = 1, Kc = 1, Ksh = 1 }) {
const Wl = W * 3600 / LB, P = P1 / PSI, notes = [];
let Kn = 1;
if (P > 1500) { if (P > 3200) throw new Error('steam above 3200 psia (22 MPa) is outside the equation'); Kn = (0.1906 * P - 1000) / (0.2292 * P - 1061); notes.push('Napier correction Kn = ' + Kn.toFixed(4) + ' applied above 1500 psia'); }
if (Ksh === 1) notes.push('Ksh = 1: right for saturated steam; for superheated steam read Ksh from the standard\'s table and give it');
const Ain = Wl / (51.5 * P * Kd * Kb * Kc * Kn * Ksh);
return { A: Ain * IN2, Kn, Ksh, notes, source: 'API 520 Part I steam: A = W / (51.5 P1 Kd Kb Kc Kn Ksh), US units converted' };
}The thrust on a valve outlet discharging a sonic gas jet to the atmosphere.
F = W·v + (P_e − P_a)·A_e; v = √(2k·Z·R·T ÷ ((k + 1)·M))
64-bit closed form; atmospheric pressure defaults to 101325 Pa.
/** Reaction force (N) on the outlet of a valve discharging a sonic gas jet to the atmosphere: F = W v + (Pe − Pa) A, v the sonic velocity at the exit temperature 2T/(k+1). */
export function reactionForce({ W, k, M, T, Z = 1, P1, Ae, Pa = 101325 }) {
const v = Math.sqrt(2 * k * Z * R * T / ((k + 1) * M)), Pe = P1 * Math.pow(2 / (k + 1), k / (k - 1)), F = W * v + (Pe - Pa) * Ae;
return { F, v, Pe, notes: ['open discharge to atmosphere with sonic flow at the exit area; a closed discharge header has a different thrust (API 520 Part II, 4.2.2)'], source: 'momentum balance on the discharge: F = W v + (Pe − Pa) Ae' };
}The pressure at the valve outlet for a gas flow through the discharge pipe, and the exit Mach number.
P₁² − P₂² = (G²·Z·R·T ÷ M) × (f·L ÷ D + K + 2·ln(P₁ ÷ P₂))
P₁ is found by repeating the formula with the last value: at most 200 passes, stopped when two successive values differ by less than 1e-6 relative. Nothing reports a run that ends at 200 passes. The set-pressure fraction subtracts the constant 101325 Pa.
/**
* Outlet line back pressure: gas of mass flow W through a pipe (inside diameter D m, length L m, fittings K) to a downstream pressure Pdown (Pa abs), isothermal compressible flow:
* P1² − P2² = (G² Z R T / M) (f L/D + K + 2 ln(P1/P2)). Returns the pressure at the valve outlet, its fraction of the set pressure, and the exit Mach number.
*/
export function outletLine({ W, k = 1.4, M, T, Z = 1, D, L, K = 0, Pdown, Pset = null, rough = 4.6e-5, mu = 1.5e-5 }) {
const A = Math.PI * D * D / 4, G = W / A, Re = G * D / mu, f = friction(Re, rough / D);
let P1 = Pdown * 1.05, notes = [];
for (let i = 0; i < 200; i++) {
const rhs = G * G * Z * R * T / M * (f * L / D + K + 2 * Math.log(P1 / Pdown)), Pn = Math.sqrt(Pdown * Pdown + rhs);
if (Math.abs(Pn - P1) < 1e-6 * P1) { P1 = Pn; break; } P1 = Pn;
}
const rho2 = Pdown * M / (Z * R * T), v2 = G / rho2, c = Math.sqrt(k * Z * R * T / M), mach = v2 / c;
if (mach >= 1) notes.push('the exit would be sonic: the downstream pressure is not what the pipe can deliver (choked outlet); enlarge the pipe or lower the flow');
const frac = Pset ? (P1 - 101325) / Pset : null;
return { P1, mach, f, Re, v2, fractionOfSet: frac, ok: mach < 1 && (frac == null || frac <= 0.1), notes, source: 'isothermal compressible pipe flow (Darcy), friction by Colebrook as hydraulics.js' };
}Tests tools/process-test/test-chem4.mjs
Pressure, temperature and vent flow of a gas-filled vessel over time as it is vented, and the time to reach a target pressure.
dm/dt = −ṁ; m·cv·dT/dt = Q − ṁ·Rs·T; P = m·Rs·T ÷ V
Classical fourth-order Runge-Kutta with a fixed step of tEnd ÷ 20000 (at least 1e-4 s); the step is not adapted. Pressure and temperature are stored every tenth step and the minimum temperature is taken from those stored samples. The time to reach the target is the end of the first step at or below it, not interpolated. Time is a running sum of the step.
export function depressuring(spec) {
const { V, P0, T0, M, k, Cd = 0.8, A, Pb = 101325, tEnd, Ptarget = null } = spec, Rs = R / M, cv = Rs / (k - 1);
const Qf = typeof spec.Q === 'function' ? spec.Q : () => spec.Q || 0;
let m = P0 * V / (Rs * T0), Tg = T0, t = 0; const out = { t: [0], P: [P0], T: [T0], mdot: [], tReach: null }, crit = Math.pow(2 / (k + 1), (k + 1) / (2 * (k - 1)));
const flow = (P, T) => {
const r = Pb / P; if (r >= 1) return 0;
if (r <= Math.pow(2 / (k + 1), k / (k - 1))) return Cd * A * P * Math.sqrt(k * M / (R * T)) * crit;
return Cd * A * P * Math.sqrt(2 * M / (R * T) * k / (k - 1) * (Math.pow(r, 2 / k) - Math.pow(r, (k + 1) / k)));
};
const f = (mm, T, tt) => { const P = mm * Rs * T / V, md = flow(P, T); return [-md, (Qf(tt) - md * Rs * T) / (mm * cv), md]; };
const dt0 = Math.max(tEnd / 20000, 1e-4); let dt = dt0, nstep = 0;
while (t < tEnd) {
const [a1, b1, c1] = f(m, Tg, t), [a2, b2] = f(m + dt / 2 * a1, Tg + dt / 2 * b1, t + dt / 2), [a3, b3] = f(m + dt / 2 * a2, Tg + dt / 2 * b2, t + dt / 2), [a4, b4] = f(m + dt * a3, Tg + dt * b3, t + dt);
void c1; m += dt / 6 * (a1 + 2 * a2 + 2 * a3 + a4); Tg += dt / 6 * (b1 + 2 * b2 + 2 * b3 + b4); t += dt;
const P = m * Rs * Tg / V; out.mdot.push(-a1);
if (++nstep % 10 === 0 || t >= tEnd) { out.t.push(t); out.P.push(P); out.T.push(Tg); }
if (Ptarget != null && out.tReach == null && P <= Ptarget) out.tReach = t;
if (P <= Pb * 1.0001) break;
}
const problems = [];
if (Ptarget != null && out.tReach == null) problems.push(problem('depressuring.target', 'warning', `The vessel does not reach ${(Ptarget / 1e5).toFixed(2)} bar within ${tEnd} s.`, [{ id: 'depressuring' }]));
return { ...out, Tmin: Math.min(...out.T), mdotInitial: out.mdot[0], problems, notes: ['ideal gas, Z = 1, no wall heat capacity: the real metal temperature and the two-phase case need the maker/engineer'] };
}The relief load when a control valve fails open, less the normal flow.
liquid: Q[m³/h] = Kv × √(ΔP[bar] ÷ SG), ṁ = Q ÷ 3600 × SG × ρ_w; load = max(flow − normal, 0)
64-bit closed forms. The water density defaults to 1000 kg/m³. The gas branch is the critical-flow expression at the relieving pressure.
export function controlValveFailure(spec) {
const problems = [];
let flow, unit;
if (spec.fluid === 'liquid') {
const { Kv, dP, SG, rho = 1000 } = spec; if (!(Kv > 0 && dP > 0 && SG > 0)) throw new Error('liquid: give Kv, dP (bar) and the specific gravity');
const Qm3h = Kv * Math.sqrt(dP / SG); flow = Qm3h / 3600 * SG * rho; unit = 'kg/s';
} else {
const { A, P1, T, M, k, Z = 1, Cd = 0.975 } = spec; if (!(A > 0 && P1 > 0)) throw new Error('gas: give the effective area A and the relieving pressure P1');
flow = Cd * A * P1 * Math.sqrt(k * M / (Z * R * T)) * Math.pow(2 / (k + 1), (k + 1) / (2 * (k - 1))); unit = 'kg/s';
}
const normal = spec.normal || 0, load = Math.max(flow - normal, 0);
if (load <= 0) problems.push(problem('cvfail.none', 'info', 'The failed-open valve passes no more than the normal flow: no relief load from this scenario.', [{ id: spec.id || 'valve' }]));
return { flow, normal, load, unit, problems, source: 'valve Kv definition (liquid) / critical orifice flow (gas); the maker\'s area or Kv at full opening is the input' };
}The total load on each flare header for each case, and the governing case.
load = Σ device loads of the scenarios in the case; governing = largest load
Loads are added in the order given, in 64-bit. A load that is not a number at or above 0 is left out and reported.
export function headerSummary({ scenarios, cases, devices = null }) {
const byId = Object.fromEntries(scenarios.map((s) => [s.id, s])), headers = new Map(), problems = [];
for (const c of cases) {
const per = new Map();
for (const id of c.scenarios) {
const s = byId[id]; if (!s) { problems.push(problem('header.scenario', 'error', `Case "${c.name}" names a scenario "${id}" that does not exist.`, [{ id }])); continue; }
const seen = per.get(s.header) || { load: 0, devices: {} };
for (const [dev, kgs] of Object.entries(s.loads || {})) { if (!(kgs >= 0)) { problems.push(problem('header.load', 'error', `Load of ${dev} in ${id} is not a number.`, [{ id: dev }])); continue; } seen.load += kgs; seen.devices[dev] = (seen.devices[dev] || 0) + kgs; }
per.set(s.header, seen);
}
for (const [hd, v] of per) { const h = headers.get(hd) || { header: hd, cases: [] }; h.cases.push({ name: c.name, load: v.load, devices: v.devices }); headers.set(hd, h); }
}
const out = [...headers.values()].map((h) => { const g = h.cases.reduce((a, b) => (b.load > a.load ? b : a)); return { header: h.header, governing: g.name, load: g.load, cases: h.cases }; });
for (const s of scenarios) if (!headers.has(s.header)) problems.push(problem('header.unused', 'info', `Scenario ${s.id} is not in any case: it does not load header ${s.header}.`, [{ id: s.id }]));
return { headers: out, problems, toFlareSources: (hd) => { const h = out.find((x) => x.header === hd); if (!h) return {}; return h.cases.find((c) => c.name === h.governing).devices; } };
}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.