Tutorial 4 – Dynamic Flowsheets

Introduction

Up until now, we have only looked at steady-state flowsheets, so in this tutorial we will look at how to build a dynamic flowsheet. To do this, we will modify the process from Tutorial 1 to simulate dynamic behavior by including:

  • accumulation and holdup,
  • a Mixer for the different feed components,
  • pressure driven flow leaving each Tank, and
  • a step-change in feed rates.

This tutorial will teach you how to:

  • Set up a dynamic flowsheet,
  • Add steady-state models to dynamic flowsheets,
  • Add constraints to existing unit models,
  • Transform the time domain using Pyomo.dae,
  • set steady-state initial conditions, and
  • plot some results.

First Steps

Firstly, we need to import the necessary libraries we will use for our Flowsheet. In addition to the imports from Tutorial 1, we need to import Constraint, Var and TransformationFactory.

from pyomo.environ import ConcreteModel, Constraint, SolverFactory, TransformationFactory, Var

from idaes.core import FlowsheetBlock, Stream

import tutorial_1_properties as properties_rxn
from idaes.models import CSTR, Mixer

We also need to import a library to use to plot our results, for which we will use matplotlib:

import matplotlib.pyplot as plt

Dynamic Flowsheets

In order to create a dynamic model for a flowsheet, we start with a ConcreteModel as before.

m = ConcreteModel()

When we add our FlowsheetBlock however, we need to indicate that this will be a dynamic flowsheet. We have already seen that FlowsheetBlock, which we have so far set to False (this is the default setting by the way, so it is not necessary to specify dynamic = False for steady-state flowsheets). For a dynamic Flowsheet, we simply need to specify dynamic = True when creating our FlowsheetBlock.

FlowsheetBlock also accepts another argument named time_set that can be used to specify a set of points within the time domain. The time domain needs to be set with a start and end point at a minimum, and the user can specify additional intermediate points if desired. time_set expects a list of points, and has the following defaults:

  • if dynamic = True, time_set = [0.0, 1.0]
  • if dynamic = False, time_set = [0.0]

For our flowsheet, let us add the following points to our time domain:

  • Start point = 0.0 s
  • End point = 5e5 s
  • Internal point at 1.0 s

Thus, we write the following to create our FlowsheetBlock:

m.fs = FlowsheetBlock(dynamic=True, time_set=[0, 1, 500000])

Adding Property Packages and Models

Dynamic FlowsheetBlocks behave the same way as steady-state FlowsheetBlocks, so we can add our property packages and unit models in the same manner as before. First, we need to create the property package, and we will also set it as the default property package for our flowsheet.

m.fs.properties_rxn = properties_rxn.PropertyParameterBlock()
m.fs.config.default_property_package = m.fs.properties_rxn

Whilst we have declared our flowsheet to be a dynamic process, there are many circumstances where we would have unit operations which react sufficiently quickly that we would like to consider them to be steady-state operations. The IDAES framework supports this by allowing us to declare Unit Models as steady-state model by specifying dynamic = False when declaring Unit Models (note that the reverse is not supported, i.e. dynamic models in steady-state flowsheets). If the dynamic argument is not provided when a Unit Model is added to a Flowsheet, it takes this argument from its parent object. This can also be applied to sub-flowsheets, allowing for steady-state sub-Flowsheets within dynamic Flowsheets.

For this tutorial, let us add a Mixer unit at the beginning of our process to mix two feed streams prior to being fed to Tank1. We will assume this is an ideal mixer, which means that we can assume it can be represented by a steady-state model.

m.fs.Mix = Mixer(dynamic=False)

Next, we can add our two Tank models as before. As we have set a default property package and we want our CSTRs to be dynamic models, we do not need to provide any arguments when adding these to our Flowsheet (the dynamic argument will be inherited from the Flowsheet).

m.fs.Tank1 = CSTR()
m.fs.Tank2 = CSTR()

Adding Variables and Constraints to Existing Models

For our dynamic process, we want to model the flow of material leaving each tank as being pressure driven flow; that is flowrate leaving each tank is proportional to the depth of fluid in each tank. However, the core CSTR model does not include the calculations necessary for modeling this behavior. Thus, we will need to add these Constraints ourselves.

First, we will need to add some additional Variables to our Tank models that we will use in our calculations. We can do this by simply adding Pyomo Var objects to the CSTR objects already constructed in our Flowsheet. We need to add the following Variables to Tank1 (Note: volume_flow is a state variable and thus technically belongs in a Property Block. As we have not covered these yet, we will place it in the Unit Model for now).

