Skip to content

Commit 2040e7b

Browse files
ronald-jaepeljbreue16
ronald-jaepel
authored and
jbreue16
committed
Update Freundlich isotherm linearization and tests
Extend linearization range into negative numbers to improve stability Extend tests Update documentation
1 parent fe8146d commit 2040e7b

File tree

3 files changed

+42
-42
lines changed

3 files changed

+42
-42
lines changed

‎doc/modelling/binding/freundlich_ldf.rst

+4-4
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,10 @@ This variant of the model is based on the linear driving force approximation (se
1313
1414
No interaction between the components is considered when the model has multiple components.
1515
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`.
16-
To address this issue an approximation of isotherm is considered near the origin.
17-
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`.
18-
The form of approximation and its derivative is given below for :math:`c_p < \varepsilon` and :math:`n>1`:
16+
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.
17+
To address these issues an approximation of the isotherm is considered below a threshold concentration :math:`c_p < \varepsilon`.
18+
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`.
19+
The form of approximation and its derivative is given below for :math:`c_p < \varepsilon`:
1920

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

‎src/libcadet/model/binding/FreundlichLDFBinding.cpp

+4-4
Original file line numberDiff line numberDiff line change
@@ -122,15 +122,15 @@ namespace cadet
122122
const ParamType kkin = static_cast<ParamType>(p->kkin[i]);
123123

124124
// Residual
125-
if ((n_param > 1) && (abs(yCp[i]) < _threshold))
125+
if (((n_param > 1) && (yCp[i] < _threshold)) || yCp[i] < 0)
126126
{
127127
const ParamType alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param);
128128
const ParamType alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param);
129129
res[bndIdx] = kkin * (y[bndIdx] - yCp[i] * (alpha_1 + alpha_2 * yCp[i]));
130130
}
131131
else
132132
{
133-
res[bndIdx] = kkin * (y[bndIdx] - kF * pow(abs(yCp[i]), 1.0 / n_param));
133+
res[bndIdx] = kkin * (y[bndIdx] - kF * pow(yCp[i], 1.0 / n_param));
134134
}
135135

136136
// Next bound component
@@ -163,15 +163,15 @@ namespace cadet
163163
jac[0] = kkin;
164164

165165
// dres / dc_{p,i}
166-
if ((n_param > 1) && (abs(yCp[i]) < _threshold))
166+
if (((n_param > 1) && (yCp[i] < _threshold)) || yCp[i] < 0)
167167
{
168168
double const alpha_1 = ((2.0 * n_param - 1.0) / n_param) * kF * pow(_threshold, (1.0 - n_param) / n_param);
169169
double const alpha_2 = ((1.0 - n_param) / n_param) * kF * pow(_threshold, (1.0 - 2.0 * n_param) / n_param);
170170
jac[i - bndIdx - offsetCp] = -kkin * (alpha_1 + 2.0 * alpha_2 * yCp[i]);
171171
}
172172
else
173173
{
174-
jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(abs(yCp[i]), (1.0 - n_param) / n_param);
174+
jac[i - bndIdx - offsetCp] = -(1.0 / n_param) * kkin * kF * pow(yCp[i], (1.0 - n_param) / n_param);
175175
}
176176

177177
++bndIdx;

‎test/BindingModels.cpp

+34-34
Original file line numberDiff line numberDiff line change
@@ -42,40 +42,40 @@ CADET_BINDINGTEST("LINEAR", "EXT_LINEAR", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0),
4242
)json", \
4343
1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING)
4444

45-
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), \
46-
R"json( "FLDF_KKIN": [1.0, 2.0],
47-
"FLDF_KF": [0.1, 0.2],
48-
"FLDF_N": [0.5, 1.2]
49-
)json", \
50-
R"json( "FLDF_KKIN": [1.0, 0.5, 2.0],
51-
"FLDF_KF": [0.1, 0.3, 0.2],
52-
"FLDF_N": [0.5, 2.2, 1.2]
53-
)json", \
54-
R"json( "EXT_FLDF_KKIN": [0.0, 0.0],
55-
"EXT_FLDF_KKIN_T": [1.0, 2.0],
56-
"EXT_FLDF_KKIN_TT": [0.0, 0.0],
57-
"EXT_FLDF_KKIN_TTT": [0.0, 0.0],
58-
"EXT_FLDF_KF": [0.0, 0.0],
59-
"EXT_FLDF_KF_T": [0.1, 0.2],
60-
"EXT_FLDF_KF_TT": [0.0, 0.0],
61-
"EXT_FLDF_KF_TTT": [0.0, 0.0],
62-
"EXT_FLDF_N": [0.0, 0.0],
63-
"EXT_FLDF_N_T": [0.5, 1.2],
64-
"EXT_FLDF_N_TT": [0.0, 0.0],
65-
"EXT_FLDF_N_TTT": [0.0, 0.0]
66-
)json", \
67-
R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0],
68-
"EXT_FLDF_KKIN_T": [1.0, 0.5, 2.0],
69-
"EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0],
70-
"EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0],
71-
"EXT_FLDF_KF": [0.0, 0.0, 0.0],
72-
"EXT_FLDF_KF_T": [0.1, 0.3, 0.2],
73-
"EXT_FLDF_KF_TT": [0.0, 0.0, 0.0],
74-
"EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0],
75-
"EXT_FLDF_N": [0.0, 0.0, 0.0],
76-
"EXT_FLDF_N_T": [0.5, 2.0, 1.2],
77-
"EXT_FLDF_N_TT": [0.0, 0.0, 0.0],
78-
"EXT_FLDF_N_TTT": [0.0, 0.0, 0.0]
45+
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), \
46+
R"json( "FLDF_KKIN": [1.0, 1.0, 1.2, 2.0],
47+
"FLDF_KF": [0.1, 0.3, 0.3, 0.2],
48+
"FLDF_N": [0.5, 1.0, 1.2, 0.8]
49+
)json", \
50+
R"json( "FLDF_KKIN": [1.0, 0.5, 1.0, 1.2, 2.0],
51+
"FLDF_KF": [0.1, 0.3, 0.3, 0.3, 0.2],
52+
"FLDF_N": [0.5, 2.2, 1.0, 1.2, 0.8]
53+
)json", \
54+
R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0, 0.0],
55+
"EXT_FLDF_KKIN_T": [1.0, 1.0, 1.2, 2.0],
56+
"EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0, 0.0],
57+
"EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0, 0.0],
58+
"EXT_FLDF_KF": [0.0, 0.0, 0.0, 0.0],
59+
"EXT_FLDF_KF_T": [0.1, 0.3, 0.3, 0.2],
60+
"EXT_FLDF_KF_TT": [0.0, 0.0, 0.0, 0.0],
61+
"EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0, 0.0],
62+
"EXT_FLDF_N": [0.0, 0.0, 0.0, 0.0],
63+
"EXT_FLDF_N_T": [0.5, 1.0, 1.2, 0.8],
64+
"EXT_FLDF_N_TT": [0.0, 0.0, 0.0, 0.0],
65+
"EXT_FLDF_N_TTT": [0.0, 0.0, 0.0, 0.0]
66+
)json", \
67+
R"json( "EXT_FLDF_KKIN": [0.0, 0.0, 0.0, 0.0, 0.0],
68+
"EXT_FLDF_KKIN_T": [1.0, 0.5, 1.0, 1.2, 2.0],
69+
"EXT_FLDF_KKIN_TT": [0.0, 0.0, 0.0, 0.0, 0.0],
70+
"EXT_FLDF_KKIN_TTT": [0.0, 0.0, 0.0, 0.0, 0.0],
71+
"EXT_FLDF_KF": [0.0, 0.0, 0.0, 0.0, 0.0],
72+
"EXT_FLDF_KF_T": [0.1, 0.3, 0.3, 0.3, 0.2],
73+
"EXT_FLDF_KF_TT": [0.0, 0.0, 0.0, 0.0, 0.0],
74+
"EXT_FLDF_KF_TTT": [0.0, 0.0, 0.0, 0.0, 0.0],
75+
"EXT_FLDF_N": [0.0, 0.0, 0.0, 0.0, 0.0],
76+
"EXT_FLDF_N_T": [0.5, 2.0, 1.0, 1.2, 0.8],
77+
"EXT_FLDF_N_TT": [0.0, 0.0, 0.0, 0.0, 0.0],
78+
"EXT_FLDF_N_TTT": [0.0, 0.0, 0.0, 0.0, 0.0]
7979
)json", \
8080
1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING)
8181

0 commit comments

Comments
 (0)