Line sizing and pressure drop

/** Inside diameter in metres: { nps } / "4" / '100mm' / a number of millimetres. null when it cannot be read. */
export function insideDiameter(size) {
  if (size == null || size === '') return null;
  if (typeof size === 'number') return Number.isFinite(size) && size > 0 ? size / 1000 : null;
  const mm = /^\s*(\d+(?:\.\d+)?)\s*mm\s*$/i.exec(String(size)); if (mm) return Number(mm[1]) / 1000;
  const k = npsKey(size); if (!k) return null;
  const [od, w] = NPS_SCH40[k]; return (od - 2 * w) * 0.0254;
}

/** Colebrook–White: 1/√f = −2 log10( ε/(3.7 D) + 2.51/(Re √f) ). */
export function colebrook(Re, rr) {
  let x = 1 / Math.sqrt(0.02);
  for (let i = 0; i < 100; i++) { const xn = -2 * Math.log10(rr / 3.7 + 2.51 * x / Re); if (Math.abs(xn - x) < 1e-13) { x = xn; break; } x = xn; }
  return 1 / (x * x);
}

/** Swamee–Jain: f = 0.25 / [log10( ε/(3.7 D) + 5.74/Re^0.9 )]². */
export const swameeJain = (Re, rr) => 0.25 / Math.pow(Math.log10(rr / 3.7 + 5.74 / Math.pow(Re, 0.9)), 2);

export function friction(Re, rr = 0) {
  if (!(Re > 0)) return 0;
  const lam = 64 / Re;
  if (Re < 2300) return lam;
  const turb = (r) => swameeJain(r, rr);
  if (Re >= 4000) return turb(Re);
  const w = (Re - 2300) / 1700; return (1 - w) * 64 / 2300 + w * turb(4000);
}

/** v_e in m/s, API RP 14E with C = 100 (rho in kg/m3). */
export const erosionalVelocity = (rho, C = 100) => (C * 0.3048) / Math.sqrt(rho * 0.062427961);

/**
 * One line. fluid = { rho (kg/m3), mu? (Pa s), mdot (kg/s), phase: 'liquid'|'gas'|'two-phase' }.
 * line = { D (m), L (m), eps?, K?, fittings?: { name: count }, dz? (m rise, + up) }.
 */
export function lineDrop(fluid, line) {
  const D = line.D, A = Math.PI * D * D / 4, notes = [];
  const rho = fluid.rho, Q = fluid.mdot / rho, v = Q / A;
  let mu = fluid.mu;
  if (!(mu > 0)) { mu = fluid.phase === 'gas' ? MU_ASSUMED.gas : MU_ASSUMED.liquid; notes.push('viscosity assumed: ' + (mu * 1000).toPrecision(2) + ' mPa s'); }
  if (fluid.phase === 'two-phase') notes.push('two-phase: homogeneous mixture, screening estimate only');
  const Re = rho * v * D / mu, eps = line.eps ?? EPS_COMMERCIAL_STEEL, f = friction(Re, eps / D);
  let K = line.K || 0;
  for (const [name, n] of Object.entries(line.fittings || {})) { const k = K_FITTINGS[name]; if (k == null) notes.push('unknown fitting ' + name + ' ignored'); else K += k * n; }
  const dpFriction = f * (line.L || 0) / D * rho * v * v / 2, dpFittings = K * rho * v * v / 2, dpStatic = rho * G * (line.dz || 0);
  return { v, Q, Re, f, K, dp: dpFriction + dpFittings + dpStatic, dpFriction, dpFittings, dpStatic, rhoV2: rho * v * v, erosional: erosionalVelocity(rho), mu, notes };
}

/** The fluid of a solved stream, in the terms lineDrop wants. */
export function fluidOf(th, s, mu) {
  const mdot = massFlow(th, s), q = volFlow(th, s), rho = q > 0 ? mdot / q : 0;
  return { rho, mdot, mu, phase: s.vf >= 0.999 ? 'gas' : s.vf <= 0.001 ? 'liquid' : 'two-phase' };
}

