Issue with consistent initialization for linear isotherm model

[schmoelder]

During the workshop, we realized something odd about using the rapid equilibrium in adsorption models.

When setting the following initial conditions in a CSTR (without any in- or outgoing flow)

linear_model.root.input.model.unit_000.init_c = [1]
linear_model.root.input.model.unit_000.init_q = [0]

and setting is_kinetic to True

## Adsorption
linear_model.root.input.model.unit_000.nbound = [1]
linear_model.root.input.model.unit_000.adsorption_model = 'LINEAR'

linear_model.root.input.model.unit_000.adsorption.is_kinetic = True
linear_model.root.input.model.unit_000.adsorption.lin_ka = [2]
linear_model.root.input.model.unit_000.adsorption.lin_kd = [1]

We get a graph that looks as expected when plotting c and q over time:

However, when setting is_kinetic = False, it seems to ignore the initial conditions. Mass is created from nowhere (expected would be 0.66 in solid phase and 0.33 in bulk phase):

[schmoelder]

After some debugging, @Immudzen and I can report some new findings:

The issue also affects column models such as the LumpedRateModelWithoutPores. Concerningly, we could observe that there is also an issue when setting is_kinetic = True.

This indicated that it was an issue related to the isotherm, and not the CSTR.

However, when testing another isotherm model (i.e. MultiComponentLangmuir), we could not observe any issues. This led us to conclude that the issue is probably inherent to the Linear model.

In case of the CSTR, the problem seem to start here:
https://github.com/modsim/CADET/blob/a84506ea07213412af132626954448d309c7d6b7/src/libcadet/model/StirredTankModel.cpp#L604

For the linear model, this steps into the following function:
https://github.com/modsim/CADET/blob/aaa25379bc6393e8559728f49550b94c593c6c3a/src/libcadet/model/binding/LinearBinding.cpp#L587
In here, q is simply calculated from c without conserving the overall mass:
https://github.com/modsim/CADET/blob/aaa25379bc6393e8559728f49550b94c593c6c3a/src/libcadet/model/binding/LinearBinding.cpp#L609
It then returns false, causing the previous loop (in CSTR implementation) to continue without the initialization part which conserves "moieties" (interesting terminology, btw.; I've never heard it in this context).

In contrast, for Langmuir, the BindingModelBase class's function is used, which simply returns true:
https://github.com/modsim/CADET/blob/2cc210e2221fb02dfa1cd2b4ee0415c9eb070a8a/src/libcadet/model/binding/BindingModelBase.hpp#L81

---

To fix the issue, we would suggest to simply treat the linear model like all the others. @sleweke Could there be any side effects? Also, we have not tested further any other isotherm models yet and don't know whether this also affects other implementations.

[sleweke-bayer]

I kind of understand why the problem appears for IS_KINETIC = false. But I do not understand why there should be a problem when IS_KINETIC = true.

In the latter case, no consistent initialization has to be performed except for interstitial-particle-flux and network connections. Do you have some insights here?

[schmoelder]

Without having done any further debugging, here are some plots for the following example with just one column and no other units in the system (velocity = 0). Thinking about it, it's probably not an issue of consistent initialization any more but something else (probably also in the linear model). I will test more.

# Column
model.root.input.model.unit_000.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES'
model.root.input.model.unit_000.ncomp = 1

model.root.input.model.unit_000.col_length = 0.25
model.root.input.model.unit_000.cross_section_area = 1e-4
model.root.input.model.unit_000.col_porosity = 0.37
model.root.input.model.unit_000.par_porosity = 0.33
model.root.input.model.unit_000.par_radius = 4.5e-5

model.root.input.model.unit_000.col_dispersion = 2.0e-7
model.root.input.model.unit_000.film_diffusion = [0.0,]

With

model.root.input.model.unit_000.adsorption.is_kinetic = False
model.root.input.model.unit_000.adsorption_model = 'LINEAR'
model.root.input.model.unit_000.adsorption.lin_ka = [2]
model.root.input.model.unit_000.adsorption.lin_kd = [1]

we still get the following solution over time (at any position, here at column inlet)

And over space (at all time points, here at t=0):

Now, if we set

model.root.input.model.unit_000.adsorption.is_kinetic = True

we get the following over time (again at all positions, here at column inlet)

And over space at t=0:

And at later times, consequently:

[schmoelder]

My bad. There is no issue! 🙈
Here is a number of problems with my demonstration:

• We should consider particle liquid phase and not bulk phase
• We can consider bulk phase, but then, film diffusion should not be 0
• Initial conditions for particle liquid phase also need to be considered.

Then, we get the expected result:

Sorry for the confusion!

[sleweke-bayer]

So it works for kinetic and quasi-stationary mode or only for kinetic?

[schmoelder]

Only for kinetic. Same problem with rapid equilibrium as the CSTR.

(Both particle and bulk at c=1)