Radial flow for chromatography column models

CMakeLists.txt

@@ -458,7 +458,7 @@ endif()
 target_compile_options(CADET::CompileOptions INTERFACE $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>,$<CXX_COMPILER_ID:GNU>>:
 		-Wall -pedantic-errors -Wextra -Wno-unused-parameter -Wno-unused-function> #-Wconversion -Wsign-conversion
 	$<$<CXX_COMPILER_ID:MSVC>:
-		/W4>
+		/W4 /wd4100>
 )
 
 add_library(CADET::AD INTERFACE IMPORTED)

doc/interface/input_group.rst

@@ -15,5 +15,6 @@ Input Group
     return_data
     sensitivities
     solver
+    parameter_dependencies
 
 

doc/interface/parameter_dependencies.rst

@@ -0,0 +1,79 @@
+.. _parameter_dependencies:
+
+Parameter Dependencies
+======================
+
+Some parameters depend on other parameters (parameter-parameter dependency) or the solution variables (parameter-state dependency).
+Parameter dependencies are defined in the unit operation scope.
+
+Parameter-Parameter Dependencies
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Group /input/model/unit_XXX
+---------------------------
+
+``COL_DISPERSION_DEP``
+
+   Parameter dependence of column dispersion on the interstitial velocity. Available for the LRM, LRMP and GRM units (with FV discretization only at the moment)
+
+   ================  =====================================  =============
+   **Type:** string  **Range:** :math:`\texttt{POWER_LAW}`  **Length:** 1
+   ================  =====================================  =============
+
+``FILM_DIFFUSION_DEP``
+
+   Parameter dependence of film diffusion on the interstitial velocity. Available for the LRMP unit (with FV discretization only at the moment)
+
+   ================  =====================================  =============
+   **Type:** string  **Range:** :math:`\texttt{POWER_LAW}`  **Length:** 1
+   ================  =====================================  =============
+
+
+**Correlations**
+""""""""""""""""
+
+Different types of parameter correlations are can be applied.
+The following correlations can be used for all parameter-parameter dependencies, but we specify the required input fields only for ``COL_DISPERSION_DEP``, for the sake of conciseness.
+
+**Power Law**
+
+.. math::
+
+    \begin{aligned}
+        p_{dep} &= p_{dep} \cdot b \ |p_{on}^x|
+    \end{aligned}
+
+Here, :math:`p_{dep}` is the dependent parameter and :math:`p_{on}` is the parameter it depends on.
+
+``COL_DISPERSION_DEP_BASE``
+
+   Base :math:`b` of the power law parameter dependence. Optional, defaults to :math:`1.0`
+
+   ================  =============================  =============
+   **Type:** double  **Range:** :math:`\mathbb{R}`  **Length:** 1
+   ================  =============================  =============
+
+``COL_DISPERSION_DEP_EXPONENT``
+
+   Exponent :math:`x` of the power law parameter dependence
+
+   ================  =============================  =============
+   **Type:** double  **Range:** :math:`\mathbb{R}`  **Length:** 1
+   ================  =============================  =============
+
+``COL_DISPERSION_DEP_ABS``
+
+   Specifies whether or not the absolute value should be computed. Optional, defaults to :math:`1`
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
+
+
+Parameter-State Dependencies
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Group /input/model/unit_XXX
+---------------------------
+
+Parameter-State Dependencies are not fully implemented yet.

doc/interface/unit_operations/index.rst

@@ -13,4 +13,5 @@ Unit Operations
     lumped_rate_model_without_pores
     2d_general_rate_model
     cstr
+    radial_flow_models
 

doc/interface/unit_operations/radial_flow_models.rst

@@ -0,0 +1,109 @@
+.. _radial_flow_models_config:
+
+Radial FLow Models
+==================
+
+Radial flow models are available for the LRM, LRMP and GRM.
+The configurations specified here complement the description for the respective model, i.e. please additionally refer to :ref:`lumped_rate_model_without_pores_config` or :ref:`lumped_rate_model_with_pores_config` or :ref:`general_rate_model_config`, respectively.
+If input variables are described in both files, then only the description provided here applies for radial flow models.
+
+The unit type input must be specified with the prefix :math:`\texttt{RADIAL_}` followed by the respective transport model name.
+In this document, we specify the unit type as the radial GRM, but the LRM and LRMP are also available and the changes to input variables are the same for all three models.
+
+Group /input/model/unit_XXX - UNIT_TYPE - RADIAL_GENERAL_RATE_MODEL
+-------------------------------------------------------------------
+
+For information on model equations, refer to :ref:`lumped_rate_model_without_pores` or :ref:`lumped_rate_model_with_pores` or :ref:`general_rate_model`, respectively.
+
+
+``UNIT_TYPE``
+
+   Specifies the type of unit operation model
+
+   ================  =====================================================  =============
+   **Type:** string  **Range:** :math:`\texttt{RADIAL_GENERAL_RATE_MODEL}`  **Length:** 1
+   ================  =====================================================  =============
+
+``COL_DISPERSION``
+
+   �Radial dispersion coefficient
+
+   **Unit:** :math:`\mathrm{m}_{\mathrm{IV}}^{2}\,\mathrm{s}^{-1}`
+
+   ================  =========================  =========================================================
+   **Type:** double  **Range:** :math:`\geq 0`  **Length:** see :math:`\texttt{COL_DISPERSION_MULTIPLEX}`
+   ================  =========================  =========================================================
+
+	In addition to the multiplex specification (e.g. component dependency, see :ref:`general_rate_model_model`), the dispersion coefficient for radial flow model usually depends on other parameters.
+	Parameter dependencies are described here :ref:`parameter_dependencies`.
+
+
+``COL_RADIUS_INNER``
+
+   Inner column radius 
+
+   **Unit:** :math:`\mathrm{m}`
+
+   ================  ======================  =============
+   **Type:** double  **Range:** :math:`> 0`  **Length:** 1
+   ================  ======================  =============
+
+``COL_RADIUS_OUTER``
+
+   Outer column radius 
+
+   **Unit:** :math:`\mathrm{m}`
+
+   ================  ======================  =============
+   **Type:** double  **Range:** :math:`> 0`  **Length:** 1
+   ================  ======================  =============
+
+``CROSS_SECTION_AREA``
+
+   Is not explicitly specified. Both the inner and outer cross section areas are implicitly given by the volumetric flow rates and either the velocity coefficient or column length
+
+``COL_LENGTH``
+
+   Column length/height (optional if :math:`\texttt{VELOCITY_COEFF}` is present, see Section :ref:`MUOPGRMflow`) 
+
+   **Unit:** :math:`\mathrm{m}`
+
+   ================  ======================  =============
+   **Type:** double  **Range:** :math:`> 0`  **Length:** 1
+   ================  ======================  =============
+
+``VELOCITY_COEFF``
+
+   Interstitial velocity coefficient of the mobile phase (optional :math:`\texttt{COL_LENGTH}` is present, see Section :ref:`MUOPGRMflow`).
+   This input replaces the ``VELOCITY`` field, which is used for axial flow models. The distinction is made to emphasize that radial flow models do not incorporate a global velocity but a variable velocity field that depends on the spatial position, for details see Section :ref:`MUOPGRMradialFlow`.
+
+   **Unit:** :math:`\mathrm{m}\,\mathrm{s}^{-1}`
+
+   ================  =============================  =======================================
+   **Type:** double  **Range:** :math:`\mathbb{R}`  **Length:** :math:`1 / \texttt{NSEC}`
+   ================  =============================  =======================================
+
+
+Group /input/model/unit_XXX/discretization - UNIT_TYPE - RADIAL_GENERAL_RATE_MODEL
+----------------------------------------------------------------------------------------
+
+``NCOL``
+
+   Number of radial column discretization points
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+Currently, there is only a first order FV spatial discretization available. Higher order spatial discretizations are planned for the future.
+Accordingly, the following specifications can be left out for radial flow models.
+
+``RECONSTRUCTION``
+
+   Type of reconstruction method for fluxes
+
+   ================  ================================  =============
+   **Type:** string  **Range:** :math:`\texttt{NONE}`  **Length:** 1
+   ================  ================================  =============
+
+Parameters specified under :ref:`flux_restruction_methods` can also be ignored.

doc/modelling/unit_operations/general_rate_model.rst

@@ -10,7 +10,7 @@ The main assumptions are:
 - The cross sections of the column are homogenous in terms of interstitial volume, fluid flow, and distribution of components.
   Thus, only one spatial coordinate in axial direction is needed and radial transport is neglected in the column bulk volume.
 
-- The bead radii :math:`r_{p}` are much smaller than the column radius :math:`r_c` and the column length :math:`L`.
+- The bead radii :math:`r_{p}` are much smaller than the column radius :math:`\mathrm{P}` and the column length :math:`L`.
   Therefore, the beads can be seen as continuously distributed inside the column (i.e., at each point there is interstitial and bead volume).
 
 .. _table_features:
@@ -276,12 +276,72 @@ The direction of flow inside the unit operation is governed by the sign of the i
 A positive sign results in (standard) forward flow, whereas a negative sign reverses the flow direction.
 Note that in case of reversed flow, the chromatogram is returned at the unit operation’s `INLET`, which may not be returned from simulation by default.
 
-The final behavior is controlled by the interplay of cross section area and interstitial velocity:
+The final behavior for axial flow models is controlled by the interplay of cross section area and interstitial velocity:
 
 - If cross section area :math:`A` is given and :math:`u` is not, :math:`u` is inferred from the volumetric flow rate.
 
 - If :math:`u` is given and :math:`A` is not, the volumetric flow rate is ignored and the provided interstitial velocity is used.
 
 - If both cross section area :math:`A` and interstitial velocity :math:`u` are given, the magnitude of the actual interstitial velocity :math:`u` is inferred from the volumetric flow rate and the flow direction is given by the sign of the provided :math:`u`.
 
+The final behavior for radial flow models is controlled by the interplay of column length/height and interstitial velocity coefficient:
+
+- If :math:`L` is given, the interstitial velocity field is inferred from the volumetric flow rate.
+
+- If :math:`u` is given and :math:`L` is not, the provided interstitial velocity coefficient is used to calculate the interstitial velocity field.
+
+
 For information on model parameters see :ref:`general_rate_model_config`.
+
+.. _MUOPGRMradialFlow:
+
+Radial flow GRM
+^^^^^^^^^^^^^^^
+
+The radial flow GRM describes transport of solute molecules through the interstitial column volume by radial convective flow, band broadening caused by radial dispersion, mass transfer resistance through a stagnant film around the beads, pore (and surface) diffusion in the porous beads :cite:`Ma1996,Schneider1968a,Miyabe2007`, and adsorption to the inner bead surfaces.
+
+The main assumptions are:
+
+- The shells of the column are homogenous in terms of interstitial volume, fluid flow, and distribution of components.
+  Thus, only one spatial coordinate in radial direction :math:`\rho` is needed and axial transport is neglected in the column bulk volume.
+
+- The bead radii :math:`r_{p}` are much smaller than the column radius :math:`\mathrm{P}-\mathrm{P}_c`, with :math:`\mathrm{P}` and :math:`\mathrm{P}_c` being the inner and outer column radius respectively, and the column length :math:`L`.
+  Therefore, the beads can be seen as continuously distributed inside the column (i.e., at each point there is interstitial and bead volume).
+
+- The fluids are incompressible, i.e. the velocity field :math:`\mathrm{V} \colon \mathbb{R}^3 \to \mathbb{R}^3` submits to :math:`\operatorname{div}\left( \mathrm{V} \right) \equiv 0`.
+  That is, the volumetric flow rate at the inner and outer column radius are the same.
+
+Consider a hollow (double walled) column with inner column diameter :math:`\mathrm{P}_c>0` and outer diameter :math:`\mathrm{P}>\mathrm{P}_c`, filled with spherical beads of (possibly) multiple types with radius :math:`r_{p,j} \ll L` (see :numref:`ModelGRMColumn`), where :math:`j` is the particle type index. The mass balance in the interstitial column volume is described by
+
+.. math::
+    :label: ModelRadialColumn
+
+    \begin{aligned}
+        \frac{\partial c^l_i}{\partial t} = -\frac{u}{\rho} \frac{\partial c^l_i}{\partial \rho} + D_{\text{rad},i} \frac{1}{\rho} \frac{\partial}{\partial \rho} \left(\rho \frac{\partial c^l_i}{\partial \rho} \right) &- \frac{1}{\beta_c} \sum_j d_j \frac{3}{r_{p,j}} k_{f,j,i} \left[ c^l_i - c^p_{j,i}(\cdot, \cdot, r_{p,j}) \right] \\
+        &+ f_{\text{react},i}^l\left(c^l\right). 
+    \end{aligned}
+
+Here, :math:`c^l_i\colon \left[0, T_{\text{end}}\right] \times [\mathrm{P}_c, \mathrm{P}] \rightarrow \mathbb{R}^{\geq 0}` denotes the concentration in the interstitial column volume, :math:`c^p_{j,i}\colon \left[0, T_{\text{end}}\right] \times [P_c, P] \times [r_{c,j}, r_{p,j}] \rightarrow \mathbb{R}^{\geq 0}` the liquid phase concentration in the beads, :math:`k_{f,j,i}` the film diffusion coefficient, :math:`D_{\text{rad},i}` the dispersion coefficient, :math:`u` the interstitial velocity, :math:`d_j` the volume fraction of particle type :math:`j`, and :math:`\beta_c = \varepsilon_c / (1 - \varepsilon_c)` the column phase ratio, where :math:`\varepsilon_c` is the column porosity (ratio of interstitial volume to total column volume).
+If reactions are considered, the term :math:`f_{\text{react},i}^l\left(c^l\right)` represents the net change of concentration :math:`c_i` due to reactions involving component :math:`i`.
+
+Danckwerts boundary conditions :cite:`Danckwerts1953` are applied to inlet and outlet of the column:
+
+.. math::
+    :label: BCOutlet
+
+    \begin{aligned}
+        u c_{\text{in},i}(t) &= u c^l_i(t,0) - D_{\text{rad},i} \frac{\partial c^l_i}{\partial \rho}(t, 0) & \forall t > 0,
+    \end{aligned}
+
+.. math::
+    :label: BCInlet
+
+    \begin{aligned}
+        \frac{\partial c^l_i}{\partial \rho}(t, \mathrm{P}) &= 0 & \forall t > 0. 
+    \end{aligned}
+
+Note that the outlet boundary condition Eq. :eq:`BCOutlet` is also known as “do nothing” or natural outflow condition.
+
+The complementing mass transport and binding equations for the liquid and solid phases of the porous beads are described by the same equations as for the axial GRM.
+
+For information on model parameters see :ref:`radial_flow_models_config` in addition to :ref:`general_rate_model_config`.

doc/modelling/unit_operations/index.rst

@@ -50,6 +50,7 @@ is given in :numref:`table_features_unit_operations`.
 
 
 Moreover, the pseudo unit operations :ref:`inlet_model`, and :ref:`outlet_model` act as sources and sinks for the system. 
+We further note that radial flow model variants are available for the LRM, LRMP and GRM.
 
 
 .. toctree::

doc/modelling/unit_operations/lumped_rate_model_with_pores.rst

@@ -14,6 +14,7 @@ Hence, the model equations are given by
     \end{aligned}
 
 .. math::
+    :label: ModelParticleLRMP
 
     \begin{aligned}
         \frac{\partial c^p_{j,i}}{\partial t} + \frac{1 - \varepsilon_{p,j}}{F_{\text{acc},j,i} \varepsilon_{p,j}} \frac{\partial}{\partial t} \sum_{m_{j,i}} c^s_{j,i,m_{j,i}} &= \frac{3}{F_{\text{acc},j,i} \varepsilon_{p,j} r_{p,j}}k_{f,j,i}\left[ c^l_i - c^p_{j,i} \right] \\
@@ -52,3 +53,40 @@ This model can also be used to simulate :ref:`MUOPGRMSizeExclusion`.
 For the specification of flow rate and direction, the same holds as for the general rate model (see Section :ref:`MUOPGRMflow`).
 
 For information on model parameters see :ref:`lumped_rate_model_with_pores_config`.
+
+Radial flow LRMP
+^^^^^^^^^^^^^^^^
+
+The radial flow LRMP describes transport of solute molecules through the interstitial column volume by radial convective flow, band broadening caused by radial dispersion, mass transfer resistance through a stagnant film around the beads, and adsorption to the inner bead surfaces.
+
+The main assumptions are:
+
+- The shells of the column are homogenous in terms of interstitial volume, fluid flow, and distribution of components.
+  Thus, only one spatial coordinate in radial direction :math:`\rho` is needed and axial transport is neglected in the column bulk volume.
+
+- The bead radii :math:`r_{p}` are much smaller than the column radius :math:`\mathrm{P}-\mathrm{P}_c`, with :math:`\mathrm{P}` and :math:`\mathrm{P}_c` being the inner and outer column radius respectively, and the column length :math:`L`.
+  Therefore, the beads can be seen as continuously distributed inside the column (i.e., at each point there is interstitial and bead volume).
+
+- The fluids are incompressible, i.e. the velocity field :math:`\mathrm{V} \colon \mathbb{R}^3 \to \mathbb{R}^3` submits to :math:`\operatorname{div}\left( \mathrm{V} \right) \equiv 0`.
+  That is, the volumetric flow rate at the inner and outer column radius are the same.
+
+Consider a hollow (double walled) column with inner column diameter :math:`\mathrm{P}_c>0` and outer diameter :math:`\mathrm{P}>\mathrm{P}_c`, filled with spherical beads of (possibly) multiple types with radius :math:`r_{p,j} \ll L` (see :numref:`ModelGRMColumn`), where :math:`j` is the particle type index. The mass balance in the interstitial column volume is described by
+
+.. math::
+
+    \begin{aligned}
+        \frac{\partial c^l_i}{\partial t} &= -\frac{u}{\rho} \frac{\partial c^l_i}{\partial \rho} + D_{\text{rad},i} \frac{1}{\rho} \frac{\partial}{\partial \rho} \left( \rho \frac{\partial c^l_i}{\partial \rho} \right) - \frac{1}{\beta_c} \sum_{j} d_j \frac{3}{r_{p,j}} k_{f,j,i}\left[ c^l_i - c^p_{j,i} \right] + f_{\text{react},i}^l\left(c^l\right),
+    \end{aligned}
+
+The equations are complemented by Eq. :ref:`ModelParticleLRMP` and the Danckwerts boundary conditions :cite:`Danckwerts1953`
+
+.. math::
+
+    \begin{aligned}
+        u c_{\text{in},i}(t) &= u c^l_i(t,0) - D_{\text{rad},i} \frac{\partial c^l_i}{\partial \rho}(t, 0) & \forall t > 0,\\
+        \frac{\partial c^l_i}{\partial \rho}(t, \mathrm{P}) &= 0 & \forall t > 0.
+    \end{aligned}
+
+The complementing binding equations are described by the same equations as for the axial LRMP.
+
+For information on model parameters see :ref:`radial_flow_models_config` in addition to :ref:`lumped_rate_model_with_pores_config`.