Chemical reactor with UV light project in python

Question

Introduction

I'm new to python and made this code which computes the decay of sulphur dioxide (SO2) pollutant in a chemical reactor by ozone (O3) when it is hit by UV light.

Code description

Short explanation of the code: using the Cantera python package I define a set of properties for the gas (air + SO2 + O3) and the reactor. The reactions are stored in an external file (photolysis.cti) which I included after the code. Then the reactions are allowed to play out but at each time instant the rate of the reaction controlled by UV light is updated. This is done based on the instantaneous concentration of the species present and the UV light distribution. This can be modelled using several models I implemented. Some models (MPSS, MSSS) are more accurate but also more complex so the simulation will take longer if they are selected (~45s vs ~1s for the rest). After the simulation finished, a plot with the decay of SO2 with time is shown:

Text

Aim of this review

My aim is to improve both the python aspect (syntax, making it easier to read, more idiomatic, object oriented...) as well as the performance of the code. Also having the whole code in one single file feels not ideal. Any tips on how to structure the project are very welcome. I know the code is long but I tried to comment all functions to improve readability and made a very simple log to know what the code is doing.

Python code

The code I wrote is the following:

"""
Code to compute a CSTR UV reactor with Cantera
"""
import re
import time
import pandas as pd
import numpy as np
import cantera as ct
import discretize.CylMesh
import matplotlib.pyplot as plt
from scipy.constants import N_A, h, c
def reactor_volume(L, r_outer, r_inner):
 """
 Computes the volume of an annular reactor (hollow cylinder).
 Parameters [units]
 ------------------
 L : int, float
 Length of the reactor/lamp [m]
 r_outer : int, float
 Outer radius [m]
 r_outer : int, float
 Inner radius [m]
 Returns [units]
 ---------------
 V_r : float
 The reactor volume [m**2]
 """
 return np.pi*L*(r_outer**2-r_inner**2)
def make_cyl_mesh(L, r, ncr, ncz):
 """
 Makes a cylindrical mesh object with one cell in the azimuthal direction.
 Parameters [units]
 ------------------
 L : int, float
 Length of the reactor/lamp [m]
 r : int, float
 Cylinder radius [m]
 ncr : int
 Number of cells in the radial direction
 ncz : int
 Number of cells in the axial direction
 Returns [units]
 ---------------
 mesh : obj
 The mesh cylinder object
 """
 dr = r/ncr # cell width r
 dz = L/ncz # cell width z
 hr = dr*np.ones(ncr)
 hz = dz*np.ones(ncz)
 offset = [0, 0, -L/2] # offset to center z on half the lamp length
 return discretize.CylMesh([hr, 2*np.pi, hz], offset)
def sigma_to_eps(sigma):
 """
 Converts the molar absorbtion cross section of a gas to the molar
 absorbtion coefficient in S.I. units.
 Parameters [units]
 ------------------
 sigma : int, float
 Absorbtion cross section (log base e) [cm**2]
 Returns [units]
 ---------------
 eps : float
 molar absorbtion coefficient (log base 10) [m**2 mol**-1]
 """
 return (sigma*N_A)/(1e4*np.log(10))
def ppm_to_C(ppm):
 """
 Converts ppm (parts per million) concentration of a gas to
 concentration in S.I. units.
 Parameters [units]
 ------------------
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 C : float
 concentration [mol m**-3]
 """
 rho_air = 1.2 # [kg m**-3]
 Mr_air = 28.96 # [g mol**-1] Note: this is for DRY air
 Vm_air = molar_volume(rho_air, Mr_air) # [m**3 mol**-1]
 C = ppm/(Vm_air*1e6)
 return C
def molar_volume(rho, Mr):
 """
 Computes the molar volume of a gas.
 Parameters [units]
 ------------------
 rho : int, float
 density [kg m**-3]
 Mr : int, float
 molar mass [g mol**-1]
 Returns [units]
 ---------------
 Vm : float
 molar volume [m**3 mol**-1]
 """
 return (Mr/rho)*1e-3
def exp_medium_absorption(r, ppm=300):
 """
 Computes the absorption of the medium considered in the fluence 
 rate models. Initially, I will treat only SO2 as absorbing species.
 Parameters [units]
 ------------------
 r : int, float, array
 Radial position (i.e. of a mesh element) [m]
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 ema : int, float, array
 Exponential medium absorption [-]
 """
 sigma = 1.551258e-19 # Value for SO2 [cm**2]
 sigma_air = sigma
 eps_air = sigma_to_eps(sigma_air)
 C_air = ppm_to_C(ppm)
 return np.exp(-np.log(10)*eps_air*C_air*r)
