Merge pull request #27445 from tanoret/interface_area_two_phase

Adding interface area to two phase

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

@@ -456,3 +456,14 @@ @article{lindsay2023moose
   year={2023},
   publisher={Elsevier}
 }
+
+@article{hibiki2000interface,
+  title={One-group interfacial area transport of bubbly flows in vertical round tubes},
+  author={Hibiki, T and Ishii, M},
+  journal={International Journal of Heat and Mass Transfer},
+  volume={43},
+  number={15},
+  pages={2711--2726},
+  year={2000},
+  publisher={Elsevier}
+}

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

@@ -0,0 +1,114 @@
+# WCNSFV2PInterfaceAreaSourceSink
+
+The interfacial area concentration is defined as the interface area between
+two phases per unit volume, i.e.,$[\xi_p] = \frac{m^2}{m^3}$. This parameter is important
+for predicting mass, momentum, and energy transfer at interfaces in two-phase flows.
+
+The general equation for interfacial area concentration transport via the mixture model reads as follows:
+
+\begin{equation}
+\frac{\partial (\rho_d \xi_p)}{\partial t} + \nabla \cdot \left( \rho_d \vec{u} \xi_p \right) - \nabla \left( D_b \nabla \xi_p \right) =
+-\frac{1}{3} \frac{D \rho_d}{Dt} + S_T + \rho_d (S_C + S_B)\,,
+\end{equation}
+
+where:
+
+- $\rho_d$ is the density of the dispersed phase $d$, e.g., the density of the gas in bubbles,
+- $t$ is time,
+- $\vec{u}$ is the velocity vector,
+- $D_b$ is a diffusion coefficient for the interface area concentration, which may be assumed to be `0` if unknown,
+- $\frac{D (\bullet)}{Dt} = \frac{\partial (\bullet)}{\partial t} + \vec{u} \cdot \nabla (\bullet)$ is the material derivative,
+- $S_T$, $S_C$, and $S_B$ are the interface area concentration sources due to mass transfer, coalescence, and breakage, respectively.
+
+The terms on the left-hand side of this equation are modeled via standard kernels for the mixture model.
+For example, in an open flow case, the time derivative may be modeled using [FVFunctorTimeKernel.md],
+the advection term using [INSFVScalarFieldAdvection.md], and the diffusion term using [FVDiffusion.md].
+
+The terms on the right-hand side are modeled using a multidimensional version of
+Hibiki and Ishii's model [!citep](hibiki2000interface).
+In this one, the sources are approximated as follows:
+
+
+## Interface area concentration source due to mass transfer
+
+The interface area concentration source due to mass transfer is modeled as follows:
+
+\begin{equation}
+S_T =
+\begin{cases}
+  \frac{6 \alpha_d}{d_p}, & \text{if } \alpha_d < \alpha_d^{co} \\
+  \frac{2}{3} \cdot h^{b \rightarrow d} \left( \frac{1}{\alpha_d} - 2.0 \right), & \text{otherwise}
+\end{cases}\,,
+\end{equation}
+
+where:
+
+- $\alpha_d$ is the volumetric fraction of the dispersed phase, e.g., the void fraction if the dispersed phase is a gas,
+- $\alpha_d^{co}$ is a cutoff fraction for mass transfer model selection,
+- $d_p$ is the best estimate for the dispersed phase particle diameter,
+- $h^{b \rightarrow d}$ is the mass exchange coefficient from the bulk to the dispersed phase,
+
+The cutoff volumetric fraction $\alpha_d^{co}$ is included because the mass transfer term
+in Hibiki and Ishii's model is not physical for low volumetric fractions of the dispersed phase.
+Below the cutoff limit, the dispersed phase is modeled as spherical particles.
+
+!alert note
+The user should select the cutoff volumetric fraction $\alpha_d^{co}$ as the limit at which
+the modeled flows transition away from bubbly flow, i.e., below the cutoff limit there is an
+implicit assumption that the flow behaves as a bubbly flow.
+
+
+## Interface area concentration source due to coalescence
+
+The reduction in interface area concentration due to coalescence of the dispersed phase is modeled as follows:
+
+\begin{equation}
+S_C = -\left( \frac{\alpha_d}{\xi_p} \right)^2 \frac{\Gamma_C \alpha_d^2 u_{\epsilon}}{\tilde{d_p}^{11/3} (\alpha_d^{max}- \alpha_d)}
+\operatorname{exp} \left( -K_C \frac{\tilde{d_p}^{5/3} \rho_b^{1/2} u_{\epsilon}^{1/3}}{\sigma^{1/2}} \right)\,,
+\end{equation}
+
+where:
+
+- $\alpha_d$ is the volumetric fraction of the dispersed phase, e.g., the void fraction if the dispersed phase is a gas,
+- $u_{\epsilon} = \left( \| \vec{u} \| \ell \| \nabla p \| / \rho_m \right)^{1/3}$ is the friction velocity due to pressure gradients, with $\| \vec{u} \|$ being the norm of the velocity vector, $\ell$ a characteristic mixing length, $\| \nabla p \|$ the norm of the pressure gradient, and $\rho_m$ the mixture density,
+- $\tilde{d_p} = \psi \frac{\alpha_d}{\xi_p}$ is the model estimate for the dispersed phase particle diameter, with $\psi$ being a shape factor, which is, for example, $\psi=6$ for spherical particles,
+- $\alpha_d^{max}$ is the maximum volumetric fraction admitted by the model, in absence of data we recommend taking $\alpha_d^{max}=1$,
+- $\rho_b$ is the bulk phase density, e.g., for air bubbles in a water flow it would be the density of water,
+- $\sigma$ is the surface tension between the two phases,
+- $\Gamma_C = 0.188$ and $K_C = 0.129$ are closure coefficients of the model.
+
+!alert note
+Many of the parameters in the coalescence model are provided as functors,
+which means that spatially dependent fields may be specified for these parameters.
+However, the validation of this model has only been performed using constant parameters.
+
+## Interface area concentration source due to breakage
+
+The increase in interface area concentration due to breakage of the dispersed phase is modeled as follows:
+
+\begin{equation}
+S_B = \left( \frac{\alpha_d}{\xi_p} \right)^2 \frac{\Gamma_B \alpha_d (1 - \alpha_d) u_{\epsilon}}{\tilde{d_p}^{11/3} (\alpha_d^{max}- \alpha_d)}
+\operatorname{exp} \left( -K_B \frac{\sigma}{\rho_b \tilde{d_p}^{5/3} u_{\epsilon}^{2/3}} \right)\,,
+\end{equation}
+
+where:
+
+- $\alpha_d$ is the volumetric fraction of the dispersed phase, e.g., the void fraction if the dispersed phase is a gas,
+- $u_{\epsilon} = \left( \| \vec{u} \| \ell \| \nabla p \| / \rho_m \right)^{1/3}$ is the friction velocity due to pressure gradients, with $\| \vec{u} \|$ being the norm of the velocity vector, $\ell$ a characteristic mixing length, $\| \nabla p \|$ the norm of the pressure gradient, and $\rho_m$ the mixture density,
+- $\tilde{d_p} = \psi \frac{\alpha_d}{\xi_p}$ is the model estimate for the dispersed phase particle diameter, with $\psi$ being a shape factor, which is, for example, $\psi=6$ for spherical particles,
+- $\alpha_d^{max}$ is the maximum volumetric fraction admitted by the model, in absence of data we recommend taking $\alpha_d^{max}=1$,
+- $\rho_b$ is the bulk phase density, e.g., for air bubbles in a water flow it would be the density of water,
+- $\sigma$ is the surface tension between the two phases,
+- $\Gamma_B = 0.264$ and $K_B = 1.370$ are closure coefficients of the model.
+
+!alert note
+Many of the parameters in the breakage model are provided as functors,
+which means that spatially dependent fields may be specified for these parameters.
+However, the validation of this model has only been performed using constant parameters.
+
+
+!syntax parameters /FVKernels/WCNSFV2PInterfaceAreaSourceSink
+
+!syntax inputs /FVKernels/WCNSFV2PInterfaceAreaSourceSink
+
+!syntax children /FVKernels/WCNSFV2PInterfaceAreaSourceSink