m.fs.Tank1.height = Var(m.fs.Tank1.time,
                        initialize=1.0,
                        doc="Depth of fluid in tank [m]")
m.fs.Tank1.area = Var(initialize=1.0,
                      doc="Cross-sectional area of tank [m^2]")
m.fs.Tank1.volume_flow = Var(m.fs.Tank1.time,
                             initialize=4.2e5,
                             doc="Volumetric flow leaving tank")
m.fs.Tank1.flow_coeff = Var(m.fs.Tank1.time,
                            initialize=5e-5,
                            doc="Tank outlet flow coefficient")

Next, we can add the extra Constraints we need. This is done using Pyomo Constraint objects, which again can be added directly to our CSTR objects. To create a Constraint, we first need to write a rule describing the equation we wish to write, and to then create a Constraint using this rule. Let us start with a Constraint relating the volume of reacting fluid, \(V_t\), to the cross-sectional area of the tank, \(A\), and the depth of the fluid, \(h_t\), at a given point in time \(t\) (we will assume \(A\) is constant with time):

\[V_t = A \times h_t\]

To do this, we first write the following rule:

def geometry(b, t):
    return b.volume[t] == b.area*b.height[t]

This is a Python method that returns an expression involving Pyomo objects (volume, area and height). In this rule we have used b to represent the object to which we are adding the Constraint (and from which we will gather these variables). When we create the Constraint in the next step, Pyomo automatically passes the object that the Constraint is being added to (in this case Tank1) to the rule, allowing us to write the rule in a general form.

We then use this rule to construct the actual Constraint in our model.

m.fs.Tank1.geometry = Constraint(m.fs.Tank1.time, rule=geometry)

We also need to write the following Constraints for volume_flow:

\[F_{vol,t} = F_{mol,t} / \rho_{mol,t}\]
\[F_{vol,t} = C_t \times h_t\]

where \(F_{vol,t}\) is volume_flow, \(F_{mol,t}\) is the total molar flowrate of material leaving the tank, \(\rho_{mol,t}\) is the total density of the fluid in the tank and \(C_t\) is flow_coeff. \(F_{mol,t}\) and \(\rho_{mol,t}\) both come from the Property Block associated with the fluid in the tank (and thus the tank outlet), which is called Tank1.holdup.properties_out. To write these Constraints, we use the following code:

def volume_flow_calculation(b, t):
    return b.volume_flow[t] == (
                b.holdup.properties_out[t].flow_mol /
                b.holdup.properties_out[t].dens_mol_mol['Liq'])
m.fs.Tank1.volume_flow_calculation = Constraint(
        m.fs.Tank1.time,
        rule=volume_flow_calculation)

def outlet_flowrate(b, t):
    return b.volume_flow[t] == b.flow_coeff[t]*b.height[t]
m.fs.Tank1.outlet_flowrate = Constraint(m.fs.Tank1.time,
                                        rule=outlet_flowrate)

We then need to repeat this for Tank2.

m.fs.Tank2.height = Var(m.fs.Tank2.time,
                        initialize=1.0,
                        doc="Depth of fluid in tank [m]")
m.fs.Tank2.area = Var(initialize=1.0,
                      doc="Cross-sectional area of tank [m^2]")
m.fs.Tank2.volume_flow = Var(m.fs.Tank2.time,
                             initialize=4.2e5,
                             doc="Volumetric flow leaving tank")
m.fs.Tank2.flow_coeff = Var(m.fs.Tank2.time,
                            initialize=5e-5,
                            doc="Tank outlet flow coefficient")

When it comes to adding the Constraints to Tank2, we can reuse the rules we wrote for Tank1, as we wrote these in a general form. When we create the Constraints in Tank2 now, the Variables in Tank2 will be used instead.

m.fs.Tank2.geometry = Constraint(m.fs.Tank2.time, rule=geometry)
m.fs.Tank2.volume_flow_calculation = Constraint(
        m.fs.Tank2.time,
        rule=volume_flow_calculation)
m.fs.Tank2.outlet_flowrate = Constraint(m.fs.Tank2.time,
                                        rule=outlet_flowrate)

Transforming the Time Domain

