Batch Reactor and Initial Condition

As a minimal example for batch reactors and initial conditions, we consider a batch reactor for ozonation. In the reactor, in the beginning a certain concentration of ozone is assumed, and then bacteria are spiked. We will here investigate the behavior of this simple system for different initial conditions.

But first, the system needs be built. For this, 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
kO3 = 100.0 / 24  # Reaction rate for the decay of ozone
kd = 200.0 / 24  # Reaction rate for the decay of the bacteria

# Initial concentrations of bacteria and ozone
IC_O3 = .5
IC_B = 1

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

# Generate the reactor
@named batch_reactor = BatchReactor(;
	initial_states = [states(disinfection).S_O3 => IC_O3,
	states(disinfection).X_B => IC_B]
	)

# Build the ODESystem
@named system = ODESystem([], t, systems = [
	batch_reactor
])
system_simplified = structural_simplify(system) # simplify

## Simulate the system (5 days)
# First create a ODEProblem and then solve it
prob = ODEProblem(system_simplified, [], (0, 5), [])
sol = solve(prob)

@assert sol[states(batch_reactor).X_B, end] > 0.1 # There not all bacteria are removed!

# Try with higher initial ozone concentration
## Simulate the system with a higher ozone concentration
prob = ODEProblem(system_simplified, [states(batch_reactor).S_O3 => 5], (0, 5), [])
sol = solve(prob)

## Plot the output
plot(sol, title = "Disinfection Batch Test: Better")

# Or adapting the parameter for decomposition rate of the bacteria with ozone
prob = ODEProblem(system_simplified, [], (0, 5), [parameters(system_simplified)[1] => 200])
sol = solve(prob)

## Plot the output
plot(sol, title = "Disinfection Batch Test: Faster rate")

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
kO3 = 100.0 / 24  # Reaction rate for the decay of ozone
kd = 200.0 / 24  # Reaction rate for the decay of the bacteria

# Initial concentrations of bacteria and ozone
IC_O3 = .5
IC_B = 1

System Creation

Then, all needed systems (The reactor) is created.

We start with setting the process for the disinfection dynamics. This process will then be automatically used for all subsequent reactors or flow elements, thus it is important to start with this step.

@set_process disinfection = OzoneDisinfection(;kO3, kd)

We continue with the batch reactor. It has an input for the provided rates (from the reaction rate system just created) and an output for the states, which are needed by the disinfection reaction rates system for computation. To create the model we provide it with the states to have in the flow, the volume is not needed as the batch process has no inflow. The compounds can be extracted as states from the OzoneDisinfection system above. Here we use the states of the corresponding rates, which is convenient. Note, however, that this is not necessary and the compounds can as well be specified manually if needed.

@named batch_reactor = BatchReactor(; 
	initial_states = [states(disinfection).S_O3 => IC_O3, 
	states(disinfection).X_B => IC_B]
	)

Batch Reactor with OzoneDisinfection 'batch_reactor':
States (2): see states(batch_reactor)
  S_O3(t) [defaults to 0.5]: S_O3
  X_B(t) [defaults to 1.0]: X_B
Parameters (2): see parameters(batch_reactor)
  batch_reactor_disinfection₊kO3 [defaults to 4.16667]
  batch_reactor_disinfection₊kd [defaults to 8.33333]

Connection

This example is a simple batch reactor. Thus it has no in- or outflows and the overall model can be directly built:

@named system = ODESystem([], t, systems = [
	batch_reactor
])

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{batch\_reactor.batch\_reactor\_reactor.states.S\_O3}\left( t \right)}{\mathrm{d}t} &= - \mathtt{batch\_reactor.batch\_reactor\_disinfection.kO3} \mathtt{batch\_reactor.batch\_reactor\_reactor.states.S\_O3}\left( t \right) \\ \frac{\mathrm{d} \mathtt{batch\_reactor.batch\_reactor\_reactor.states.X\_B}\left( t \right)}{\mathrm{d}t} &= - \mathtt{batch\_reactor.batch\_reactor\_disinfection.kd} \mathtt{batch\_reactor.batch\_reactor\_reactor.states.X\_B}\left( t \right) \mathtt{batch\_reactor.batch\_reactor\_reactor.states.S\_O3}\left( t \right) \end{align} \]

Simulation and Plotting

Finally, we simulate and plot the result.

## Simulate the system (5 days)
# First create a ODEProblem and then solve it
prob = ODEProblem(system_simplified, [], (0, 5), [])
sol = solve(prob)

## Plot the output
plot(sol, title = "Disinfection Batch Test")

Adaptation of initial conditions

From the plot above, we see that there is not enough ozone in the tank and bacteria remain in it. Thus we would like to increase its initial condition and try with 10 times the value. But how do we do this? It is necessary to rebuild the whole system now?

No, luckily there is another option, and we can as well provide the initial conditions directly when generating the problem:

## Simulate the system with a higher ozone concentration
prob = ODEProblem(system_simplified, [states(batch_reactor).S_O3 => 5], (0, 5), [])
sol = solve(prob)

## Plot the output
plot(sol, title = "Disinfection Batch Test")

And as we see here, the initial condition of the ozone equation is indeed changed, and it is high enough for eliminating all bacteria.

Adaptation of parameters

Alternatively to changing the initial condiditions, we could (of course only mathematically and not in reality) as well just increase the rate at which the bacteria decay with the ozone.

Adapting this parameter is as well possible when generating the problem, similar to adapting the initial conditions:

## Simulate the system with a higher ozone concentration
prob = ODEProblem(system_simplified, [], (0, 5), [parameters(system_simplified)[1] => 200])
sol = solve(prob)

## Plot the output
plot(sol, title = "Disinfection Batch Test")