Stirring Tank instant bulk solution initial concentration change

Konstantin · July 19, 2022, 8:33pm

I’d like to ask how can I change the initial concentration of the bulk solution after some time but keep solid and pore concentration profiles from the previous test?

I have a stirring tank with a feed rate equal to 0 but I want to change solution consecration after some time. I see that there is an option to do cubic polynomial on feed concertation but in my case, it is always zero.

I tried to run independent consecutive tests but I can not set specific distribution on solid and pore concentrations among the particle radius, only the average.

Thanks!

j.schmoelder · July 20, 2022, 6:57am

Hi Konstantin,

the initial concentration can only be set at the start of the simulation. It’s not possible to directly modify concentrations later on. The cubic polynomial is to modify the boundary conditions of unit operations. However, this requires flow into the unit.

But there are still some hacks to achieve the desired behaviour.

You could use a large flow rate to essentially flush the CSTR with new substrate. But to be honest, I’d have to check if that also works if you only want to change individual component concentrations.
Alternatively, I recommend running a first simulation until your desired point. Then use the final state of the simulation as initial state of the next simulation, only changing the state of the CSTR.

Hope that helps.

If you’re interested, we also have our first Community Call later today. Feel free to join!

Konstantin · July 20, 2022, 3:23pm

j.schmoelder thanks!
How can I use “final state of the simulation as initial state of the next simulation”? Or how can I set up the concertation profile among the particle radius in the code?
Should I just put:

model.root.input.model.unit_001.init_c = new_init_c

and run it again like that:

run_simulation(model)

Will it keep concentration profile among the particle radius from my pervious run? or it create a new fresh instance of model?
Also, have GENERAL_RATE_MODEL in case if this makes the difference.
Thanks!

j.schmoelder · July 24, 2022, 9:26am

First, you have to tell CADET to write the full state of the CSTR unit at the end of the simulation by adding the return flag (see: Return data — CADET):

/input/return/unit_xxx/WRITE_SOLUTION_LAST_UNIT = 1

After simulation, you can now read the last state in the output group (see Output Group — CADET)

/output/solution/unit_XXX/LAST_STATE_Y

Then you can set overwrite the corresponding entries and set the initial state for the next simulation (see: Continuous stirred tank reactor model — CADET)

/input/model/unit_XXX/INIT_STATE

You have to be careful with the ordering of the data though and I’m not entirely sure about the internal structure. Maybe @s.leweke can comment here.

Konstantin · July 25, 2022, 9:24pm

Johannes thanks for the response!

So, I tried do stuff like that:

model.root.input['return'].unit_000.WRITE_SOLUTION_LAST_UNIT =1
utils.run_simulation(model) #-first simulation
vector = model.root.output.solution.unit_001.last_state_y
model.root.input.model.unit_001.INIT_STATE = vector
model.root.input.model.unit_001.init_c = [new_balk_concentration,]
utils.run_simulation(model) #-second simulation

But I don’t see bulk concentration is equal to “new_balk_concentration” at the beginning of second run, it actually goes down to zero.
Should I look what is inside of that “vector” array? It is hard to figure out what to change there.
Thanks,
K

s.leweke · July 25, 2022, 9:39pm

If you specify INIT_STATE, the array provided there has precedence over INIT_C. This is the reason why the line

model.root.input.model.unit_001.init_c = [new_balk_concentration,]

is ignored.

For the CSTR, the layout of the state is

concentration liquid phase components, concentration bound states, volume.

Hence, you should do something like

model.root.input['return'].unit_000.WRITE_SOLUTION_LAST_UNIT =1
utils.run_simulation(model) #-first simulation
vector = model.root.output.solution.unit_001.last_state_y
vector[:len(new_balk_concentration)] = new_balk_concentration
model.root.input.model.unit_001.INIT_STATE = vector
utils.run_simulation(model) #-second simulation

Konstantin · July 26, 2022, 3:48pm

Samuel,
It worked! Even for General Rate Model (GRM) with flow rate zero that I have, not the CSTR.
Thanks!
K

Konstantin · October 5, 2022, 8:17am

Samuel,
I have a problem when do the same trick for the third simulation, it says:

File “h5py\h5o.pyx”, line 202, in h5py.h5o.link.
OSError: Unable to create link (name already exists).

