diff --git a/doc/modelling/binding/freundlich_ldf.rst b/doc/modelling/binding/freundlich_ldf.rst index e225f32bd..2b50cf28c 100644 --- a/doc/modelling/binding/freundlich_ldf.rst +++ b/doc/modelling/binding/freundlich_ldf.rst @@ -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:: @@ -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`. diff --git a/src/libcadet/model/binding/FreundlichLDFBinding.cpp b/src/libcadet/model/binding/FreundlichLDFBinding.cpp index 1b76304f8..910aab16e 100644 --- a/src/libcadet/model/binding/FreundlichLDFBinding.cpp +++ b/src/libcadet/model/binding/FreundlichLDFBinding.cpp @@ -122,7 +122,7 @@ namespace cadet const ParamType kkin = static_cast(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); @@ -130,7 +130,7 @@ namespace cadet } 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 @@ -163,7 +163,7 @@ 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); @@ -171,7 +171,7 @@ namespace cadet } 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; diff --git a/test/BindingModels.cpp b/test/BindingModels.cpp index d0aae3c63..de55e9f57 100644 --- a/test/BindingModels.cpp +++ b/test/BindingModels.cpp @@ -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)