Pressure vessel thickness and support formulas

/** Cylindrical shell: R inside radius (or D inside diameter). → { tCirc, tLong, t (governing), tNominalMin (+ ca), governs } */
export function shell({ P, R, D, S, E = 1, ca = 0 }) {
  R = R ?? D / 2; SC(S, 'S'); SC(E, 'E');
  const tCirc = P * R / (S * E - 0.6 * P), tLong = P * R / (2 * S * E + 0.4 * P), t = Math.max(tCirc, tLong);
  return { tCirc, tLong, t, tMin: t + ca, governs: tCirc >= tLong ? 'circumferential' : 'longitudinal', notes: P > 0.385 * S * E ? ['P exceeds 0.385 S E: UG-27 thin-wall equations do not apply (Appendix 1-2)'] : [], source: SOURCES.shell };
}

/** Formed head. type: 'ellipsoidal' (2:1; other ratios use K = (2 + (D/2h)²)/6), 'torispherical' (L, r given), 'hemispherical' (R inside radius). */
export function head({ type = 'ellipsoidal', P, D, S, E = 1, ca = 0, h = null, L = null, r = null }) {
  let t, notes = [];
  if (type === 'ellipsoidal') { const K = h ? (2 + Math.pow(D / (2 * h), 2)) / 6 : 1; t = P * D * K / (2 * S * E - 0.2 * P); }
  else if (type === 'torispherical') {
    L = L ?? D; r = r ?? 0.06 * D; const M = 0.25 * (3 + Math.sqrt(L / r)); t = P * L * M / (2 * S * E - 0.2 * P);
    if (r < 0.06 * L) notes.push('knuckle radius below 6 % of L: outside UG-32 limits'); if (r < 3 * t) notes.push('knuckle radius under 3t');
  } else if (type === 'hemispherical') t = P * (D / 2) / (2 * S * E - 0.2 * P);
  else throw new Error('head type: ellipsoidal | torispherical | hemispherical');
  return { t, tMin: t + ca, notes, source: SOURCES.heads };
}

export function cone({ P, D, alpha, S, E = 1, ca = 0 }) {
  const a = alpha * Math.PI / 180, notes = alpha > 30 ? ['half-apex angle above 30 degrees: needs a special analysis (UG-33 / Appendix 1-5)'] : [];
  const t = P * D / (2 * Math.cos(a) * (S * E - 0.6 * P)); return { t, tMin: t + ca, notes, source: SOURCES.cone };
}

/** Maximum allowable working pressure of a shell of corroded thickness t (inverse of UG-27). */
export function mawp({ t, R, S, E = 1 }) { const Pc = S * E * t / (R + 0.6 * t), Pl = 2 * S * E * t / (R - 0.4 * t); return { P: Math.min(Pc, Pl), governs: Pc <= Pl ? 'circumferential' : 'longitudinal' }; }

/** Hydrostatic test pressure (UG-99(b)): 1.3 × MAWP × (S at test temperature / S at design temperature). */
export const hydrotest = ({ mawp: m, Sratio = 1 }) => 1.3 * m * Sratio;

/** External pressure on a plain cylinder (no stiffening rings): critical P by the smaller of Windenburg's short-cylinder and Bresse's long-cylinder relations, divided by 4. */
export function externalCylinder({ Do, t, L, E = 2.0e5, nu = 0.3, Sy = 250 }) {
  const rr = t / Do, long = 2 * E * Math.pow(rr, 3) / (1 - nu * nu);
  const short = 2.42 * E * Math.pow(rr, 2.5) / (Math.pow(1 - nu * nu, 0.75) * (L / Do - 0.45 * Math.sqrt(rr)));
  const Pcr = Math.max(Math.min(long, short > 0 ? short : long), 0), Pa = Pcr / 4, sigma = Pcr * Do / (2 * t), notes = [];
  if (sigma > 0.5 * Sy) notes.push('the buckling stress exceeds half the yield strength: the elastic formula over-predicts; the Section II-D chart (not built) governs. Do not use this Pa.');
  return { Pcr, Pa, sigmaCr: sigma, notes, source: SOURCES.external };
}

