Water hammer and pump trip

/** Wave speed of a liquid in an elastic pipe (Korteweg–Joukowsky): a = sqrt( (K/rho) / (1 + (K D / (E e)) c1) ), c1 = 1 − nu²/2... anchored: 1 − nu² (both ends fixed), else 1 − nu/2 (anchored at the upstream end only) — the usual published pair. */
export function waveSpeed({ K = 2.19e9, rho = 998, D, e, E = 2.0e11, nu = 0.3, anchored = 'both' }) {
  const c1 = anchored === 'joints' ? 1 : anchored === 'upstream' ? 1 - nu / 2 : 1 - nu * nu;
  if (!(e > 0)) return Math.sqrt(K / rho);
  return Math.sqrt((K / rho) / (1 + (K * D / (E * e)) * c1));
}

/** Joukowsky: the head or pressure rise for an instantaneous velocity change dV: dP = rho a dV. */
export const joukowsky = (a, dV, rho) => rho * a * dV;

/**
 * Water hammer by the method of characteristics in one pipe: a reservoir upstream, a valve downstream (closing or opening by tau(t) = relative opening 0..1).
 *   spec = { L, D, a (wave speed), f? (Darcy friction, quasi-steady), Hres (m, reservoir head above the valve outlet), V0 (initial velocity m/s, steady-state before the event),
 *            tau: (t) => 0..1 or [[t, tau], …], tEnd (s), N? (reaches), Hv? (vapour head in the same datum, m: below it the run says column separation may occur), rho? }
 *   → { t[], Hvalve[], Qvalve[], Hmid[], Hmax, Hmin, tMax, period (4L/a), problems }
 * Heads are measured from the valve outlet; the steady head at the valve is Hres − friction. Not modelled: column separation (flagged), pump trip, unsteady friction.
 */
export function hammer(spec) {
  const { L, D, a, f = 0, Hres, V0, tEnd, N = 50, Hv = null, rho = 998 } = spec;
  const dx = L / N, dt = dx / a, A = Math.PI * D * D / 4, B = a / (G * A), Rr = f * dx / (2 * G * D * A * A);
  const tau = typeof spec.tau === 'function' ? spec.tau : ((pts) => (t) => { if (t <= pts[0][0]) return pts[0][1]; for (let i = 1; i < pts.length; i++) if (t <= pts[i][0]) { const w = (t - pts[i - 1][0]) / (pts[i][0] - pts[i - 1][0]); return pts[i - 1][1] + w * (pts[i][1] - pts[i - 1][1]); } return pts[pts.length - 1][1]; })(spec.tau);
  const Q0 = V0 * A;
  let Q = new Array(N + 1).fill(Q0), Hh = new Array(N + 1);
  for (let i = 0; i <= N; i++) Hh[i] = Hres - Rr * Q0 * Q0 * (i / 1) * 1;   // head falls linearly along the pipe by friction
  const Hv0 = Hh[N];   // steady head at the valve (reservoir head less friction over N reaches)
  const CvOpen = Q0 * Q0 / Hv0;   // Q² = Cv τ² H at the valve with τ = 1 giving Q0 at the steady head
  const T = [], Hval = [], Qval = [], Hmid = [], problems = [];
  let Hmax = -Infinity, Hmin = Infinity, tMax = 0; const steps = Math.ceil(tEnd / dt);
  for (let s = 0; s <= steps; s++) {
    const t = s * dt;
    if (s > 0) {
      const Qn = new Array(N + 1), Hn = new Array(N + 1);
      for (let i = 1; i < N; i++) {
        const CP = Hh[i - 1] + (B - 0) * Q[i - 1] - Rr * Q[i - 1] * Math.abs(Q[i - 1]) * 0 - Rr * Q[i - 1] * Math.abs(Q[i - 1]);
        const CM = Hh[i + 1] - B * Q[i + 1] + Rr * Q[i + 1] * Math.abs(Q[i + 1]);
        // Wylie–Streeter: HP = 0.5 (CP + CM), QP = (CP − CM) / (2B) with CP = H(i−1) + B Q(i−1) − R Q|Q|, CM = H(i+1) − B Q(i+1) + R Q|Q|
        Hn[i] = 0.5 * (CP + CM); Qn[i] = (CP - CM) / (2 * B);
      }
      // upstream reservoir: C− from node 1
      { const CM = Hh[1] - B * Q[1] + Rr * Q[1] * Math.abs(Q[1]); Hn[0] = Hres; Qn[0] = (Hres - CM) / B; }
      // downstream valve: C+ from node N−1;  Q² = Cv τ² HP  with HP = CP − B Q
      { const CP = Hh[N - 1] + B * Q[N - 1] - Rr * Q[N - 1] * Math.abs(Q[N - 1]), tv = tau(t), c = CvOpen * tv * tv;
        if (c <= 0) { Qn[N] = 0; Hn[N] = CP; } else { const q = (-B * c + Math.sqrt(B * B * c * c + 4 * c * CP)) / 2; Qn[N] = q; Hn[N] = CP - B * q; } }
      Q = Qn; Hh = Hn;
    }
    T.push(t); Hval.push(Hh[N]); Qval.push(Q[N]); Hmid.push(Hh[N >> 1]);
    if (Hh[N] > Hmax) { Hmax = Hh[N]; tMax = t; } for (const hv of Hh) if (hv < Hmin) Hmin = hv;
  }
  if (Hv != null && Hmin < Hv) problems.push(prob('hammer.cavitation', 'error', `The head falls to ${Hmin.toFixed(1)} m, below the vapour head ${Hv.toFixed(1)} m: the liquid column can separate and this model (no column separation) is not valid after that moment.`, 'hammer'));
  return { t: T, Hvalve: Hval, Qvalve: Qval, Hmid, Hmax, Hmin, tMax, period: 4 * L / a, steady: Hv0, dtStep: dt, pressureMax: rho * G * Hmax, problems };
}