def coord_ls(L, nls):
 """
 Computes the coordinates of the number light source elements 
 which have been chosen to model the UV lamp.
 Parameters [units]
 ------------------
 L : int, float
 Length of the reactor/lamp [m]
 nls : int
 Number of light source elements
 Returns [units]
 ---------------
 coord_ls : array
 Coordinates of the number light source elements [m]
 """
 d_ls = L/nls
 return np.arange(-L/2, L/2+d_ls, d_ls)
def optimize_nls(case_variables, tol=1E-2):
 """
 Computes the optimum number of light sources (nls) for which to
 model a UV lamp (for the MPSS or MSSS models). This optimum is 
 the lowest nls that makes the fluence rate computed by the model
 (without medium absorbtion) match the fluence rate computed by 
 the LSI model for a given tolerance. Based on Powell and Lawryshyn 
 2015.
 Note: this routine takes some time, proceed with caution.
 Parameters [units]
 ------------------
 case_variables : array
 Array contining the following case variables:
 case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz] 
 P = lamp power [W]
 L = lamp length [m]
 eff = lamp efficiency at 254 nm wavelength
 r_outer = reactor radius [m]
 r_inner = quartz sleeve radius [m]
 ncr = number of mesh cells in r
 ncz = number of mesh cells in z
 tol : int, float
 Tolerance for the optimization routine
 Returns [units]
 ---------------
 nls : int
 (Optimum) number of light source elements
 """
 P, eff, L, r_outer, r_inner, ncr, ncz = case_variables
 V_r = reactor_volume(L, r_outer, r_inner)
 mesh_outer = make_cyl_mesh(L, r_outer, ncr, ncz)
 mesh_inner = make_cyl_mesh(L, r_inner, ncr, ncz)
 err, nls = 100, 10 # Initialize error and smallest possible number of light sources
 # Compute reference LSI vafr
 I_LSI = compute_vafr(case_variables, model='LSI')
 # Loop to find the optimum value for the number of light sources
 while err > tol*I_LSI:
 nls += 1
 I_i_outer = per_cell_fluence_rate_MPSS(
 mesh_outer, P, eff, L, nls=nls, no_abs=True)
 I_i_inner = per_cell_fluence_rate_MPSS(
 mesh_inner, P, eff, L, nls=nls, no_abs=True)
 I_MPSS = volume_avg_fluence_rate(
 I_i_outer, mesh_outer, I_i_inner, mesh_inner, V_r)
 err = np.abs(I_MPSS-I_LSI)
 return nls
def per_cell_fluence_rate_RAD(mesh, P, eff, L, no_abs=False, ppm=300):
 """
 Computes the fluence rate for each cell in a given mesh following
 the RAD model (Coenen 2013).
 Parameters [units]
 ------------------
 mesh : obj
 Mesh object
 P : int, float
 Lamp power [W]
 eff : float
 Lamp efficiency at 254 nm wavelength
 L : int, float
 Length of the reactor/lamp [m]
 no_abs : bool
 Flag to enable/disable medium absorption
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 I_i : array
 Fluence rate for each cell in the mesh [W m**-2]
 """
 R = mesh.gridCC[:, 0]
 if no_abs == True:
 I_i = ((P*eff)/(2*np.pi*R*L))
 else:
 I_i = ((P*eff)/(2*np.pi*R*L))*exp_medium_absorption(R, ppm=ppm)
 return I_i
def per_cell_fluence_rate_LSI(mesh, P, eff, L):
 """
 Computes the fluence rate for each cell in a given mesh following
 the LSI model (Blatchley 1998). It assumes no medium absorption
 Parameters [units]
 ------------------
 mesh : obj
 Mesh object
 P : int, float
 Lamp power [W]
 eff : float
 Lamp efficiency at 254 nm wavelength
 L : int, float
 Length of the reactor/lamp [m]
 Returns [units]
 ---------------
 I_i : array
 Fluence rate for each cell in the mesh [W m**-2]
 """
 R = mesh.gridCC[:, 0]
 H = mesh.gridCC[:, 2]
 I_i = ((P*eff)/(4*np.pi*L*R))*(np.arctan((L/2+H)/R)+np.arctan((L/2-H)/R))
 return I_i
def per_cell_fluence_rate_RAD_LSI(mesh, P, eff, L):
 """
 Computes the fluence rate for each cell in a given mesh following
 the RAD-LSI model (Liu 2004). It assumes no medium absorption
 Parameters [units]
 ------------------
 mesh : obj
 Mesh object
 P : int, float
 Lamp power [W]
 eff : float
 Lamp efficiency at 254 nm wavelength
 L : int, float
 Length of the reactor/lamp [m]
 Returns [units]
 ---------------
 I_i : array
 Fluence rate for each cell in the mesh [W m**-2]
 """
 R = mesh.gridCC[:, 0]
 H = mesh.gridCC[:, 2]
 I_i_LSI = ((P*eff)/(4*np.pi*L*R)) * \
 (np.arctan((L/2+H)/R)+np.arctan((L/2-H)/R))
 I_i_RAD = ((P*eff)/(2*np.pi*R*L))
 I_i = np.minimum(I_i_RAD, I_i_LSI)
 return I_i
