Adding new Reactors

This section explains how to add new reactors to the framework.

In this section, first a general note about reactors in this package is given, followed by a detailed description on how to add a new reactor.

About Reactors

First, before we go into the process of defining new reactors, we will have a look into what we consider to be a Reactor in this framework and what possible elements it comprises:

Thus, a reactor contains only information about the mass balance over it and is independent from the processes taking place within it. Instead, it gets the already computed reaction rates from a Process system, to which it is to be connected via the states and the rates ports. This has the advantage, that the reactor is independent of the processes in it and thus it allows to define the reactor only once for all processes and even to combine different processes in a single reactor by adding up the individual process rates before connecting them to the reactor.

For it's computations, each Reactor has the following properties (some are optional):

• A name
• (Multiple) inflow ports (can be empty)
• (Multiple) outflow ports (can be empty)
• A rates port
• A states port
• equations describing the differential relation between the inflow and outflow depending on the rates
• t: The independent variable for time
• parameters for the reactor (can be empty)

Adding a new Reactor

Adding a new reactor has to be done by writing Julia code which generates a Reactor struct which has all neccessary fields defined.

To do so, four steps are required:

1. Creating connectors for the inflow and states inputs and for the outflow and rates outputs
2. Generating symbols for the parameters
3. Specifying the differential equations
4. Generating a Reactor struct with all this information

To be able to reuse it and to make the usage in building a complex model more clear, this should be done in a (properly-named) function which takes the name, desired states, possible initial values for the states and values for the parameters as inputs and returns the resulting Reactor struct.

In the following, the process of adding a new Reactor is explained on the example of a continuously strirred tank reactor (CSTR). A CSTR is defined by the following set of equations:

\[\begin{align*} q_{in} &= q_{out}\\ \frac{dC_i}{dt} &= q_{in} * C_{i, in} / V - q_{out} * C_i / V + r_{C_i} \qquad \text{for concentrations } C_i \text{ of every desired compound}\\ C_{i, out} &= C_i \qquad \text{for concentrations } C_i \text{ of every desired compound} \end{align*}\]

Hereby, $V$ is the reactor volume, $q_{in}$ and $q_{out}$ are the incoming, resp. outflowing, flowrates, $C_{i, in}$ is the concentration of compound $i$ in the inflow and $C_{i}$ the concentration in the reactor. Further $r_{C_i}$ is the reaction rate for the compound $i$, which is obtained externally.

In the process of adding the new reactor, the biggest challenge will be to keep exactly this independence of the actual compounds in the flow. However, this is key to using the same reactor for arbitrary processes and combinations therof. How this can be achieved is explained in the following, procedure, in the step where the differential equations are specified.

The function header

First, the function header has to be built. It should take the volume and name of the reactor as well as the states for the compounds and optionally the initial values for the states.

@reactor function CSTR(V, states; initial_states = nothing, name)

end 

Here the volume V and the states are used as required positional arguments, the initial_states and the name are keyword arguments, with the initial_states having a default of nothing, which means that a default set of initial values will be picked.

The @reactor macro has been used on the function to mark it as reactor, which leads to automatic generation of a function for calling the reactor with one (or multiple) process (e.g. CSTR(1000, ASM1(); name=:react)) without caring about the connections in-between. For this, it is important that the reactor argument for supplying the states is called states and that the name is provided by a keyword argument called name (as in the here-presented example).

Creating the connectors

Then, we need to create the connectors for inflow and outflow as well as rates and states. To do so, we use the InflowPort, OutflowPort, ReactionInputPort and StateOutputPort functions. First we start with the inflow and call the InflowPort function:

 @named inflow = InflowPort(states)

Here we provide the states as input, just as we got them as input. This allows that the handling and parsing of it is unified in the ports, which allow for multiple inputs. However, the most convenient is to get it from the process that is used with using the unknowns function. Further, this function would also need a name but instead of explicitly specifying it, we used the convenient @named macro here, which automatically assignes the variable name (here influent) as name.

Then, we do it equivalently for the outflow and rates ports:

@named outflow = OutflowPort(states)
@named rates = ReactionInputPort(states)

For the remaining states port we use the StateOutputPort function:

@named states = StateOutputPort(states; initial_states, set_guess = false)

This function is additionally provided with keyword arguments for the initial_states (which are forwarded from the input of this function) and a set_guess argument, which is set to false here. This additional argument tells the StateOutputPort how to specify the initial conditions, the possibilities are either as strict initial conditions, or only as guesses which do not have to be followed strictly. By default it uses the setting as guess method, as providing too many fixed initial conditions can result in contradicting initial conditions and also if they match it would lead to an overdetermined systems of initial values. However, here we want the initial_states to be followed strictly (and as they are associated with differential equations which need an initial condition, it is save), thus we pass set_guess = false here.

