Error when using quasi stationary binding

Hey everyone,

I am runnning into an issue as soon as I change the “is_kinetic” parameter to False or 0.
Is there any condition I am missing in order to use quasi stationary binding?
The function I am using and the error are pasted below.

Thank you in advance for your help!

def simulation_framework(n_comp, isotherm_name, component_name, L_c, D_c, A_c, V_c,  
                         por_total, Q_ms,V_inj, t_inj, t_elute,c_inj_mM, column_col_disp,
                         isotherm_parameters, n_col, n_par, plotting, inlet_profile):
    
    t_start = time.time()

    model = get_cadet_template(n_units=3, split_components_data=False)
    
    ## Sections
    section_intervals = [t_inj,  t_elute]
    section_times = [0]
    for index, section in enumerate(section_intervals): 
        section_times.append(section_times[index] + section) # s

    model.root.input.solver.sections.nsec = len(section_intervals)
    model.root.input.solver.sections.section_times = section_times

    ## Inlet Profile
    model.root.input.model.unit_000.sec_000.const_coeff = c_inj_mM
    model.root.input.model.unit_000.sec_001.const_coeff = [0]*n_comp
    
    # set the times that the simulator writes out data for
    model.root.input.solver.user_solution_times = np.linspace(0,  t_inj+  t_elute, 10001)

    ## Flowsheet
    # INLET 
    model.root.input.model.unit_000.unit_type = 'INLET'
    model.root.input.model.unit_000.ncomp = n_comp
    model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
    
    # Column
    model.root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
    model.root.input.model.unit_001.ncomp = n_comp
    model.root.input.model.unit_001.col_length = L_c
    model.root.input.model.unit_001.total_porosity = por_total
    model.root.input.model.unit_001.col_dispersion = column_col_disp
    model.root.input.model.unit_001.cross_section_area = A_c

    kA=isotherm_parameters[0]
    kD=isotherm_parameters[1]
    q_max=isotherm_parameters[2]
    
    model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
    model.root.input.model.unit_001.adsorption.is_kinetic = 0
    model.root.input.model.unit_001.adsorption.mcl_kA = kA
    model.root.input.model.unit_001.adsorption.mcl_kD = kD
    model.root.input.model.unit_001.adsorption.mcl_qmax = q_max 
    model.root.input.model.unit_001.init_c = [0]*n_comp
    model.root.input.model.unit_001.init_q = [0]*n_comp

    ## Outlet
    model.root.input.model.unit_002.unit_type = 'OUTLET'
    model.root.input.model.unit_002.ncomp = n_comp

    ## Switches and connections
    model.root.input.model.connections.nswitches = 1
    model.root.input.model.connections.switch_000.section = 0
    model.root.input.model.connections.switch_000.connections = [0, 1, -1, -1, Q_ms,
                                                                 1, 2, -1, -1, Q_ms
                                                                  ]

    ## Run simulation
    now2 = datetime.now()
    time_now = now2.strftime("%H:%M:%S")
    print(time_now+': Simulation started')
    set_discretization(model=model ,n_bound=n_comp*[1], n_col=n_col, n_par=n_par)
    run_simulation(model)
        
    ## Plots
    added_input_solutions = [0]* len(model.root.output.solution.unit_000.solution_outlet[:,0])
    added_output_solutions = [0]* len(model.root.output.solution.unit_002.solution_outlet[:,0])
    for i3 in range(n_comp):
        added_input_solutions = added_input_solutions+model.root.output.solution.unit_000.solution_outlet[:,i3]*M[i3]/1000
        added_output_solutions = added_output_solutions+model.root.output.solution.unit_002.solution_outlet[:,i3]*M[i3]/1000
           
                
    if plotting ==True:
        fontsize=12
        figsize=[25,12]
        fig3, ax3 = plt.subplots(figsize=cm2inch(figsize[0], figsize[1]), dpi=120)
        ax3.tick_params(axis='x', labelsize=fontsize)
        ax3.tick_params(axis='y', labelsize=fontsize)  
        if inlet_profile == True:
            ax3.plot(model.root.output.solution.solution_times/60, added_input_solutions, 'k', label=component_name +'combined simulated inlet' )        
        
        for i2 in range(n_comp):
            ax3.plot(model.root.output.solution.solution_times/60, model.root.output.solution.unit_002.solution_outlet[:,i2]*M[i2]/1000, color=colors[i2])
        
        ax3.plot(model.root.output.solution.solution_times/60, added_output_solutions, 'k--', label=component_name+' 15-30 combined simulated outlet')
        ax3.set_xlabel('Time (min)',fontsize=fontsize)
        ax3.set_ylabel(component_name + ' concentration (g/L)',fontsize=fontsize)

        lines, labels = ax3.get_legend_handles_labels()

        #ax3.legend(lines, labels,bbox_to_anchor=(1.13, 1),loc='upper left', borderaxespad=0)
        ax3.legend(lines, labels,fontsize=fontsize)
        #plt.subplots_adjust(right=0.7)
        plt.tight_layout()

    t_end = time.time()
    t_sim = format(round(t_end - t_start,1))
    print('Simulation duration =',t_sim,'seconds')
    
    return t_sim, model.root.output.solution.solution_times/60, added_output_solutions