def per_cell_fluence_rate_MPSS(mesh, P, eff, L, nls, no_abs=False, ppm=300):
 """
 Computes the fluence rate for each cell in a given mesh following
 the MPSS model (Bolton 2000).
 Parameters [units]
 ------------------
 mesh : obj
 Mesh object
 P : int, float
 Lamp power [W]
 eff : float
 Lamp efficiency at 254 nm wavelength
 L : int, float
 Length of the reactor/lamp [m]
 nls : int
 Number of light source elements
 no_abs : bool
 Flag to enable/disable medium absorption
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 I_i : array
 Fluence rate for each cell in the mesh [W m**-2]
 """
 x_ls = coord_ls(L, nls)
 R = mesh.gridCC[:, 0]
 H = mesh.gridCC[:, 2]
 I_i = np.zeros(len(R))
 for i in range(len(R)):
 r = R[i]*np.ones(len(x_ls))
 h = H[i]*np.ones(len(x_ls))
 rho = np.sqrt(r**2+(h-x_ls)**2)
 if no_abs == True:
 I_ls = (((P*eff)/nls)/(4*np.pi*rho**2))
 else:
 I_ls = (((P*eff)/nls)/(4*np.pi*rho**2)) * \
 exp_medium_absorption(r, ppm)
 I_i[i] = sum(I_ls)
 return I_i
def per_cell_fluence_rate_MSSS(mesh, P, eff, L, nls, no_abs=False, ppm=300):
 """
 Computes the fluence rate for each cell in a given mesh following
 the MSSS model (Liu 2004).
 Parameters [units]
 ------------------
 mesh : obj
 Mesh object
 P : int, float
 Lamp power [W]
 eff : float
 Lamp efficiency at 254 nm wavelength
 L : int, float
 Length of the reactor/lamp [m]
 nls : int
 Number of light source elements
 no_abs : bool
 Flag to enable/disable medium absorption
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 I_i : array
 Fluence rate for each cell in the mesh [W m**-2]
 """
 x_ls = coord_ls(L, nls)
 R = mesh.gridCC[:, 0]
 H = mesh.gridCC[:, 2]
 I_i = np.zeros(len(R))
 for i in range(len(R)):
 r = R[i]*np.ones(len(x_ls))
 h = H[i]*np.ones(len(x_ls))
 rho = np.sqrt(r**2+(h-x_ls)**2)
 if no_abs == True:
 I_ls = (((P*eff)/nls)/(4*np.pi*rho**2))
 else:
 I_ls = (((P*eff)/nls)/(4*np.pi*rho**2)) * \
 exp_medium_absorption(r, ppm)
 I_ls = I_ls*np.cos(np.arctan(r/h))
 I_i[i] = sum(I_ls)
 return I_i
def per_cell_fluence_rate(mesh, P, eff, L, model='RAD', nls=40, ppm=300):
 """
 Generic function that computes the fluence rate for each cell 
 in a given mesh for various model.
 Parameters [units]
 ------------------
 mesh : obj
 Mesh object
 P : int, float
 Lamp power [W]
 eff : float
 Lamp efficiency at 254 nm wavelength
 L : int, float
 Length of the reactor/lamp [m]
 model : str
 name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
 nls : int
 Number of light source elements
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 I_i : array
 Fluence rate for each cell in the mesh [W m**-2]
 """
 if isinstance(model, str) == 0:
 raise TypeError("The variable 'model' has to be a string")
 if model == 'RAD':
 I_i = per_cell_fluence_rate_RAD(mesh, P, eff, L, ppm=ppm)
 elif model == 'LSI':
 I_i = per_cell_fluence_rate_LSI(mesh, P, eff, L)
 elif model == 'MPSS':
 I_i = per_cell_fluence_rate_MPSS(mesh, P, eff, L, nls=nls, ppm=ppm)
 elif model == 'MSSS':
 I_i = per_cell_fluence_rate_MSSS(mesh, P, eff, L, nls=nls, ppm=ppm)
 elif model == 'RAD_LSI':
 I_i = per_cell_fluence_rate_RAD_LSI(mesh, P, eff, L)
 else:
 raise ValueError(
 "Only 'RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI' models are implemented")
 return I_i
