Pipe stress checks

/** Pipe section properties from NPS text (schedule 40) or {D, t} in mm. */
export function pipeSection(size) {
  let D, t;
  if (size && typeof size === 'object') { D = size.D; t = size.t; }
  else { const k = npsKey(size); if (!k) throw new Error('pipe size not understood: ' + size); D = NPS_SCH40[k][0] * 25.4; t = NPS_SCH40[k][1] * 25.4; }
  const d = D - 2 * t, A = Math.PI / 4 * (D * D - d * d), I = Math.PI / 64 * (D ** 4 - d ** 4);
  return { D, t, d, A, I, J: 2 * I, Z: I / (D / 2) };
}

/** Appendix D stress intensification of a tee at the run: { h, ii, io }. kind 'welding' (h = 4.4 T/r), 'unreinforced' (h = T/r), 'reinforced' (h = (T + Te/2)^2.5 / (T^1.5 r)). Typed from the code's table: check. */
export function teeSIF({ D, t, kind = 'welding', Te = 0 }) {
  const r = (D - t) / 2, T = t;
  const h = kind === 'welding' ? 4.4 * T / r : kind === 'unreinforced' ? T / r : kind === 'reinforced' ? Math.pow(T + Te / 2, 2.5) / (Math.pow(T, 1.5) * r) : null;
  if (h == null) throw new Error('tee kind: welding | unreinforced | reinforced');
  const io = Math.max(0.9 / Math.pow(h, 2 / 3), 1), ii = Math.max(0.75 * io + 0.25, 1);
  return { h, io, ii };
}

/** Appendix D stress intensification for a bend: { h, ii, io, k (flexibility, NOT applied) }. R1 bend radius (mm). */
export function bendSIF({ D, t, R1 }) {
  const r = (D - t) / 2, h = t * R1 / (r * r); return { h, ii: Math.max(0.9 / Math.pow(h, 2 / 3), 1), io: Math.max(0.75 / Math.pow(h, 2 / 3), 1), k: Math.max(1.65 / h, 1) };
}

/** Replace each bend of the polyline by its two tangent points (radius R1): a chord between them is the bend. Returns the new points, the bend ranges and a map from the old distance to the new. */
function roundBends(pts, R1) {
  const out = [pts[0]], ranges = [], knots = [[0, 0]], notes = []; let cumOld = 0, cumNew = 0;
  const push = (p, oldS) => { cumNew += len(sub(p, out[out.length - 1])); out.push(p); knots.push([oldS, cumNew]); return cumNew; };
  for (let i = 1; i < pts.length - 1; i++) {
    const l1 = len(sub(pts[i], pts[i - 1])), l2 = len(sub(pts[i + 1], pts[i])); cumOld += l1;
    const u1 = unit(sub(pts[i], pts[i - 1])), u2 = unit(sub(pts[i + 1], pts[i])), th = Math.acos(Math.max(-1, Math.min(1, dot(u1, u2))));
    if (th < 1e-3 || th > Math.PI - 1e-3) { push(pts[i], cumOld); continue; }
    const a = R1 * Math.tan(th / 2);
    if (a > 0.5 * Math.min(l1, l2)) { notes.push('the bend at point ' + (i + 1) + ' is too tight for its neighbours at radius ' + R1.toFixed(0) + ' mm: left as a stiff corner'); push(pts[i], cumOld); continue; }
    const P1 = pts[i].map((v, k) => v - u1[k] * a), P2 = pts[i].map((v, k) => v + u2[k] * a);
    const s1 = push(P1, cumOld - a), s2 = push(P2, cumOld + a); ranges.push([s1, s2]);
  }
  cumOld += len(sub(pts[pts.length - 1], pts[pts.length - 2]));
  // the last point: offset between old and new distances at the end
  push(pts[pts.length - 1], cumOld);
  return { pts: out, ranges, notes, mapS: (s) => { for (let k = 0; k + 1 < knots.length; k++) if (s <= knots[k + 1][0] + 1e-9) { const [a0, b0] = knots[k], [a1, b1] = knots[k + 1]; return a1 === a0 ? b0 : b0 + (b1 - b0) * (s - a0) / (a1 - a0); } return knots[knots.length - 1][1]; } };
}

