Maintainer triage (June 29, 2026)
- Status: Verified bug against the current
master source in this fork.
- Severity: High.
- Area labels: area:reactors.
- Tackle batch: P4 - high-impact numerical/model 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
Part of the model-correctness audit — see #17.
Severity: high | Category: units | Location: PharmaPy/Reactors.py lines 1567
What's wrong
In PlugFlowReactor.energy_steady the specific heat-transfer area is computed as a_prime = self.diam / 4, whereas the transient energy balance (energy_balances, line 1686) correctly uses a_prime = 4 / self.diam. self.diam is the tube inner diameter (set at line 1472 from diam_in).
Why it's incorrect
For a cylindrical tube of diameter D, the heat-transfer area per unit reactor volume is a' = (perimeterL)/(cross-sectionL) = (piDL)/((pi*D^2/4)*L) = 4/D, with units m^2/m^3. Writing D/4 gives the reciprocal scaled wrong; the ratio of the two expressions is (4/D)/(D/4) = 16/D^2, so for any D != 4 m the steady-state non-adiabatic heat-transfer term is off by a factor of 16/D^2 (e.g. for D = 0.05 m it under/over-estimates the wall heat duty by ~6400x). This makes the non-adiabatic steady-state PFR temperature profile physically wrong and inconsistent with the transient solver, which uses the correct 4/D.
Suggested regression test
Solve a non-adiabatic PFR to steady state with solve_steady(adiabatic=False) and compare the temperature profile to the long-time limit of solve_unit (transient) for the same geometry and U; they should converge to the same profile. They will differ unless line 1567 is changed to 4 / self.diam.
Verification: hand-confirmed by direct code reading.
Maintainer triage (June 29, 2026)
mastersource in this fork.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
Part of the model-correctness audit — see #17.
Severity: high | Category: units | Location:
PharmaPy/Reactors.pylines 1567What's wrong
In PlugFlowReactor.energy_steady the specific heat-transfer area is computed as a_prime = self.diam / 4, whereas the transient energy balance (energy_balances, line 1686) correctly uses a_prime = 4 / self.diam. self.diam is the tube inner diameter (set at line 1472 from diam_in).
Why it's incorrect
For a cylindrical tube of diameter D, the heat-transfer area per unit reactor volume is a' = (perimeterL)/(cross-sectionL) = (piDL)/((pi*D^2/4)*L) = 4/D, with units m^2/m^3. Writing D/4 gives the reciprocal scaled wrong; the ratio of the two expressions is (4/D)/(D/4) = 16/D^2, so for any D != 4 m the steady-state non-adiabatic heat-transfer term is off by a factor of 16/D^2 (e.g. for D = 0.05 m it under/over-estimates the wall heat duty by ~6400x). This makes the non-adiabatic steady-state PFR temperature profile physically wrong and inconsistent with the transient solver, which uses the correct 4/D.
Suggested regression test
Solve a non-adiabatic PFR to steady state with solve_steady(adiabatic=False) and compare the temperature profile to the long-time limit of solve_unit (transient) for the same geometry and U; they should converge to the same profile. They will differ unless line 1567 is changed to 4 / self.diam.
Verification: hand-confirmed by direct code reading.