Storm drainage: rational method, pipe sizing and grade lines
Works out the runoff reaching each pipe, the pipe size that carries it, and the hydraulic and energy grade lines through the network.
How the precision is kept
The arithmetic is 64-bit. The accuracy is set by the empirical formulas (rational method, Manning) and by constants written to three or four digits: Manning's factor 1.486 and gravity 9.81 or 32.174. Normal depth is found by bisection run to the limit of 64-bit numbers.
Open item. Use the exact Manning factor 1 ÷ 0.3048^(1/3) and standard gravity 9.80665 m/s² (32.17404… ft/s²) instead of the rounded 1.486, 9.81 and 32.174.
Calculations under this policy
9 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.
The area of a catchment polygon in acres or hectares.
A = |Σ (x_j + x_i)(y_j − y_i)| ÷ 2; ÷ 43560 for acres, ÷ 10000 for hectares
The trapezoid form of the shoelace sum on raw coordinates, with no shift to a local origin, so large coordinates cost digits in the sum. The divisors 43560 (square feet in an acre) and 10000 (square meters in a hectare) are exact.
hydro.js
lines 65 to 70
CivilCraft 0.8.1
export function catchmentArea(pts, units = 'ft') {
let a = 0;
for (let i = 0, j = pts.length - 1; i < pts.length; j = i++) a += (pts[j].x + pts[i].x) * (pts[j].y - pts[i].y);
a = Math.abs(a) / 2;
return units === 'm' ? a / 10000 : a / 43560;
}
Rainfall intensity
64-bit
The design rainfall intensity for a time of concentration.
i = a ÷ (tc + b)^c
One 64-bit expression with Math.pow. The coefficients a, b and c come from the project; the defaults in HYDRO_DEFAULTS are placeholders.
The unit factor in Manning's equation and the acceleration of gravity used for velocity head.
k = 1 (m) or 1.486 (ft); g = 9.81 (m) or 32.174 (ft)
All three non-unit constants are rounded: 1.486 stands for 1 ÷ 0.3048^(1/3) = 1.48592…, 9.81 for 9.80665, and 32.174 for 9.80665 ÷ 0.3048 = 32.17405…. The rounded values differ from the exact ones by about 5 parts in 100,000 (1.486), 3 parts in 10,000 (9.81) and 2 parts in a million (32.174), and every flow, velocity head and junction loss carries that difference.
The flow and velocity of a circular pipe running full.
Q = (k ÷ n) × A × R^(2/3) × √S; A = πD² ÷ 4; R = D ÷ 4
One 64-bit expression using Math.PI, Math.pow with the exponent 2 ÷ 3 and Math.sqrt. A zero or negative diameter, slope or roughness returns zero flow instead of dividing.
hydro.js
lines 77 to 82
CivilCraft 0.8.1
export function manningFull(D, S, n, units = 'ft') {
if (!(D > 0) || !(S > 0) || !(n > 0)) return { Q: 0, V: 0, A: (Math.PI * D * D) / 4 };
const A = (Math.PI * D * D) / 4, R = D / 4;
const Q = (K(units) / n) * A * Math.pow(R, 2 / 3) * Math.sqrt(S);
return { Q, V: Q / A, A };
}
Part-full pipe section
64-bit
The flow area, wetted perimeter and top width of a circular pipe flowing to a given depth.
θ = 2·acos(1 − 2y ÷ D); A = D² ÷ 8 × (θ − sin θ); P = Dθ ÷ 2; T = D sin(θ ÷ 2)
The depth is clamped to 0…D and the acos argument to −1…1. For a very shallow depth θ − sin θ subtracts two nearly equal numbers and the area loses digits; the area concerned is itself very small.
The depth at which a pipe carries a given flow at its own slope.
solve (k ÷ n) × A(y) × (A(y) ÷ P(y))^(2/3) × √S = Q for y
The depth of greatest flow is taken from 201 evenly spaced samples, then 80 bisections between zero and that depth narrow the bracket below the spacing of 64-bit numbers. Because the maximum is the best sample and not the true peak, a flow between the two is reported as more than the pipe carries.
hydro.js
lines 89 to 99
CivilCraft 0.8.1
export function normalDepth(Q, D, S, n, units = 'ft') {
if (!(Q > 0)) return 0;
if (!(S > 0)) return null;
// the most a circular pipe carries is a little below full (about 0.938 D)
let hi = D, best = 0;
for (let k = 0; k <= 200; k++) { const y = (D * k) / 200, q = qAt(y, D, S, n, units); if (q > best) { best = q; hi = y; } }
if (Q > best) return null;
let lo = 0;
for (let it = 0; it < 80; it++) { const m = (lo + hi) / 2; if (qAt(m, D, S, n, units) < Q) lo = m; else hi = m; }
return (lo + hi) / 2;
}
Design flow in each pipe
64-bit
Accumulates C × A downstream, takes the longest time of concentration, and gives each pipe its flow.
ΣCA is a running sum down the network in pipe order. The 360 is exact for mm/h and hectares; the feet form leaves out the 1.008 factor. Pipe length is hypot of the plan coordinates of its two structures, and travel time uses the full-flow velocity or an assumed 3 ft/s or 1 m/s when the pipe has no slope.
hydro.js
lines 117 to 142
CivilCraft 0.8.1
export function rationalFlows(net, units = 'ft', { velocity = null } = {}) {
const h = cleanHydro(net?.hydro, units), problems = [];
const structs = new Map((net?.structures || []).map((s) => [s.id, { CA: 0, tc: 0 }]));
for (const c of net?.catchments || []) {
const s = structs.get(c.to);
if (!s) { problems.push({ id: c.id, kind: 'catchment', msg: t('hy.no.outlet', { name: c.name }) }); continue; }
s.CA += c.C * catchmentArea(c.pts, units);
s.tc = Math.max(s.tc, c.tc);
}
const pipes = new Map(), byId = new Map((net?.structures || []).map((x) => [x.id, x]));
const ord = order(net || { pipes: [] });
if (!ord) { problems.push({ id: null, kind: 'network', msg: t('hy.loop') }); return { structs, pipes, problems, hydro: h }; }
for (const p of ord) {
const up = structs.get(p.from), dn = structs.get(p.to);
if (!up || !dn) continue;
const tc = Math.max(up.tc, h.minTc), i = idf(tc, h), Q = units === 'm' ? (up.CA * i) / 360 : up.CA * i;
const a = byId.get(p.from), b = byId.get(p.to), len = Math.hypot(b.x - a.x, b.y - a.y);
const S = fin(p.invUp) && fin(p.invDn) && len > 0 ? (p.invUp - p.invDn) / len : 0;
const v = velocity?.(p) ?? (S > 0 ? manningFull(diaOf(p, units), S, h.n, units).V : 0);
const travel = len / Math.max(v > 0 ? v : h.vAssumed, 1e-6) / 60; // minutes
pipes.set(p.id, { CA: up.CA, tc, i, Q, travel, len, S });
dn.CA += up.CA;
dn.tc = Math.max(dn.tc, tc + travel);
}
return { structs, pipes, problems, hydro: h };
}
Pipe size
64-bit
The smallest standard diameter whose full-flow capacity is at least the design flow.
Each standard size in inches or millimeters is divided by the units-per-length factor from pipes.js and passed to the full-pipe formula; the comparison capacity ≥ flow has no tolerance.
hydro.js
lines 144 to 157
CivilCraft 0.8.1
export function sizePipes(net, flows, units = 'ft', { sizes = PIPE_SIZES[units === 'm' ? 'm' : 'ft'] } = {}) {
const per = DIA_UNIT[units === 'm' ? 'm' : 'ft'].perUnit, h = flows.hydro || cleanHydro(net.hydro, units);
const out = [];
for (const p of net.pipes) {
const f = flows.pipes.get(p.id);
if (!f) continue;
if (!(f.S > 0)) { out.push({ id: p.id, dia: p.dia, Q: f.Q, cap: 0, ok: false, flat: true }); continue; }
let pick = null;
for (const d of sizes) { const cap = manningFull(d / per, f.S, h.n, units).Q; if (cap >= f.Q) { pick = { dia: d, cap }; break; } }
if (!pick) { const d = sizes[sizes.length - 1]; pick = { dia: d, cap: manningFull(d / per, f.S, h.n, units).Q, none: true }; }
out.push({ id: p.id, dia: pick.dia, Q: f.Q, cap: pick.cap, ok: !pick.none });
}
return out;
}
Hydraulic and energy grade lines
64-bit
The water level and energy level at each end of each pipe and at each structure, from the outlet upstream.
Elevations are sums of inverts, depths and head terms in 64-bit. The pipe counts as full when the water below is more than 1e-9 above the crown (absolute). Velocity head and junction loss use the rounded gravity constant. The flooding message rounds the height above the rim to 2 decimals for the text only.
hydro.js
lines 159 to 192
CivilCraft 0.8.1
export function hydraulicGrade(net, flows, units = 'ft', surface = null) {
const h = flows.hydro || cleanHydro(net.hydro, units), g = G(units), problems = [];
const pipes = new Map(), structs = new Map();
const ord = order(net);
if (!ord) return { pipes, structs, problems: [{ id: null, kind: 'network', msg: t('hy.loop') }] };
const sById = new Map(net.structures.map((s) => [s.id, s]));
const rimOf = (s) => (fin(s.rim) ? s.rim : surface ? surface.elevationAt(s.x, s.y, 1e-6) : null);
for (const p of ord.slice().reverse()) { // outlet first
const f = flows.pipes.get(p.id);
if (!f || !fin(p.invUp) || !fin(p.invDn)) continue;
const D = diaOf(p, units), full = manningFull(D, f.S, h.n, units);
const below = structs.get(p.to);
const outlet = !net.pipes.some((q) => q.from === p.to);
const yn = f.S > 0 ? normalDepth(f.Q, D, f.S, h.n, units) : null;
const surch = yn === null;
const startDn = Math.max(p.invDn + (surch ? D : yn), outlet ? (h.tailwater ?? -Infinity) : below?.hgl ?? -Infinity);
// friction slope with the pipe full (pressure flow), from Manning
const A = full.A, R = D / 4, Sf = A > 0 ? Math.pow((f.Q * h.n) / (K(units) * A * Math.pow(R, 2 / 3)), 2) : 0;
const fullFlow = surch || startDn > p.invDn + D + 1e-9;
let hglDn = startDn, hglUp;
if (fullFlow) hglUp = Math.max(hglDn + Sf * f.len, p.invUp + D);
else hglUp = p.invUp + yn;
const aFlow = fullFlow ? A : partFull(yn, D).A, V = aFlow > 0 ? f.Q / aFlow : 0, hv = (V * V) / (2 * g);
pipes.set(p.id, { hglUp, hglDn, eglUp: hglUp + hv, eglDn: hglDn + hv, yn, full: fullFlow, V, Sf });
if (fullFlow) problems.push({ id: p.id, kind: 'pipe', msg: t('hy.full', { name: p.name }) });
// the structure above: the outgoing pipe's upper HGL + the junction loss
const sh = hglUp + h.kj * hv;
const prev = structs.get(p.from);
if (!prev || sh > prev.hgl) structs.set(p.from, { hgl: sh, egl: sh + hv, rim: rimOf(sById.get(p.from)) });
if (!structs.has(p.to) && outlet) structs.set(p.to, { hgl: hglDn, egl: hglDn + hv, rim: rimOf(sById.get(p.to)) });
}
for (const [id, s] of structs) if (fin(s.rim) && s.hgl > s.rim + 1e-9) problems.push({ id, kind: 'structure', msg: t('hy.flood', { name: sById.get(id)?.name || id, by: (s.hgl - s.rim).toFixed(2) }) });
return { pipes, structs, problems };
}
In detail
Method
Rational method Q = C·i·A with i = a ÷ (tc + b)^c; Manning's equation for full and part-full circular pipes; hydraulic grade line from the outlet upstream with friction slope and a junction loss k·V² ÷ 2g.
Number types
64-bit numbers.
Tolerances
Normal depth: 201 samples to find the depth of greatest flow, then 80 bisections. Full-flow test: water more than 1e-9 above the crown. Flooding test: grade line more than 1e-9 above the rim.
Limits
Manning's factor is 1.486 (the exact value is 1 ÷ 0.3048^(1/3) = 1.48592…). Gravity is 9.81 m/s² or 32.174 ft/s², not 9.80665 and its exact foot equivalent. In feet the rational formula leaves out the 1.008 factor, as the file header states. The default rainfall curve values are placeholders, as the file header states.
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.