1.4 Reactor design problem (Hart et al., 2017; Bequette, 2003):

In this example, we will consider a chemical reactor designed to produce product B from reactant A using a reaction scheme known as the Van de Vusse reaction:

\[A^{\underrightarrow{k_{1}}}B^{\underrightarrow{k_{2}}}C\]

\[2A^{\underrightarrow{k_{3}}}D\]

Under appropriate assumptions, \(F\) is the volumetric flowrate through the tank. The concentation of component \(A\) in the feed is \(c_{Af}\) , and the concentrations in the reactor are equivalent to the concentrations of each component flowing out of the reactor, given by \(c_{A}, c_{B}, c_{C}, c_{D}\)

If the reactor is too small, we will not produce sufficient quantity of \(B\), and if the reactor is too large, much of \(B\) will be further reacted to form the undesired product \(C\). Therefore, our goal is to solve for the reactor volume that maximizes the outlet concentration for product \(B\).

The steady-state mole balances for each of the four components are given by,

\[0 = \frac{F}{V}c_{Af} - \frac{F}{V}c_{A} - k_{1}c_{A} - 2k_{3}c_{A}^{2}\]

\[0 = -\frac{F}{V}c_{B} + k_{1}c_{A} - k_{2}c_{B}\]

\[0 = -\frac{F}{V}c_{C} + k_{2}c_{B}\]

\[0 = -\frac{F}{V}c_{D} + k_{3}c_{A}^{2}\]

The known parameters for the system are,

\[c_{Af} = 10 \frac{gmol}{m^3} \;\;\;\; k_{1} = \frac{5}{6} min^{-1} \;\;\;\; k_{2} = \frac{5}{3} min^{-1} \;\;\;\; k_{3} = \frac{1}{6000} \frac{m^{3}}{mol\;min}\]

Below we formulate and solve this optimization problem using Pyomo. Since the volumetric flowrate \(F\) always appears as the numerator over the reactor volume \(V\) , it is common to consider this ratio as a single variable, called the space-velocity \(SV\).

import pyomo.environ as pyo

# create the concrete model
model = pyo.ConcreteModel()

# set the data (native python data)
k1 = 5.0/6.0     # min^-1
k2 = 5.0/3.0     # min^-1
k3 = 1.0/6000.0  # m^3/(gmol min)
caf = 10000.0    # gmol/m^3

# create the variables
model.sv = pyo.Var(initialize = 1.0, within=pyo.PositiveReals)
model.ca = pyo.Var(initialize = 5000.0, within=pyo.PositiveReals)
model.cb = pyo.Var(initialize = 2000.0, within=pyo.PositiveReals)
model.cc = pyo.Var(initialize = 2000.0, within=pyo.PositiveReals)
model.cd = pyo.Var(initialize = 1000.0, within=pyo.PositiveReals)

# create the objective
model.obj = pyo.Objective(expr = model.cb, sense=pyo.maximize)

# create the constraints
model.ca_bal = pyo.Constraint(expr = (0 == model.sv * caf \
                 - model.sv * model.ca - k1 * model.ca \
                 -  2.0 * k3 * model.ca ** 2.0))

model.cb_bal = pyo.Constraint(expr=(0 == -model.sv * model.cb \
                 + k1 * model.ca - k2 * model.cb))

model.cc_bal = pyo.Constraint(expr=(0 == -model.sv * model.cc \
                 + k2 * model.cb))

model.cd_bal = pyo.Constraint(expr=(0 == -model.sv * model.cd \
                 + k3 * model.ca ** 2.0))

solver = pyo.SolverFactory('ipopt')

solver.solve(model)
model.pprint()

5 Var Declarations
    ca : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 3874.2588672317133 :  None : False : False : PositiveReals
    cb : Size=1, Index=None
        Key  : Lower : Value             : Upper : Fixed : Stale : Domain
        None :     0 : 1072.437200108632 :  None : False : False : PositiveReals
    cc : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 1330.0935334088806 :  None : False : False : PositiveReals
    cd : Size=1, Index=None
        Key  : Lower : Value             : Upper : Fixed : Stale : Domain
        None :     0 : 1861.605199625387 :  None : False : False : PositiveReals
    sv : Size=1, Index=None
        Key  : Lower : Value              : Upper : Fixed : Stale : Domain
        None :     0 : 1.3438117610672782 :  None : False : False : PositiveReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize :         cb

4 Constraint Declarations
    ca_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                                       : Upper : Active
        None :   0.0 : 10000.0*sv - sv*ca - 0.8333333333333334*ca - 0.0003333333333333333*ca**2.0 :   0.0 :   True
    cb_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                    : Upper : Active
        None :   0.0 : - sv*cb + 0.8333333333333334*ca - 1.6666666666666667*cb :   0.0 :   True
    cc_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                            : Upper : Active
        None :   0.0 : - sv*cc + 1.6666666666666667*cb :   0.0 :   True
    cd_bal : Size=1, Index=None, Active=True
        Key  : Lower : Body                                     : Upper : Active
        None :   0.0 : - sv*cd + 0.00016666666666666666*ca**2.0 :   0.0 :   True

10 Declarations: sv ca cb cc cd obj ca_bal cb_bal cc_bal cd_bal