Error when completing simulation with polynomial flow ramp

Hi there,

I am having some issues with polynomial flow terms. I am trying to run a simple bind and elute experiment using the multistate steric mass action model in CADET-Process (dev). As part of this I am trying to incorporate some pump flow ramps. Currently I have just implemented a step ramp using the code at the bottom of this message and the simulation runs fine. However when I change one line of code

from :
process.add_event('load_pump_ramp', 'flow_sheet.inlet.flow_rate', (5/(60*1e6)), 0)
to
process.add_event('load_pump_ramp', 'flow_sheet.inlet.flow_rate', [0,5e-9], 0)

to incorporate a linear ramp instead of a step, my events look good in that it appears to have had the desired effect but when i come to run the actual; simulation it doesnt complete and the kernel gets stuck on ‘busy’.

I tried updating the CADET dev branch by running ‘pip install git+https://github.com/fau-advanced-separations/cadet-process.git@dev’ , but this didn’t work either.

import numpy as np
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import MultistateStericMassAction
from CADETProcess.processModel import Inlet,Cstr,TubularReactor,  LumpedRateModelWithoutPores, Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.processModel import ComponentSystem

# Component System
component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('A')
component_system.add_component('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, 100,100,100,100]
binding_model.desorption_rate = [0.0, 30,30,3,3]
binding_model.conversion_rate = [0,0,10,100,10,0,0,0,0]
binding_model.characteristic_charge = [0.0, 10,10,25,25]
binding_model.steric_factor = [0.0, 2,2,2,2]
binding_model.capacity = 300.0
binding_model.reference_liquid_phase_conc = 1500
binding_model.reference_solid_phase_conc = 300

# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet_2 = Inlet(component_system, name='inlet_2')
tubing= TubularReactor(component_system, 'tubing')
tubing.length=1
tubing.diameter=0.001
tubing.axial_dispersion=0.004
tubing.c=[150,0,0]

probe= Cstr(component_system, name='probe')
probe.c=[150,0,0]
probe.V=2e-07

connector= TubularReactor(component_system, name='connector')
connector.length=3e-3
connector.diameter=0.02
connector.axial_dispersion= 0
connector.c=[150,0,0]

mixer= Cstr(component_system, name='mixer')
mixer.V=1e-06
mixer.c=[150,0,0]

column = LumpedRateModelWithoutPores(component_system, name='column')
column.binding_model = binding_model
column.cross_section_area = 3.5e-4
column.length = (1.005e-6)/(column.cross_section_area)
column.total_porosity = 0.8
column.axial_dispersion = 1e-8
column.c = [150, 0,0]
column.q = [binding_model.capacity, 0,0,0,0]

conductivity= Cstr(component_system, name='conductivity')
conductivity.V=2E-07
conductivity.c=[150,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(probe)
flow_sheet.add_unit(tubing)
flow_sheet.add_unit(connector)
flow_sheet.add_unit(mixer)
flow_sheet.add_unit(conductivity)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet, product_outlet=True)

flow_sheet.add_connection(inlet, tubing)
flow_sheet.add_connection(inlet_2, tubing)
flow_sheet.add_connection(tubing,probe)
flow_sheet.add_connection(probe,connector)
flow_sheet.add_connection(connector,mixer)
flow_sheet.add_connection(mixer,column)
flow_sheet.add_connection(column, conductivity)
flow_sheet.add_connection(conductivity, outlet)

# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 14*60

load_duration = 2.2*60
t_gradient_start = 4.5*60
grad_end=(10*60)
gradient_duration = grad_end - t_gradient_start

c_load = np.array([150, 0.00125,0.00125])
c_wash = np.array([150.0, 0.0,0])
c_elute = np.array([1350.0, 0.0,0])
gradient_slope = (c_elute - c_wash)/gradient_duration
c_gradient_poly = np.array(list(zip(c_wash, gradient_slope)))

process.add_event('load_pump_ramp', 'flow_sheet.inlet.flow_rate',  (5/(60*1e6)), 0)
process.add_event('load', 'flow_sheet.inlet.c', c_load,0)
process.add_event('load_flow_on', 'flow_sheet.inlet.flow_rate',  (65/(60*1e6)), 40)
process.add_event('buffer_flow_off', 'flow_sheet.inlet_2.flow_rate',  (0/(60*1e6)), dependencies=['load_pump_ramp'])
process.add_event('load_flow_off', 'flow_sheet.inlet.flow_rate',  0, load_duration )
process.add_event('buffer_flow_on', 'flow_sheet.inlet_2.flow_rate',  (10/(60*1e6)), 150)

process.add_event('wash', 'flow_sheet.inlet_2.c',  c_wash, 0)

process.add_event('grad_start', 'flow_sheet.inlet_2.c', c_gradient_poly, t_gradient_start)
process.add_event('grad_end', 'flow_sheet.inlet_2.c',c_elute, grad_end)



_ = process.plot_events()

if __name__ == '__main__':
    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 = '$c_{salt}$'

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

Hey George,
thanks for reporting. I can reproduce the bug and will investigate further.
Cheers

Edit: bug is now tracked on Github.

1 Like

Hey George,

we found the issue but it’s not trivial to fix. It has to do with consistent initialization when using polynomial flow rates. Currently, CADET-Core seems to struggle when the initial flow rate is 0. So until we have a proper fix, maybe you could start with a very small constant value for the ramp. It shouldn’t change the outlet profile by any substantial amount but it will make the simulation run. This way, you could at least continue until we have a proper solution.

process.add_event('load_pump_ramp', 'flow_sheet.inlet.flow_rate', [1e-12,5e-9], 0)
1 Like

Amazing, thanks Johannes!

1 Like