def volume_avg_fluence_rate(I_i_outer, mesh_outer, I_i_inner, mesh_inner, V_r):
 """
 Computes the volume averaged fluence rate (vafr) for the
 annular reactor (hollow cylinder) by subtracting the
 average of the small cylinder (_inner) from the large 
 cylinder (_outer).
 Parameters [units]
 ------------------
 I_i_outer : array
 Fluence rate for each cell in the outer cylinder mesh [W m**-2]
 mesh_outer : obj
 outer cylinder mesh object
 I_i_inner : array
 Fluence rate for each cell in the inner cylinder mesh [W m**-2]
 mesh_inner : obj
 Inner cylinder mesh object
 V_r : float
 The reactor volume [m**2]
 Returns [units]
 ---------------
 I_vol_avg : array
 Volume average fluence rate for the annular reactor [W m**-2]
 """
 # volume averaged fluence rate (vafr)
 I_outer = np.sum(I_i_outer*mesh_outer.vol)
 I_inner = np.sum(I_i_inner*mesh_inner.vol)
 return (I_outer - I_inner)/V_r
def compute_vafr(case_variables, model='RAD', nls=40, optimize=False, ppm=300):
 """
 Generic function to compute the volume average fluence rate.
 Parameters [units]
 ------------------
 case_variables : array
 Array contining the following case variables:
 case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz] 
 P = lamp power [W]
 L = lamp length [m]
 eff = lamp efficiency at 254 nm wavelength
 r_outer = reactor radius [m]
 r_inner = quartz sleeve radius [m]
 ncr = number of mesh cells in r
 ncz = number of mesh cells in z
 model : str
 name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
 nls : int
 Number of light source elements
 no_abs : bool
 Flag to enable/disable nls optimization for the MPSS or MSSS models
 ppm : int, float
 ppm concentration of a gas [-]
 Returns [units]
 ---------------
 I_vol_avg : array
 Volume average fluence rate for the annular reactor [W m**-2]
 """
 # Unpack the case variables
 P, eff, L, r_outer, r_inner, ncr, ncz = case_variables
 # Compute reactor volume
 V_r = reactor_volume(L, r_outer, r_inner)
 # Make the outer and inner mesh objects
 mesh_outer = make_cyl_mesh(L, r_outer, ncr, ncz)
 mesh_inner = make_cyl_mesh(L, r_inner, ncr, ncz)
 # Find optimum nls if required
 if (model == 'MPSS' or model == 'MSSS') and optimize == True:
 nls = optimize_nls(case_variables, tol=1E-2)
 # Calculate the per-cell fluence rate for the inner and outer meshes
 I_i_outer = per_cell_fluence_rate(
 mesh_outer, P, eff, L, model=model, nls=nls, ppm=ppm)
 I_i_inner = per_cell_fluence_rate(
 mesh_inner, P, eff, L, model=model, nls=nls, ppm=ppm)
 # Compute volume average LSI fluence rate
 I_vol_avg = volume_avg_fluence_rate(
 I_i_outer, mesh_outer, I_i_inner, mesh_inner, V_r)
 return I_vol_avg
def set_gas(mechanism, T, P, C):
 """
 Define the properties of the cantera gas object.
 Parameters [units]
 ------------------
 mechanism : string
 Path to the CTI mechanism to be used
 T : int, float
 Temperature [K]
 P : int, float
 Pressure [in atm]
 C : dict
 Dictionary containing the concentrations of the different gas species [molar fraction]
 Returns [units]
 ---------------
 gas : obj
 cantera gas object
 """
 gas = ct.Solution('data/photolysis.cti')
 gas.TPX = T, P, C
 return gas
def set_reactor(gas, V_r, residence_t, p_valve_coeff):
 """
 Set-up the reactor network using Cantera functions and syntax.
 Parameters [units]
 ------------------
 gas : obj
 Cantera gas object
 V_r : float
 The reactor volume [m**2]
 residence_t : int, float
 Time that the gas will stay in the reactor [s]
 p_valve_coeff : int, float
 This is the "conductance" of the pressure valve and will 
 determine its efficiency in holding the reactor pressure 
 to the desired conditions. Set to 0.01 for default
 Returns [units]
 ---------------
 stirred_reactor : obj
 Cantera reactor object
 reactor_network : obj
 Cantera network object
 """
 fuel_air_mixture_tank = ct.Reservoir(gas)
 exhaust = ct.Reservoir(gas)
 stirred_reactor = ct.IdealGasReactor(gas, energy='off', volume=V_r)
 mass_flow_controller = ct.MassFlowController(
 upstream=fuel_air_mixture_tank, downstream=stirred_reactor, mdot=stirred_reactor.mass/residence_t)
 pressure_regulator = ct.Valve(
 upstream=stirred_reactor, downstream=exhaust, K=p_valve_coeff)
 reactor_network = ct.ReactorNet([stirred_reactor])
 return [stirred_reactor, reactor_network]
def intialize_results(stirred_reactor):
 """
 Initializes the DataFrame where results will be stored.
 Parameters [units]
 ------------------
 stirred_reactor : obj
 Cantera reactor object
 Returns [units]
 ---------------
 time_history : DataFrame
 Pandas DataFrame to store time evolution of reasults. It
 includes the gas's mass, volume, temperature and species
 concentration
 """
 column_names = [stirred_reactor.component_name(
 item) for item in range(stirred_reactor.n_vars)]
 time_history = pd.DataFrame(columns=column_names)
 return time_history
