Simulation not running with 'LumpedRateModelWithoutPores' in latest CADET Process Dev

Hi There,

I recently installed the latest branch of CADET Process using the following commands from an environment with CADET-Process 8.0 previously installed.

pip uninstall cadet-process
pip install git+https://github.com/fau-advanced-separations/CADET-Process.git@dev

When I run a simulation for the MSSMA using the general rate model , the simulation runs fine. (First Code below). However if I simply change the column model to lumped rate without pores the simulation doesn’t complete (second code below). Are people able to replicate this error as I can’t tell if this is an installation issue on my end. I have an environment with an older branch of CADET that seems to run fine but I was unable to run multi simulation optimisations using that so had to upgrade, as reported here.

GeneralRateModel CODE

import numpy as np

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import MultistateStericMassAction
from CADETProcess.processModel import Inlet, GeneralRateModel,TubularReactor, Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process

# Component System
component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein A')
component_system.add_component('Protein B')

# Binding Model
binding_model = MultistateStericMassAction(component_system,bound_states=[1,2,2], name='MSSMA')
binding_model.is_kinetic = True
binding_model.adsorption_rate = [0.0, 20,30,1,10]
binding_model.desorption_rate = [0.0, 3,30,3,3]
binding_model.conversion_rate = [0,100,10,100,100,100,100,100,100]
binding_model.characteristic_charge = [0.0, 2,4,6,8]
binding_model.steric_factor = [0.0, 2,2,2,2]
binding_model.capacity = 300.0
binding_model.reference_liquid_phase_conc = 1350
binding_model.reference_solid_phase_conc = 300

# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = 0/(60*1e6)
inlet_2 = Inlet(component_system, name='inlet_2')
inlet_2.flow_rate = 0

V_column_void_in= TubularReactor(component_system, name='V_column_void_in')
V_column_void_in.length=0.003545772724442985
V_column_void_in.diameter=0.0215
V_column_void_in.axial_dispersion= 0
V_column_void_in.c=[150,0,0]

column = GeneralRateModel(component_system, name='column')
column.binding_model = binding_model

column.length = 0.25
column.diameter = 0.0115
column.bed_porosity = 0.37
column.particle_radius = 4.5e-5
column.particle_porosity = 0.33
column.axial_dispersion = 2.0e-7
column.film_diffusion = [2.0e-5, 2.0e-7, 2.0e-7]
column.pore_diffusion = [7e-5, 1e-9, 1e-9]
column.surface_diffusion = [0.0, 0.0, 0.0,0,0]

outlet = Outlet(component_system, name='outlet')

# Flow Sheet
flow_sheet = FlowSheet(component_system)

flow_sheet.add_unit(inlet)
flow_sheet.add_unit(inlet_2)
flow_sheet.add_unit(V_column_void_in)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet, product_outlet=True)

flow_sheet.add_connection(inlet, column)
flow_sheet.add_connection(inlet_2, V_column_void_in)
flow_sheet.add_connection(V_column_void_in, column)
flow_sheet.add_connection(column, outlet)

# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 19000.0
gradient_end = 15000.0
## Create Events and Durations
wash_start = 7500.0
gradient_start = 9500.0
concentration_difference = np.array([1350.0, 0.0]) - np.array([150.0, 0.0])
gradient_duration = gradient_end - gradient_start
gradient_slope = concentration_difference/gradient_duration

_ = process.add_event('load', 'flow_sheet.inlet.c', [150.0,((0.000016275/4.17)*2.36), ((0.000016275/4.17)*(4.17-2.36))],0)
_ = process.add_event('load_flow_on', 'flow_sheet.inlet.flow_rate',10/(60*1e6),0 )
_ = process.add_event('wash_flow_set', 'flow_sheet.inlet_2.flow_rate',0/(60*1e6),0 )
_ = process.add_event('load_flow_off', 'flow_sheet.inlet.flow_rate',0,wash_start,)
_ = process.add_event('wash', 'flow_sheet.inlet_2.c', [150.0, 0.0, 0.0], 0)
_ = process.add_event('wash_flow', 'flow_sheet.inlet_2.flow_rate',10/(60*1e6) , wash_start,dependencies=['load_flow_off'] )
_ = process.add_event(
    'grad_start',
    'flow_sheet.inlet_2.c',
    [[150.0, gradient_slope[0]], [0, gradient_slope[1]], [0, gradient_slope[1]]],
    gradient_start
)
_ = process.add_event('gradient_end', 'flow_sheet.inlet_2.c', [1350.0, 0.0, 0.0], gradient_end)

column.c = [150, 0, 0]
column.q = [binding_model.capacity, 0, 0,0,0]
_ = process.plot_events()

