Skip to content

Commit

Permalink
Update Freundlich isotherm linearization and tests
Browse files Browse the repository at this point in the history 
Extend linearization range into negative numbers to improve stability
Extend tests
Update documentation
  • Loading branch information
@ronald-jaepel
ronald-jaepel committed Jan 15, 2025
1 parent 3e8c128 commit 86d5c42
Show file tree
Hide file tree
Showing 3 changed files with 42 additions and 42 deletions.
8 changes: 4 additions & 4 deletions doc/modelling/binding/freundlich_ldf.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,10 @@ This variant of the model is based on the linear driving force approximation (se
No interaction between the components is considered when the model has multiple components.
One of the limitation of this isotherm is the first order Jacobian :math:`\left(\frac{dq^*}{dc_p}\right)` tends to infinity as :math:`c_{p} \rightarrow 0` for :math:`n>1`.
To address this issue an approximation of isotherm is considered near the origin.
This approximation matches the isotherm in such a way that :math:`q=0` at :math:`c_p=0` and also matches the first derivative of the isotherm at :math:`c_p = \varepsilon`, where :math:`\varepsilon` is a very small number, for example :math:`1e-14`.
The form of approximation and its derivative is given below for :math:`c_p < \varepsilon` and :math:`n>1`:
Additionally, the isotherm is undefined for :math:`c_{p} < 0` if :math:`\frac{1}{n_i}` can be expressed as :math:`\frac{p}{q}` with :math:`p,q \in \mathbb{N}` where :math:`q` is an even number.
To address these issues an approximation of the isotherm is considered below a threshold concentration :math:`c_p < \varepsilon`.
This approximation matches the isotherm in such a way that :math:`q=0` at :math:`c_p=0` and also matches the value and the first derivative of the isotherm at :math:`c_p = \varepsilon`, where :math:`\varepsilon` is a very small number, for example :math:`1e-14`.
The form of approximation and its derivative is given below for :math:`c_p < \varepsilon`:

.. math::
Expand All @@ -39,6 +40,5 @@ where :math:`\alpha_0=0` and :math:`\alpha_1` and :math:`\alpha_2` are determine
\end{aligned}
This approximation can be used for any pore phase concentration :math:`c_p < \varepsilon` given :math:`n>1`.
For the case, when :math:`n \le 1` no special treatment near the origin is required.
For more information on model parameters required to define in CADET file format, see :ref:`freundlich_ldf_config`.

8 changes: 4 additions & 4 deletions src/libcadet/model/binding/FreundlichLDFBinding.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,15 +122,15 @@ namespace cadet
const ParamType kkin = static_cast<ParamType>(p->kkin[i]);

// Residual
if ((n_param > 1) && (abs(yCp[i]) < _threshold))
if (((n_param > 1) && (yCp[i] < _threshold)) || yCp[i] < 0)
{
const ParamType alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param);
const ParamType alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param);
res[bndIdx] = kkin * (y[bndIdx] - yCp[i] * (alpha_1 + alpha_2 * yCp[i]));
}
else
{
res[bndIdx] = kkin * (y[bndIdx] - kF * pow(abs(yCp[i]), 1.0 / n_param));
res[bndIdx] = kkin * (y[bndIdx] - kF * pow(yCp[i], 1.0 / n_param));
}

// Next bound component
Expand Down Expand Up @@ -163,15 +163,15 @@ namespace cadet
jac[0] = kkin;

// dres / dc_{p,i}
if ((n_param > 1) && (abs(yCp[i]) < _threshold))
if (((n_param > 1) && (yCp[i] < _threshold)) || yCp[i] < 0)
{
double const alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param);
double const alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param);
jac[i - bndIdx - offsetCp] = -kkin * (alpha_1 + 2.0 * alpha_2 * yCp[i]);
}
else
{
jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(abs(yCp[i]), (1.0 - n_param) / n_param);
jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(yCp[i], (1.0 - n_param) / n_param);
}

++bndIdx;
Expand Down
68 changes: 34 additions & 34 deletions test/BindingModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,40 +42,40 @@ CADET_BINDINGTEST("LINEAR", "EXT_LINEAR", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0),
)json", \
1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING)

