Disinfection with Ozone

As a small example, we consider the following code, which implements a simple connection of two continuously stirred tank reactors (CSTR) as commonly used in water treatment. It serves as a quick example to start with. For more detailed explanations, we referred to the API of the ProcessSimulator and the more advanced examples.

This example implements two reactors used for ozonation of water. In the first reactor, the concentration of ozone is kept constant, to simulate a supply of ozone in the beginning. The stoichiometric matrix for the ozonation process is defined in the library and accessible under the name OzoneDisinfection. This system is depicted as follows:

For building up this system, it is proceeded as follows:

Plug and Play Example

using BioChemicalTreatment # Reactors etc.
using ModelingToolkit # Modeling Framework
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Define the needed constants
v = 1000.0/2    # Reactor volume
q = 10000.0/24  # Flow rate through the reactor
kO3 = 10.0/24  # Reaction rate for the decay of ozone
kd = 1500.0/24  # Reaction rate for the decay of the bacteria

# Set default process
@set_process disinfection = OzoneDisinfection(;kO3, kd)

# The reactors
@named CSTR1 = CSTR(v, OzoneDisinfection(;kO3=0, kd, name=:disin_supply); initial_states = [0.5, 0]) # Override the default process!
@named CSTR2 = CSTR(v; initial_states = [0.5, 0])

# The influent
@named influent = Influent([0.5, 1], flowrate=q)

# The connections
connection_eqns = [
	# Connect the flows
	connect(outflows(influent)[1], inflows(CSTR1)[1]), # Influent -> CSTR1
	connect(outflows(CSTR1)[1], inflows(CSTR2)[1]), # CSTR1 -> CSTR2
]

# Building the overall system
@named system = ODESystem(connection_eqns, t, systems = [
	influent,
	CSTR1,
	CSTR2,
])
system_simplified = structural_simplify(system)

## Simulate the system (for 1 day)
# First create a ODEProblem for the time span (0, 1) and then solve it
prob = ODEProblem(system_simplified, [], (0, 1), [])
sol = solve(prob)

## Plot the output
plot(sol,
	idxs = [states(CSTR1).X_B ,states(CSTR2).X_B], # Plot bacteria in tank
	legend = :right,
	title = "Remaining fraction of bacteria in each tank")

Step-by-Step Explanations

Setup

First, the needed packages are imported and constants defined:

using BioChemicalTreatment # Reactors etc.
using ModelingToolkit # Modeling Framework
using DifferentialEquations # Solve the equations
using Plots # Plotting

# Define the needed constants
v = 1000.0/2    # Reactor volume
q = 10000.0/24  # Flow rate through the reactor
kO3 = 10.0/24  # Reaction rate for the decay of ozone
kd = 1500.0/24  # Reaction rate for the decay of the bacteria

System Creation

Then, all needed systems (here inflow and the reactors) are created.

We start with setting the default reaction rates system for the disinfection dynamics. This system defines all compounds needed for this specific reaction, which are needed afterwards for defining the CSTRs.

@set_process disinfection = OzoneDisinfection(;kO3, kd)

We continue with the two CSTRs. But wait: The first reactor is supposed to have zero ozone decomposition rate (to simulate the supply, as explained above), but the default ozonation process set above has it fixed to 10/24! Thus this process must be overridden. This can be done by supplying a new process directly to the reactor and setting use_default_processes=false:

@named CSTR1 = CSTR(v, OzoneDisinfection(;kO3=0, kd, name=:disin_supply); initial_states = [0.5, 0], use_default_processes=false)
@named CSTR2 = CSTR(v; initial_states = [0.5, 0])

Finally, the influent is to be created: It needs the variables to provide an input for, and the corresponding values. We take the values directly from the default state, which are:

default_states()

\[ \begin{equation} \left[ \begin{array}{c} \mathtt{S\_O3}\left( t \right) \\ \mathtt{X\_B}\left( t \right) \\ \end{array} \right] \end{equation} \]

