Reactors

/** Reaction enthalpy J per mol of extent at T (ideal-gas formation basis). */
export function reactionDH(th, rxn, T) { const nu = nuVector(th, rxn); let h = 0; for (let i = 0; i < th.n; i++) if (nu[i]) h += nu[i] * th.hIg(i, T); return h; }

/** ln K(T) of the reaction (ideal-gas standard state, 1 bar): SVN7 eq. 13.18 with the Cp polynomial integrated analytically. */
export function lnKeq(th, rxn, T) {
  const nu = nuVector(th, rxn); let dA = 0, dB = 0, dC = 0, dD = 0, dH0 = 0, dG0 = 0;
  for (let i = 0; i < th.n; i++) {
    if (!nu[i]) continue; const c = th.c[i]; if (c.gf == null) throw new Error('no Gibbs energy of formation for ' + th.ids[i]);
    dA += nu[i] * c.cp.A; dB += nu[i] * c.cp.B; dC += nu[i] * c.cp.C; dD += nu[i] * c.cp.D; dH0 += nu[i] * c.hf; dG0 += nu[i] * c.gf;
  }
  const J = dH0 / R - dA * T0 - dB / 2 * T0 * T0 - dC / 3 * T0 ** 3 + dD / T0;
  const gT = (t) => J - dA * t * Math.log(t) - dB / 2 * t * t - dC / 6 * t ** 3 - dD / (2 * t);       // (dG/R) without the integration-constant term
  const Ic = (dG0 / R - gT(T0)) / T0;
  return -(gT(T) + Ic * T) / T;
}

  const nu = nuVector(th, r), k = th.index(r.key), xi = r.X * n[k] / Math.abs(nu[k]); extents.push(xi); n = applyExtents(n, [nu], [xi]);
}
if (n.some((v) => v < -1e-9 * s.F)) notes.push('a conversion asked for more of a reactant than there is');
n = n.map((v) => Math.max(v, 0)); const F = sum(n), z = n.map((v) => v / F), P = s.P - (p.dP || 0);

const res = (xi) => {
  const n = applyExtents(n0, nus, xi), F = sum(n), lp = lnphi(n.map((v) => Math.max(v, floor))), out = [];
  for (let j = 0; j < nus.length; j++) { let s2 = 0; for (let i = 0; i < n.length; i++) if (nus[j][i]) s2 += nus[j][i] * (Math.log(Math.max(n[i], floor) / F) + lp[i] + Math.log(P / P0)); out.push(s2 - lnK[j]); }
  return out;
};

const nu = nus[0]; let lo = -Infinity, hi = Infinity;
for (let i = 0; i < nu.length; i++) { if (nu[i] > 0) lo = Math.max(lo, -n0[i] / nu[i]); else if (nu[i] < 0) hi = Math.min(hi, -n0[i] / nu[i]); }
const eps = 1e-12 * Math.max(1, Ftot0); lo = lo === -Infinity ? -1e6 * Ftot0 : lo + eps * Math.abs(lo > -Infinity ? 1 : 1); hi = hi === Infinity ? 1e6 * Ftot0 : hi - eps;
const f = (x) => res([x])[0]; const r0 = illinois(f, lo, hi, { tol: 1e-13 }); xi = [r0 == null ? (f(lo) > 0 ? lo : hi) : r0];

const target = s.F * s.H, hOut = (t) => { const x = solveAt(t), nn = applyExtents(n0, nus, x); const F = sum(nn); return F * th.flashPT(t, P, nn.map((v) => v / F)).H - target; };
T = illinois(hOut, 250, 3500, { tol: 1e-9 }); if (T == null) throw new Error('equilibrium reactor: adiabatic temperature not bracketed');

function rates(th, p, n, T, P, vol) {
  // mol/(m3 s) of each reaction; concentrations from the ideal gas (phase 'gas') or the inlet volumetric flow (phase 'liquid', constant density)
  const F = sum(n), out = [];
  for (const r of p.reactions) {
    let rate = r.k0 * Math.exp(-(r.Ea || 0) / (R * T));
    for (const [id, a] of Object.entries(r.orders || Object.fromEntries(Object.entries(r.stoich).filter(([, v]) => v < 0).map(([k, v]) => [k, -v])))) {
      const i = th.index(id), C = p.phase === 'gas' ? (n[i] / F) * P / (R * T) : n[i] / vol;
      rate *= Math.pow(Math.max(C, 0), a);
    }
    out.push(rate);
  }
  return out;
}

function cpMix(th, n, T, P, phase) {
  const F = sum(n), x = n.map((v) => v / F);
  if (phase === 'gas') return F * th.cpIgMix(x, T);
  const d = 1e-3 * T; return F * (th.phase(T + d, P, x, 'L').H - th.phase(T - d, P, x, 'L').H) / (2 * d);
}

const f = (v, y) => {
  const n = y.slice(0, th.n), T = adiabatic ? y[th.n] : isoT, rt = rates(th, p, n, T, P, vol0), dn = new Array(th.n).fill(0); let q = 0;
  for (let j = 0; j < nus.length; j++) { for (let i = 0; i < th.n; i++) dn[i] += nus[j][i] * rt[j]; q += -reactionDHPhase(th, p.reactions[j], T, P, phs) * rt[j]; }
  return adiabatic ? [...dn, q / cpMix(th, n, T, P, phase)] : dn;
};
const y0 = adiabatic ? [...n0, s.T] : n0, sol = rk45(f, y0, 0, p.V, { rtol: 1e-10, atol: 1e-14 * s.F });

const Fin = H(n0, s.T), solveV = (V, guess) => {
  const res = (x) => {
    const n = x.slice(0, th.n), T = adiabatic ? x[th.n] : isoT, rt = rates(th, p, n, T, P, vol0), r = [];
    for (let i = 0; i < th.n; i++) { let g = 0; for (let j = 0; j < nus.length; j++) g += nus[j][i] * rt[j]; r.push((n[i] - n0[i] - V * g) / (s.F)); }
    if (adiabatic) r.push((H(n, T) - Fin) / (s.F * 1e5));
    return r;
  };
  return newtonFD(res, guess, { feasible: (x) => x.slice(0, th.n).every((v) => v >= -1e-12 * s.F) && (!adiabatic || (x[th.n] > 100 && x[th.n] < 5000)), tol: 1e-12, maxIter: 80 });
};
let g = adiabatic ? [...n0, s.T] : n0.slice(), ok = true;
const steps = 24; for (let k = 1; k <= steps; k++) { const sol = solveV(p.V * (k / steps) ** 2, g); g = sol.x; ok = sol.ok; }