Merge pull request idaholab#21230 from lindsayad/friction-doc

Change the way we do Darcy/Forchheimer

modules/navier_stokes/doc/content/bib/navier_stokes.bib

@@ -336,3 +336,44 @@ @phdthesis{jasak1996error
   school={Imperial College London (University of London)}
 }
 
+@InProceedings{suikkanen,
+      author = {H. Suikkanen and V. Rintala and R. Kyrki-Rajamaki},
+       title = {{Development of a Coupled Multi-Physics Code System for Pebble Bed Reactor Core Modeling}},
+  booktitle  = {{Proceedings of the HTR 2014}},
+  year       = 2014
+}
+
+@InProceedings{y_li,
+      author = {Y. Li and W. Ji},
+       title = {{Thermal Analysis of Pebble-Bed Reactors Based on a Tightly Coupled Mechanical-Thermal Model}},
+  booktitle  = {{Proceedings of NURETH}},
+  year       = 2016
+}
+
+@Article{gunn1987_htc,
+   author = {D.J. Gunn},
+    title = {{Transfer of Heat or Mass to Particles in Fixed and Fluidised Beds}},
+  journal = {International Journal of Heat and Mass Transfer},
+     year = 1978,
+   volume = 21,
+    pages = {467--476},
+      DOI = {10.1016/0017-9310(78)90080-7}
+}
+
+@Article{littman,
+   author = {H. Littman and R.G. Barile and A.H. Pulsifer},
+    title = {{Gas-Particle Heat Transfer Coefficients in Packed Beds at Low {Reynolds} Numbers}},
+  journal = {I \& EC Technology},
+     year = 1968,
+   volume = 7,
+    pages = {554--561},
+      DOI = {10.1021/i160028a005}
+}
+
+@InProceedings{becker,
+   author = {S. Becker and E. Laurien},
+    title = {{Three-Dimensional Numerical Simulation of Flow and Heat Transport in High-Temperature Nuclear Reactors}},
+  booktitle  = {{High Performance Computing in Science and Engineering}},
+  year       = 2002,
+        DOI = {10.1016/S0029-5493(03)00011-6}
+}

modules/navier_stokes/doc/content/modules/navier_stokes/pinsfv.md

@@ -1,20 +1,24 @@
 # Finite Volume Incompressible Porous media Navier Stokes
 
+!alert note
+The weakly compressible formulation is also implemented for porous flow finite volume.
+
 ## Equations
 
 This module implements the porous media Navier Stokes equations. They are expressed in terms of the superficial
-velocity $\vec{v}_d = \epsilon \vec{V}$ where $\epsilon$ is the porosity and $\vec{V}$ the interstitial velocity. The
+velocity $\vec{v}_D = \epsilon \vec{v}$ where $\epsilon$ is the porosity and $\vec{v}$ the interstitial velocity. The
 superficial velocity is also known as the extrinsic or Darcy velocity. The other non-linear variables used are
 pressure and temperature. This is known as the primitive superficial set of variables.
 
 Mass equation:
 \begin{equation}
-\nabla \cdot \rho \vec{v}_d = 0
+\nabla \cdot \rho \vec{v}_D = 0
 \end{equation}
 
 Momentum equation, with friction and gravity force as example forces:
 \begin{equation}
 \dfrac{\partial \rho \mathbf{v}_D}{\partial t} + \nabla \cdot (\dfrac{\rho}{\epsilon} \mathbf{v}_D \otimes \mathbf{v}_D) = \nabla \cdot (\mu \nabla \dfrac{\mathbf{v}_D}{\epsilon}) - \epsilon \nabla p + \epsilon (\mathbf{F}_g + \mathbf{F}_f)
+\label{eq:pinsfv_mom}
 \end{equation}
 
 Fluid phase energy equation, with a convective heat transfer term:
