psr_solver

Steady perfectly stirred reactor (PSR) extinction-curve solver.

Traces the burning (upper) branch of the steady PSR temperature-vs-residence-time response curve through the extinction turning point, using pseudo-arclength continuation so the fold (saddle-node bifurcation) can be passed without the singular Jacobian or basin-of-attraction problems that defeat naive time-integration and residence-time-parametrized Newton solves.

The reactor is modeled as adiabatic and constant pressure, so the steady governing equations are

\[(Y_k - Y_{k,\mathrm{in}}) - \tau\, \frac{\dot\omega_k W_k}{\rho} = 0 \qquad (k = 1, \dots, K)\]
\[h(T, Y) - h_{\mathrm{in}} = 0\]

for the species mass fractions \(Y_k\) and temperature \(T\), with the residence time \(\tau\) as the continuation parameter. Here \(\dot\omega_k\) are the net molar production rates, \(W_k\) the molecular weights, \(\rho\) the mass density, and \(h\) the mass-specific mixture enthalpy. The corrector is scipy.optimize.root() (hybr) supplied with an analytic Jacobian assembled from Cantera’s kinetics derivatives, falling back to a finite-difference damped Newton when needed.

pymars.psr_solver.DS_MAX_MIDDLE = 0.05

Maximum pseudo-arclength step (scaled) on the unstable middle branch. A coarser step there lets the corrector jump off onto a spurious solution; this cap is applied only after the extinction fold, during a full stop_at_extinction=False march toward the ignition turning point.

pymars.psr_solver.PENALTY = 1000.0

Residual returned for an unphysical state so the corrector backtracks.

pymars.psr_solver.TREF = 1000.0

Temperature scale (K) used to nondimensionalize the unknown vector.

pymars.psr_solver.trace_extinction_curve(gas, stop_at_extinction=True, ds0=0.1, ds_min=0.001, ds_max=0.4, tau_start=1.0, tau_second=0.5, max_steps=6000)

Trace the PSR response curve through its turning point(s).

Parameters:

  • gas (cantera.Solution) – Gas object already set to the inlet state (inlet temperature, pressure, and composition). It is used as scratch space for kinetics/thermo evaluations and its state is modified during the solve.

  • stop_at_extinction (bool, optional) – If True (the default, used by model reduction), stop as soon as the extinction turning point is found, returning only the upper (burning) branch. If False, continue past the extinction fold along the unstable middle branch and, for a sufficiently high inlet temperature, around the lower (ignition) turning point onto the weakly-reacting branch, returning the complete S-curve. For a cold inlet the ignition turning point lies at an impractically large residence time (any reactor below the ignition delay sits near the inlet temperature), so the full march instead terminates at max_steps.

  • ds0 (float, optional) – Initial, minimum, and maximum pseudo-arclength step size (in the scaled unknown space).

  • ds_min (float, optional) – Initial, minimum, and maximum pseudo-arclength step size (in the scaled unknown space).

  • ds_max (float, optional) – Initial, minimum, and maximum pseudo-arclength step size (in the scaled unknown space).

  • tau_start (float, optional) – Residence times (s) of the two seed points used to start continuation.

  • tau_second (float, optional) – Residence times (s) of the two seed points used to start continuation.

  • max_steps (int, optional) – Maximum number of continuation steps.

Returns:

A mapping with two keys:

  • branch – a numpy.ndarray of shape (N, 2) whose columns are the residence time and temperature (tau, T) along the traced curve. With stop_at_extinction=True it ends at the extinction turning point; with stop_at_extinction=False it continues onto the middle branch (and, for a high enough inlet temperature, around the ignition turning point onto the weakly-reacting branch).

  • points – a dict with keys "extinction", "near_0.1s", and "log_mid", each mapping to a tuple (tau, T, Y) of the residence time, temperature, and mass-fraction vector at that sample point (all on the upper branch).

Return type:

dict

Raises:

RuntimeError – If the burning branch cannot be seeded, or the continuation never reaches the extinction turning point.