Hey Stefanos,
the problem was, that the SMA model doesn’t work well with salt concentrations of exactly 0. My theory is that due to numerical noise concentrations can sometimes dip below 0 and would then get raised to a non-integer, which doesn’t work. To solve this, please set up your model to contain a bit of salt before starting the simulation.
The two changes are: adding the salt concentration to your column.c and column.cp initialization states and also initialize the column with all binding sites saturated with salt:
starting_salt_conc = 1 # or any other non-zero number
column.c = [starting_salt_conc, 0]
column.cp = [starting_salt_conc, 0]
column.q = [binding_model.capacity, 0]
and changing the inlet concentrations to also use non-zero concentrations:
_ = process.add_event('Load', 'flow_sheet.inlet.c', [starting_salt_conc, 1], 0)
_ = process.add_event('Wash', 'flow_sheet.inlet.c', [starting_salt_conc, 0], wash_start)
_ = process.add_event(
'grad_start',
'flow_sheet.inlet.c',
[[starting_salt_conc, gradient_slope[0]], [0, gradient_slope[1]]],
gradient_start
)
The complete working code is:
from CADETProcess.processModel import ComponentSystem, LumpedRateModelWithPores
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet, GeneralRateModel, Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
import numpy as np
Q_ml_min = 0.5 # ml/min
Q_m3_s = Q_ml_min / (60 * 1e6) # m3/s
V_tracer = 5e-8 # m3
starting_salt_conc = 1
component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein')
binding_model = StericMassAction(component_system, name='SMA')
binding_model.is_kinetic = True
binding_model.adsorption_rate = [0, 1]
binding_model.desorption_rate = [0, 1]
binding_model.steric_factor = [0, 1]
binding_model.capacity = 800
# binding_model.reference_solid_phase_conc = 800
# binding_model.reference_liquid_phase_conc = 1000
binding_model.characteristic_charge = [0, 1.6]
# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = Q_m3_s
column = GeneralRateModel(component_system, name='column')
column.binding_model = binding_model
column.length = 0.1
column.diameter = 0.0077
column.bed_porosity = 0.3
column.particle_radius = 3.4e-5
column.particle_porosity = 0.857
column.axial_dispersion = 2e-7
column.film_diffusion = [2.0e-3, 2.0e-5]
column.pore_diffusion = [7e-3, 7e-5]
column.c = [starting_salt_conc, 0]
column.cp = [starting_salt_conc, 0]
column.q = [starting_salt_conc, 0]
outlet = Outlet(component_system, name='outlet')
# Flow Sheet
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(inlet)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet)
flow_sheet.add_connection(inlet, column)
flow_sheet.add_connection(column, outlet)
# Process
process = Process(flow_sheet, 'lwe')
wash_start = V_tracer / Q_m3_s
gradient_start = (1200 + wash_start)
concentration_difference = np.array([1000, 0]) - np.array([0, 0])
gradient_volume = 6.11e-5
process.cycle_time = gradient_start + (gradient_volume / Q_m3_s)
gradient_duration = process.cycle_time - gradient_start
gradient_slope = concentration_difference / gradient_duration
_ = process.add_event('Load', 'flow_sheet.inlet.c', [starting_salt_conc, 1], 0)
_ = process.add_event('Wash', 'flow_sheet.inlet.c', [starting_salt_conc, 0], wash_start)
_ = process.add_event(
'grad_start',
'flow_sheet.inlet.c',
[[starting_salt_conc, gradient_slope[0]], [0, gradient_slope[1]]],
gradient_start
)
from CADETProcess.simulator import Cadet
from CADETProcess.plotting import SecondaryAxis
sec = SecondaryAxis()
sec.components = ['Salt']
sec.y_label = '$c_{salt}$'
simulator = Cadet()
sim_res = simulator.simulate(process)
sim_res.solution.outlet.outlet.plot(secondary_axis=sec)