from CADETProcess.simulator import Cadet
process_simulator = Cadet()

simulation_results = process_simulator.simulate(process)

from CADETProcess.plotting import SecondaryAxis
sec = SecondaryAxis()
sec.components = ["Salt"]
sec.y_label = 'csalt'

_ = simulation_results.solution.column.outlet.plot(secondary_axis=sec)

Lumped rate without pores code

import numpy as np

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import MultistateStericMassAction
from CADETProcess.processModel import Inlet, LumpedRateModelWithoutPores,TubularReactor, Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process

# Component System
component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein A')
component_system.add_component('Protein B')

# Binding Model
binding_model = MultistateStericMassAction(component_system,bound_states=[1,2,2], name='MSSMA')
binding_model.is_kinetic = True
binding_model.adsorption_rate = [0.0, 20,30,1,10]
binding_model.desorption_rate = [0.0, 3,30,3,3]
binding_model.conversion_rate = [0,100,10,100,100,100,100,100,100]
binding_model.characteristic_charge = [0.0, 2,4,6,8]
binding_model.steric_factor = [0.0, 2,2,2,2]
binding_model.capacity = 300.0
binding_model.reference_liquid_phase_conc = 1350
binding_model.reference_solid_phase_conc = 300

# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = 0/(60*1e6)
inlet_2 = Inlet(component_system, name='inlet_2')
inlet_2.flow_rate = 0

V_column_void_in= TubularReactor(component_system, name='V_column_void_in')
V_column_void_in.length=0.003545772724442985
V_column_void_in.diameter=0.0215
V_column_void_in.axial_dispersion= 0
V_column_void_in.c=[150,0,0]

column = LumpedRateModelWithoutPores(component_system, name='column')
column.binding_model = binding_model

column.length = 0.25
column.diameter = 0.0115
column.total_porosity = 0.33
column.axial_dispersion = 2.0e-7


outlet = Outlet(component_system, name='outlet')

# Flow Sheet
flow_sheet = FlowSheet(component_system)

flow_sheet.add_unit(inlet)
flow_sheet.add_unit(inlet_2)
flow_sheet.add_unit(V_column_void_in)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet, product_outlet=True)

flow_sheet.add_connection(inlet, column)
flow_sheet.add_connection(inlet_2, V_column_void_in)
flow_sheet.add_connection(V_column_void_in, column)
flow_sheet.add_connection(column, outlet)

# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 19000.0
gradient_end = 15000.0
## Create Events and Durations
wash_start = 7500.0
gradient_start = 9500.0
concentration_difference = np.array([1350.0, 0.0]) - np.array([150.0, 0.0])
gradient_duration = gradient_end - gradient_start
gradient_slope = concentration_difference/gradient_duration

_ = process.add_event('load', 'flow_sheet.inlet.c', [150.0,((0.000016275/4.17)*2.36), ((0.000016275/4.17)*(4.17-2.36))],0)
_ = process.add_event('load_flow_on', 'flow_sheet.inlet.flow_rate',10/(60*1e6),0 )
_ = process.add_event('wash_flow_set', 'flow_sheet.inlet_2.flow_rate',0/(60*1e6),0 )
_ = process.add_event('load_flow_off', 'flow_sheet.inlet.flow_rate',0,wash_start,)
_ = process.add_event('wash', 'flow_sheet.inlet_2.c', [150.0, 0.0, 0.0], 0)
_ = process.add_event('wash_flow', 'flow_sheet.inlet_2.flow_rate',10/(60*1e6) , wash_start,dependencies=['load_flow_off'] )
_ = process.add_event(
    'grad_start',
    'flow_sheet.inlet_2.c',
    [[150.0, gradient_slope[0]], [0, gradient_slope[1]], [0, gradient_slope[1]]],
    gradient_start
)
_ = process.add_event('gradient_end', 'flow_sheet.inlet_2.c', [1350.0, 0.0, 0.0], gradient_end)

column.c = [150, 0, 0]
column.q = [binding_model.capacity, 0, 0,0,0]
_ = process.plot_events()

from CADETProcess.simulator import Cadet
process_simulator = Cadet()

simulation_results = process_simulator.simulate(process)

from CADETProcess.plotting import SecondaryAxis
sec = SecondaryAxis()
sec.components = ["Salt"]
sec.y_label = 'csalt'

_ = simulation_results.solution.column.outlet.plot(secondary_axis=sec)

Hi George,

cannot reproduce, runs fine for me on dev, seems like an installation issue. I do realize that we need to make a new release for CADET-Process which contains a lot of the recent bug fixes and I hope we can do that rather sooner than later but it takes some time…

Cheers

Jo

1 Like

Thanks Jo,

That’s good will have a look at the installation. Thank you !