It happens when it does model.save()
I have to run four cycles, but it gives me this exception after second one.
Please advise!
Thanks!
K

s.leweke · October 7, 2022, 8:24am

Hey Konstantin,

I’m not sure what’s the cause of this problem. If the file already exists, it should be overwritten / edited. Maybe the file is still open in a different place in the code / there’s another handle to the file?

As @j.schmoelder and @w.heymann are more familiar with the Python side of things, I hope they can help here.

Konstantin · October 7, 2022, 8:05pm

Samuel,
Thanks for your response!
I put my script here/below, in case if I did something wrong.
It crashes on last model.save(). I tried to do new .h5 file before but it does not help.
K

import numpy as np
from cadet import Cadet
Cadet.cadet_path = 'C:\\Users\\konstanting\\Anaconda3\\pkgs\\cadet-4.3.0-hbe28382_2\\bin'
import utils
total_time =300*60 # s
time_frame = 1000         



col_length =  0.01 # m
cross_section_area = 0.1 # m2 

col_porosity = 0.996168582 

par_porosity = 0.44 # see 
 
par_radius = 3.25E-04 #m 
  
col_dispersion = 3.789E-6 #m2/s 

film_diffusion = [6.034E-02] #m/s 

par_diffusion = [(1.62E-09),] # m2/s 

par_surfdiffusion =   [5e-13,] # m^2 / s (solid phase)

k_factor = 1
mcl_ka = [10.58735191/k_factor,]    
mcl_kd = [1/k_factor,] # was just 1
mcl_qmax = [22471.33561/(mcl_ka[0]*k_factor),] #mol / m^3
is_kinetic = False


concentration_ini = 1.96E-03 #initial concentration in the solution mol/m3
init_q = 0 #initial concentration in solids mol/m3

const_coeff = [concentration_ini,] # concentration mol/m3
Q = 0   # Flow rate in m^3*s^-1


#Model
model = Cadet()
model.root.input.model.nunits = 3

#Inlet Model
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 1
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'

#General Rate Model
model.root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL'
model.root.input.model.unit_001.ncomp = 1

## Geometry
model.root.input.model.unit_001.col_length = col_length                  # m,  
model.root.input.model.unit_001.cross_section_area = cross_section_area  #  m^2,  

model.root.input.model.unit_001.col_porosity = col_porosity     # -
model.root.input.model.unit_001.par_porosity = par_porosity     # -
model.root.input.model.unit_001.par_radius = par_radius         # m  

                                                            
## Transport
model.root.input.model.unit_001.col_dispersion = col_dispersion # m^2 / s (interstitial volume) 
model.root.input.model.unit_001.film_diffusion = film_diffusion # m / s 
model.root.input.model.unit_001.par_diffusion = par_diffusion   # m^2 / s (mobile phase) 
model.root.input.model.unit_001.par_surfdiffusion = par_surfdiffusion #  m^2 / s (solid phase)


# # Adsorption Model
#model.root.input.model.unit_001.adsorption_model = 'NONE' 
model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
model.root.input.model.unit_001.adsorption.is_kinetic = is_kinetic   # not the Kinetic binding
model.root.input.model.unit_001.adsorption.mcl_ka = mcl_ka    # m^3 / (mol * s)   (mobile phase)
model.root.input.model.unit_001.adsorption.mcl_kd = mcl_kd       # 1 / s (desorption)
model.root.input.model.unit_001.adsorption.mcl_qmax = mcl_qmax  # mol / m^3   (solid phase)


# Initial Conditions
model.root.input.model.unit_001.init_c = [concentration_ini,]
model.root.input.model.unit_001.init_q = [init_q,]

#Setting up the Discretization

### Grid cells
model.root.input.model.unit_001.discretization.ncol = 1 
model.root.input.model.unit_001.discretization.npar = 500 

### Bound states
model.root.input.model.unit_001.discretization.nbound = [1] 

### Other options
model.root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR'    
model.root.input.model.unit_001.discretization.use_analytic_jacobian = 1
model.root.input.model.unit_001.discretization.reconstruction = 'WENO'
model.root.input.model.unit_001.discretization.gs_type = 1
model.root.input.model.unit_001.discretization.max_krylov = 0
model.root.input.model.unit_001.discretization.max_restarts = 10
model.root.input.model.unit_001.discretization.schur_safety = 1.0e-8

