Piping network, pumps and compressors

/** Quadratic pump curve H = c0 + c1 Q + c2 Q² through 3 points, or H = H0 − k Q² through 2 points (the first is taken as the shut-off point if Q = 0). */
export function fitPumpCurve(pts) {
  if (!Array.isArray(pts) || pts.length < 2) throw new Error('a pump curve needs at least two (Q, H) points');
  const p = pts.slice().sort((x, y) => x[0] - y[0]);
  if (p.length === 2) { const [[q1, h1], [q2, h2]] = p; const k = (h1 - h2) / (q2 * q2 - q1 * q1); return [h1 + k * q1 * q1, 0, -k]; }
  const [[q1, h1], [q2, h2], [q3, h3]] = p;
  const s = solveLinear([[1, q1, q1 * q1], [1, q2, q2 * q2], [1, q3, q3 * q3]], [h1, h2, h3]);
  if (!s) throw new Error('pump curve points are not distinct'); return s;
}

/** Net positive suction head available (m): (P − Pv)/(rho g) + z − hLoss. P, Pv absolute (Pa), z the liquid level above the pump centreline, hLoss the suction line loss (m). */
export const npsha = ({ P, Pv, rho, z = 0, hLoss = 0 }) => (P - Pv) / (rho * G) + z - hLoss;

/** Polytropic head (m of gas column, J/kg / g given as J/kg too): Hp = n/(n−1) Z R T1 / M [ (P2/P1)^((n−1)/n) − 1 ]. Returns { head (J/kg), power (W per kg/s) }. */
export function polytropicHead({ M, T, Z = 1, n, P1, P2, eta = 1 }) {
  const e = (n - 1) / n, head = (n / (n - 1)) * Z * R * T / M * (Math.pow(P2 / P1, e) - 1); return { head, power: head / eta };
}

/** Surge check against a surge line map = [[Q, H], …] ascending in Q (the left edge of the compressor map): margin = (Q − Qsurge(H)) / Qsurge(H). */
export function compressorPoint(surgeLine, Q, H) {
  const s = surgeLine.slice().sort((a, b) => a[1] - b[1]); let q = null;
  for (let i = 1; i < s.length; i++) if (H >= s[i - 1][1] && H <= s[i][1]) { const w = (H - s[i - 1][1]) / (s[i][1] - s[i - 1][1]); q = s[i - 1][0] + w * (s[i][0] - s[i - 1][0]); }
  if (q == null) return { qSurge: null, margin: null, surge: false, notes: ['the head is outside the surge line'] };
  const margin = (Q - q) / q; return { qSurge: q, margin, surge: margin <= 0, notes: margin < 0.1 ? ['less than 10 % flow margin to surge (a common design margin; the project sets its own)'] : [] };
}

/** Greitzer's B parameter (1976): B = U/(2 as) sqrt(Vp/(Ac Lc)). B above about 0.7–0.8 tends to surge, below to rotating stall (published rule of thumb). */
export const greitzerB = ({ U, as, Vp, Ac, Lc }) => U / (2 * as) * Math.sqrt(Vp / (Ac * Lc));

/** Control valve: Kv in m3/h of water at 1 bar drop. Equal-percentage Kv(x) = Kvs R^(x−1), linear Kv(x) = Kvs x, quick opening Kvs sqrt(x). */
export function cvalveKv(l, x = l.opening ?? 1) {
  const Kvs = l.Kv, c = l.charact || 'equal'; x = Math.min(1, Math.max(0, x));
  if (x <= 0) return 0; return c === 'linear' ? Kvs * x : c === 'quick' ? Kvs * Math.sqrt(x) : Kvs * Math.pow(l.R || 50, x - 1);
}

