Modelling a reaction with variable stoichiometric coefficients

[tk3016]

Hi,

I am trying to model a reactor for Fischer-Tropsch synthesis using a PFR (ControlVolume1D) unit.
The reaction package that I am trying to implement does not have fixed stoichiometric coefficients, i.e., they are also variables.
Can you suggest a workaround for modelling this?

From what I can understand by looking at the source code of ControlVolume1D (c.f. Lines 751-812) is that I can only define stoichiometric coefficients as parameters.

Hence, I wanted to ask about a workaround without having to make too many changes to the original code myself.

Looking forward to hearing back from you.

Best wishes,
Tanuj

[dallan-keylogic]

You say the stoichiometric coefficients are variables---are they considered parameters (but not Param objects) that need to be identified, or are they functions of local properties like temperature, pressure, and composition? The fix may look different depending on which it is.

A quick fix is to manually delete the constraint associated with the reaction in the control volume, then replace it with a modified version with variable coefficients. For a rate reaction, the constraint you need to replace is

@self.Constraint(
    self.flowsheet().time,
    self.length_domain,
    pc_set,
    doc="Kinetic reaction stoichiometry constraint",
)
def rate_reaction_stoichiometry_constraint(b, t, x, p, j):
    if (
        b.config.transformation_scheme != "FORWARD"
        and x == b.length_domain.first()
    ) or (
        b.config.transformation_scheme == "FORWARD"
        and x == b.length_domain.last()
    ):
        return Constraint.Skip
    elif (p, j) in pc_set:
        rparam = rblock[t, x].config.parameters
        return b.rate_reaction_generation[t, x, p, j] == (
            sum(
                rparam.rate_reaction_stoichiometry[r, p, j]
                * b.rate_reaction_extent[t, x, r]
                for r in b.config.reaction_package.rate_reaction_idx
            )
        )
    else:
        return Constraint.Skip

You delete that by writing path.to.CV.control_volume.del_component("rate_reaction_stoichiometry_constraint"), then create a new constraint:

@path.to.CV.control_volume.Constraint(
     path.to.CV.control_volume.flowsheet().time,
     path.to.CV.control_volume.length_domain,
     pc_set,
     doc="whatever"
)
def rate_reaction_stoichiometry_constraint(b, t, x, p, j):
    if (
        b.config.transformation_scheme != "FORWARD"
        and x == b.length_domain.first()
    ) or (
        b.config.transformation_scheme == "FORWARD"
        and x == b.length_domain.last()
    ):
        return Constraint.Skip
    elif (p, j) in pc_set:
        return b.rate_reaction_generation[t, x, p, j] == whatever
    else:
        return Constraint.Skip

Note that using existing initialization methods will require the stoichiometric coefficient Var to be fixed and whatever constraint is used in its correlation to be deactivated.

[tk3016]

Hi @dallan-keylogic,

Many thanks for your prompt and comprehensive response - much appreciated.
The stoichiometric coefficients are variables and are a function of temperature inside the reactor.
If you need more context, please see this publication.

Here, the stoichiometric coefficient is a function of the chain length probability (the amount of Cn+ compounds produced in the reaction), which in turn is a function of the temperature inside the reactor.

Is the workaround you suggested still valid for this case?

Best wishes,
Tanuj

[dallan-keylogic]

Thank you for the reference. Yes, I would still use the workaround I suggested. There was a simpler method available if the stoichiometric coefficient was an uncertain global parameter that needed to be estimated by creating a Var on the reaction parameter block and unfixing it, but since it's dependent on local conditions that method doesn't work.

[tk3016]

Thanks for the clarification @dallan-keylogic.
I will try this and may get back to you if I need further assistance (let's hope not!)

Best wishes,
Tanuj