def update_results(time_history, stirred_reactor, t):
 """
 Updates the results at a given time instant.
 Parameters [units]
 ------------------
 time_history : DataFrame
 Pandas DataFrame to store time evolution of reasults. It
 includes the gas's mass, volume, temperature and species
 concentration
 stirred_reactor : obj
 Cantera reactor object
 t : int, float
 Current simulation time [s]
 Returns [units]
 ---------------
 time_history : DataFrame
 Pandas DataFrame to store time evolution of reasults. It
 includes the gas's mass, volume, temperature and species
 concentration
 """
 state = np.hstack([stirred_reactor.mass, stirred_reactor.volume,
 stirred_reactor.T, stirred_reactor.thermo.X])
 time_history.loc[t] = state
 return time_history
def uv_update_reactions(stirred_reactor, case_variables, model='RAD', nls=40):
 """
 Modifies the rate constant of the photolysis reaction of ozone. This
 has to be done for each loop iteration as the concentrations of
 the different gas species are changing, which will affect the absorption
 of the medium.
 Parameters [units]
 ------------------
 stirred_reactor : obj
 Cantera reactor object
 case_variables : array
 Array contining the following case variables:
 case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz] 
 P = lamp power [W]
 L = lamp length [m]
 eff = lamp efficiency at 254 nm wavelength
 r_outer = reactor radius [m]
 r_inner = quartz sleeve radius [m]
 ncr = number of mesh cells in r
 ncz = number of mesh cells in z
 model : str
 name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
 nls : int
 Number of light source elements
 Returns [units]
 ---------------
 None
 """
 ppm = stirred_reactor.thermo.X[stirred_reactor.component_index(
 'SO2')-3]*1e6
 k_uv = (0.9) * ((254E-9)/(N_A*h*c)) * (np.log(10)) * sigma_to_eps(1.132935E-17) * \
 compute_vafr(case_variables, model=model, nls=nls, ppm=ppm)
 ID = np.size(gas.reactions())-1
 reaction = gas.reactions()[ID]
 reaction.rate = ct.Arrhenius(A=k_uv, b=0, E=0)
 gas.modify_reaction(ID, reaction)
 return None
def time_loop(reactor_network, stirred_reactor, case_variables, time_history, max_simulation_t, model='RAD'):
 """
 Computes the time loop iterations for the chemical reactor
 until the maximum simulation time is reached.
 Parameters [units]
 ------------------
 stirred_reactor : obj
 Cantera reactor object
 reactor_network : obj
 Cantera network object
 case_variables : array
 Array contining the following case variables:
 case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz] 
 P = lamp power [W]
 L = lamp length [m]
 eff = lamp efficiency at 254 nm wavelength
 r_outer = reactor radius [m]
 r_inner = quartz sleeve radius [m]
 ncr = number of mesh cells in r
 ncz = number of mesh cells in z
 time_history : DataFrame
 Pandas DataFrame to store time evolution of reasults. It
 includes the gas's mass, volume, temperature and species
 concentration
 max_simulation_t : int
 Maximum simulation time [s]
 model : str
 name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
 Returns [units]
 ---------------
 time_history : DataFrame
 Pandas DataFrame to store time evolution of reasults. It
 includes the gas's mass, volume, temperature and species
 concentration
 """
 counter, t = 1, 0
 while t < max_simulation_t:
 if (t == 0) or (counter % 10 == 0): # Update every 10 iterations
 if (model == 'MPSS') or (model == 'MSSS'):
 if t == 0:
 nls = optimize_nls(case_variables, tol=1E-2)
 uv_update_reactions(
 stirred_reactor, case_variables, model=model, nls=nls)
 else:
 uv_update_reactions(stirred_reactor, case_variables, model=model)
 time_history = update_results(time_history, stirred_reactor, t)
 t = reactor_network.step()
 counter += 1
 return time_history