/** The `pipe` unit. Gas is solved with the density at the mean pressure (ideal-gas scaling, isothermal), so a long gas line gets its density right. */
export function pipeUnit(th, [s], p = {}) {
  if (isZero(s)) return { outs: [emptyStream(th)], Q: 0, W: 0, notes: [] };
  const D = p.D ?? insideDiameter(p.nps ?? p.size);
  if (!(D > 0)) throw new Error('pipe needs D (m) or a nominal size nps');
  const base = fluidOf(th, s, p.mu), line = { D, L: p.L ?? 0, eps: p.eps, K: p.K, fittings: p.fittings, dz: p.dz };
  let r = lineDrop(base, line), Pout = s.P - r.dp;
  if (base.phase !== 'liquid') {
    for (let i = 0; i < 40; i++) {
      const Pm = Math.max(1e3, (s.P + Pout) / 2), rho = base.rho * Pm / s.P, rr = lineDrop({ ...base, rho }, line), Pn = s.P - rr.dp;
      r = rr; if (Math.abs(Pn - Pout) < 1e-6 * s.P) { Pout = Pn; break; } Pout = Pn;
    }
  }
  const notes = [...r.notes];
  if (!(Pout > 0)) { notes.push('the line cannot pass this flow: the drop is larger than the inlet pressure'); throw new Error(notes[notes.length - 1]); }
  const out = fromFlash(s.F, th.flashPT(s.T, Pout, s.z));
  return { outs: [out], Q: 0, W: 0, notes, dP: s.P - Pout, v: r.v, Re: r.Re, f: r.f, erosional: r.erosional, D };
}

/** A pipe run's length, rise and fittings (a corner of 60 degrees or more is a 90 degree elbow, a smaller turn a 45; mm in, metres out). */
export function runGeometry(run) {
  const p = run.points || []; let L = 0; const fittings = {};
  for (let i = 1; i < p.length; i++) L += Math.hypot(p[i][0] - p[i - 1][0], p[i][1] - p[i - 1][1], (p[i][2] || 0) - (p[i - 1][2] || 0));
  let bends = 0;
  for (let i = 1; i + 1 < p.length; i++) {
    const a = [p[i][0] - p[i - 1][0], p[i][1] - p[i - 1][1], (p[i][2] || 0) - (p[i - 1][2] || 0)], b = [p[i + 1][0] - p[i][0], p[i + 1][1] - p[i][1], (p[i + 1][2] || 0) - (p[i][2] || 0)];
    const la = Math.hypot(...a), lb = Math.hypot(...b); if (!la || !lb) continue;
    const ang = Math.acos(Math.max(-1, Math.min(1, (a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) / (la * lb)))) * 180 / Math.PI;
    if (ang >= 60) { fittings.elbow90 = (fittings.elbow90 || 0) + 1; bends++; } else if (ang >= 15) { fittings.elbow45 = (fittings.elbow45 || 0) + 1; bends++; }
  }
  for (const v of run.valves || []) { const k = VALVE_K[v.kind]; if (k) fittings[k] = (fittings[k] || 0) + 1; }
  const dz = p.length > 1 ? ((p[p.length - 1][2] || 0) - (p[0][2] || 0)) / 1000 : 0;
  return { L: L / 1000, dz, fittings, bends };
}

