Heads, pressures, flows and velocities in a steady network of pipes, pumps and fittings.
How the precision is kept
The network is solved by Newton iteration in 64-bit until the total flow change is below a relative 1e-8 and the largest continuity error is below 1e-6 of the total flow. Friction comes from empirical formulas (Hazen–Williams, Swamee–Jain), and fluid properties are 4-digit values. Node pressure always uses a density of 998.2 kg/m³.
Open item. Use the fluid's own density, not a fixed 998.2 kg/m³, when a node's pressure is worked out from its head.
Calculations under this policy
9 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.
fluid.jsShared shell, every app64-bit8 calculationsReviewed
Standard gravity and the built-in density and viscosity of water, air and methane.
g = 9.80665 m/s²; water ρ = 998.2 kg/m³, μ = 1.002 × 10⁻³ Pa·s
g is the exact standard value. The densities and viscosities are values of 3 or 4 digits for one temperature each; a result that scales with them has no more digits than they do.
The Darcy friction factor from the Reynolds number and relative roughness.
Re < 2000: f = 64 ÷ Re; Re > 4000: f = 0.25 ÷ log₁₀(ε ÷ 3.7D + 5.74 ÷ Re^0.9)²; a linear blend between
Swamee–Jain is an explicit approximation of the Colebrook–White equation, so f carries that formula's error, not rounding error. Between Re 2000 and 4000 the value is a straight-line blend. Re is floored at 1e-6 to avoid dividing by zero.
The implicit friction factor, offered for comparison.
1 ÷ √f = −2 log₁₀(ε ÷ 3.7D + 2.51 ÷ (Re √f))
Fixed-point iteration from the Swamee–Jain value, stopped when f changes by less than 1e-14 (absolute) or after 60 passes. The network solver itself calls friction(), not this function.
fluid.js
lines 37 to 42
/** Colebrook–White by fixed-point iteration (the implicit equation Swamee–Jain approximates). */
export function colebrook(Re, relRough = 0) {
let f = friction(Re, relRough);
for (let i = 0; i < 60; i++) { const n = Math.pow(-2 * Math.log10(relRough / 3.7 + 2.51 / (Re * Math.sqrt(f))), -2); if (Math.abs(n - f) < 1e-14) return n; f = n; }
return f;
}
Head loss of a link
64-bit
The loss of a pipe, fitting, pump or gas pipe as a function of flow.
Hazen–Williams: h = 10.67 L Q^1.852 ÷ (C^1.852 D^4.871); Darcy–Weisbach: h = (f L ÷ D + K) Q|Q| ÷ (2 g A²); pump: H = ω² h₀ − a Q²
Hazen–Williams is an empirical law with a 4-digit coefficient. The Darcy–Weisbach slope is a central difference with step max(1e-6 |Q|, 1e-12); the power-law slopes add 1e-9 to |Q|. Slopes steer the iteration and do not enter the equations the converged flows satisfy. A pump's check valve is a resistance of 1e6 for backflow.
fluid.js
lines 47 to 79
/** The law of one link: returns { loss(Q), slope(Q) } or null for a closed link. Throws on missing parameters. */
export function linkLaw(l, kind = 'water') {
if (l.status === 'closed') return null;
const need = (k) => { if (!(l[k] > 0)) throw new Error(`link ${l.id}: ${k} must be > 0`); return l[k]; };
switch (l.type) {
case 'pipe': {
const D = need('D'), L = need('L');
if (l.law === 'dw' || l.eps != null || l.fluid) {
const fl = { ...FLUIDS.water, ...(typeof l.fluid === 'string' ? FLUIDS[l.fluid] : l.fluid) }, A = area(D), rr = (l.eps || 0) / D, K = l.K || 0;
const loss = (Q) => { const Re = 4 * fl.rho * Math.abs(Q) / (Math.PI * D * fl.mu); const f = friction(Re, rr); return (f * L / D + K) / (2 * G * A * A) * Q * Math.abs(Q); };
return { loss, slope: (Q) => Math.max(slopeOf(loss, Q), MIN_SLOPE) };
}
const C = need('C'), k = 10.67 * L / (Math.pow(C, 1.852) * Math.pow(D, 4.871));
return { loss: (Q) => k * Math.pow(Math.abs(Q), 0.852) * Q, slope: (Q) => Math.max(1.852 * k * Math.pow(Math.abs(Q) + 1e-9, 0.852), MIN_SLOPE) };
}
case 'minor': { const k = need('K') / (2 * G * area(need('D')) ** 2); return { loss: (Q) => k * Q * Math.abs(Q), slope: (Q) => Math.max(2 * k * Math.abs(Q), MIN_SLOPE) }; }
case 'resistance': { const r = need('r'), n = l.n || 2; return { loss: (Q) => r * Math.pow(Math.abs(Q), n - 1) * Q, slope: (Q) => Math.max(n * r * Math.pow(Math.abs(Q) + 1e-9, n - 1), MIN_SLOPE) }; }
case 'pump': {
const w = l.speed ?? 1; let { h0, a } = l;
if (l.curve) { h0 = l.curve[0]; a = (l.curve[0] - l.curve[2]) / (l.curve[1] ** 2); } // curve = [H at zero flow, design Q, H at design Q]
if (!(h0 > 0) || !(a > 0)) throw new Error(`link ${l.id}: a pump needs h0 > 0 and a > 0 (or curve [H0, Qdesign, Hdesign])`);
const Rcv = 1e6; // check valve: continuous at Q = 0, very steep for backflow
return { loss: (Q) => (Q >= 0 ? -w * w * h0 + a * Q * Q : -w * w * h0 + Rcv * Q), slope: (Q) => (Q >= 0 ? Math.max(2 * a * Q, MIN_SLOPE) : Rcv) };
}
case 'gaspipe': {
const D = need('D'), L = need('L'), A = area(D), fl = { ...FLUIDS.methane, ...(typeof l.fluid === 'string' ? FLUIDS[l.fluid] : l.fluid) };
const T = l.T || 288.15, Z = l.Z || 1, R = fl.R || 518.3, rr = (l.eps || 0) / D;
const loss = (m) => { const Re = 4 * Math.abs(m) / (Math.PI * D * fl.mu); return friction(Re, rr) * L * R * T * Z / (D * A * A) * m * Math.abs(m); };
return { loss, slope: (m) => Math.max(slopeOf(loss, m), MIN_SLOPE) };
}
default: throw new Error(`link ${l.id}: unknown type ${l.type}`);
}
}
Network iteration
64-bit
Updates every node head and link flow until they balance.
The nodal matrix is dense and rebuilt each pass. Convergence needs the summed flow change below 1e-8 of the summed flow and the largest node imbalance below 1e-6 × max(1, Σ|Q|). The imbalance test is absolute at 1e-6 whenever the summed flow is below 1, so it does not tighten for small networks.
fluid.js
lines 107 to 137
for (let it = 0; it < (opts.maxIter || 100) && !converged; it++) {
iterations = it + 1;
const N = Array.from({ length: nu }, () => new Array(nu).fill(0)), rhs = new Array(nu).fill(0);
// F1_i = loss_i(Q_i) − (H_a − H_b); Dinv_i = 1/slope_i; A21 rows: −1 at a (outflow), +1 at b (inflow)
links.forEach((k, i) => {
const h = k.law.loss(Q[i]), d = k.law.slope(Q[i]), inv = 1 / d, F1 = h - (H[k.a] - H[k.b]);
const ia = freeIdx.get(k.a), ib = freeIdx.get(k.b);
// N = A21 Dinv A12 with A12[i,a] = −1, A12[i,b] = +1
if (ia != null) { N[ia][ia] += inv; rhs[ia] += inv * F1; } // −A21 Dinv F1 : A21[a,i] = −1 → +inv·F1 … sign handled below
if (ib != null) { N[ib][ib] += inv; rhs[ib] -= inv * F1; }
if (ia != null && ib != null) { N[ia][ib] -= inv; N[ib][ia] -= inv; }
});
// continuity residual of the current Q: (q − A21 Q) with q = demand, A21 Q = inflow − outflow
const net_in = new Array(nodes.length).fill(0);
links.forEach((k, i) => { net_in[k.b] += Q[i]; net_in[k.a] -= Q[i]; });
// (A21 Dinv A12) dH = −A21 Dinv F1 − (q − A21 Q); with q = −demand? demand is outflow: inflow − outflow = demand (consumed at the node)
const cont = free.map((ni) => net_in[ni] - demand[ni]); // A21 Q − q
const b = free.map((_, r) => rhs[r] + cont[r]);
const dH = nu ? solveLinear(N, b) : [];
if (!dH) { problems.push(problem('network.singular', 'error', 'The nodal equations are singular: part of the network is not connected to a reservoir, or a link is closed so that some nodes are cut off.', [{ kind: 'network', id: 'all' }])); return { ok: false, problems, nodes: {}, links: {}, iterations, residual: NaN, notes }; }
free.forEach((ni, r) => { H[ni] += dH[r]; });
let dq = 0, sq = 0;
links.forEach((k, i) => {
const h = k.law.loss(Q[i]), d = k.law.slope(Q[i]), F1 = h - (H[k.a] - H[k.b]); // after the head update: dQ = −Dinv (F1 + A12 dH) is the same as −Dinv·(loss_i − ΔH_new)
const dQ = -F1 / d; Q[i] += dQ; dq += Math.abs(dQ); sq += Math.abs(Q[i]);
});
const ni2 = new Array(nodes.length).fill(0); links.forEach((k, i) => { ni2[k.b] += Q[i]; ni2[k.a] -= Q[i]; });
resid = Math.max(0, ...free.map((ni) => Math.abs(ni2[ni] - demand[ni])));
if (!Number.isFinite(sq)) { problems.push(problem('network.diverged', 'error', 'The iteration diverged (flows are not finite). Check pump curves and closed links.', [{ kind: 'network', id: 'all' }])); return { ok: false, problems, nodes: {}, links: {}, iterations, residual: NaN, notes }; }
if (dq / Math.max(sq, 1e-30) < (opts.tol || 1e-8) && resid < 1e-6 * Math.max(1, sq)) converged = true;
}
Node pressure
64-bit
The pressure reported at each node.
water: p = (H − z) × 998.2 × g; gas: p = √(max(p², 0))
The pressure of a node in a water-type network is its pressure head times a fixed 998.2 kg/m³, whatever fluid the links were given, so a network of another liquid or of air gets a pressure scaled for water. For gas the unknown is p²; a negative value is clamped to 0 before the root and reported as a problem.
A second, independent way to balance flows in loops given by hand.
ΔQ = −Σ dir × h(Q) ÷ Σ h′(Q) per loop
Stops when the largest correction in a pass is below 1e-10, an absolute flow, or after 500 passes.
fluid.js
lines 168 to 185
/**
* Hardy Cross loop corrections (1936) for links with a power law, from continuity-consistent starting flows.
* loops: [[{ link, dir: +1 | −1 }, …], …] (dir = +1 when the link's positive flow runs along the loop's direction). → { Q: { id: Q }, iterations, ok }
*/
export function hardyCross(net, loops, Q0, { tol = 1e-10, maxIter = 500 } = {}) {
const law = new Map(net.links.map((l) => [l.id, linkLaw(l, net.kind)])), Q = { ...Q0 };
for (let it = 1; it <= maxIter; it++) {
let worst = 0;
for (const loop of loops) {
let num = 0, den = 0;
for (const { link, dir } of loop) { const lw = law.get(link), q = Q[link]; num += dir * lw.loss(q); den += lw.slope(q); }
const dq = -num / den; worst = Math.max(worst, Math.abs(dq));
for (const { link, dir } of loop) Q[link] += dir * dq;
}
if (worst < tol) return { Q, iterations: it, ok: true };
}
return { Q, iterations: maxIter, ok: false };
}
numeric.jsShared shell, every app64-bit1 calculationReviewed
Solves the nodal equations of the pipe network and the load-flow Jacobian.
Gaussian elimination with partial pivoting
Works on a copy. Singular when a pivot is at most 1e-14 of the largest matrix entry, so the test scales with the matrix. No refinement step.
numeric.js
lines 7 to 21
/** Solve A x = b (A: n×n row arrays). Returns null when singular (pivot below `tiny` of the matrix scale). */
export function solveLinear(A, b) {
const n = b.length, M = A.map((r, i) => [...r, b[i]]);
let scale = 0; for (const r of A) for (const v of r) scale = Math.max(scale, Math.abs(v));
const tiny = Math.max(scale, 1e-300) * 1e-14;
for (let c = 0; c < n; c++) {
let p = c; for (let r = c + 1; r < n; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
if (Math.abs(M[p][c]) <= tiny) return null;
[M[c], M[p]] = [M[p], M[c]];
for (let r = c + 1; r < n; r++) { const f = M[r][c] / M[c][c]; if (f) for (let k = c; k <= n; k++) M[r][k] -= f * M[c][k]; }
}
const x = new Array(n);
for (let i = n - 1; i >= 0; i--) { let s = M[i][n]; for (let k = i + 1; k < n; k++) s -= M[i][k] * x[k]; x[i] = s / M[i][i]; }
return x;
}
In detail
Method
Global gradient algorithm (the file cites Todini and Pilati 1988); Hardy Cross loop corrections as a second method; Darcy–Weisbach with laminar 64/Re and Swamee–Jain (1976); Colebrook–White by fixed-point iteration.
Number types
64-bit numbers in plain arrays; a dense nodal matrix.
Tolerances
Converged when Σ|ΔQ| ÷ Σ|Q| < 1e-8 and the largest continuity error < 1e-6 × max(1, Σ|Q|); 100 iterations. Slopes have a floor of 1e-9. Numerical slope step max(1e-6 |Q|, 1e-12). Colebrook: |Δf| < 1e-14, 60 iterations. Hardy Cross: largest correction < 1e-10 (absolute), 500 iterations. Linear solve: a pivot at most 1e-14 of the largest matrix entry is singular.
Checked by
The engines test runs a one-pipe example and compares the junction head with 100 minus the Hazen–Williams loss within 1e-4, and the flow with the demand.
Limits
g = 9.80665. Water 998.2 kg/m³ and 1.002e-3 Pa·s, air 1.204 and 1.825e-5, methane 0.668, 1.1e-5 and R = 518.3 are rounded data. Hazen–Williams uses the coefficient 10.67. Gas pipes neglect the kinetic term.
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.