def run_simulation(gas, V_r, residence_t, p_valve_coeff, case_variables, max_simulation_t, model='RAD'):
 """
 Generic function to run the chemical reactor with UV light
 modeling.
 Parameters [units]
 ------------------
 gas : obj
 Cantera gas object
 V_r : float
 The reactor volume [m**2]
 residence_t : int, float
 Time that the gas will stay in the reactor [s]
 p_valve_coeff : int, float
 This is the "conductance" of the pressure valve and will 
 determine its efficiency in holding the reactor pressure 
 to the desired conditions. Set to 0.01 for default
 case_variables : array
 Array contining the following case variables:
 case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz] 
 P = lamp power [W]
 L = lamp length [m]
 eff = lamp efficiency at 254 nm wavelength
 r_outer = reactor radius [m]
 r_inner = quartz sleeve radius [m]
 ncr = number of mesh cells in r
 ncz = number of mesh cells in z
 max_simulation_t : int
 Maximum simulation time [s]
 model : str
 name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
 Returns [units]
 ---------------
 time_history : DataFrame
 Pandas DataFrame to store time evolution of reasults. It
 includes the gas's mass, volume, temperature and species
 concentration
 """
 # Define the reactor object
 [stirred_reactor, reactor_network] = set_reactor(
 gas, V_r, residence_t, p_valve_coeff)
 # Initialize the results dataFrame
 time_history = intialize_results(stirred_reactor)
 # Run the time loop
 time_history = time_loop(reactor_network, stirred_reactor, case_variables, time_history,
 max_simulation_t, model=model)
 return time_history
def plot_results(specie, log=False):
 """
 Plots the results from the reactor simulation.
 Parameters [units]
 ------------------
 specie : str
 Specie for which to plot the time evolution
 log : bool
 Flag to enable/disable the x axis as logarithmic
 Returns [units]
 ---------------
 None
 """
 plt.style.use('ggplot')
 plt.style.use('seaborn-pastel')
 plt.rcParams['axes.labelsize'] = 18
 plt.rcParams['xtick.labelsize'] = 14
 plt.rcParams['ytick.labelsize'] = 14
 plt.rcParams['figure.autolayout'] = True
 plt.figure()
 if log == True:
 plt.semilogx(time_history.index, time_history[specie]*1e6, '-')
 else:
 plt.plot(time_history.index, time_history[specie]*1e6, '-')
 plt.xlabel('Time (s)')
 plt.ylabel(re.sub(r"([0-9]+(\.[0-9]+)?)", r"_1円",
 '$'+specie+'$' + r' Mole Fraction : $ppmv$'))
 plt.show()
 return None
if __name__ == "__main__":
 start = time.time()
 print('Defining the gas properties')
 gas = set_gas('photolysis.cti', 295, ct.one_atm,
 {'N2': 0.78014, 'O2': 0.20946, 'O3': 0.0001, 'SO2': 0.0003, 'H2O': 0.01})
 print('Defining the geometry variables')
 case_variables = P, eff, L, r_outer, r_inner, ncr, ncz = [
 17, 0.33, 0.28, 0.0425, 0.0115, 20, 20]
 V_r = reactor_volume(L, r_outer, r_inner)
 print('Defining the simulation variables')
 p_valve_coeff, max_p_rise, residence_t, max_simulation_t = 0.01, 0.01, 10, 100
 print('Running the simulation')
 time_history = run_simulation(gas, V_r, residence_t, p_valve_coeff, case_variables,
 max_simulation_t, model='RAD')
 end = time.time()
 print('Plotting the simulation results')
 plot_results('SO2', log=True)
 print(f'Elapsed time: {end - start:.3f}s')

Kinetic mechanism

If you wish to run the code you will require the kinetic mechanism. Save it as photolysis.cti:

"""
-----------------------------------------------------------------------|
 |
 Photolysis mechanism |
 |
-----------------------------------------------------------------------|
"""
units(length='cm', time='s', quantity='molec', act_energy='K')
ideal_gas(name='gas',
 elements="O H N S",
 species="""O O(1D) O2 O3 H2O HO2 OH N2
 SO2 SO3 HSO3 H2SO4 H""",
 reactions='all',
 initial_state=state(temperature=300.0, pressure=OneAtm))
#-------------------------------------------------------------------------------
# Species data
#-------------------------------------------------------------------------------
species(name='O',
 atoms='O:1',
 thermo=(NASA([200.00, 1000.00],
 [ 2.94642900E+00, -1.63816600E-03, 2.42103200E-06,
 -1.60284300E-09, 3.89069600E-13, 2.91476400E+04,
 2.96399500E+00]),
 NASA([1000.00, 5000.00],
 [ 2.54206000E+00, -2.75506200E-05, -3.10280300E-09,
 4.55106700E-12, -4.36805200E-16, 2.92308000E+04,
 4.92030800E+00])),
 note='120186')
species(name='O(1D)',
 atoms='O:1',
 thermo=(NASA([200.00, 1000.00],
 [ 2.49993786E+00, 1.71935346E-07, -3.45215267E-10,
 3.71342028E-13, -1.70964494E-16, 5.19965317E+04,
 4.61684555E+00]),
 NASA([1000.00, 6000.00],
 [ 2.49368475E+00, 1.37617903E-05, -1.00401058E-08,
 2.76012182E-12, -2.01597513E-16, 5.19986304E+04,
 4.65050950E+00])),
 note='121286')
