High: [Reactors] Semibatch energy balance subtracts fictitious outflow enthalpy

[sergey-gusev94]

Maintainer triage (June 29, 2026)

• Status: Verified bug against the current master source in this fork.
• Severity: High.
• Area labels: area:reactors.
• Tackle batch: P2 - shared thermo, phase, and basis correctness.

The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.

---

Original audit details

Severity: High · Part of #67 (round-2 audit batch)

Convective enthalpy term weights the second enthalpy by holdup instead of inlet concentration, contradicting the no-outflow material balance

Location: PharmaPy/Reactors.py:1046-1048 (CSTR.energy_balances flow_term), reached for SemibatchReactor via lines 1289-1294 and 1013-1069

SemibatchReactor reuses CSTR.energy_balances. flow_term = inlet_flow*(h_in - h_temp) with h_in = sum(inlet_conch_inj) (inlet composition) and h_temp = sum(mole_conch_tempj) (HOLDUP composition) (lines 1046-1048). This is the correct in-minus-out form for a CSTR, but a semibatch reactor has NO outflow: its material balance (1289-1294, dN_i/dt = FC_in,i + r_iV; dvol_dt = inlet_flow at 1292) is pure accumulation. Deriving the consistent temperature-form balance gives NCpdT/dt = Fsum(C_in,i*(h_in,i(T_in) - h_i(T))) + source - ht, i.e. BOTH enthalpy terms must be weighted by inlet C_in,i. The code instead subtracts Fsum(C_ih_i(T)) using holdup C_i, an enthalpy of a nonexistent outflow, leaving a spurious error term F*sum((C_in,i - C_i)*h_i(T)). Since h_i is sensible enthalpy referenced to 298.15 K, the error does not vanish even when T_in ~ T.

• Trigger: SemibatchReactor(vol_tank=..., isothermal=False) (isothermal=False is documented, lines 1233-1235) with a feed composition different from the initial charge (the normal case), then solve_unit(). Exercised in tests/Flowsheet/flowsheet_tests.py:131 (holdup B+solvent, feed pure A).
• Wrong result: dtemp_dt (line 1069) uses the wrong convective enthalpy, so the temperature trajectory is wrong; through Arrhenius temp_term this also corrupts reaction rates and hence the concentration profiles of every non-isothermal semibatch run with feed composition != charge. No crash; silently wrong.
• Suggested fix: Override energy_balances for SemibatchReactor (or branch in CSTR.energy_balances on 'vol' in states_uo): replace h_temp with the inlet-weighted holdup enthalpy, flow_term = inlet_flow*(inlet_conc*(h_inj - h_tempj)).sum(axis=1)1000 = Fsum(C_in,i*(h_in,i(T_in) - h_i(T))), removing the spurious outflow-at-holdup-composition term.
• Net-new (not a duplicate): Not High: [Reactors] CSTR equilibrium branch calls getHeatOfRxn with the wrong signature #34 (CSTR getHeatOfRxn signature) nor Medium: [Containers] Mixer energy balance uses inconsistent inlet and outlet enthalpy references #38/High: [Containers] Mixer slurry energy balance weights liquid enthalpy by total slurry Cp #25/High: [Drying] Condensed-phase energy balance doubles latent heat #24/High: [Drying] Energy terms multiply molar drying rates by mass-basis heat quantities #48/Critical: [Crystallizers] Energy balances crash in adiabatic/jacket modes and mix enthalpy bases #68/Medium: [DynamicExtraction] Energy balance assigns inlet temperatures by feed/solvent names, not phase roles #80 (enthalpy-reference issues in other units). No filed item addresses the semibatch convective enthalpy basis / material-energy inconsistency.

Related existing issues (referenced, not modified): #34 (CSTR equilibrium getHeatOfRxn signature), #38/#25/#24/#48/#68/#80 (Mixer/crystallizer/drying/extraction enthalpy-reference issues) -- reference only; none address the semibatch convective enthalpy term.

---

Found by a multi-agent source audit (round 2); confirmed by re-reading the source and cross-checked against all 80 existing open issues. Each defect was verified by an independent code-truth skeptic and a novelty skeptic.

[parkyr]

Triage: reproduced ✅

Independently confirmed the root cause already stated in the audit above: SemibatchReactor inherits CSTR.energy_balances (no override — SemibatchReactor(CSTR) at PharmaPy/Reactors.py:1213), whose convective term weights the second enthalpy by the holdup composition (h_temp at Reactors.py:1047), i.e. an outflow that a semibatch does not have. Nothing to add to the diagnosis; this comment adds the executable evidence.

Environment

• Base commit: 1c3a046 (origin/master, unmodified)
• Interpreter: envs/pharmapy (PharmaPy imported from the repo checkout; assimulo present but not needed — the test calls energy_balances directly)
• Run: python -m pytest test_semibatch_flow_term.py -v