/** Limits of reinforcement, UG-40: half-width along the vessel wall = the larger of the finished opening diameter d and (Rn + tn + t) (Rn = d/2), height normal to the wall = the smaller of 2.5 t and 2.5 tn + te. mm. */
export function reinforcementLimits({ d, t, tn, te = 0, Rn = null }) {
  Rn = Rn ?? d / 2;
  return { alongWall: Math.max(d, Rn + tn + t), normal: Math.min(2.5 * t, 2.5 * tn + te), source: 'ASME VIII-1 UG-40(b), (c)' };
}

/**
 * Opening reinforcement, UG-37. Areas in mm². vessel = { t (shell nominal, corroded), tr (required shell thickness), S, E1: 1 }, nozzle = { d (finished opening diameter, corroded),
 * tn (nozzle wall, corroded), trn (required nozzle wall), Sn, fr1?, leg41?, pad?: {Dp, te, fr4} }. F = 1.
 */
export function nozzleReinforcement({ vessel: v, nozzle: n }) {
  const fr1 = n.fr1 ?? Math.min(1, n.Sn / v.S), fr2 = fr1, E1 = v.E1 ?? 1, d = n.d;
  const A = d * v.tr + 2 * n.tn * v.tr * (1 - fr1);
  const x = E1 * v.t - v.tr, A1 = Math.max(d * x - 2 * n.tn * x * (1 - fr1), 2 * (v.t + n.tn) * x - 2 * n.tn * x * (1 - fr1));
  const A2 = Math.min(5 * (n.tn - n.trn) * fr2 * v.t, 5 * (n.tn - n.trn) * fr2 * n.tn);
  const lim = reinforcementLimits({ d, t: v.t, tn: n.tn, te: n.pad ? n.pad.te : 0 }), notes = [];
  // UG-40: material counts only inside the limits; a pad wider than the limit is cut to it
  const padW = n.pad ? Math.min(n.pad.Dp, 2 * lim.alongWall) : 0;
  if (n.pad && n.pad.Dp > 2 * lim.alongWall) notes.push('the pad is wider than the UG-40 limit of ' + (2 * lim.alongWall).toFixed(0) + ' mm: only that width is counted');
  const A41 = (n.leg41 ?? 0) ** 2 * fr2, A5 = n.pad ? (padW - d - 2 * n.tn) * n.pad.te * (n.pad.fr4 ?? fr2) : 0;
  const avail = Math.max(A1, 0) + Math.max(A2, 0) + A41 + Math.max(A5, 0);
  notes.push('limits of reinforcement (UG-40): along the wall ' + lim.alongWall.toFixed(0) + ' mm each side of the centreline, normal to it ' + lim.normal.toFixed(1) + ' mm; the areas above use them (A1 and A2 as UG-37 prints them, the pad cut to the limit)');
  return { A, A1: Math.max(A1, 0), A2: Math.max(A2, 0), A41, A5: Math.max(A5, 0), available: avail, ok: avail >= A, margin: avail / A - 1, limits: lim, source: SOURCES.nozzle, notes };
}

/** Horizontal vessel on two saddles (Zick). W total weight (N), L tangent-to-tangent length, R mean radius, H head depth, a saddle distance from the tangent line, t shell thickness. → stresses in MPa. */
export function saddle({ W, L, R, H, a, t }) {
  const Q = W / 2, k = (R * R - H * H) / (2 * a * L);
  const M1 = Q * a * (1 - (1 - a / L + k) / (1 + 4 * H / (3 * L)));
  const M2 = Q * L / 4 * ((1 + 2 * (R * R - H * H) / (L * L)) / (1 + 4 * H / (3 * L)) - 4 * a / L);
  const s1 = M1 / (Math.PI * R * R * t), s2 = M2 / (Math.PI * R * R * t);
  return { M1, M2, sigmaSaddle: s1, sigmaMid: s2, notes: ['a / R should be at most 0.5 (head stiffening at the saddle); circumferential and web stresses are not included'], source: SOURCES.saddle };
}

/** Wind velocity pressure (N/m²). */
export const windPressure = ({ V, Kz = 1, Kzt = 1, Kd = 0.95, I = 1 }) => 0.613 * Kz * Kzt * Kd * V * V * I;