species(name='O2',
 atoms='O:2',
 thermo=(NASA([200.00, 1000.00],
 [ 3.21293600E+00, 1.12748600E-03, -5.75615000E-07,
 1.31387700E-09, -8.76855400E-13, -1.00524900E+03,
 6.03473800E+00]),
 NASA([1000.00, 5000.00],
 [ 3.69757800E+00, 6.13519700E-04, -1.25884200E-07,
 1.77528100E-11, -1.13643500E-15, -1.23393000E+03,
 3.18916600E+00])),
 note='121386')
species(name='O3',
 atoms='O:3',
 thermo=(NASA([200.00, 1000.00],
 [ 2.46260900E+00, 9.58278100E-03, -7.08735900E-06,
 1.36336800E-09, 2.96964700E-13, 1.60615200E+04,
 1.21418700E+01]),
 NASA([1000.00, 5000.00],
 [ 5.42937100E+00, 1.82038000E-03, -7.70560700E-07,
 1.49929300E-10, -1.07556300E-14, 1.52352700E+04,
 -3.26638700E+00])),
 note='121286')
species(name='H2O',
 atoms='H:2 O:1',
 thermo=(NASA([200.00, 1000.00],
 [ 3.38684200E+00, 3.47498200E-03, -6.35469600E-06,
 6.96858100E-09, -2.50658800E-12, -3.02081100E+04,
 2.59023300E+00]),
 NASA([1000.00, 5000.00],
 [ 2.67214600E+00, 3.05629300E-03, -8.73026000E-07,
 1.20099600E-10, -6.39161800E-15, -2.98992100E+04,
 6.86281700E+00])),
 note='20387')
species(name='HO2',
 atoms='H:1 O:2',
 thermo=(NASA([200.00, 1000.00],
 [ 2.97996300E+00, 4.99669700E-03, -3.79099700E-06,
 2.35419200E-09, -8.08902400E-13, 1.76227400E+02,
 9.22272400E+00]),
 NASA([1000.00, 5000.00],
 [ 4.07219100E+00, 2.13129600E-03, -5.30814500E-07,
 6.11226900E-11, -2.84116500E-15, -1.57972700E+02,
 3.47602900E+00])),
 note='20387')
species(name='OH',
 atoms='H:1 O:1',
 thermo=(NASA([200.00, 1000.00],
 [ 3.63726600E+00, 1.85091000E-04, -1.67616500E-06,
 2.38720300E-09, -8.43144200E-13, 3.60678200E+03,
 1.35886000E+00]),
 NASA([1000.00, 5000.00],
 [ 2.88273000E+00, 1.01397400E-03, -2.27687700E-07,
 2.17468400E-11, -5.12630500E-16, 3.88688800E+03,
 5.59571200E+00])),
 note='121286')
species(name='N2',
 atoms='N:2',
 thermo=(NASA([100.00, 1000.00],
 [ 3.29867700E+00, 1.40824000E-03, -3.96322200E-06,
 5.64151500E-09, -2.44485500E-12, -1.02090000E+03,
 3.95037200E+00]),
 NASA([1000.00, 5000.00],
 [ 2.92664000E+00, 1.48797700E-03, -5.68476100E-07,
 1.00970400E-10, -6.75335100E-15, -9.22797700E+02,
 5.98052800E+00])),
 note='121286')
species(name='SO2',
 atoms='O:2 S:1',
 thermo=(NASA([200.00, 1000.00],
 [ 2.91143900E+00, 8.10302200E-03, -6.90671000E-06,
 3.32901600E-09, -8.77712100E-13, -3.68788200E+04,
 1.11174000E+01]),
 NASA([1000.00, 5000.00],
 [ 5.25449800E+00, 1.97854500E-03, -8.20422600E-07,
 1.57638300E-10, -1.12045100E-14, -3.75688600E+04,
 -1.14605600E+00])),
 note='121286')
species(name='SO3',
 atoms='O:3 S:1',
 thermo=(NASA([200.00, 1000.00],
 [ 2.57528300E+00, 1.51509200E-02, -1.22987200E-05,
 4.24025700E-09, -5.26681200E-13, -4.89441100E+04,
 1.21951200E+01]),
 NASA([1000.00, 5000.00],
 [ 7.05066800E+00, 3.24656000E-03, -1.40889700E-06,
 2.72153500E-10, -1.94236500E-14, -5.02066800E+04,
 -1.10644300E+01])),
 note='121286')
species(name='HSO3',
 atoms='H:1 O:3 S:1',
 thermo=(NASA([200.00, 1000.00],
 [ 3.18471474E+00, 1.92277645E-02, -1.85975873E-05,
 7.95502482E-09, -9.42254910E-13, -4.42156599E+04,
 1.30496771E+01]),
 NASA([1000.00, 6000.00],
 [ 8.19420431E+00, 3.77828016E-03, -1.34903219E-06,
 2.17197023E-10, -1.29874848E-14, -4.55013223E+04,
 -1.23824851E+01])),
 note='HO-SO2T10/10')