CADET_BINDINGTEST("FREUNDLICH_LDF", "EXT_FREUNDLICH_LDF", (1, 1), (1, 0, 1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \
R"json( "FLDF_KKIN": [1.0, 2.0],
"FLDF_KF": [0.1, 0.2],
"FLDF_N": [0.5, 1.2]
)json", \
R"json( "FLDF_KKIN": [1.0, 0.5, 2.0],
"FLDF_KF": [0.1, 0.3, 0.2],
"FLDF_N": [0.5, 2.2, 1.2]
)json", \
R"json( "EXT_FLDF_KKIN": [0.0, 0.0],
"EXT_FLDF_KKIN_T": [1.0, 2.0],
"EXT_FLDF_KKIN_TT": [0.0, 0.0],
"EXT_FLDF_KKIN_TTT": [0.0, 0.0],
"EXT_FLDF_KF": [0.0, 0.0],
"EXT_FLDF_KF_T": [0.1, 0.2],
"EXT_FLDF_KF_TT": [0.0, 0.0],
"EXT_FLDF_KF_TTT": [0.0, 0.0],
"EXT_FLDF_N": [0.0, 0.0],
"EXT_FLDF_N_T": [0.5, 1.2],
"EXT_FLDF_N_TT": [0.0, 0.0],
"EXT_FLDF_N_TTT": [0.0, 0.0]
)json", \
R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0],
"EXT_FLDF_KKIN_T": [1.0, 0.5, 2.0],
"EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0],
"EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0],
"EXT_FLDF_KF": [0.0, 0.0, 0.0],
"EXT_FLDF_KF_T": [0.1, 0.3, 0.2],
"EXT_FLDF_KF_TT": [0.0, 0.0, 0.0],
"EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0],
"EXT_FLDF_N": [0.0, 0.0, 0.0],
"EXT_FLDF_N_T": [0.5, 2.0, 1.2],
"EXT_FLDF_N_TT": [0.0, 0.0, 0.0],
"EXT_FLDF_N_TTT": [0.0, 0.0, 0.0]
CADET_BINDINGTEST("FREUNDLICH_LDF", "EXT_FREUNDLICH_LDF", (1, 1, 1, 1), (1, 0, 1, 1, 1), (1.0, 0.0, 2.0, -1e-4, 0.0, 0.0, 0.0, 0.0), (1.0, 3.0, 0.0, 2.0, -1e-4, 0.0, 0.0, 0.0, 0.0), \
R"json( "FLDF_KKIN": [1.0, 1.0, 1.2, 2.0],
"FLDF_KF": [0.1, 0.3, 0.3, 0.2],
"FLDF_N": [0.5, 1.0, 1.2, 0.8]
)json", \
R"json( "FLDF_KKIN": [1.0, 0.5, 1.0, 1.2, 2.0],
"FLDF_KF": [0.1, 0.3, 0.3, 0.3, 0.2],
"FLDF_N": [0.5, 2.2, 1.0, 1.2, 0.8]
)json", \
R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KKIN_T": [1.0, 1.0, 1.2, 2.0],
"EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KF": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KF_T": [0.1, 0.3, 0.3, 0.2],
"EXT_FLDF_KF_TT": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_N": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_N_T": [0.5, 1.0, 1.2, 0.8],
"EXT_FLDF_N_TT": [0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_N_TTT": [0.0, 0.0, 0.0, 0.0]
)json", \
R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KKIN_T": [1.0, 0.5, 1.0, 1.2, 2.0],
"EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KF": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KF_T": [0.1, 0.3, 0.3, 0.3, 0.2],
"EXT_FLDF_KF_TT": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_N": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_N_T": [0.5, 2.0, 1.0, 1.2, 0.8],
"EXT_FLDF_N_TT": [0.0, 0.0, 0.0, 0.0, 0.0],
"EXT_FLDF_N_TTT": [0.0, 0.0, 0.0, 0.0, 0.0]
)json", \
1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING)

Expand Down

0 comments on commit 86d5c42

Please sign in to comment.