/** Wind overturning of a vertical vessel: F = q Cf D_eff H, moment at the base F H/2 (uniform pressure). Units N, mm. */
export const windOnVessel = ({ q, Cf = 0.7, D, H }) => { const F = q * 1e-6 * Cf * D * H; return { F, M: F * H / 2 }; };

/** Static seismic base shear (N) = Cs W; overturning at 2/3 of the height as a first estimate. */
export const seismicBase = ({ Cs, W, H }) => ({ V: Cs * W, M: Cs * W * 2 * H / 3 });

/** Skirt-supported vertical vessel: compressive stress at the skirt-to-shell junction from weight + overturning, and the tension side. D outside mean diameter (mm), t skirt thickness. */
export function skirt({ W, M, D, t }) {
  const A = Math.PI * D * t, Z = Math.PI * D * D * t / 4;
  return { sigmaComp: W / A + M / Z, sigmaTens: M / Z - W / A, notes: ['buckling allowable of the skirt (ASME / Brownell & Young) is the engineer\'s: not computed'] };
}

/** Flat head or cover, UG-34(c)(2) (circular, attachment factor C from Figure UG-34, the caller's): t = d sqrt(C P / (S E)) + corrosion allowance. d diameter mm, P MPa, S MPa. */
export function flatHead({ P, d, S, E = 1, C = 0.33, ca = 0 }) {
  SC(S, 'S'); SC(E, 'E'); SC(C, 'C');
  const t = d * Math.sqrt(C * P / (S * E));
  return { t, tMin: t + ca, C, notes: C === 0.33 ? ['C = 0.33 is the plain-plate value; the attachment detail of Figure UG-34 sets the real C (0.17 to 0.33 and others): give it'] : [], source: 'ASME VIII-1 UG-34(c)(2), t = d sqrt(C P / S E)' };
}

/**
 * Minimum design metal temperature (MDMT), UCS-66, with the figures left to the person. The base temperature of the material curve (Figure UCS-66, curves A to D, by thickness) and the
 * temperature reduction against the stress ratio (Figure UCS-66.1) are READ FROM THE BOOK by the caller: `baseMDMT` in °C, `reduction` = [[ratio, deltaT °C], …] ascending in ratio.
 * This function forms the ratio  tr E* / (tn − c)  (UCS-66(b)(1)(b)), interpolates the reduction and subtracts it. A ratio at or below 0.35 is the UCS-66(b)(3) exemption
 * (−105 °C; typed from memory: check it). Nothing is looked up for you and a missing figure is said, never guessed.
 */
export function mdmt({ baseMDMT = null, tr, tn, E = 1, ca = 0, reduction = null }) {
  const ratio = tr * E / (tn - ca), notes = [], src = 'ASME VIII-1 UCS-66(b)';
  if (!(tn - ca > 0) || !(tr > 0)) throw new Error('give the required thickness tr and the nominal thickness tn (> corrosion allowance)');
  if (ratio <= 0.35) return { ratio, mdmt: -105, exempt: true, notes: ['stress ratio ≤ 0.35: UCS-66(b)(3) lets the MDMT go to −105 °C without impact testing (typed from memory: check against the book)'], source: src };
  if (!Array.isArray(reduction) || reduction.length < 2 || baseMDMT == null) return { ratio, mdmt: null, exempt: false, notes: ['type the base temperature of the material curve (Figure UCS-66) and the reduction points of Figure UCS-66.1 as [[ratio, deltaT °C], …]; none are built in'], source: src };
  const pts = reduction.slice().sort((a, b) => a[0] - b[0]);
  if (ratio > pts[pts.length - 1][0]) return { ratio, mdmt: baseMDMT, exempt: false, notes: ['the ratio is above the last reduction point: no reduction is taken'], source: src };
  let dT = pts[0][1];
  if (ratio >= pts[0][0]) { for (let k = 0; k + 1 < pts.length; k++) if (ratio >= pts[k][0] && ratio <= pts[k + 1][0]) { const f = (ratio - pts[k][0]) / ((pts[k + 1][0] - pts[k][0]) || 1); dT = pts[k][1] + f * (pts[k + 1][1] - pts[k][1]); break; } }
  else notes.push('the ratio is below the first reduction point: the first reduction is used');
  return { ratio, mdmt: baseMDMT - dT, reduction: dT, exempt: false, notes, source: src };
}