species(name='H2SO4',
 atoms='H:2 O:4 S:1',
 thermo=(NASA([200.00, 1000.00],
 [ 4.53388173E+00, 3.10347679E-02, -4.10421795E-05,
 2.95752341E-08, -8.81459071E-12, -9.05459072E+04,
 3.93961412E+00]),
 NASA([1000.00, 6000.00],
 [ 1.13355392E+01, 5.60829109E-03, -1.94574192E-06,
 3.07136054E-10, -1.81109544E-14, -9.21087435E+04,
 -2.96094003E+01])),
 note='T8/03')
species(name='H',
 atoms='H:1',
 thermo=(NASA([200.00, 1000.00],
 [ 2.50000000E+00, 0.00000000E+00, 0.00000000E+00,
 0.00000000E+00, 0.00000000E+00, 2.54716300E+04,
 -4.60117600E-01]),
 NASA([1000.00, 5000.00],
 [ 2.50000000E+00, 0.00000000E+00, 0.00000000E+00,
 0.00000000E+00, 0.00000000E+00, 2.54716300E+04,
 -4.60117600E-01])),
 note='120186')
#-------------------------------------------------------------------------------
# Reaction data
#-------------------------------------------------------------------------------
# Reaction 1
reaction('O + O3 => 2 O2', [8.000000e+12, 0.0, 2060.0])
# Reaction 2
reaction('O(1D) + H2O => 2 OH', [1.630000e+10, 0.0, 65.0])
# Reaction 3
reaction('O + OH => O2 + H', [2.400000e+11, 0.0, 110.0])
# Reaction 4
reaction('O + HO2 => OH + O2', [2.700000e+11, 0.0, 224.0])
# Reaction 5
three_body_reaction('H + O2 + M => HO2 + M', [2.170000e+29, 1.1, 0.0])
# Reaction 6
reaction('OH + O3 => HO2 + O2', [1.700000e+12, 0.0, 940.0])
# Reaction 7
reaction('OH + HO2 => H2O + O2', [4.800000e+11, 0.0, 250.0])
# Reaction 8
reaction('OH + SO2 => HSO3', [7.070000e+11, 0.7, 0.0])
# Reaction 9
three_body_reaction('O + SO2 + M => SO3 + M', [4.000000e+32, 0.0, 1000.0])
# Reaction 10
reaction('SO3 + H2O => H2SO4', [1.200000e+15, 0.0, 0.0])
# Reaction 11
reaction('HSO3 + O2 => HO2 + SO3', [1.300000e+12, 0.0, 330.0])
#-------------------------------------------------------------------------------
# Photolysis reaction
#-------------------------------------------------------------------------------
# Reaction 12
reaction('O3 => O(1D) + O2', [1.000000e+00, 0.0, 0.0])

Answer 1

Overall

This code is fairly well-structured and very well-documented. Excellent job.

Type hints

It's good that you specified types here:

 L : int, float
 Length of the reactor/lamp [m]
 r_outer : int, float
 Outer radius [m]
 r_outer : int, float

Better would be to codify these in type hints:

def reactor_volume(
 L: Union[int, float], 
 r_outer: Union[int, float], 
 r_inner: Union[int, float],
) -> float:

Some alternatives to the union are explained here:

https://stackoverflow.com/questions/50928592/mypy-type-hint-unionfloat-int-is-there-a-number-type

Order of operations

(sigma*N_A)/(1e4*np.log(10))

can be

sigma * N_A * 1e-4 / np.log(10)

Simplification

The exp will cancel out the log here:

np.exp(-np.log(10)*eps_air*C_air*r)

to become

np.power(10, -eps_air*C_air*r)

Boolean comparison

if no_abs == True should simply be if no_abs.

Enums

model should capture RAD, LSI etc. in an Enum.

set_reactor return value

You're forming a list, which is a little odd. Typically, multi-return values are implicit tuples, i.e.

return stirred_reactor, reactor_network

Spelling

reasults -> results

Redundant return

return None

can be deleted; that's what happens by default.

Dictionary update

plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 14
plt.rcParams['ytick.labelsize'] = 14
plt.rcParams['figure.autolayout'] = True

should use

plt.rcParams.update({ ...

Graph bounds

The lower bound for your time axis should be 1e-3 or maybe 1e-2. Remove the long, flat part to focus on the interesting bits.

• Marc: Thanks got the feedback! I'll apply all the changes but I just wanted some clarification on the Enums. You are suggesting to make a model a class and then use the Enum package to define the different models within the class, right? I don't have much experience with OOP programming, could you give me a pointer to a good resource that I can use as an example?
• Reinderien: docs.python.org/3/library/enum.html has good examples and is the only "authoritative source"