At this point in previous tutorials, we would now call post_transform_build to finish constructing our models. However, in dynamic flowsheets, we first need to transform the differential equations into a form that can be handled by our solver. When the time domain is created for a dynamic flowsheet, it is created as a Pyomo ContinuousSet object and the associated derivatives as DerivativeVars. However, most solvers do not understand these types of objects so they need to be transformed prior to solving. This is done using Pyomo’s TransformationFactory; for more details on what happens during DAE transformation, see the Pyomo documentation.

Pyomo provides support for a number of types of DAE transformation:

  • finite difference methods - backward, forward and central difference methods, and
  • orthogonal collocation methods - Lagrange-Radau and Lagrange-Legendre roots.

For time domains, we generally use a 1st order backward finite difference method.

To use Pyomo’s DAE Transformation, we first need to create a Transformation object, and to then apply it to our model object.

discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m.fs,
                     nfe=200,
                     wrt=m.fs.time,
                     scheme='BACKWARD')

Once the DAE transformation has been applied to the time domain, we can call post_transform_build to finish the model construction.

m.fs.post_transform_build()

Connecting Units

As before, we cannot start connecting units together until post_transform_build has been called. Now that that is done, we can add Streams as necessary:

m.fs.stream_1 = Stream(source=m.fs.Mix.outlet,
                       destination=m.fs.Tank1.inlet)

m.fs.stream_2 = Stream(source=m.fs.Tank1.outlet,
                       destination=m.fs.Tank2.inlet)

Setting Design and Operating Constraints

For this tutorial, let us split our feed into two parts (with the same temperature and pressure):

  • a stream with 10.0 mol/s of “a” and 1.0 mol/s of “c”
  • a stream with 20.0 mol/s of “b”

As we did not add Feed Blocks to this flowsheet, we will fix these through the inlet Port objects of the Mixer. Note that for a Mixer, there is a single inlet object which is indexed by time and inlet name. Thus, when we fix the feed conditions, we need to do so at all points in time for each inlet. We can do this using slice notation as shown below:

m.fs.Mix.inlet[:, "1"].vars["flow_mol_comp"]["a"].fix(10.0)
m.fs.Mix.inlet[:, "1"].vars["flow_mol_comp"]["b"].fix(0.0)
m.fs.Mix.inlet[:, "1"].vars["flow_mol_comp"]["c"].fix(1.0)
m.fs.Mix.inlet[:, "1"].vars["flow_mol_comp"]["d"].fix(0.0)
m.fs.Mix.inlet[:, "1"].vars["flow_mol_comp"]["e"].fix(0.0)
m.fs.Mix.inlet[:, "1"].vars["flow_mol_comp"]["f"].fix(0.0)
m.fs.Mix.inlet[:, "1"].vars["temperature"].fix(303.15)
m.fs.Mix.inlet[:, "1"].vars["pressure"].fix(101325.0)

m.fs.Mix.inlet[:, "2"].vars["flow_mol_comp"]["a"].fix(0.0)
m.fs.Mix.inlet[:, "2"].vars["flow_mol_comp"]["b"].fix(20.0)
m.fs.Mix.inlet[:, "2"].vars["flow_mol_comp"]["c"].fix(0.0)
m.fs.Mix.inlet[:, "2"].vars["flow_mol_comp"]["d"].fix(0.0)
m.fs.Mix.inlet[:, "2"].vars["flow_mol_comp"]["e"].fix(0.0)
m.fs.Mix.inlet[:, "2"].vars["flow_mol_comp"]["f"].fix(0.0)
m.fs.Mix.inlet[:, "2"].vars["temperature"].fix(303.15)
m.fs.Mix.inlet[:, "2"].vars["pressure"].fix(101325.0)

We also need to fix values for area, flow_coeff and heat duty in both Tanks.

m.fs.Tank1.area.fix(0.5)
m.fs.Tank1.flow_coeff.fix(5e-6)
m.fs.Tank1.heat.fix(0.0)

m.fs.Tank2.area.fix(0.5)
m.fs.Tank2.flow_coeff.fix(5e-6)
m.fs.Tank2.heat.fix(0.0)

Setting Initial Conditions

When setting up dynamic models, we also need to provide initial conditions for the system. One approach is to specify that the system is at steady-state initially (i.e. all time derivatives are equal to zero at the starting time). In order to simplify specifying this, IDAES FlowsheetBlocks have the following method:

m.fs.fix_initial_conditions('steady-state')

Steady-state initial conditions are often only useful for modeling simple systems under ideal circumstances, and users should choose a set of initial conditions that suit their actual circumstances.

Initializing the Flowsheet

