rocketry/src/engine/engineDesignCalcs.js

// Pure calculation functions for engine design (no React)

/**
 * Combustion chamber geometry from thermodynamic results and design inputs.
 * Returns null if required thermodynamic values are not yet available.
 */
export function calcChamber(thermo, chamber) {
  const At = thermo.At
  if (!At || !isFinite(At)) return null

  const { Lstar, contractionRatio, convAngleDeg } = chamber
  if (!Lstar || !contractionRatio || !convAngleDeg) return null

  const Ac = contractionRatio * At
  const rc = Math.sqrt(Ac / Math.PI)
  const rt = Math.sqrt(At / Math.PI)
  const Dc = 2 * rc
  const Dt = 2 * rt

  const Vc = Lstar * At
  const theta = (convAngleDeg * Math.PI) / 180
  const L_conv = (rc - rt) / Math.tan(theta)
  const V_conv = (Math.PI / 3) * L_conv * (rc * rc + rc * rt + rt * rt)
  const L_cyl = Math.max(0, (Vc - V_conv) / Ac)
  const Lc = L_cyl + L_conv

  return { Dc, Dt, rc, rt, Ac, At, Lc, L_cyl, L_conv, Vc, V_conv, contractionRatio }
}

/**
 * Nozzle geometry from thermodynamic results and nozzle design inputs.
 * Returns null if required thermodynamic values are not yet available.
 */
export function calcNozzle(thermo, nozzle) {
  const At = thermo.At
  const Ae = thermo.Ae
  if (!At || !Ae || !isFinite(At) || !isFinite(Ae)) return null

  const rt = Math.sqrt(At / Math.PI)
  const re = Math.sqrt(Ae / Math.PI)
  const Dt = 2 * rt
  const De = 2 * re

  const { type, divAngleDeg } = nozzle
  let Ln
  if (type === 'bell') {
    // 80% bell nozzle length relative to equivalent conical at 15°
    Ln = 0.8 * (re - rt) / Math.tan((15 * Math.PI) / 180)
  } else {
    // conical
    const theta = (divAngleDeg * Math.PI) / 180
    Ln = (re - rt) / Math.tan(theta)
  }

  return { Dt, De, rt, re, Ln, type }
}

/**
 * Injector orifice sizing for a given element type.
 * Returns null if required thermodynamic values are not yet available.
 */
export function calcInjector(thermo, injector) {
  const { mdot_f, mdot_ox, p0 } = thermo
  if (!mdot_f || !mdot_ox || !p0 || !isFinite(mdot_f) || !isFinite(mdot_ox) || !isFinite(p0)) return null

  const { type, N, dpFraction, Cd, rhoFuel, rhoOx } = injector
  if (!N || !dpFraction || !Cd || !rhoFuel || !rhoOx) return null

  const deltaP = dpFraction * p0

  const v_f = Cd * Math.sqrt(2 * deltaP / rhoFuel)
  const v_ox = Cd * Math.sqrt(2 * deltaP / rhoOx)

  // N elements, each element has one fuel and one oxidiser orifice
  const A_f_each = mdot_f / (N * rhoFuel * v_f)
  const A_ox_each = mdot_ox / (N * rhoOx * v_ox)

  const d_f = Math.sqrt(4 * A_f_each / Math.PI)
  const d_ox = Math.sqrt(4 * A_ox_each / Math.PI)

  return { deltaP, v_f, v_ox, A_f_each, A_ox_each, d_f, d_ox, N, type }
}

/**
 * Cooling analysis based on selected method.
 * For regenerative cooling: simplified Bartz heat flux estimate.
 * For film cooling: propellant fraction and Isp penalty estimate.
 * For ablative/uncooled: informational only.
 */
export function calcCooling(thermo, cooling, chamberGeom) {
  const { p0, T0, cstar, mdot } = thermo
  const { method } = cooling

  if (method === 'regenerative') {
    if (!p0 || !T0 || !cstar || !chamberGeom || !isFinite(p0)) return { method }

    const { Dt, Dc, Lc } = chamberGeom
    // Simplified Bartz heat flux [W/m²] using typical exhaust gas properties
    const { mu = 6e-5, cp = 2000, Pr = 0.7, T_wall = 800 } = cooling

    const q_est = (0.026 / Math.pow(Dt, 0.2)) *
      (Math.pow(mu, 0.2) * cp / Math.pow(Pr, 0.6)) *
      Math.pow(p0 / cstar, 0.8) *
      (T0 - T_wall)

    const chamberArea = Math.PI * Dc * Lc
    const q_total = q_est * chamberArea

    const { channelCount } = cooling
    const q_perChannel = channelCount > 0 ? q_total / channelCount : 0

    return { method, q_est, q_total, channelCount, q_perChannel }
  }

  if (method === 'film') {
    const { filmFraction } = cooling
    const mdot_film = (mdot || 0) * filmFraction
    // Film cooling reduces effective propellant utilisation — rough Isp penalty
    const ispPenalty = filmFraction * 100 // % estimate
    return { method, filmFraction, mdot_film, ispPenalty }
  }

  const notes = {
    ablative: 'Ablative liner — consult manufacturer data for material thickness and char rate.',
    uncooled: 'Uncooled — confirm combustion gas temperature is within material thermal limits for the burn duration.',
  }
  return { method, note: notes[method] ?? '' }
}

/**
 * Feed system sizing.
 * Pressure-fed: pressurant mass estimate.
 * Pump-fed: pump head and turbine power estimate.
 */
export function calcFeedSystem(thermo, feedSystem, burnTime) {
  const { p0, mdot_f, mdot_ox } = thermo
  if (!p0 || !mdot_f || !mdot_ox || !isFinite(p0)) return null

  const { type, feedFactor, rhoFuel, rhoOx, pressurantR, pressurantT } = feedSystem
  const tb = burnTime && burnTime > 0 ? burnTime : 30

  const V_fuel = (tb * mdot_f) / (rhoFuel || 800)
  const V_ox  = (tb * mdot_ox) / (rhoOx || 1140)
  const V_prop = V_fuel + V_ox

  if (type === 'pressure_fed') {
    const p_tank = p0 * (feedFactor || 1.3)
    const R_press = pressurantR || 2077  // J/(kg·K) — helium
    const T_press = pressurantT || 300   // K — ambient temperature
    const m_press = (p_tank * V_prop) / (R_press * T_press)
    return { type, p_tank, V_fuel, V_ox, V_prop, m_press }
  }

  // pump-fed
  const p_tank = 0.5e6 // 0.5 MPa typical inlet tank pressure
  const dP_pump = p0 - p_tank
  // Power = flow-rate × pressure rise / density (simplified, single-stage)
  const P_turbine = (mdot_f / (rhoFuel || 800) + mdot_ox / (rhoOx || 1140)) * dP_pump
  return { type, p_tank, V_fuel, V_ox, V_prop, dP_pump, P_turbine }
}