modules/navier_stokes/include/fvkernels/WCNSFV2PInterfaceAreaSourceSink.h

@@ -0,0 +1,75 @@
+//* 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 "FVElementalKernel.h"
+#include "INSFVVelocityVariable.h"
+
+/**
+ * Computes source the sink terms for the interface area in the mixture model of two-phase flows.
+ */
+class WCNSFV2PInterfaceAreaSourceSink : public FVElementalKernel
+{
+public:
+  static InputParameters validParams();
+
+  WCNSFV2PInterfaceAreaSourceSink(const InputParameters & parameters);
+
+protected:
+  ADReal computeQpResidual() override;
+
+protected:
+  /// The dimension of the domain
+  const unsigned int _dim;
+
+  /// x-velocity
+  const Moose::Functor<ADReal> & _u_var;
+  /// y-velocity
+  const Moose::Functor<ADReal> * _v_var;
+  /// z-velocity
+  const Moose::Functor<ADReal> * _w_var;
+
+  /// Characterisitc Length
+  const Moose::Functor<ADReal> & _characheristic_length;
+
+  /// Mixture Density
+  const Moose::Functor<ADReal> & _rho_mixture;
+
+  /// Dispersed Phase Density
+  const Moose::Functor<ADReal> & _rho_d;
+
+  /// Pressure field
+  const Moose::Functor<ADReal> & _pressure;
+
+  /// Interface Mass Exchange Coeefficient
+  const Moose::Functor<ADReal> & _mass_exchange_coefficient;
+
+  /// Void Fraction
+  const Moose::Functor<ADReal> & _f_d;
+
+  /// Maximum Void Fraction
+  const Real _f_d_max;
+
+  /// Surface Tension
+  const Moose::Functor<ADReal> & _sigma;
+
+  /// Particle Diameter
+  const Moose::Functor<ADReal> & _particle_diameter;
+
+  /// Cutoff fraction at which the full mass transfer model is activated
+  const Real _cutoff_fraction;
+
+  /// Internal closure coefficients
+  static constexpr Real _gamma_c = 0.188;
+  static constexpr Real _Kc = 0.129;
+  static constexpr Real _gamma_b = 0.264;
+  static constexpr Real _Kb = 1.37;
+  static constexpr Real _shape_factor = 6.0;
+};