model.root.input.model.unit_001.discretization.weno.boundary_model = 0
model.root.input.model.unit_001.discretization.weno.weno_eps = 1e-10
model.root.input.model.unit_001.discretization.weno.weno_order = 3

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

#Time Sections

model.root.input.solver.sections.nsec = 1
model.root.input.solver.sections.section_times = [0.0,total_time,]   # s
model.root.input.solver.sections.section_continuity = []

model.root.input.model.unit_000.sec_000.const_coeff = const_coeff # mol / m^3
model.root.input.model.unit_000.sec_000.lin_coeff = [0.0,]
model.root.input.model.unit_000.sec_000.quad_coeff = [0.0,]
model.root.input.model.unit_000.sec_000.cube_coeff = [0.0,]


#System Connectivity
                                            
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,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, Q]  # [unit_001, unit_002, all components, all components, Q/ m^3*s^-1 

#Setting Up the Simulator and Running the Simulation

model.root.input.model.solver.gs_type = 1
model.root.input.model.solver.max_krylov = 0
model.root.input.model.solver.max_restarts = 10
model.root.input.model.solver.schur_safety = 1e-8

# Number of cores for parallel simulation
model.root.input.solver.nthreads = 1

# Tolerances for the time integrator
model.root.input.solver.time_integrator.abstol = 1e-6
model.root.input.solver.time_integrator.algtol = 1e-10
model.root.input.solver.time_integrator.reltol = 1e-6
model.root.input.solver.time_integrator.init_step_size = 1e-6
model.root.input.solver.time_integrator.max_steps = 1000000


# Return data
model.root.input['return'].split_components_data = 0
model.root.input['return'].split_ports_data = 0
model.root.input['return'].unit_000.write_solution_bulk = 1
model.root.input['return'].unit_000.write_solution_inlet = 1
model.root.input['return'].unit_000.write_solution_outlet = 1
model.root.input['return'].unit_000.write_coordinates = 1 
model.root.input['return'].unit_000.write_solution_solid = 1 
model.root.input['return'].unit_000.write_solution_particle =1 
model.root.input['return'].unit_000.write_solution_volume  = 1 

model.root.input['return'].unit_000.WRITE_SOLUTION_LAST_UNIT = 1

# Copy settings to the other unit operations
model.root.input['return'].unit_001 = model.root.input['return'].unit_000
model.root.input['return'].unit_002 = model.root.input['return'].unit_001


# Solution times
model.root.input.solver.user_solution_times = np.linspace(0, total_time, 1001)


model.filename = 'C:\\Users\\konstanting\\modelA.h5'
model.save()

# Simulation A
data = model.run()
model.load()

stateAfterA = model.root.output.solution.unit_001.last_state_y

# Simulation B
new_balk_concentration = [2.103E-02,]
stateAfterA[:len(new_balk_concentration)] = new_balk_concentration # inject new_balk_solution in index 1
model.root.input.model.unit_001.INIT_STATE = stateAfterA

model.filename = 'C:\\Users\\konstanting\\modelB.h5'
model.save()
data = model.run()
model.load()
stateAfterB = model.root.output.solution.unit_001.last_state_y
c_out = model.root.output.solution.unit_001.solution_outlet
solidConcentration =     model.root.output.solution.unit_001.solution_solid[-1, 0, :, 0]

# Simulation C
new_balk_concentration = [1.769E-01,]
stateAfterB[:len(new_balk_concentration)] = new_balk_concentration # inject new_balk_solution in index 1
model.root.input.model.unit_001.INIT_STATE = stateAfterB

model.filename = 'C:\\Users\\konstanting\\modelC1.h5'
model.save() # Crashes here!!!, gives "OSError: Unable to create link (name already exists)"
data = model.run()
model.load()

l.thiel · October 10, 2022, 9:21am

Hey Konstantin,

the error arises from the initial state specification for the simulations B and C. Johannes opened an issue on Github related to that topic and we will work on that.

To get your model working, you can not use capital letters for that. It needs to be:

# Simulation B
model.root.input.model.unit_001.init_state= stateAfterA

# Simulation C
model.root.input.model.unit_001.init_state = stateAfterB

For more information, see also here.

Konstantin · October 10, 2022, 3:23pm

Lukas,
It worked!
Thanks for your help!
K