Numeric solver runaway using ReactorNet.advance_to_steady_state

[philipp1337]

I want to simulate burning methane in an IdealGasConstPressureReactor. I use different mass-flows and equivalence-ratios using two vectors you can see in my example code below. However, using ReactorNet.advance_to_steady_state breaks the solver at certain points, leading to huge reactor masses and volumes and therefore CanteraErrors. If I varying the equivalence-ratio leads to different combustors failing, varying the mass-flows to each reactor has no effect on this. Using a ConstPressureReactor also leads to the same results, using an IdealGasReactor eliminates the error while the pressure before and after the reactor stays about the same (0.1% differences).

Steps to reproduce

My example code:

import cantera as ct
import numpy as np

# Initial constants
inlet_temperature = 712 # K
inlet_pressure = 2171756.0738435777 # Pa
inlet_density = 10.529 # kg/m³
phi_mean = 0.3
phi_std = 0.25
mdot = 350 # kg/s

# Vectors for phi and mass flow ratio
x=[0.1125, 0.1375, 0.1625, 0.1875, 0.2125, 0.2375, 0.2625, 0.2875, 0.3125, 0.3375,
 0.3625, 0.3875, 0.4125, 0.4375, 0.4625, 0.4875]
y_norm=[0.05223705, 0.05602466, 0.05948904, 0.0625391,  0.06509137, 0.0670737,
 0.06842868, 0.0691164,  0.0691164,  0.06842868, 0.0670737,  0.06509137,
 0.0625391,  0.05948904, 0.05602466, 0.05223705]

psr_exhaust_gases = []


# Run Combustors for each Massflow/Phi
n = 0
for phi, ratio in zip(x,y_norm):
    gas = ct.Solution("gri30.xml")
    gas.TP = inlet_temperature, inlet_pressure
    print(f"{n} Combustor using phi={phi} and mass flow {ratio*mdot}")
    gas.set_equivalence_ratio(phi, 'CH4:1.0', 'O2:1.0, N2:3.76')
    gas.equilibrate("HP")
    inlet = ct.Reservoir(gas)
    exhaust = ct.Reservoir(gas)
    reactor = ct.IdealGasConstPressureReactor(gas)
    reactor.volume = 1
    inlet_mfc = ct.MassFlowController(inlet, reactor, mdot=mdot*ratio)
    outlet_mfc = ct.PressureController(reactor, exhaust, master=inlet_mfc, K=0.98)  
    
    # Perform simulation
    sim = ct.ReactorNet([reactor])
    try:
        sim.advance_to_steady_state()
    except Exception:
        print(f"Reactor {n} failed")
    print(f"Reactor Mass: {reactor.mass}")
    psr_exhaust_gases.append(gas)
    n+=1

Behavior

0 Combustor using phi=0.1125 and mass flow 18.2829675
Reactor Mass: 7.509287323526787
1 Combustor using phi=0.1375 and mass flow 19.608631
Reactor 1 failed
Reactor Mass: -31989052.24073112
2 Combustor using phi=0.1625 and mass flow 20.821164
Reactor Mass: 6.698140629935816
3 Combustor using phi=0.1875 and mass flow 21.888685
Reactor Mass: 6.363181017677751
4 Combustor using phi=0.2125 and mass flow 22.7819795
Reactor 4 failed
Reactor Mass: 5706937067178.77
5 Combustor using phi=0.2375 and mass flow 23.475795
Reactor 5 failed
Reactor Mass: -503360543931.0006
6 Combustor using phi=0.2625 and mass flow 23.950038000000003
Reactor Mass: 5.5565265127696755
7 Combustor using phi=0.2875 and mass flow 24.190739999999998
Reactor 7 failed
Reactor Mass: -1743850606.9494267
8 Combustor using phi=0.3125 and mass flow 24.190739999999998
Reactor Mass: 5.139242033695134
9 Combustor using phi=0.3375 and mass flow 23.950038000000003
Reactor Mass: 4.957399144199033
10 Combustor using phi=0.3625 and mass flow 23.475795
Reactor 10 failed
Reactor Mass: -1583462983.9286568
11 Combustor using phi=0.3875 and mass flow 22.7819795
Reactor 11 failed
Reactor Mass: -1327952461.4260874
12 Combustor using phi=0.4125 and mass flow 21.888685
Reactor Mass: 4.4945519221887915
13 Combustor using phi=0.4375 and mass flow 20.821164
Reactor 13 failed
Reactor Mass: 1025354594.860802
14 Combustor using phi=0.4625 and mass flow 19.608631
Reactor 14 failed
Reactor Mass: -6298273.947422755
15 Combustor using phi=0.4875 and mass flow 18.2829675
Reactor 15 failed
Reactor Mass: -2568280.108493603

System information

• Cantera version: 2.5
• OS: Windows 10
• Python software version: 3.8.5

[bryanwweber]

Hi! Thank you for reporting this. Although the errors and failures here are suboptimal, they aren't totally unexpected. The *ConstPressureReactor classes adjust the volume of the reactor to keep the pressure constant. On the other hand, the *Reactor classes are constant volume, so the outlet mass flow rate adjusts to keep the pressure constant. In general, you want to use the latter class when you are modelling flow through a device. The most robust setup that I've found is: Reservoir->MassFlowController->Reactor->PressureController->Reservoir, with the master attribute of the PressureController set to the MassFlowController instance. You can use an IdealGasReactor instead of Reactor if the ideal gas model is appropriate. Hope that helps!

[speth]

Should we consider using a constant-pressure reactor with both inlets and outlets an error? I'm trying to think of a case where a user would want the behavior where the reactor volume is varying wildly, as opposed to the behavior from the setup using a PressureController.

[bryanwweber]

I'm not in favor of an error, or at least, not in favor of an error that can't be disabled by users. Just because we can't think of a case where this would be useful doesn't mean there aren't any. It does seem like this comes up quite frequently though, so we should have some better user interface here.

[ischoegl]

My 2 cents are that a ConstPressureReactor in combination with a PressureController combines two approaches where each has a mechanism to hold the pressure constant (the former adjusts volume, the latter adjusts mass flow rates). Using both at the same time is asking for trouble (not just numerically, where negative mass flow rates are indeed problematic; I also cannot think of an experimental setup where this would work). For the constant pressure reactor scenario, you'd have to imagine a cylinder with a frictionless piston and two ports: the inlet has a mass flow controller and the outlet has a pressure controller. This is imho an ill defined problem, and in practice you'll have two controllers 'fighting' with each other. So myself, I'd be in favor of not allowing for the configuration that led to this issue report.

[bryanwweber]

I also cannot think of an experimental setup where this would work

What I'm saying is that it doesn't matter if there's an experimental setup where this would work. IMO what matters is that there is an experimental setup or a thought experiment where a ConstPressureReactor and a PressureController would be an appropriate model, which doesn't necessarily correspond to any physical components. In that case, IMO, users should have an escape hatch to say, no I really know what I'm doing and I'm willing to deal with the consequences (e.g., numerical instability).

ETA: That said, clearly the number of messages we get about this specific problem here and on the Users' Group warrants some sort of response. I just want to temper it down from a "thou shalt not do this!" response 😄