Hi Jan-Christoph,

for quasi stationary bindings you should make sure that the consistent initialization mode is set to 1, i.e. full initialization (model.root.input.solver.CONSISTENT_INIT_MODE = 1).
Have you played around with adsorption parameters? If your simulation doesnt run for moderate parameters, something in the setup is most certainly wrong.
Can you upload an H5 file of a simulation that does not run so I can have a look at it?

Cheers
Jan

Hey Jan,
thanks for your reply!

I will include this line of code concerning the consistent initialization mode, but it doesn’t seem to fix the issue.
Basically I have determined an isotherm for my system and since I only have the equilibrium adsorption constant I want to use quasi stationary binding.
The determined parameters are:

q_max:
[45.54340912 42.76794019 40.31132087 38.12159098 36.15749867 34.385877
32.77975578 31.31697895 29.97917628 28.75098795 27.61947209 26.57364705
25.60413379 24.70287385 23.86290479 23.07818002]

kA:
[0.09103335 0.1239764 0.1684097 0.22574546 0.29968001 0.3991841
0.52960163 0.70465885 0.93170029 1.23003143 1.61832123 2.12644535
2.81337946 3.7184139 4.88759507 6.12738621]

kD:
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

For dynamic binding and the bold assumption kD=1 my simulation works just fine (yielding peaks that are way too broad). But I think that quasi-stationary binding is the way to go for my case.

The .h5 files for the non-working quasi-stationary case is attached.
(I noticed the .h5 file for the dynamic simulation is enormous (~1GB), maybe due to 16 components and a 200 minute experiment)

model_quasi_stationary.h5 (663.3 KB)

Cheers
Jan

It turns out that using a “lean initialization”, model.root.input.solver.CONSISTENT_INIT_MODE = 3, solves the issue. I will still look into why full consistent initialization doesnt work for this case.
However, since every concentration is zero inside the column at the beginning of your simulation, its sufficient to use this “lean” initialization (which does not initialize algebraic equations of rapid-equilibrium bindings).

Yes for this new initialisation mode the model works.
Thanks a lot for your help!

Best regards
Jan

It is worth mentioning that while setting k_d = 1 is seemingly nonphysical, it is convenient from a parametric standpoint. This way your k_a is equivalent to k_{eq}. When estimating these parameters from standard chromatographic experiments (e.g., breakthrough curves and/or elution) these parameters are inextricably correlated, so it is beneficial to fix k_d and fit k_a. Of course, it is unlikely that k_d = 1. These parameters can be determined using alternative approaches, such as the biophysical technique of surface plasmon resonance, which provides separable measurements of the kinetic on/off rates. However, this is an entirely separate task and I am only mentioning it for argument’s sake.

Another way to look at it is that you divide the terms in the adsorption isotherm’s kinetic form by k_d, which recovers a term on the equation’s left-hand-side; often called k_{kin} which is equal to \frac{1}{k_d}. This becomes a separate isotherm parameter that governs the peak width and can be regressed, or also set to 1, effectively ignoring it.

2 Likes