const cum = [0]; for (let i = 1; i < pts.length; i++) cum.push(cum[i - 1] + len(sub(pts[i], pts[i - 1])));
const marks = new Map(); pts.forEach((p, i) => marks.set(+cum[i].toFixed(3), { point: p, corner: i > 0 && i < pts.length - 1 }));
const place = (sNeed) => {
  const s = Math.min(Math.max(sNeed, 0), cum[cum.length - 1]); const key = +s.toFixed(3);
  if (!marks.has(key)) { let i = 1; while (i < cum.length - 1 && cum[i] < s) i++; const a = pts[i - 1], b = pts[i], f = (s - cum[i - 1]) / (cum[i] - cum[i - 1]); marks.set(key, { point: a.map((v, k) => v + f * (b[k] - v)), corner: false }); }
  return key;
};
for (const sp of supportsIn) place(sp.s);
for (const tt of teesIn) { const key = place(tt.s); marks.get(key).tee = teeSIF({ D: sec.D, t: sec.t, kind: tt.kind, Te: tt.Te }); }
const order = [...marks.keys()].sort((x, y) => x - y), nodes = order.map((s, k) => ({ id: 'N' + k, s, point: marks.get(s).point, corner: marks.get(s).corner, tee: marks.get(s).tee }));

const dI = sec.d, wContents = (spec.contentsRho || 0) * 1e-9 * Math.PI / 4 * dI * dI * G, wIns = (spec.insulationKgPerM || 0) / 1000 * G;
const W = { id: 'W', name: 'Weight', kind: 'dead', selfWeight: 1, loads: fm.members.map((e) => ({ type: 'udl', member: e.id, dir: 'Z', w: -(wContents + wIns) })).filter((l) => l.w !== 0) };
// thermal: equivalent nodal loads of the free expansion along each member
const dT = (spec.Thot ?? 20) - (spec.Tinstall ?? 20), eps = m.alpha * dT, Fax = m.E * sec.A * eps, nodeF = new Map();
fm.members.forEach((e, k) => { const u = unit(sub(nodes[k + 1].point, nodes[k].point)); const add = (id, s) => { const f = nodeF.get(id) || [0, 0, 0, 0, 0, 0]; for (let q = 0; q < 3; q++) f[q] += s * Fax * u[q]; nodeF.set(id, f); }; add(nodes[k].id, -1); add(nodes[k + 1].id, 1); });
const T = { id: 'T', name: 'Thermal expansion', kind: 'other', selfWeight: 0, loads: [...nodeF].map(([node, F]) => ({ type: 'node', node, F })) };

const rowsFor = (R, kind) => {
  const rows = [];
  for (let k = 0; k < nn; k++) {
    const ms = getMoment(R, k); if (!ms.length) continue;
    let n = null; if (b.nodes[k].corner && ms.length === 2) { const c = cross(ms[0].ex, ms[1].ex); if (len(c) > 1e-6) n = unit(c); }
    let worst = null;
    for (const a of ms) {
      const Mb = len(a.M), Mo = n ? Math.abs(dot(a.M, n)) : 0, Mi = Math.sqrt(Math.max(Mb * Mb - Mo * Mo, 0));
      const tee = b.nodes[k].tee, ii = tee ? tee.ii : (n ? sif.ii : 1), io = tee ? tee.io : (n ? sif.io : 1);
      const tors = a.T;
      let val;
      if (kind === 'S') { const Isi = Math.max(0.75 * ii, 1), Iso = Math.max(0.75 * io, 1); val = Math.sqrt((Isi * Mi) ** 2 + (Iso * Mo) ** 2) / sec.Z + Math.abs(a.N) / sec.A; }
      else if (kind === 'O') val = Math.sqrt((Math.max(0.75 * ii, 1) * Mi) ** 2 + (Math.max(0.75 * io, 1) * Mo) ** 2) / sec.Z + Math.abs(a.N) / sec.A;
      else { const Sb = Math.sqrt((ii * Mi) ** 2 + (io * Mo) ** 2) / sec.Z, St = tors / (2 * sec.Z); val = Math.sqrt(Sb * Sb + 4 * St * St); }
      if (!worst || val > worst.val) worst = { val, Mi, Mo, tors };
    }
    rows.push({ node: b.nodes[k].id, s: b.nodes[k].s, point: b.nodes[k].point, val: worst.val, Mi: worst.Mi, Mo: worst.Mo, T: worst.tors });
  }
  return rows;
};

