Resources Numerical accuracy Calculation policy
Reads a reaction mechanism file and follows the species concentrations and temperature of a closed, constant-volume vessel over time.
Rate constants, equilibrium constants and species thermodynamics are closed forms in 64-bit from the coefficients in the person's file. The gas constant here is written to all its digits. The concentrations over time come from the stiff integrator with a relative tolerance of 1e-8 and an absolute tolerance of 1e-12 by default.
8 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests tools/process-test/test-chem4.mjs
R as used for rate constants, equilibrium constants and species energies in this file.
R = 8.31446261815324 J/(mol K)
All fifteen digits of the defined value. The rest of the process engine imports a ten-digit R from data.js, so the two parts of the engine do not use the same number.
const RGAS = 8.31446261815324;The factor that turns the file's activation-energy unit into J/mol.
cal/mol = 4.184; kcal/mol = 4184; K = R; eV = 96485.33212
The calorie factors are the exact thermochemical definition. The electron-volt factor is typed to ten significant digits.
const EA = { 'J/kmol': 1e-3, 'J/mol': 1, 'kJ/mol': 1e3, 'kJ/kmol': 1, 'cal/mol': 4.184, 'kcal/mol': 4184, 'K': RGAS, 'eV': 96485.33212 };C_p ÷ R, H ÷ RT and S ÷ R of a species from its NASA 7 coefficients.
C_p⁄R = a₀ + a₁T + a₂T² + a₃T³ + a₄T⁴; H⁄RT = a₀ + a₁T⁄2 + a₂T²⁄3 + a₃T³⁄4 + a₄T⁴⁄5 + a₅⁄T; S⁄R = a₀ ln T + a₁T + a₂T²⁄2 + a₃T³⁄3 + a₄T⁴⁄4 + a₆
The polynomials are written term by term with powers, not in nested form. The coefficient set for the temperature is picked by coeffsAt.
export function speciesThermo(sp, T) {
const a = coeffsAt(sp, T); if (!a) throw new Error('species ' + sp.name + ' has no NASA7 thermo');
const cpR = a[0] + a[1] * T + a[2] * T * T + a[3] * T ** 3 + a[4] * T ** 4;
const hRT = a[0] + a[1] * T / 2 + a[2] * T * T / 3 + a[3] * T ** 3 / 4 + a[4] * T ** 4 / 5 + a[5] / T;
const sR = a[0] * Math.log(T) + a[1] * T + a[2] * T * T / 2 + a[3] * T ** 3 / 3 + a[4] * T ** 4 / 4 + a[6];
return { cpR, hRT, sR, gRT: hRT - sR };
}The rate constant's pre-exponential factor in m³, mol and s.
A_SI = A × (L³ ÷ Q)^(order − 1); E_a (J/mol) = E_a × unit factor
The reaction order is the sum of the reactant coefficients, plus 1 for a third body. One power and one multiplication; the length and quantity factors are powers of ten.
const order = Object.values(reactants).reduce((a, b) => a + b, 0) + (type === 'three-body' ? 1 : 0);
const A = rate.A * Math.pow(L * L * L / Q, order - 1), Ea = (rate.Ea || 0) * eaFactor; // m3, mol, s; J/molThe forward rate constant of a reaction at a temperature.
k = A T^b exp(−E_a ÷ RT)
One 64-bit expression: a power of T and an exponential, multiplied by the pre-exponential factor.
export function rateConstant(rx, T) { return rx.A * Math.pow(T, rx.b) * Math.exp(-rx.Ea / (RGAS * T)); }K_c of a reaction, used to get the reverse rate from the forward one.
K_c = exp(−Σ νᵢ gᵢ⁄RT) × (P₀ ÷ RT)^Δν
The Gibbs sum adds products and subtracts reactants in the order of the file; the terms are large and of opposite sign, so the exponent has fewer relative digits than its terms.
/** Concentration equilibrium constant Kc of a reaction (mol/m3 units). */
export function kcOf(mech, rx, T) {
let dG = 0, dn = 0;
for (const [n, nu] of Object.entries(rx.products)) { dG += nu * speciesThermo(mech.species[mech.index[n]], T).gRT; dn += nu; }
for (const [n, nu] of Object.entries(rx.reactants)) { dG -= nu * speciesThermo(mech.species[mech.index[n]], T).gRT; dn -= nu; }
return Math.exp(-dG) * Math.pow(mech.p0 / (RGAS * T), dn);
}The rate at which each species is made or used by all reactions together.
rate = (k_f Π C^ν − k_f ÷ K_c Π C^ν) × third-body sum; ωᵢ = Σ νᵢ × rate
Near equilibrium the forward and reverse rates are nearly equal and their difference loses digits. Negative concentrations are taken as zero. Reactions are added in the order of the file.
/** Net molar production rates (mol/m3/s) at concentrations C (mol/m3) and temperature T. */
export function production(mech, C, T) {
const w = new Array(mech.names.length).fill(0), M = C.reduce((a, b) => a + Math.max(b, 0), 0);
for (const rx of mech.reactions) {
let kf = rateConstant(rx, T), kr = 0;
if (rx.reversible) kr = kf / kcOf(mech, rx, T);
let rf = kf, rr = kr;
for (const [n, nu] of Object.entries(rx.reactants)) rf *= Math.pow(Math.max(C[mech.index[n]], 0), (rx.orders && rx.orders[n]) ?? nu);
if (rx.reversible) for (const [n, nu] of Object.entries(rx.products)) rr *= Math.pow(Math.max(C[mech.index[n]], 0), nu);
let mult = 1;
if (rx.efficiencies) { mult = 0; mech.names.forEach((n, i) => { mult += (rx.efficiencies[n] ?? rx.efficiencies.default) * Math.max(C[i], 0); }); }
const rate = (rf - rr) * mult;
for (const [n, nu] of Object.entries(rx.reactants)) w[mech.index[n]] -= nu * rate;
for (const [n, nu] of Object.entries(rx.products)) w[mech.index[n]] += nu * rate;
}
return w;
}Concentrations, and temperature when adiabatic, over time.
dCᵢ/dt = ωᵢ; adiabatic: dT/dt = −Σ uᵢ ωᵢ ÷ Σ Cᵢ c_v,ᵢ, u = RT (H⁄RT − 1), c_v = R (C_p⁄R − 1)
Integrated by the stiff solver with rtol 1e-8 and atol 1e-12 by default. The points returned are the solver's accepted steps.
/** spec: { T (K), C0 { name: mol/m3 }, t (s), adiabatic?, samples? } → { t[], T[], C[][], names, problems } */
export function batchReactor(mech, spec) {
const n = mech.names.length, C0 = mech.names.map((nm) => spec.C0[nm] || 0), adiabatic = !!spec.adiabatic;
const f = (t, y) => {
const C = y.slice(0, n), T = adiabatic ? y[n] : spec.T, w = production(mech, C, T), d = w.slice();
if (adiabatic) { // constant volume: sum(C cv) dT/dt = -sum(u_i w_i), u = h - RT
let num = 0, den = 0; mech.species.forEach((sp, i) => { const th = speciesThermo(sp, T), u = RGAS * T * (th.hRT - 1), cv = RGAS * (th.cpR - 1); num += u * w[i]; den += Math.max(C[i], 0) * cv; });
d.push(-num / Math.max(den, 1e-30));
}
return d;
};
const y0 = adiabatic ? [...C0, spec.T] : C0;
const r = stiff(f, 0, y0, spec.t, { rtol: spec.rtol || 1e-8, atol: spec.atol || 1e-12 });
const pts = r.points;
return { names: mech.names, t: pts.map((p) => p[0]), C: pts.map((p) => p[1].slice(0, n)), T: pts.map((p) => (adiabatic ? p[1][n] : spec.T)), final: r.y.slice(0, n), problems: mech.problems.slice(), steps: r.steps };
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.