/** Link head-loss function: returns { h (loss a→b in the network's head unit), p (dh/dQ, > 0) } for flow Q (m3/s liquid, kg/s gas). */
function linkLoss(net, st, l, Q) {
  const f = net.fluid, liquid = f.kind !== 'gas';
  const A = Math.PI * l.D * l.D / 4, floor = 1e-9;
  if (l.kind === 'pump') {
    const c = l.coef; const Qa = Math.abs(Q);
    if (Q < 0) { const r = 1e3; return { h: r * Q * Math.abs(Q) - c[0], p: Math.max(2 * r * Qa, floor) }; }   // no reverse flow through a pump: a large resistance
    return { h: -pumpHead(c, Q), p: Math.max(-(c[1] + 2 * c[2] * Q), 1e-6 * (Math.abs(c[2]) + 1)) };
  }
  if (l.kind === 'check' && Q < 0) return { h: 1e6 * Q, p: 1e6 };
  if (liquid) {
    const rho = f.rho, mu = f.mu > 0 ? f.mu : MU_ASSUMED.liquid, Qa = Math.max(Math.abs(Q), 1e-12), v = Qa / A, Re = rho * v * l.D / mu;
    if (l.kind === 'cvalve') { const Kv = cvalveKv(l); if (Kv <= 0) return { h: 1e9 * Q, p: 1e9 }; const k = rho / 1000 * 1e5 / (G * rho) * Math.pow(3600 / Kv, 2); return { h: k * Q * Qa, p: Math.max(2 * k * Qa, floor) }; }
    const ff = (q) => { const vv = Math.max(Math.abs(q), 1e-12) / A; return friction(rho * vv * l.D / mu, (l.eps ?? EPS_COMMERCIAL_STEEL) / l.D); };
    let K = l.K || 0; for (const [n, c] of Object.entries(l.fittings || {})) K += (K_FITTINGS[n] || 0) * c;
    const Kf = fitK(net, l, Re); K += Kf; l._K = K; l._Kfit = Kf; l._Re = Re;
    const kk = (ff(Qa) * (l.L || 0) / l.D + K) / (2 * G * A * A), dz = 0;
    const h = kk * Q * Qa;
    // d/dQ of f(Q) Q|Q| neglected beyond the quadratic part: the turbulent f varies slowly, the iteration absorbs the rest
    const lam = Re < 2300 ? 1 : 0; const p = Math.max(2 * kk * Qa * (1 - 0.5 * lam), floor);
    return { h: h + dz, p: lam ? Math.max(p, 32 * mu * (l.L || 0) / (rho * G * l.D * l.D * A)) : p };
  }
  // gas: relation in P² for the isothermal compressible line  P1² − P2² = (R T / M) (W/A)² [ f L / D + K + 2 ln(P1/P2) ]
  const W = Q, Wa = Math.max(Math.abs(W), 1e-12), RT = R * f.T * (f.Z || 1) / f.M, mu = f.mu > 0 ? f.mu : MU_ASSUMED.gas;
  const Re = 4 * Wa / (Math.PI * l.D * mu);
  if (l.kind === 'cvalve') throw new Error('control valves on gas lines: give the valve as a pipe with K for now (gas Cv sizing is not built)');
  let K = l.K || 0; for (const [n, c] of Object.entries(l.fittings || {})) K += (K_FITTINGS[n] || 0) * c;
  const Kf = fitK(net, l, Re); K += Kf; l._K = K; l._Kfit = Kf; l._Re = Re;
  const fk = friction(Re, (l.eps ?? EPS_COMMERCIAL_STEEL) / l.D) * (l.L || 0) / l.D + K;
  let ln = 0; if (st && st.P1 > 0 && st.P2 > 0) ln = 2 * Math.log(st.P1 / st.P2);
  const k = RT * (fk + ln) / (A * A);
  return { h: k * W * Wa, p: Math.max(2 * k * Wa, floor) };
}

for (; iter < maxIter; iter++) {
  // node i (free):  Σ_in Q − Σ_out Q = demand,  Q = c0 + g (Ha − Hb)  →  Σ g (Hi − Hn) = −demand + Σ_in c0 − Σ_out c0 + Σ g H(fixed neighbour)
  const A2 = Array.from({ length: nf }, () => new Array(nf).fill(0)), rhs = free.map((ni) => -dem[ni]), rec = [];
  links.forEach((l, k) => {
    const st = gas ? { P1: pOf(net.nodes[l.ia], H[l.ia]), P2: pOf(net.nodes[l.ib], H[l.ib]) } : null;
    const { h, p } = linkLoss(net, st, l, Q[k]), g = 1 / p, c0 = Q[k] - h / p; rec.push({ g, c0 });
    const ra = fi[l.ia], rb = fi[l.ib];
    if (rb >= 0) { A2[rb][rb] += g; rhs[rb] += c0; if (ra >= 0) A2[rb][ra] -= g; else rhs[rb] += g * H[l.ia]; }
    if (ra >= 0) { A2[ra][ra] += g; rhs[ra] -= c0; if (rb >= 0) A2[ra][rb] -= g; else rhs[ra] += g * H[l.ib]; }
  });
  const Hn = nf ? solveLinear(A2, rhs) : [];
  if (nf && !Hn) throw new Error('the network is not connected to a fixed-pressure node, or a link has no resistance');
  let dmax = 0; free.forEach((ni, r) => { const hv = Hn[r]; dmax = Math.max(dmax, Math.abs(hv - H[ni]) / (Math.abs(hv) + 1)); H[ni] = gas ? Math.max(hv, 1) : hv; });
  let dq = 0; links.forEach((l, k) => { const qn = rec[k].c0 + rec[k].g * (H[l.ia] - H[l.ib]); dq = Math.max(dq, Math.abs(qn - Q[k]) / (Math.abs(qn) + 1e-12)); Q[k] = qn; });
  err = Math.max(dmax, dq); if (err < tol) { ok = true; iter++; break; }
}