To this, the value vector has to be aligned and additionally the flowrate supplied. The ozone concentration is set to 0.5 to keep this level in the first reactor (Consider the ozone to be supplied before the first reactor) and the bacteria concentration is set to 1, such that one can interpret the bacteria in the effluent as the percentage of bacteria that are remaining.

@named influent = Influent([0.5, 1], flowrate=q)

Connection

Then, the connection equations for the systems are created. For this, simply the inflows and outflows can be connected as desired:

connection_eqns = [
	# Connect the flows
	connect(outflows(influent)[1], inflows(CSTR1)[1]), # Influent -> CSTR1
	connect(outflows(CSTR1)[1], inflows(CSTR2)[1]), # CSTR1 -> CSTR2
]

With the connections, we can then build the overall model:

@named system = ODESystem(connection_eqns, t, systems = [
	influent,
	CSTR1,
	CSTR2,
])

And simplify it to get the minimum set of equations:

system_simplified = structural_simplify(system)
full_equations(system_simplified) # Display the equations

\[ \begin{align} \frac{\mathrm{d} \mathtt{CSTR1.CSTR1\_reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{CSTR1.CSTR1\_reactor.states.S\_O3}\left( t \right)}{\mathtt{CSTR1.CSTR1\_reactor.V}} + \frac{\mathtt{influent.S\_O3\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{CSTR1.CSTR1\_reactor.V}} - \mathtt{CSTR1.disin\_supply.kO3} \mathtt{CSTR1.CSTR1\_reactor.states.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CSTR1.CSTR1\_reactor.states.X\_B}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.X\_B\_const.k} \mathtt{influent.q\_const.k}}{\mathtt{CSTR1.CSTR1\_reactor.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{CSTR1.CSTR1\_reactor.states.X\_B}\left( t \right)}{\mathtt{CSTR1.CSTR1\_reactor.V}} - \mathtt{CSTR1.disin\_supply.kd} \mathtt{CSTR1.CSTR1\_reactor.states.S\_O3}\left( t \right) \mathtt{CSTR1.CSTR1\_reactor.states.X\_B}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CSTR2.CSTR2\_reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{influent.q\_const.k} \mathtt{CSTR1.CSTR1\_reactor.states.S\_O3}\left( t \right)}{\mathtt{CSTR2.CSTR2\_reactor.V}} + \frac{ - \mathtt{influent.q\_const.k} \mathtt{CSTR2.CSTR2\_reactor.states.S\_O3}\left( t \right)}{\mathtt{CSTR2.CSTR2\_reactor.V}} - \mathtt{CSTR2.CSTR2\_disinfection.kO3} \mathtt{CSTR2.CSTR2\_reactor.states.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{CSTR2.CSTR2\_reactor.states.X\_B}\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{influent.q\_const.k} \mathtt{CSTR2.CSTR2\_reactor.states.X\_B}\left( t \right)}{\mathtt{CSTR2.CSTR2\_reactor.V}} + \frac{\mathtt{influent.q\_const.k} \mathtt{CSTR1.CSTR1\_reactor.states.X\_B}\left( t \right)}{\mathtt{CSTR2.CSTR2\_reactor.V}} - \mathtt{CSTR2.CSTR2\_disinfection.kd} \mathtt{CSTR2.CSTR2\_reactor.states.S\_O3}\left( t \right) \mathtt{CSTR2.CSTR2\_reactor.states.X\_B}\left( t \right) \end{align} \]

Simulation and Plotting

Finally, we simulate and plot the result.

## Simulate the system (for 1 day)
# First create a ODEProblem for the time span (0, 1) and then solve it
prob = ODEProblem(system_simplified, [], (0, 1), [])
sol = solve(prob)

## Plot the output
plot(sol,
	idxs = [CSTR1.states.X_B ,CSTR2.states.X_B], # Plot bacteria in tank
	legend = :right,
	title = "Remaining fraction of bacteria in each tank")