/** Pressure drop, velocity and erosion screening for every stream that has a line. Needs the solved result. */
export function hydraulics(model, result, ctx = {}) {
  const out = [], th = result._thermo, st = structure(model);
  if (!th) return out;
  for (const id of st.streams) {
    const link = model.links[id], s = result.streams[id];
    if (!link || link.kind !== 'stream' || !s || !(s.F > 0)) continue;
    const run = runOf(ctx.plant, link.run), pl = lineOf(ctx.pid, link.ln);
    const size = sizeFor(link, run, pl), D = insideDiameter(size);
    if (!D) { out.push({ id, ln: link.ln, size: size || '', missing: 'size' }); continue; }
    const g = run ? runGeometry(run) : { L: link.L || 0, dz: link.dz || 0, fittings: link.fittings || {}, bends: 0 };
    const geom = { D, L: link.L ?? g.L, dz: link.dz ?? g.dz, fittings: { ...g.fittings, ...(link.fittings || {}) }, K: link.K, eps: link.eps };
    const fl = fluidOf(th, s, link.mu), r = lineDrop(fl, geom);
    out.push({ id, ln: link.ln, size: nominalOf(size) || String(size), D, L: geom.L, v: r.v, Re: r.Re, f: r.f, K: r.K, dp: r.dp, rho: fl.rho, phase: fl.phase, erosional: r.erosional, notes: r.notes, dpMax: link.dpMax, mdot: fl.mdot });
  }
  return out;
}

/** Stream table for the P&ID: one row per stream, the line it runs in, its ends, state and composition. Numbers are rounded for reading; the model keeps the full values. */
export function streamTable(model, result, { hyd = null, maxComps = 6 } = {}) {
  const st = structure(model), th = result._thermo, ids = th ? th.ids : [];
  const tot = ids.map((_, i) => st.streams.reduce((a, s) => a + (result.streams[s] ? result.streams[s].F * result.streams[s].z[i] : 0), 0));
  const top = ids.map((c, i) => i).sort((a, b) => tot[b] - tot[a]).slice(0, maxComps);
  const keys = ['stream', 'line', 'from', 'to', 'T', 'P', 'vf', 'F', 'mass', 'v', 'dp', ...top.map((i) => 'x:' + ids[i])];
  const columns = ['Stream', 'Line', 'From', 'To', 'T (°C)', 'P (bar)', 'Vapour fraction', 'Flow (mol/s)', 'Mass flow (kg/h)', 'Velocity (m/s)', 'Line drop (bar)', ...top.map((i) => ids[i] + ' (mol fraction)')];
  const rows = [], round = (v, d) => (Number.isFinite(v) ? String(Number(v.toPrecision(d))) : '');
  const hmap = new Map((hyd || []).map((h) => [h.id, h]));
  for (const id of st.streams) {
    const s = result.streams[id]; if (!s) continue;
    const l = model.links[id], p = st.prod[id], c = (st.cons[id] || [])[0], h = hmap.get(id);
    rows.push([id, l?.ln || '', p ? unitTag(model, p) : '', c ? unitTag(model, c) : '', round(s.T - 273.15, 5), round(bar(s.P), 5), round(s.vf, 3), round(s.F, 5), th ? round(massFlow(th, s) * 3600, 5) : '',
      h && h.v != null ? round(h.v, 3) : '', h && h.dp != null ? round(bar(h.dp), 3) : '', ...top.map((i) => round(s.z[i], 4))]);
  }
  return { name: 'Stream table', keys, columns, rows };
}

/** The smallest nominal size whose velocity is at or under the target (m/s); the biggest when none is. */
export function pickSize(mdot, rho, vmax) {
  const Q = mdot / rho;
  for (const nps of NPS_LIST) { const D = insideDiameter(nps); if (Q / (Math.PI * D * D / 4) <= vmax) return nps; }
  return NPS_LIST[NPS_LIST.length - 1];
}

const ci = th.index(spec.scale.comp), nowKgS = s ? s.F * s.z[ci] * th.c[ci].MWkg : 0;
if (!(nowKgS > 0)) throw new Error('stream ' + sid + ' carries no ' + spec.scale.comp + ' to scale to');
const k = (spec.scale.tPerDay * 1000 / 86400) / nowKgS;
const f = def.feeds.S1;
if (f.flows) for (const c of Object.keys(f.flows)) f.flows[c] *= k; else if (f.massFlows) for (const c of Object.keys(f.massFlows)) f.massFlows[c] *= k; else f.F *= k;
result = solveModel(model);