Failing test (red)

Physical setup that isolates the term: feed enters at the reactor temperature (T_in == T = 350 K, T_ref = 298.15 K) with no reaction, charge is pure B, feed is pure A. Adding material at the current temperature cannot change temperature and there is no outflow, so the convective driving term must be exactly 0 W. Any nonzero value is the spurious F·Σ(C_in,i − C_i)·h_i(T) residual.

import numpy as np
from PharmaPy.Reactors import SemibatchReactor

class _StubLiquid:            # sensible enthalpy h_i(T) = cp_i*(T - T_ref)
    def __init__(self, cp): self.cp = np.asarray(cp, float)
    def getCpPure(self, temp): return None, self.cp
    def getEnthalpy(self, temp, temp_ref, total_h=False, basis='mole'):
        temp = np.atleast_1d(temp)
        return self.cp[None, :] * (temp[:, None] - temp_ref)
    def getHeatOfRxn(self, stoich, temp, mask, dh_ref, tref): return np.zeros((1, 1))

class _StubKinetics:
    delta_hrxn = np.zeros(1); tref_hrxn = 298.15; stoich_matrix = np.zeros((1, 2))
    def get_rxn_rates(self, conc, temp, overall_rates=False, delta_hrxn=None):
        return np.zeros((1, 1))   # no reaction

def test_semibatch_convective_term_uses_inlet_not_holdup_composition():
    sb = SemibatchReactor(vol_tank=1.0, isothermal=True, temp_ref=298.15)
    sb.mask_species = np.array([True, True])
    liq = _StubLiquid([100.0, 200.0]); sb.Liquid_1 = liq; sb.Inlet = liq
    sb._Kinetics = _StubKinetics()

    inputs = {'Inlet': {'vol_flow': 2.0,
                        'mole_conc': np.array([1.0, 0.0]),   # feed: pure A
                        'temp': 350.0}}
    heat_profile = sb.energy_balances(
        0.0, np.array([0.0, 1.0]), 1.0, 350.0, 350.0,  # holdup: pure B, T = 350 K
        inputs, heat_prof=True)
    flow_term = float(heat_profile[0, -1])
    assert flow_term == 0.0, f"expected 0 W, got {flow_term} W"

E  AssertionError: expected 0 W, got -10370000.000000004 W
E  assert -10370000.000000004 == 0.0
1 failed in 0.49s

Observed vs. expected

Quantity (feed at reactor T, no reaction, charge=pure B, feed=pure A)	Observed (base)	Expected
convective flow_term (W)	-1.037e7	0

The residual equals F·Σ(C_in,i − C_i)·h_i(T) = 2·(0−1)·200·(350−298.15)·1000 W, exactly the fictitious outflow-at-holdup-composition enthalpy.

Acceptance criteria

1. Correct value: with feed at reactor temperature and no reaction, a SemibatchReactor's convective flow_term is 0 for any feed composition (including feed ≠ charge).
2. Balance consistency: the energy balance's convective term is weighted by the same inlet composition C_in,i as the material balance's accumulation term — no term weighted by holdup composition C_i.
3. No CSTR regression: CSTR.energy_balances behavior is unchanged (a CSTR has a real outflow, so holdup weighting of the second term stays correct there).
4. Existing reactor/flowsheet tests still pass.

Scope: which reactors are affected (informational)

The bug is limited to SemibatchReactor, but the buggy code lives in CSTR.energy_balances, which two classes reach:

Reactor	energy_balances source	Convective term	Affected?
CSTR	own (Reactors.py:1013)	F·(h_in − h_temp), h_temp holdup-weighted	No — correct
SemibatchReactor	inherits CSTR's	same	Yes — the bug
BatchReactor	own (Reactors.py:661)	none (closed system)	No
PlugFlowReactor	own (Reactors.py:1651)	−flow_in·ΔT/Δvol (axial gradient)	No

Why the same code is correct for a CSTR: a CSTR runs at constant volume with outflow = inflow, and the outflow carries the holdup composition C_i, so F·Σ(C_i·h_i(T)) is a real, present outflow enthalpy — the in-minus-out form is right. The identical term is wrong for a semibatch only because a semibatch has no outflow. BatchReactor (no flow term) and PlugFlowReactor (spatial-gradient convection form, no inlet-vs-holdup weighting) use their own methods and are structurally unrelated.

This is why acceptance criterion #3 above pins "no CSTR regression": the correct fix must be scoped so CSTR is untouched — either override energy_balances in SemibatchReactor, or branch inside CSTR.energy_balances on 'vol' in self.states_uo (true only for the semibatch). The trap to avoid is unconditionally switching h_temp to inlet-weighting inside CSTR.energy_balances, which would silently break the currently-correct CSTR.

Next: gh-issue-to-pr 87 — the fix PR should reuse this failing test (promoting it into the repo's test suite) and satisfy the criteria above.