The best method for initializing dynamic problems is an open question, and it is often difficult to find a good way to do this. For this tutorial, it is sufficient to initialize the entire model at a single, steady-state operating condition and to then introduce a time disturbance later. Thus, we can just use the same approach that we have used previously to initialize the flowsheet.

m.fs.Mix.initialize()
m.fs.Tank1.initialize(state_args={
        "flow_mol_comp": {
            "a": m.fs.Mix.outlet[0].vars["flow_mol_comp"]["a"].value,
            "b": m.fs.Mix.outlet[0].vars["flow_mol_comp"]["b"].value,
            "c": m.fs.Mix.outlet[0].vars["flow_mol_comp"]["c"].value,
            "d": m.fs.Mix.outlet[0].vars["flow_mol_comp"]["d"].value,
            "e": m.fs.Mix.outlet[0].vars["flow_mol_comp"]["e"].value,
            "f": m.fs.Mix.outlet[0].vars["flow_mol_comp"]["f"].value},
        "pressure": m.fs.Tank1.outlet[0].vars["pressure"].value,
        "temperature": m.fs.Tank1.outlet[0].vars["temperature"].value})
m.fs.Tank2.initialize(state_args={
        "flow_mol_comp": {
            "a": m.fs.Tank1.outlet[0].vars["flow_mol_comp"]["a"].value,
            "b": m.fs.Tank1.outlet[0].vars["flow_mol_comp"]["b"].value,
            "c": m.fs.Tank1.outlet[0].vars["flow_mol_comp"]["c"].value,
            "d": m.fs.Tank1.outlet[0].vars["flow_mol_comp"]["d"].value,
            "e": m.fs.Tank1.outlet[0].vars["flow_mol_comp"]["e"].value,
            "f": m.fs.Tank1.outlet[0].vars["flow_mol_comp"]["f"].value},
        "pressure": m.fs.Tank1.outlet[0].vars["pressure"].value,
        "temperature": m.fs.Tank1.outlet[0].vars["temperature"].value})

Before moving on, we should also solve the entire flowsheet at our initial conditions, as the Constraints within the Streams are not yet initialized (even though the values on both sides should be the same).

solver = SolverFactory('ipopt')
results = solver.solve(m, tee=True)

Creating a Disturbance

Now that our model is initialized, let’s create a disturbance in the system so that we can observe the response. For now, let us halve the flowrate of component “b” at time 1, which is done as follows:

for t in m.fs.time:
    if t >= 1.0:
        m.fs.Mix.inlet[t, "2"].vars["flow_mol_comp"]["b"].fix(10.0)

Then, we can solve the model with the new disturbance and print the results so that we can see if the problem converged.

results = solver.solve(m, tee=True)
print(results)

Hopefully you will see that the solver found an optimal solution, and that the problem had 40588 variables and constraints.

Plotting the Results

The best way to observe the results of our disturbance is to plot the flowrates of the different components leaving Tank2. To do this, we first need to collect the time and flow information in a form that can be understood by our plotting library.

To begin with, we create an empty list to hold the values of each variable we wish to plot (and the time domain).

time = []
a = []
b = []
c = []
d = []
e = []
f = []

Next, we iterate over the time domain, and append the value of each variable to the appropriate list.

for t in m.fs.time:
    time.append(t)
    a.append(m.fs.Tank2.outlet[t].vars["flow_mol_comp"]["a"].value)
    b.append(m.fs.Tank2.outlet[t].vars["flow_mol_comp"]["b"].value)
    c.append(m.fs.Tank2.outlet[t].vars["flow_mol_comp"]["c"].value)
    d.append(m.fs.Tank2.outlet[t].vars["flow_mol_comp"]["d"].value)
    e.append(m.fs.Tank2.outlet[t].vars["flow_mol_comp"]["e"].value)
    f.append(m.fs.Tank2.outlet[t].vars["flow_mol_comp"]["f"].value)

Next, we need to create a plot object, and to add each curve that we wish to plot.

plt.figure(1)
plt.plot(time, a, label='a')
plt.plot(time, b, label='b')
plt.plot(time, c, label='c')
plt.plot(time, d, label='d')
plt.plot(time, e, label='e')
plt.plot(time, f, label='f')

We can also add legends and axis labels as follows:

plt.legend()
plt.grid()
plt.xlabel("Time [s]")
plt.ylabel("Molar Flowrate [mol/s]")

Finally, we need to display the plot object so that we can see it.

plt.show(block=True)

If all goes well, you should see something like this:

../../_images/tutorial_4_figure.png