modules/navier_stokes/src/fvkernels/WCNSFV2PInterfaceAreaSourceSink.C

@@ -0,0 +1,154 @@
+//* 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
+
+#include "WCNSFV2PInterfaceAreaSourceSink.h"
+#include "NS.h"
+#include "NonlinearSystemBase.h"
+#include "NavierStokesMethods.h"
+
+registerMooseObject("NavierStokesApp", WCNSFV2PInterfaceAreaSourceSink);
+
+InputParameters
+WCNSFV2PInterfaceAreaSourceSink::validParams()
+{
+  InputParameters params = FVElementalKernel::validParams();
+  params.addClassDescription(
+      "Source and sink of interfacial area for two-phase flow mixture model.");
+  params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
+  params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
+  params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
+
+  params.addParam<MooseFunctorName>("L", 1.0, "The characteristic dissipation length.");
+  params.addRequiredParam<MooseFunctorName>(NS::density, "Continuous phase density.");
+  params.addRequiredParam<MooseFunctorName>(NS::density + std::string("_d"),
+                                            "Dispersed phase density.");
+  params.addRequiredParam<MooseFunctorName>(NS::pressure, "Continuous phase density.");
+  params.addParam<MooseFunctorName>(
+      "k_c", 0.0, "Mass exchange coefficients from continous to dispersed phases.");
+  params.addParam<MooseFunctorName>("fd", 0.0, "Fraction dispersed phase.");
+  params.addParam<Real>("fd_max", 1.0, "Maximum dispersed phase fraction.");
+
+  params.addParam<MooseFunctorName>("sigma", 1.0, "Surface tension between phases.");
+  params.addParam<MooseFunctorName>("particle_diameter", 1.0, "Maximum particle diameter.");
+
+  params.addParam<Real>("cutoff_fraction",
+                        0.1,
+                        "Void fraction at which the interface area density mass transfer model is "
+                        "activated. Below this fraction, spherical bubbles are assumed.");
+
+  params.set<unsigned short>("ghost_layers") = 2;
+
+  return params;
+}
+
+WCNSFV2PInterfaceAreaSourceSink::WCNSFV2PInterfaceAreaSourceSink(const InputParameters & params)
+  : FVElementalKernel(params),
+    _dim(_subproblem.mesh().dimension()),
+    _u_var(getFunctor<ADReal>("u")),
+    _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
+    _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
+    _characheristic_length(getFunctor<ADReal>("L")),
+    _rho_mixture(getFunctor<ADReal>(NS::density)),
+    _rho_d(getFunctor<ADReal>(NS::density + std::string("_d"))),
+    _pressure(getFunctor<ADReal>(NS::pressure)),
+    _mass_exchange_coefficient(getFunctor<ADReal>("k_c")),
+    _f_d(getFunctor<ADReal>("fd")),
+    _f_d_max(getParam<Real>("fd_max")),
+    _sigma(getFunctor<ADReal>("sigma")),
+    _particle_diameter(getFunctor<ADReal>("particle_diameter")),
+    _cutoff_fraction(getParam<Real>("cutoff_fraction"))
+{
+  if (_dim >= 2 && !_v_var)
+    paramError("v", "In two or more dimensions, the v velocity must be supplied!");
+
+  if (_dim >= 3 && !_w_var)
+    paramError("w", "In three or more dimensions, the w velocity must be supplied!");
+}
+
+ADReal
+WCNSFV2PInterfaceAreaSourceSink::computeQpResidual()
+{
+
+  // Useful Arguments
+  const auto state = determineState();
+  const auto elem_arg = makeElemArg(_current_elem);
+  const bool is_transient = _subproblem.isTransient();
+
+  // Current Variables
+  const auto u = _u_var(elem_arg, state);
+  const auto rho_d = _rho_d(elem_arg, state);
+  const auto rho_d_grad = _rho_d.gradient(elem_arg, state);
+  const auto xi = _var(elem_arg, state);
+  const auto rho_m = _rho_mixture(elem_arg, state);
+  const auto f_d = _f_d(elem_arg, state);
+  const auto sigma = _sigma(elem_arg, state);
+  const auto rho_l = (rho_m - f_d * rho_d) / (1.0 - f_d + libMesh::TOLERANCE);
+  const auto complement_fd = std::max(_f_d_max - f_d, libMesh::TOLERANCE);
+  const auto f_d_o_xi = f_d / (_var(elem_arg, state) + libMesh::TOLERANCE) + libMesh::TOLERANCE;
+  const auto f_d_o_xi_old =
+      f_d / (raw_value(_var(elem_arg, state)) + libMesh::TOLERANCE) + libMesh::TOLERANCE;
+
+  // Adding bubble compressibility term
+  ADReal material_time_derivative_rho_d = u * rho_d_grad(0);
+  if (_dim > 1)
+  {
+    material_time_derivative_rho_d += (*_v_var)(elem_arg, state) * rho_d_grad(1);
+    if (_dim > 2)
+    {
+      material_time_derivative_rho_d += (*_w_var)(elem_arg, state) * rho_d_grad(2);
+    }
+  }
+  if (is_transient)
+    material_time_derivative_rho_d +=
+        raw_value(_rho_d(elem_arg, state) - _rho_d(elem_arg, Moose::oldState())) / _dt;
+  const auto bubble_compressibility = material_time_derivative_rho_d * xi / 3.0;
+
+  // Adding area growth due to added mass
+  ADReal bubble_added_mass;
+  if (f_d < _cutoff_fraction)
+    bubble_added_mass = raw_value(_rho_d(elem_arg, state)) *
+                        (f_d * 6.0 / _particle_diameter(elem_arg, state) - _var(elem_arg, state));
+  else
+    bubble_added_mass = 2. / 3. * _mass_exchange_coefficient(elem_arg, state) *
+                        (1.0 / (f_d + libMesh::TOLERANCE) - 2.0);
+
+  // Model parameters
+  const auto db = _shape_factor * f_d_o_xi_old + libMesh::TOLERANCE;
+
+  ADRealVectorValue velocity(u);
+  if (_v_var)
+    velocity(1) = (*_v_var)(elem_arg, state);
+  if (_w_var)
+    velocity(2) = (*_w_var)(elem_arg, state);
+
+  const ADReal velocity_norm = NS::computeSpeed(velocity);
+  const auto pressure_gradient = raw_value(_pressure.gradient(elem_arg, state));
+  const Real pressure_grad_norm =
+      MooseUtils::isZero(pressure_gradient) ? 1e-42 : pressure_gradient.norm();
+
+  const auto u_eps =
+      std::pow(velocity_norm * _characheristic_length(elem_arg, state) * pressure_grad_norm / rho_m,
+               1. / 3.);
+
+  const auto interaction_prefactor =
+      Utility::pow<2>(f_d_o_xi) * u_eps / (std::pow(db, 11. / 3.) / complement_fd);
+
+  // Adding coalescence term
+  const auto f_c = interaction_prefactor * _gamma_c * Utility::pow<2>(f_d);
+  const auto exp_c = std::exp(-_Kc * std::pow(db, 5. / 6.) * std::sqrt(rho_l / sigma) * u_eps);
+  const auto s_rc = f_c * exp_c;
+
+  // Adding breakage term
+  const auto f_b = interaction_prefactor * _gamma_b * f_d * (1. - f_d);
+  const auto exp_b =
+      std::exp(-_Kb * sigma / (rho_l * std::pow(db, 5. / 3.) * Utility::pow<2>(u_eps)));
+  const auto s_rb = f_b * exp_b;
+
+  return -bubble_added_mass + bubble_compressibility + s_rc - s_rb;
+}