Expected behavior.
Simulation runs with CSTR init_liquid_volume < 1e-4
Actual behavior
Simulation hangs with CSTR init_liquid_volume < 1e-4
How to produce bug (including a minimal reproducible example)
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# import cadet packages
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import LangmuirLDF
from CADETProcess.processModel import Inlet, LumpedRateModelWithoutPores, Outlet
from CADETProcess.processModel import TubularReactor, Cstr
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem
from CADETProcess.optimization import U_NSGA3
# setup experiment
# conversion
mLmin_to_m3s = 60*1e6
# input parameters
exp_flow = 8. # MV/min (mL/min)
Q = exp_flow/(mLmin_to_m3s) # convert to m^3/s
# load and wash
loadVol = 60.1 #42.3 #43.0 #60.1 # mL
washVol = 34.8 #33.9 #25.0 #34.8 # mL
protein_c = 3.5 #3.7 #3.6 #3.5 # g/L
# --------------------------------------
total_vol = loadVol + washVol # mL
exp_time = total_vol/exp_flow*60 # seconds
# timepoints
load_time = loadVol/(exp_flow/60.)
wash_time = washVol/(exp_flow/60.)
# device parameters
device_Dax = 5.413e-9 # m^2/s
# protein size
size_kDa = 150.
Qmax = 80 # mg/mL = g/L
Qmax_mMol = Qmax/size_kDa# -> mMolar
# --------------------
# create component system
# --------------------
component_system = ComponentSystem()
component_system.add_component('pH') # pH currently dummy variable
component_system.add_component('protein')
# --------------------
# binding model
# ---------------------
binding_model = LangmuirLDF(component_system, name='MultiComponentLangmuir_LDF')
binding_model.is_kinetic = True
binding_model.equilibrium_constant =[0,20000]
binding_model.driving_force_coefficient = [0,1.0]
binding_model.capacity = [0,Qmax_mMol]
binding_model.bound_states=[0,1]
binding_model.check_required_parameters()
# ---------------------
# feed
# ---------------------
feed = Inlet(component_system, name='feed')
feed.c = [1,protein_c/size_kDa]
# ---------------------
# washout
# ---------------------
wash = Inlet(component_system, name='wash')
wash.c = [1,0]
# ----------------------
# device
# ----------------------
device = LumpedRateModelWithoutPores(component_system, name='device')
device.length = 3.305/1000 #
device.diameter = 2.85/100 #
device.total_porosity = 0.804
device.axial_dispersion = 4.25e-07
device.binding_model = binding_model
device.c = [1,0]
#device.q = [1,0]
# ----------------------
# feed tubing
# ----------------------
load_tubing = TubularReactor(component_system, name='load_tubing')
load_tubing.c = [1,0]
load_tubing.length = 270/100
load_tubing.diameter = 0.75/1000
load_tubing.axial_dispersion = 0.005
# --------------------
# dead_vol
# --------------------
dead_vol = TubularReactor(component_system, name='dead_vol')
dead_vol.c = [1,0]
dead_vol.length = 1.34/100
dead_vol.diameter = 4/1000
dead_vol.axial_dispersion = 0.005531511
# ----------------------
# outlet
# ----------------------
outlet = Outlet(component_system, name='outlet')
#outlet.solution_recorder.write_solution_bulk = True
# -----------------------
# CSTR
# -----------------------
inlet_mixer = Cstr(component_system, name='inlet_mixer')
inlet_mixer.c = [1,0]
#inlet_mixer.V = 0.001 # units?
#inlet_mixer.V = 0.0001
#!!!! ISSUE IS HERE !!!!
inlet_mixer.init_liquid_volume = 1e-5 #m^3
inlet_mixer.const_solid_volume = 0.0
inlet_mixer.flow_rate_filter = Q # 8 mL/min
inlet_mixer.check_required_parameters()
inlet_mixer.const_solid_volume = 0
#inlet_mixer.solution_recorder.write_solution_bulk=True
# ------------------------------------------------
# now we create the Flow Sheet
# ------------------------------------------------
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(feed)
flow_sheet.add_unit(wash)
flow_sheet.add_unit(load_tubing)
flow_sheet.add_unit(device)
flow_sheet.add_unit(dead_vol)
flow_sheet.add_unit(outlet)
flow_sheet.add_unit(inlet_mixer)
flow_sheet.add_connection(feed, load_tubing)
flow_sheet.add_connection(wash, load_tubing)
flow_sheet.add_connection(load_tubing, inlet_mixer)
flow_sheet.add_connection(inlet_mixer,device)
flow_sheet.add_connection(device, dead_vol)
flow_sheet.add_connection(dead_vol, outlet)
# --------------------------------------------------------------
# now we create the process, i.e., turning feed on, swap to wash buffer, same flow rate
# --------------------------------------------------------------
process = Process(flow_sheet, 'breakthrough')
# breakthrough phase
process.add_event('feed_on','flow_sheet.feed.flow_rate',Q,time=0)
process.add_event('feed_off','flow_sheet.feed.flow_rate',0,time=load_time)
# wash
process.add_event('wash_on','flow_sheet.wash.flow_rate',Q,time=load_time)
process.add_event('wash_off','flow_sheet.wash.flow_rate',0,time=exp_time)
# cycle time
process.cycle_time = exp_time # seconds
print(Q,'m^3/s')
# -----------------------------
# simluate the process
# -----------------------------
if __name__ == '__main__':
from CADETProcess.simulator import Cadet
process_simulator = Cadet()
process_simulator.time_resolution = 0.01
simulation_results = process_simulator.simulate(process)
print('simulation finished')
# save x- and y- axis information
x_axis_sim_init = simulation_results.solution.outlet.outlet.time
x_axis_sim_min = x_axis_sim_init/60.
y_axis_sim_init = simulation_results.solution.outlet.outlet.solution
y_axis_sim_init = y_axis_sim_init[:,1]
#initial_y_max= np.max(y_axis_sim_init)
plt.plot(x_axis_sim_min,y_axis_sim_init,'-r',label='initial simulation')
plt.ylabel('mMol')
plt.xlabel('min')
plt.legend()
File produced by conda env export > environment.yml
environment.yml (7.56 KB)