const nodes = {}; net.nodes.forEach((n, i) => { nodes[n.id] = { P: pOf(n, H[i]), head: gas ? null : H[i], demand: dem[i] }; });
const out = {};
links.forEach((l, k) => {
  const A = Math.PI * l.D * l.D / 4, q = Q[k];
  const r = { Q: q, dH: H[l.ia] - H[l.ib] }; if (gas) { r.mdot = q; r.v = q / (A * Math.max(1, (nodes[l.a].P + nodes[l.b].P) / 2) * net.fluid.M / (R * net.fluid.T * (net.fluid.Z || 1))); r.dP = nodes[l.a].P - nodes[l.b].P; }
  else { r.mdot = q * rho; r.v = q / A; r.dP = nodes[l.a].P - nodes[l.b].P; if (l.kind === 'pump') { r.head = pumpHead(l.coef, q); if (q < 0) problems.push(prob('net.pump-reverse', 'error', `Pump ${l.id} would run backwards at this solution.`, l.id)); } }
  if (!gas && l.kind !== 'pump') { const lim = 3; if (Math.abs(r.v) > lim) problems.push(prob('net.velocity', 'warning', `Line ${l.id}: velocity ${Math.abs(r.v).toFixed(1)} m/s is above the ${lim} m/s screening limit for liquids.`, l.id)); }
  if (gas) { const cs = Math.sqrt((net.fluid.k || 1.3) * R * net.fluid.T * (net.fluid.Z || 1) / net.fluid.M); if (Math.abs(r.v) > 0.7 * cs) problems.push(prob('net.mach', 'warning', `Line ${l.id}: Mach number ${(Math.abs(r.v) / cs).toFixed(2)} (isothermal model is not valid close to choking).`, l.id)); }
  if (l._K != null) { r.K = l._K; r.Kfit = l._Kfit; r.Re = l._Re; }
  out[l.id] = r;
});

export function compressorHead({ P1, P2, T1, M, k, Z = 1, eta, kind = 'polytropic', mdot = null } = {}) {
  const bad = (what) => { throw new Error(what); };
  if (!(P1 > 0 && P2 > 0)) bad('P1 and P2 are absolute pressures in Pa');
  if (!(P2 > P1)) bad('this is a compressor: the discharge pressure must be above the suction pressure');
  if (!(T1 > 0)) bad('T1 is the suction temperature in K');
  if (!(M > 0)) bad('M is the molar mass in kg/mol (air 0.02897)');
  if (!(k > 1)) bad('k = cp/cv must be above 1');
  if (!(Z > 0)) bad('Z must be positive');
  if (!(eta > 0 && eta <= 1)) bad('give the efficiency as a number above 0 and up to 1');
  if (kind !== 'polytropic' && kind !== 'isentropic') bad('kind is polytropic or isentropic');
  const r = P2 / P1, a = (k - 1) / k, RT = Z * R * T1 / M, rs = Math.pow(r, a);
  const Hs = RT / a * (rs - 1);
  let etaP, etaS, m;
  if (kind === 'polytropic') { etaP = eta; m = a / etaP; etaS = (rs - 1) / (Math.pow(r, m) - 1); }
  else { etaS = eta; m = Math.log(1 + (rs - 1) / etaS) / Math.log(r); etaP = a / m; }
  const Hp = RT / m * (Math.pow(r, m) - 1), work = kind === 'polytropic' ? Hp / etaP : Hs / etaS, n = 1 / (1 - m), T2 = T1 * Math.pow(r, m);
  const notes = ['Ideal gas with constant k and Z, adiabatic, no inter-cooling. A real gas near its critical point, or a multistage machine with inter-coolers, needs an equation of state and a stage-by-stage calculation.'];
  if (Z !== 1) notes.push(`Z = ${Z} is used as given and held constant from suction to discharge.`);
  const out = { r, m, n, etaS, etaP, Hs, Hp, work, T2, headM: Hp / G, notes };
  if (mdot != null) { if (!(mdot >= 0)) bad('mdot is a mass flow in kg/s'); out.power = mdot * work; }
  return out;
}

/** The compressor's operating point against a surge line [[Q, H], …] drawn in the unit of `headUnit` ('J/kg' or 'm'); the polytropic head of the run is converted to that unit. */
export function compressorSurge(line, Q, run, headUnit = 'J/kg') {
  const H = headUnit === 'm' ? run.Hp / G : run.Hp;
  return { ...compressorPoint(line, Q, H), H, headUnit };
}