const pressureLong = (spec.P || 0) * sec.D / (4 * sec.t), Sh = m.Sh, Sc = m.Sc ?? m.Sh, f = spec.f ?? 1, SAbase = f * (1.25 * Sc + 0.25 * Sh);
const sust = rowsFor(res.results['case:W'], 'S').map((r) => ({ ...r, S_L: r.val + pressureLong }));
const exp = rowsFor(res.results['case:T'], 'E');
const occ = res.results['case:Wd'] ? rowsFor(res.results['case:Wd'], 'O') : null;
const stresses = b.nodes.map((n, k) => {
  const sL = sust.find((r) => r.node === n.id), sE = exp.find((r) => r.node === n.id), sO = occ?.find((r) => r.node === n.id);
  const SL = sL?.S_L ?? pressureLong, SE = sE?.val ?? 0, SA = SAbase + Math.max(0, Sh - SL), occTotal = sO ? SL + sO.val : null;
  return { node: n.id, s: n.s, point: n.point, S_L: SL, ratioSustained: SL / Sh, S_E: SE, SA, ratioExpansion: SE / SA, ...(occTotal != null ? { S_occ: occTotal, ratioOccasional: occTotal / (1.33 * Sh) } : {}) };
});

const growth = b.dT * m.alpha * (b.nodes[nn - 1].s - b.nodes[0].s);

/** Lift-off of rests and Coulomb friction (opt-in through supports of type 'rest' and `mu`). The operating case is weight + thermal; a rest with a negative reaction there is released and the run solved again. */
function solveNonlinear(b, spec) {
  const mu = spec.mu ?? 0, released = [], friction = [];
  const restNodes = (b.model.supports || []).filter((sp) => sp.fix[2] && !sp.fix[0] && !sp.fix[1] && !sp.fix[3] && !sp.fix[4] && !sp.fix[5]).map((sp) => sp.node);
  let off = new Set(), fric = new Map(), model = b.model, res = null;
  const make = () => {
    const mdl = JSON.parse(JSON.stringify(b.model)); mdl.supports = mdl.supports.filter((sp) => !off.has(sp.node));
    if (fric.size) { const T = mdl.cases.find((c) => c.id === 'T'); for (const [node, F] of fric) T.loads.push({ type: 'node', node, F: [F[0], F[1], 0, 0, 0, 0] }); }
    return mdl;
  };
  for (let it = 0; it < 12; it++) {
    model = make(); res = solveFrame(model);
    if (!res.ok) { if (off.size) { off = new Set([...off].slice(0, -1)); model = make(); res = solveFrame(model); } break; }
    const idx = (id) => model.nodes.findIndex((n) => n.id === id), Rw = res.results['case:W'].reactions, Rt = res.results['case:T'].reactions;
    let changed = false;
    for (const node of restNodes) {
      if (off.has(node)) continue;
      const i = idx(node), op = Rw[6 * i + 2] + Rt[6 * i + 2];
      if (op < -1e-6) { off.add(node); changed = true; }
    }
    if (changed) { fric = new Map(); continue; }
    if (mu > 0) {
      const u = res.results['case:T'].u; let moved = false; const next = new Map();
      for (const node of restNodes) {
        if (off.has(node)) continue;
        const i = idx(node), N = Rw[6 * i + 2] + Rt[6 * i + 2], ux = u[6 * i], uy = u[6 * i + 1], h = Math.hypot(ux, uy);
        if (h < 1e-6 || N <= 0) continue;
        const F = [-mu * N * ux / h, -mu * N * uy / h], prev = fric.get(node) || [0, 0], relaxed = [prev[0] + 0.6 * (F[0] - prev[0]), prev[1] + 0.6 * (F[1] - prev[1])];
        if (Math.hypot(relaxed[0] - prev[0], relaxed[1] - prev[1]) > 1e-3 * Math.max(1, Math.hypot(...F))) moved = true;
        next.set(node, relaxed);
      }
      fric = next; if (moved) continue;
    }
    break;
  }
  for (const node of off) { const nd = b.nodes.find((n) => n.id === node); released.push({ node, s: nd ? nd.s : 0 }); }
  for (const [node, F] of fric) friction.push({ node, F: Math.hypot(F[0], F[1]), mu });
  return { res, model, released, friction };
}

/** Loads a nozzle/anchor passes to the equipment (CH.9): the reaction at a support node, per case, in N and N·mm, and the same against a limit the caller gets from the vessel check. */
export function nozzleLoads(analysis, node, limit = null) {
  const out = {}; for (const [c, r] of Object.entries(analysis.reactions)) out[c] = r[node] || null;
  if (limit && out.thermal) { const F = Math.hypot(...out.thermal.slice(0, 3)), M = Math.hypot(...out.thermal.slice(3)); out.check = { F, M, ok: (!limit.F || F <= limit.F) && (!limit.M || M <= limit.M) }; }
  return out;
}