@@ -27,7 +31,7 @@ Solid phase energy equation, with convective heat transfer and an energy source
 \dfrac{\partial (1-\epsilon) \rho c_{ps} T_s}{\partial t} = \nabla \cdot (\kappa_s \nabla T_s) + \alpha (T_f - T_s) + (1-\epsilon) \dot{Q}
 \end{equation}
 
-where $\rho$ is the fluid density, $\mu$ the viscosity, $c_p$ the specific heat capacity $\alpha$ the convective heat transfer coefficient.
+where $\rho$ is the fluid density, $\mu$ the dynamic viscosity, $c_p$ the specific heat capacity $\alpha$ the convective heat transfer coefficient. The effective diffusivities $\kappa$ include the effect of porosity, which is often approximated as $\epsilon k$ with $k$ the thermal diffusivity.
 
 ## Implementation for a straight channel
 

modules/navier_stokes/doc/content/source/functormaterials/FunctorErgunDragCoefficients.md

@@ -0,0 +1,33 @@
+# FunctorErgunDragCoefficients
+
+!syntax description /Materials/FunctorErgunDragCoefficients
+
+This class implements the Darcy and Forchheimer coefficients for the Ergun drag
+force model which is discussed in detail in
+[PINSFVMomentumFriction.md#friction_example]. We also give details on the
+definition of Darcy and Forchheimer coefficients there.
+
+To summarize, this class implements the Darcy coefficient as
+
+\begin{equation}
+\frac{150}{D_h^2}
+\end{equation}
+
+where $D_h$ is the hydraulic diameter defined as $\frac{\epsilon d_p}{1 - \epsilon}$
+where $d_p$ is the diameter of the pebbles in the bed. The Forchheimer
+coefficient is given as
+
+\begin{equation}
+\frac{2 \cdot 1.75}{D_h}
+\end{equation}
+
+where we have made the $2(1.75)$ multiplication explicit, instead of writing
+$3.5$, to make the 1.75 factor from the
+[Ergun wikipedia page](https://en.wikipedia.org/wiki/Ergun_equation) more
+recognizable.
+
+!syntax parameters /Materials/FunctorErgunDragCoefficients
+
+!syntax inputs /Materials/FunctorErgunDragCoefficients
+
+!syntax children /Materials/FunctorErgunDragCoefficients

modules/navier_stokes/doc/content/source/functormaterials/FunctorKappaFluid.md

@@ -0,0 +1,23 @@
+# FunctorKappaFluid
+
+!syntax description /Materials/FunctorKappaFluid
+
+## Description
+
+Most macroscale models neglect thermal dispersion [!cite](suikkanen,y_li), in which case $\kappa_f$ is given as
+
+\begin{equation}
+\label{eq-KappaFluid}
+\kappa_f=\epsilon k_f\ .
+\end{equation}
+
+Neglecting thermal dispersion is expected to be a reasonable approximation for high Reynolds numbers [!cite](gunn1987_htc,littman),
+but for low Reynolds numbers more sophisticated models should be used [!cite](becker).
+Because thermal dispersion acts to increase the diffusive effects, neglecting thermal dispersion is
+(thermally) conservative in the sense that peak temperatures are usually higher [!cite](becker).
+
+!syntax parameters /Materials/FunctorKappaFluid
+
+!syntax inputs /Materials/FunctorKappaFluid
+
+!syntax children /Materials/FunctorKappaFluid

modules/navier_stokes/doc/content/source/fvkernels/PINSFVMomentumFriction.md

@@ -2,19 +2,115 @@
 
 This kernel adds the friction term to the porous media Navier Stokes momentum
 equations. This kernel must be used with the canonical PINSFV variable set,
-e.g. pressure and superficial velocity.  A variety of models are supported for
-the friction force:
+e.g. pressure and superficial velocity, and supports Darcy and
+Forchheimer friction models:
 
 Darcy drag model
 \begin{equation}
-F_i = - \dfrac{f_i}{\epsilon} \rho v_d
+\label{darcy}
+\epsilon F_i = - f_i \mu \epsilon \frac{v_{D,i}}{\epsilon} = -f_i \mu v_{D,i}
 \end{equation}
 Forchheimer drag model
 \begin{equation}
-F_i = - \dfrac{f_i}{\epsilon} \rho v_d
+\label{forchheimer}
+\epsilon F_i = - f_i \frac{\rho}{2} \epsilon \frac{v_{D,i}}{\epsilon}\frac{|\vec{v_D}|}{\epsilon} = - f_i \frac{\rho}{2} v_{D,i} \frac{|\vec{v_D}|}{\epsilon}
 \end{equation}
-where $F_i$ is the i-th component of the friction force, f the friction factor, which may be anisotropic,
-$\epsilon$ the porosity and $\rho$ the fluid density and $v_d$ the fluid superficial velocity.
+
+where $F_i$ is the i-th component of the friction force (denoted by
+$\mathbf{F_f}$ in [!eqref](pinsfv.md#eq:pinsfv_mom)), $f_i$ the friction factor,
+which may be anisotropic, $\epsilon$ the porosity, $\mu$ the fluid dynamic
+viscosity, $\rho$ the fluid density, and $v_{D,i}$ the i-th component of the
+fluid superficial velocity. We have used a negative sign to match the notation
+used in [!eqref](pinsfv.md#eq:pinsfv_mom) where the friction force is on the
+right-hand-side of the equation. When moved to the left-hand side, which is done
+when setting up a Newton scheme, the term becomes positive which is what is
+shown in the source code itself.  Darcy and Forchheimer terms represent
+fundamentally different friction effects. Darcy is meant to represent viscous
+effects and as shown in [darcy] has a linear dependence on the fluid
+velocity. Meanwhile, Forchheimer is meant to represent inertial effects and as
+shown in [forchheimer] has a quadratic dependence on velocity.
+
+## Computation of friction factors and pre-factors id=friction_example
+
+To outline how friction factors for Darcy and Forchheimer may be calculated,
+let's consider a specific example. We'll draw from the Ergun equation, which is
+outlined [here](https://en.wikipedia.org/wiki/Ergun_equation). Let's consider
+the form:
+
+\begin{equation}
+\Delta p = \frac{150\mu L}{d_p^2} \frac{(1-\epsilon)^2}{\epsilon^3} v_D + \frac{1.75 L \rho}{d_p} \frac{(1-\epsilon)}{\epsilon^3} |v_D| v_D
+\end{equation}
+
+where $L$ is the bed length, $\mu$ is the fluid dynamic viscosity and $d_p$ is
+representative of the diameter of the pebbles in the pebble bed. We can divide
+the equation through by $L$, recognize that $\Delta p$ denotes $p_0 - p_L$ such
+that $\Delta p/L \approx -\nabla p$, multiply the equation through by
+$-\epsilon$, move all terms to the left-hand-side, and do
+some term manipulation in order to yield:
+
+\begin{equation}
+\epsilon \nabla p + 150\mu\epsilon\frac{(1-\epsilon)^2}{\epsilon^2 d_p^2} \frac{\vec{v_D}}{\epsilon} +
+1.75\epsilon\frac{1-\epsilon}{\epsilon d_p} \rho \frac{|\vec{v}_D|}{\epsilon} \frac{\vec{v_D}}{\epsilon} = 0
+\end{equation}
+
+If we define the hydraulic diameter as $D_h = \frac{\epsilon d_p}{1 - \epsilon}$,
+then the above equation can be rewritten as:
+
+\begin{equation}
+\epsilon \nabla p + \frac{150\mu}{D_h^2} \vec{v_D} +
+\frac{1.75}{D_h} \rho \vec{v_D} \frac{|\vec{v}_D|}{\epsilon} = 0
+\end{equation}
+
+Let's introduce the interstitial fluid velocity $v = \vec{v_D} / \epsilon$ to rewrite
+the above equation as:
+
+\begin{equation}
+\epsilon \nabla p + \epsilon\frac{150\mu}{D_h^2} \vec{v} +
+\epsilon\frac{1.75}{D_h} \rho \vec{v} |\vec{v}| = 0
+\end{equation}
+
+Then dividing through by $\epsilon$:
+
+\begin{equation}
+\label{derived_ergun}
+\nabla p + \frac{150\mu}{D_h^2} \vec{v} +
+\frac{1.75}{D_h} \rho \vec{v} |\vec{v}| = 0
+\end{equation}
+
+We are now very close the forms for Darcy and Forchheimer espoused by
+[Holzmann](https://holzmann-cfd.com/community/blog-and-tools/darcy-forchheimer)
+and
+[SimScale](https://www.simscale.com/knowledge-base/predict-darcy-and-forchheimer-coefficients-for-perforated-plates-using-analytical-approach/)
+which is:
+
+\begin{equation}
+\label{holzmann}
+\nabla p + \mu D \vec{v} + \frac{\rho}{2} F |\vec{v}| \vec{v}
+\end{equation}
+
+Looking at [holzmann] we can rearrange [derived_ergun]:
+
+\begin{equation}
+\nabla p + \mu \frac{150}{D_h^2} \vec{v} +
+\frac{\rho}{2} \frac{2(1.75)}{D_h} \vec{v} |\vec{v}| = 0
+\end{equation}
+
+and arrive at the Ergun expression for the Darcy coefficient:
+
+\begin{equation}
+\frac{150}{D_h^2}
+\end{equation}
+
+and the Ergun expression for the Forchheimer coefficient:
+
+\begin{equation}
+\frac{2 \cdot 1.75}{D_h}
+\end{equation}
+
+where we have made the $2 \cdot 1.75$ multiplication explicit to make the 1.75 factor
+from the [Ergun wikipedia page](https://en.wikipedia.org/wiki/Ergun_equation)
+more recognizable. We perform a similar separation in the implementation of the
+Ergun Forchheimer coefficient outlined in [FunctorErgunDragCoefficients.md].
 
 !syntax parameters /FVKernels/PINSFVMomentumFriction
 

modules/navier_stokes/include/base/NS.h

@@ -22,6 +22,7 @@ using namespace HeatConduction;
 static const std::string directions[3] = {"x", "y", "z"};
 
 // geometric quantities
+static const std::string pebble_diameter = "pebble_diameter";
 static const std::string infinite_porosity = "infinite_porosity";
 static const std::string axis = "axis";
 static const std::string center = "center";

modules/navier_stokes/include/base/NSFVBase.h

@@ -1943,7 +1943,6 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
     params.template set<MooseFunctorName>(NS::density) = _density_name;
     params.template set<UserObjectName>("rhie_chow_user_object") =
         prefix() + "pins_rhie_chow_interpolator";
-    params.template set<MooseFunctorName>(NS::porosity) = _flow_porosity_functor_name;
 
     for (unsigned int block_i = 0; block_i < num_used_blocks; ++block_i)
     {
@@ -1967,10 +1966,16 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
         {
           const auto upper_name = MooseUtils::toUpper(_friction_types[block_i][type_i]);
           if (upper_name == "DARCY")
+          {
+            params.template set<MooseFunctorName>(NS::mu) = _dynamic_viscosity_name;
             params.template set<MooseFunctorName>("Darcy_name") = _friction_coeffs[block_i][type_i];
+          }
           else if (upper_name == "FORCHHEIMER")
+          {
             params.template set<MooseFunctorName>("Forchheimer_name") =
                 _friction_coeffs[block_i][type_i];
+            params.template set<MooseFunctorName>(NS::speed) = NS::speed;
+          }
         }
 
         getProblem().addFVKernel(kernel_type,
@@ -1990,7 +1995,6 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
         corr_params.template set<MooseFunctorName>(NS::density) = _density_name;
         corr_params.template set<UserObjectName>("rhie_chow_user_object") =
             prefix() + "pins_rhie_chow_interpolator";
-        corr_params.template set<MooseFunctorName>(NS::porosity) = _flow_porosity_functor_name;
         corr_params.template set<Real>("consistent_scaling") =
             parameters().template get<Real>("consistent_scaling");
         for (unsigned int d = 0; d < _dim; ++d)
@@ -2001,11 +2005,17 @@ NSFVBase<BaseType>::addINSMomentumFrictionKernels()
           {
             const auto upper_name = MooseUtils::toUpper(_friction_types[block_i][type_i]);
             if (upper_name == "DARCY")
+            {
+              corr_params.template set<MooseFunctorName>(NS::mu) = _dynamic_viscosity_name;
               corr_params.template set<MooseFunctorName>("Darcy_name") =
                   _friction_coeffs[block_i][type_i];
+            }
             else if (upper_name == "FORCHHEIMER")
+            {
               corr_params.template set<MooseFunctorName>("Forchheimer_name") =
                   _friction_coeffs[block_i][type_i];
+              corr_params.template set<MooseFunctorName>(NS::speed) = NS::speed;
+            }
           }
 
           getProblem().addFVKernel(correction_kernel_type,
@@ -2875,14 +2885,15 @@ template <class BaseType>
 void
 NSFVBase<BaseType>::addEnthalpyMaterial()
 {
-  InputParameters params = getFactory().getValidParams("INSFVEnthalpyMaterial");
+  InputParameters params = getFactory().getValidParams("INSFVEnthalpyFunctorMaterial");
   assignBlocks(params, _blocks);
 
   params.template set<MooseFunctorName>(NS::density) = _density_name;
   params.template set<MooseFunctorName>(NS::cp) = _specific_heat_name;
   params.template set<MooseFunctorName>("temperature") = _fluid_temperature_name;
 
-  getProblem().addMaterial("INSFVEnthalpyMaterial", prefix() + "ins_enthalpy_material", params);
+  getProblem().addMaterial(
+      "INSFVEnthalpyFunctorMaterial", prefix() + "ins_enthalpy_material", params);
 }
 
 template <class BaseType>

modules/navier_stokes/include/functormaterials/FunctorDragCoefficients.h

@@ -0,0 +1,25 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorMaterial.h"
+
+/**
+ * Abstract base class material providing the drag coefficients for linear and quadratic friction
+ * models. The expression of the coefficients is highly dependent on the formulation of the kernel.
+ * The reader should consult the PINSFVMomentumFrictionKernel documentation for details
+ */
+class FunctorDragCoefficients : public FunctorMaterial
+{
+public:
+  FunctorDragCoefficients(const InputParameters & parameters);
+
+  static InputParameters validParams();
+};

modules/navier_stokes/include/functormaterials/FunctorEffectiveFluidThermalConductivity.h

@@ -0,0 +1,29 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "FunctorMaterial.h"
+/**
+ * This is a base class material to calculate the effective thermal
+ * conductivity of the fluid phase. Unlike for the solid effective thermal
+ * conductivity, a factor of porosity is removed such that all parameters in the
+ * fluid effective thermal conductivity calculation can be expressed as porosity
+ * multiplying some parameter as $\kappa_f=\epsilon K$, where this material
+ * actually computes $K$. This is performed such that the correct spatial
+ * derivative of $\kappa_f$ can be performed, since the porosity is in general a
+ * function of space, and spatial derivatives of materials cannot be computed.
+ */
+class FunctorEffectiveFluidThermalConductivity : public FunctorMaterial
+{
+public:
+  FunctorEffectiveFluidThermalConductivity(const InputParameters & parameters);
+
+  static InputParameters validParams();
+};