def get_cadet_template(n_units=3, split_components_data=True):
cadet_template = Cadet()
cadet_template.root.input.model.nunits = n_units
# Store solution
cadet_template.root.input['return'].split_components_data = split_components_data
cadet_template.root.input['return'].split_ports_data = 0
cadet_template.root.input['return'].unit_000.write_solution_inlet = 1
cadet_template.root.input['return'].unit_000.write_solution_outlet = 1
cadet_template.root.input['return'].unit_000.write_solution_bulk = 1
cadet_template.root.input['return'].unit_000.write_solution_particle = 1
cadet_template.root.input['return'].unit_000.write_solution_solid = 1
cadet_template.root.input['return'].unit_000.write_solution_flux = 1
cadet_template.root.input['return'].unit_000.write_solution_volume = 1
cadet_template.root.input['return'].unit_000.write_coordinates = 1
cadet_template.root.input['return'].unit_000.write_sens_outlet = 1
for unit in range(n_units):
cadet_template.root.input['return']['unit_{0:03d}'.format(unit)] = cadet_template.root.input['return'].unit_000
# Tolerances for the time integrator
cadet_template.root.input.solver.time_integrator.abstol = 1e-6
cadet_template.root.input.solver.time_integrator.algtol = 1e-10
cadet_template.root.input.solver.time_integrator.reltol = 1e-6
cadet_template.root.input.solver.time_integrator.init_step_size = 1e-6
cadet_template.root.input.solver.time_integrator.max_steps = 1000000
# Solver settings
cadet_template.root.input.model.solver.gs_type = 1
cadet_template.root.input.model.solver.max_krylov = 0
cadet_template.root.input.model.solver.max_restarts = 10
cadet_template.root.input.model.solver.schur_safety = 1e-8
# Run the simulation on single thread
cadet_template.root.input.solver.nthreads = 1
return cadet_template
def set_discretization(model, n_bound=None, n_col=20, n_par_types=1):
columns = {'RADIAL_LUMPED_RATE_MODEL_WITH_PORES', 'GENERAL_RATE_MODEL', 'LUMPED_RATE_MODEL_WITH_PORES', 'LUMPED_RATE_MODEL_WITHOUT_PORES'}
for unit_name, unit in model.root.input.model.items():
if 'unit_' in unit_name and unit.unit_type in columns:
unit.discretization.ncol = n_col
unit.discretization.npar = 5
unit.discretization.npartype = n_par_types
if n_bound is None:
n_bound = unit.ncomp*[0]
unit.discretization.nbound = n_bound
unit.discretization.par_disc_type = 'EQUIDISTANT_PAR'
unit.discretization.use_analytic_jacobian = 1
unit.discretization.reconstruction = 'WENO'
unit.discretization.gs_type = 1
unit.discretization.max_krylov = 0
unit.discretization.max_restarts = 10
unit.discretization.schur_safety = 1.0e-8
unit.discretization.weno.boundary_model = 0
unit.discretization.weno.weno_eps = 1e-10
unit.discretization.weno.weno_order = 3
def run_simulation(cadet, file_name=None):
if file_name is None:
f = next(tempfile._get_candidate_names())
cadet.filename = os.path.join(tempfile.tempdir, f + '.h5')
else:
cadet.filename = file_name
# save the simulation
cadet.save()
# run the simulation and load results
data = cadet.run()
cadet.load()
# Remove files
if file_name is None:
os.remove(os.path.join(tempfile.tempdir, f + '.h5'))
# Raise error if simulation fails
if data.returncode == 0:
print("Simulation completed successfully")
else:
print(data)
raise Exception("Simulation failed")
#paper 2020 Hagemann
t_cycle = 12500
langmuir_model = get_cadet_template(n_units=6)
# Sections and Switches
langmuir_model.root.input.solver.sections.nsec = 1
langmuir_model.root.input.solver.sections.section_times = [0.0, t_cycle]
langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [5.33e-3]
#set the times that the simulator writes out data for
langmuir_model.root.input.solver.user_solution_times = np.linspace(0, t_cycle, int(t_cycle) + 1)
n_comp = 1
# INLET
langmuir_model.root.input.model.unit_000.unit_type = 'INLET'
langmuir_model.root.input.model.unit_000.ncomp = n_comp
langmuir_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
# GENERAL_RATE_MODEL
langmuir_model.root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL'
langmuir_model.root.input.model.unit_001.ncomp = n_comp
langmuir_model.root.input.model.unit_001.col_length = 55e-3#59e-3
langmuir_model.root.input.model.unit_001.cross_section_area = 2.1e-6
langmuir_model.root.input.model.unit_001.col_porosity = 0.25
langmuir_model.root.input.model.unit_001.par_porosity = 0.8
langmuir_model.root.input.model.unit_001.par_radius = 46.8e-6
langmuir_model.root.input.model.unit_001.col_dispersion = 5.99e-6
langmuir_model.root.input.model.unit_001.film_diffusion = [1.96e-6]
langmuir_model.root.input.model.unit_001.par_diffusion = 5.96e-12#[6.5e-12]
langmuir_model.root.input.model.unit_001.par_surfdiffusion = [0]
langmuir_model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
langmuir_model.root.input.model.unit_001.adsorption.is_kinetic = False
langmuir_model.root.input.model.unit_001.adsorption.mcl_ka = [25574]
langmuir_model.root.input.model.unit_001.adsorption.mcl_kd = [1]
langmuir_model.root.input.model.unit_001.adsorption.mcl_qmax = [0.5433]
langmuir_model.root.input.model.unit_001.init_c = [0.0]
langmuir_model.root.input.model.unit_001.init_q = [0.0]
# OUTLET
langmuir_model.root.input.model.unit_002.unit_type = 'OUTLET'
langmuir_model.root.input.model.unit_002.ncomp = n_comp
# Mixer
langmuir_model.root.input.model.unit_003.unit_type = 'CSTR'
langmuir_model.root.input.model.unit_003.ncomp = n_comp
langmuir_model.root.input.model.unit_003.init_volume =0.95e-7
langmuir_model.root.input.model.unit_003.init_c = n_comp*[0.0]
## Tubing-inlet
langmuir_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
langmuir_model.root.input.model.unit_004.ncomp = n_comp
langmuir_model.root.input.model.unit_004.col_length = 1.8#27e-2
langmuir_model.root.input.model.unit_004.cross_section_area = 3.14*((0.75e-3)/2)**2
langmuir_model.root.input.model.unit_004.total_porosity = 1.0
langmuir_model.root.input.model.unit_004.col_dispersion = 3.5e-4
langmuir_model.root.input.model.unit_004.adsorption_model = 'NONE'
langmuir_model.root.input.model.unit_004.nbound = [1]*n_comp
langmuir_model.root.input.model.unit_004.init_c = n_comp*[0.0]
langmuir_model.root.input.model.unit_004.init_q = n_comp*[0.0]
set_discretization(langmuir_model, n_col=40)
langmuir_model.root.input['return'].unit_004 = langmuir_model.root.input['return'].unit_001
## Tubing-otlet
langmuir_model.root.input.model.unit_005.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
langmuir_model.root.input.model.unit_005.ncomp = n_comp
langmuir_model.root.input.model.unit_005.col_length = 1.8#62e-2
langmuir_model.root.input.model.unit_005.cross_section_area = 3.14*((0.75e-3)/2)**2
langmuir_model.root.input.model.unit_005.total_porosity = 1.0
langmuir_model.root.input.model.unit_005.col_dispersion = 3.5e-4
langmuir_model.root.input.model.unit_005.adsorption_model = 'NONE'
langmuir_model.root.input.model.unit_005.nbound = [1]*n_comp
langmuir_model.root.input.model.unit_005.init_c = n_comp*[0.0]
langmuir_model.root.input.model.unit_005.init_q = n_comp*[0.0]
set_discretization(langmuir_model, n_col=40)
langmuir_model.root.input['return'].unit_005 = langmuir_model.root.input['return'].unit_002
## Discretization
set_discretization(langmuir_model, n_bound=[1,1])
# Connections
langmuir_model.root.input.model.connections.nswitches = 1
langmuir_model.root.input.model.connections.switch_000.section = 0
langmuir_model.root.input.model.connections.switch_000.connections = [
0, 3, -1, -1, 3.3333e-9,
3, 4, -1, -1, 3.3333e-9,
4, 1, -1, -1, 3.3333e-9,
1, 5, -1, -1, 3.3333e-9,
5, 2, -1, -1, 3.3333e-9
]
run_simulation(langmuir_model)
#experimental data
# Load data
file_path = "/Users/kmaalim/Documents/CADET/trial/data_hagemann.csv" # Update with the correct file path if needed
df = pd.read_csv(file_path)
# Extract columns
time = df["time"]
concentration = df["c (mM)"]
# Plot breakthrough curve
plt.figure(figsize=(8,5))
plt.plot(time, concentration, marker="o", linestyle="-", color="b", label="experiment")
plt.xlabel("Time (s)")
plt.ylabel("Concentration (mM)")
plt.title("Breakthrough Curve")
plt.legend()
plt.grid(False)
# Plot the simulation
time = langmuir_model.root.output.solution.solution_times
c = langmuir_model.root.output.solution.unit_001.solution_outlet_comp_000
plt.plot(time, c, linestyle="--", color="r", label="sim")
plt.title('Break through curve')
plt.xlabel('$time~/~min$')
plt.ylabel('$concentration~/~mM$')
plt.legend()
plt.show()
Hi Kauthar,
could you please clarify a bit what your actual question is, it is not clear from your code snippet.
Also, please make sure to provide code that we can actually reproduce. Since we do not have access to your file system, we cannot load the experimental data.
Finally, I would be curious as to why you are using CADET-Python directly. Is there something missing in CADET-Process that you need?
Thanks and best regards
Jo
Hi Jo,
I am getting output, which does not seem to be right compared to experimental data
Experimental data
data_hagemann.csv (625 Bytes)
Cadet Python code
ProtA_model_cadet_python.ipynb (85.9 KB)
The paper that I am replicating
Thanks for updating your post.
Before looking deeper into the code it would be interesting to know how you selected the model and how you determined the values of the model parameters.
Also, I’m still curious as to why you are using CADET-Python directly instead of CADET-Process.
Since this is a very broad topic, I would also recommend visiting our next Office Hours where we can give you some hands-on support. Unfortunately, you just missed the last date yesterday. Next one will be May 8th.
1 Like
I use the general rate model (GRM) for the mass transfer and diffusion, and adsorption onto the solid I selected the Langmuir isotherm. The values of the model parameters I just extracted from the paper