Simulation doesn't converge with SMA and non-integer characteristic charge

Hello everyone,

I used the Yamamoto method in Cadet Process, resulting in an optimum of approximately 1.69 for the characteristic charge. Unfortunately, the model does not converge afterwards. However, if I then arbitrarily set the characteristic charge to 1, 2 or another natural number, it converges. The model consistently fails to converge when the characteristic charge is not an integer.

Could someone help me?

Good evening Stefanos,

can you share a minimum reproducible example including your code and input data?

Hello Ronald,

many thanks for the quick reply and I apologise for not having sent a minimal example before.

Here is the example:

from CADETProcess.processModel import ComponentSystem
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 

component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein')


binding_model = StericMassAction(component_system, name='SMA')
binding_model.is_kinetic = False
binding_model.adsorption_rate = [1, 1]
binding_model.desorption_rate = [1, 1]
binding_model.steric_factor = [1, 1]
binding_model.capacity = 800 
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]


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', [0, 1],0)
_ = process.add_event('Wash', 'flow_sheet.inlet.c', [0, 0], wash_start)
_ = process.add_event(
    'grad_start',
    'flow_sheet.inlet.c',
    [[0, gradient_slope[0]], [0, gradient_slope[1]]],
    gradient_start
)
from CADETProcess.simulator import Cadet

simulator = Cadet() 
sim_res = simulator.simulate(process)

As soon as I change the characeristic charge of the protein to a natural number the model converges.

Best regards
Stefanos

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)

1 Like

Hello Ronald,
thank you very much for the very quick and detailed explanation! :slight_smile:

The SMA model also is not really mechanistically valid at salt concentration c_s = 0. It is based on stoichiometric displacement where salt ions displace adsorbed proteins. There are other models such as Langmuir that can be used with c_s = 0.

3 Likes