Convergence of 1D Control Volume unit.

[tk3016]

Hi IDAES community,

I am trying to simulate a fixed-bed catalytic reactor with a pseudo-homogeneous assumption.
I am using the PFR unit in IDAES to do so.
I am defining the PFR unit as follows:
m.fs.R101 = PFR( property_package=m.fs.thermo_params, reaction_package=m.fs.reaction_params, has_equilibrium_reactions=False, has_heat_of_reaction=True, has_heat_transfer=True, has_pressure_change=False, transformation_method='dae.finite_difference', transformation_scheme='BACKWARD', finite_elements=24 )

Even though I have a pressure drop, for now, to keep the model debugging process simpler, I am assuming there is no pressure drop within the reactor.
Initially, I initialise and simulate the reactor by fixing the temperature across the control volume to be the same as the temperature at the inlet of the reactor (isothermal assumption).

The model solves correctly. I do this to improve the initialisation of the model.
Then I unfix the temperature inside the control volume, and I add the following constraint to model the external cooling:

def deltaT_expr(b, z):
    return m.fs.R101.heat_duty[0,z] == m.fs.R101.U * m.fs.R101.area * m.fs.R101.eff_area_expr * (m.fs.R101.coolant_temperature - b.properties[0,z].temperature)
m.fs.R101.control_volume.deltaT_expr = Constraint(
    m.fs.R101.control_volume.length_domain,
    rule=deltaT_expr,
    doc="deltaT")

But on doing so, my model fails to converge regardless of the solver used, and it always complains that the model is infeasible.
I have thoroughly double-checked all the model equations of my reaction package.
Although when I reduce the number of finite elements to 1, the model converges, giving me the qualitatively expected conditions at the outlet of the reactor.

Can someone kindly advise on some ways I can troubleshoot/fix this?

Looking forward to hearing back from you.

Best wishes,
Tanuj

[dallan-keylogic]

Are you initializing the heat duty to some reasonable value? What does your IPOPT outlet look like? Are you scaling any variables or constraints?

I'd suggest using the Diagnostics Toolbox to look at the model both before and after the solve that fails. If the model is infeasible, either 1) you've got at least one variable at a bound, which display_variables_near_bounds should show, or 2) your Jacobian matrix is rank-deficient, which the SVD Toolbox will show.

[tk3016]

Hi @dallan-keylogic,

Thanks for your insightful response as always.
Fortunately, right before you replied, I managed to resolve the issue.
Even then, I will make a note of the tips you provided! :)

For everyone's information (hopefully it helps someone out there in the future) - the issue was that there was a thermal runaway condition which caused the numerics to break down when I used higher number of finite elements.
I verified by solving my model with BARON (a global solver) for NFE=1 to NFE=5.
I noticed the hotspot temperature increasing quickly on increasing the NFE and for NFE > 5 the model became infeasible.

Accordingly, I finetuned my reaction package and it works well.
Regarding my previous question about variable stoichiometric coefficients, I figured out a simpler workaround.
I introduced dummy reactions where the effective rates for the dummy reactions where the variable stoichiometric coefficient times the original reaction.
I calculate the variable stoichiometric coefficients in my reaction package.

I hope this helps someone in the future.
I am very grateful for the active community here and all the quick help that I get.

Best wishes,
Tanuj