export const k2K = ({ K1, Kinf }, Re, Dinch) => K1 / Re + Kinf * (1 + 1 / Dinch);

export const k3K = ({ K1, Ki, Kd }, Re, Dinch) => K1 / Re + Ki * (1 + Kd / Math.pow(Dinch, 0.3));

export function pumpTrip(spec) {
  const { L, D, a, Hs, Hd, Q0, H0, eta, w0, I, tEnd, N = 40, Hv = null, rho = 998, tTrip = 0 } = spec;
  const h = spec.curve ? spec.curve.slice() : fitPumpCurve(spec.points);
  const A = Math.PI * D * D / 4, B = a / (G * A), dx = L / N, dt = dx / a;
  const hpump = (Q, al) => al * al * (h[0] + h[1] * Q / al + h[2] * Q * Q / (al * al));
  // calibrate friction so the steady state is exact: Hs + H(Q0) = Hd + loss
  const Hp0 = hpump(Q0, 1), loss = Hs + Hp0 - Hd;
  if (!(loss > 0)) throw new Error('the pump head at the rated flow does not exceed the static lift: no steady flow');
  const Rr = loss / (N * Q0 * Q0);                                      // per reach, head per Q|Q|
  let Q = new Array(N + 1).fill(Q0), Hh = Array.from({ length: N + 1 }, (_, i) => Hs + Hp0 - Rr * Q0 * Q0 * i);
  let w = w0, closed = false, tClose = null; const T = [], Hpump = [], Qpump = [], Spd = [], problems = [];
  let Hmin = Infinity, Hmax = -Infinity; const steps = Math.ceil(tEnd / dt);
  for (let s = 0; s <= steps; s++) {
    const t = s * dt;
    if (s > 0) {
      const Qn = new Array(N + 1), Hn = new Array(N + 1);
      for (let i = 1; i < N; i++) {
        const CP = Hh[i - 1] + B * Q[i - 1] - Rr * Q[i - 1] * Math.abs(Q[i - 1]), CM = Hh[i + 1] - B * Q[i + 1] + Rr * Q[i + 1] * Math.abs(Q[i + 1]);
        Hn[i] = 0.5 * (CP + CM); Qn[i] = (CP - CM) / (2 * B);
      }
      { const CP = Hh[N - 1] + B * Q[N - 1] - Rr * Q[N - 1] * Math.abs(Q[N - 1]); Hn[N] = Hd; Qn[N] = (CP - Hd) / B; }
      // pump end: C- from node 1
      const CM = Hh[1] - B * Q[1] + Rr * Q[1] * Math.abs(Q[1]);
      const tripped = !spec.noTrip && t >= tTrip;
      if (closed) { Qn[0] = 0; Hn[0] = CM; }
      else {
        let al = w / w0, Qp = Q[0], wNew = w;
        for (let it = 0; it < 30; it++) {
          // solve  Hs + al^2 h0 + al h1 Q + h2 Q^2 = CM + B Q
          const qa = h[2], qb = al * h[1] - B, qc = Hs + al * al * h[0] - CM;
          const disc = qb * qb - 4 * qa * qc;
          let q; if (Math.abs(qa) < 1e-14) q = -qc / qb; else if (disc < 0) q = -qb / (2 * qa); else { const r1 = (-qb + Math.sqrt(disc)) / (2 * qa), r2 = (-qb - Math.sqrt(disc)) / (2 * qa); q = Math.abs(r1 - Qp) < Math.abs(r2 - Qp) ? r1 : r2; }
          Qp = q;
          if (tripped) { const Hp = hpump(Math.max(Qp, 0), Math.max(al, 1e-6)), torque = rho * G * Math.max(Qp, 0) * Hp / (eta * Math.max(w, 1e-3)); wNew = Math.max(w - dt / I * torque, 0); al = 0.5 * (w + wNew) / w0; }
        }
        if (tripped) w = wNew;
        if (tripped && Qp <= 0) { closed = true; tClose = t; Qn[0] = 0; Hn[0] = CM; } else { Qn[0] = Qp; Hn[0] = CM + B * Qp; }
      }
      Q = Qn; Hh = Hn;
    }
    T.push(t); Hpump.push(Hh[0]); Qpump.push(Q[0]); Spd.push(w / w0);
    for (const hv of Hh) { if (hv < Hmin) Hmin = hv; if (hv > Hmax) Hmax = hv; }
  }
  if (Hv != null && Hmin < Hv) problems.push(problem('pumptrip.cavitation', 'error', `After the trip the head falls to ${Hmin.toFixed(1)} m, below the vapour head ${Hv.toFixed(1)} m: the liquid can separate and this model (no column separation) is not valid.`, [{ id: 'pump' }]));
  return { t: T, Hpump, Qpump, speed: Spd, Hmin, Hmax, tClose, steadyH: Hs + Hp0, friction: Rr, period: 4 * L / a, dtStep: dt, problems, notes: ['no reverse flow through the pump is modelled: the check valve is assumed to close at zero flow'] };
}