/** Data sheet: rows of label / value / unit + the lines every sheet must carry. */
export function dataSheet(spec) {
  const rows = [['Design pressure', spec.P, 'MPa'], ['Design temperature', spec.T, '°C'], ['Inside diameter', spec.D, 'mm'], ['Shell thickness required', spec.shellT?.toFixed(2), 'mm'], ['Shell thickness provided (nominal)', spec.tNom, 'mm'], ['Corrosion allowance', spec.ca, 'mm'], ['Allowable stress S', spec.S, 'MPa'], ['Joint efficiency E', spec.E, ''], ['Head', spec.headType, ''], ['Head thickness required', spec.headT?.toFixed(2), 'mm'], ['MAWP (corroded)', spec.mawp?.toFixed(3), 'MPa'], ['Hydrotest pressure', spec.hydro?.toFixed(3), 'MPa']];
  const text = rows.map(([a, b, c]) => `${a}: ${b ?? '-'} ${c}`).join('\n') + '\n\nDrafting and study aid, not an engineering approval. To be checked by an engineer. No representation of ASME Code compliance; no U-stamp. Allowable stress and materials are the user\'s input.';
  return { rows, text, problems: [] };
}

/** The vessel as BimCraft Plant equipment: its type and dimension parameters (mm). */
export function plantParams({ orientation = 'vertical', D, L, support }) {
  return orientation === 'horizontal' ? { type: 'vessel-h', params: { d: Math.round(D), len: Math.round(L), leg: Math.round(support || 500) } } : { type: 'vessel-v', params: { d: Math.round(D), h: Math.round(L), leg: Math.round(support || 600) } };
}

export function externalPressure({ Do, t, A, curve = null, E = null }) {
  const ratio = Do / t, notes = [];
  if (!(Do > 0 && t > 0 && A > 0)) throw new Error('give Do, t and the chart factor A');
  if (!Array.isArray(curve) || curve.length < 2) {
    if (E) return { Pa: 2 * A * E / (3 * ratio), B: null, region: 'elastic', notes: ['elastic range assumed (left of the curve) with the given E; give the material chart curve for the real answer'], needsSource: NEEDS_SOURCE };
    return { Pa: null, B: null, notes: ['give the material chart curve [[A, B], …] for the design temperature (or E for the elastic line)'], needsSource: NEEDS_SOURCE };
  }
  const pts = curve.slice().sort((a, b) => a[0] - b[0]);
  if (A <= pts[0][0]) {
    if (!E) return { Pa: null, B: null, notes: ['A is left of the curve: give E (the modulus at the temperature) for Pa = 2AE/(3 Do/t)'], needsSource: NEEDS_SOURCE };
    return { Pa: 2 * A * E / (3 * ratio), B: null, region: 'elastic', notes, source: 'ASME VIII-1 UG-28(c)(2)' };
  }
  if (A >= pts[pts.length - 1][0]) { notes.push('A is beyond the last chart point: the last B is used (flat)'); return { Pa: 4 * pts[pts.length - 1][1] / (3 * ratio), B: pts[pts.length - 1][1], region: 'chart', notes, source: 'ASME VIII-1 UG-28(c)(1)' }; }
  for (let k = 0; k + 1 < pts.length; k++) if (A >= pts[k][0] && A <= pts[k + 1][0]) {
    const f = (Math.log(A) - Math.log(pts[k][0])) / (Math.log(pts[k + 1][0]) - Math.log(pts[k][0])), B = pts[k][1] + f * (pts[k + 1][1] - pts[k][1]);
    return { Pa: 4 * B / (3 * ratio), B, region: 'chart', notes, source: 'ASME VIII-1 UG-28(c)(1)' };
  }
  return { Pa: null, B: null, notes: ['chart not readable'], needsSource: NEEDS_SOURCE };
}

export function stiffeningRing({ Do, Ls, t, As, A, Is = null }) {
  const req = Do * Do * Ls * (t + As / Ls) * A / 14;
  return { Is_required: req, ok: Is == null ? null : Is >= req, source: 'ASME VIII-1 UG-29(a), Is = Do² Ls (t + As/Ls) A / 14', notes: Is == null ? ['give the ring section moment of inertia Is to check it'] : [] };
}