Generating symbols for the parameters

To generate the symbols for the parameters we can use the @parameters macro:

parameters = @parameters V=V

Here, the parameter for the reactor volumen called V is created and the current value of the variable V is set as its default value. Further the parameter is from now on available under the variable V (overwriting the default value, as it is already set and not needed anymore after here, this is not a problem) and additionally the variable parameters now is a vector with all parameters specified in this line.

  parameters = @parameters param1=param1_value param2=param2_value

Specifying the equations

Now, as we have all ports and parameters, we can specify the equations. First we start with the equations independent of the compounds, here the one which sets the flowrates of in- and outflow to be equal:

eqs = [inflow.q ~ outflow.q]

Here, inflow.q takes the flowrate q variable of the inflow port, and similarly outflow.q designates the flowrate of the outflow. Hereby, the ~ designates the equal sign of the equation (=) and is the symbol for equality in equations within the framework (the normal = is reserved for assignements). If the reactor does not have any equations independent of the compounds, just create an empty vector of equations here (eqs = Equation[]).

Then, for the equations for each compound, we loop over the variables in the states which we get using unknowns(states). This is is then transformed to a Symbol representing the name of the variable using the symbolic_to_namestring function and then converting the resulting string to a symbol:

for sym = Symbol.(symbolic_to_namestring.(unknowns(states)))

end

Within this loop, we can then use the getproperty function to get the corresponding variable of all needed ports (inflow, rates, states, outflow):

for sym = Symbol.(symbolic_to_namestring.(unknowns(states)))
    in = getproperty(inflow, sym) # The compound in the inflow
    r = getproperty(rates, sym) # The reaction rate for the compound
    s = getproperty(states, sym) # The compound in the state
    out = getproperty(outflow, sym) # The compound in the outflow
end

And finally, we can add the equation for each compound to the eqs vector of equations:

for sym = Symbol.(symbolic_to_namestring.(unknowns(states)))
    in = getproperty(inflow, sym)
    r = getproperty(rates, sym)
    s = getproperty(states, sym)
    out = getproperty(outflow, sym)
    push!(eqs, D(s) ~ in * inflow.q / V - s * outflow.q / V + r)
    push!(eqs, out ~ s)
end

Hereby, the operator D is used to get the derivative of the state variable s with respect to time. In this framework, this operator is defined but if you use this code outside of it, you will need to use BioChemicalTreatment.ProcessSimulator.D instead, as it is generally not intended for public use. On the right hand side, exactly the equation written symbolically above is written.

Generating the Reactor

Finally, we can assemble the Reactor struct using everything above:

Reactor(eqs, t, parameters, [inflow], [outflow], rates, states, name=name)

Hereby the arguments are:

• The vector of equations
• The time (as independent variable). The library exports this as t.
• The vector of parameters (if none specify as [])
• The inflow ports (Possibly multiples, which is why it has to be a vector and thus [inflow])
• The outflow ports (Possibly multiples, which is why it has to be a vector and thus [outflow])
• The rates port
• The states port
• The name as keyword argument, here just forwarding the name provided to the function

The resulting function

Putting everything together, we get the following function:

@reactor function CSTR(V, states; initial_states = nothing, name)
    # Get all connectors
    @named inflow = InflowPort(states)
    @named rates = ReactionInputPort(states)
    @named states = StateOutputPort(states; initial_states, set_guess = false)
    @named outflow = OutflowPort(states)
    # The parameter for the Volume
    parameters = @parameters V = V

    # The CSTR equation (inflow rate = outflow rate)
    eqs = [inflow.q ~ outflow.q]
    # For each Component c add the following equations:
    #   dc/dt = c_in * q_in/V - c * q_out/V + r_c (Conservation of mass + reaction term (supplied from connector))
    #   c_out = c (The output is equal to the internal state)
    for sym = Symbol.(symbolic_to_namestring.(unknowns(states)))
        in = getproperty(inflow, sym)
        r = getproperty(rates, sym)
        s = getproperty(states, sym)
        out = getproperty(outflow, sym)
        push!(eqs, D(s) ~ in * inflow.q / V - s * outflow.q / V + r)
        push!(eqs, out ~ s)
    end

    # Compose and return the CSTR system
    Reactor(eqs, t, parameters, [inflow], [outflow], rates, states, name=name)
end

Now, this function can then be used as any other reactor that is already defined in the package.