Merge pull request #28952 from freiler/linearFV-add-ons

Linear FV Fluid Energy equation

framework/include/linearfvkernels/LinearFVSource.h

@@ -33,4 +33,7 @@ class LinearFVSource : public LinearFVElementalKernel
 protected:
   /// The functor for the source density
   const Moose::Functor<Real> & _source_density;
+
+  /// Scale factor
+  const Real _scale;
 };

framework/src/linearfvkernels/LinearFVSource.C

@@ -21,11 +21,14 @@ LinearFVSource::validParams()
       "Represents the matrix and right hand side contributions of a "
       "solution-independent source term in a partial differential equation.");
   params.addParam<MooseFunctorName>("source_density", 1.0, "The source density.");
+  params.addParam<Real>("scaling_factor", 1.0, "Coefficient to multiply the body force term with.");
   return params;
 }
 
 LinearFVSource::LinearFVSource(const InputParameters & params)
-  : LinearFVElementalKernel(params), _source_density(getFunctor<Real>("source_density"))
+  : LinearFVElementalKernel(params),
+    _source_density(getFunctor<Real>("source_density")),
+    _scale(getParam<Real>("scaling_factor"))
 {
 }
 
@@ -41,6 +44,6 @@ Real
 LinearFVSource::computeRightHandSideContribution()
 {
   // The contribution to the right hand side is s_C*V_C
-  return _source_density(makeElemArg(_current_elem_info->elem()), determineState()) *
+  return _scale * _source_density(makeElemArg(_current_elem_info->elem()), determineState()) *
          _current_elem_volume;
 }

modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVEnergyAdvection.md

@@ -0,0 +1,23 @@
+# LinearFVEnergyAdvection
+
+This kernel adds the contributions of the energy advection term to the matrix and right hand side of the energy equation system for the finite volume SIMPLE segregated solver [SIMPLE.md].
+
+This term is described by $\nabla \cdot \left(\rho\vec{u} c_p T \right)$ present in the energy equation conservation for an incompressible/weakly-compressible formulation. Currently, the kernel only supports +constant specific heat values+ for $c_p$ and solves for a temperature variable. The specific heat value needs to be prescribed, otherwise it will default to its default value of $c_p=1$.
+
+For FV, the integral of the advection term over a cell can be expressed as:
+
+\begin{equation}
+\int\limits_{V_C} \nabla \cdot \left(\rho\vec{u} c_p T \right) dV \approx \sum\limits_f (\rho \vec{u}\cdot \vec{n})_{RC} c_p T_f |S_f| \,
+\end{equation}
+
+where $T_f$ is a face temperature. The temperature acts as the advected quantity and an interpolation scheme (e.g. upwind) can be used to compute the face value. This kernel adds the face contribution for each face $f$ to the right hand side and matrix.
+
+The face mass flux $(\rho \vec{u}\cdot \vec{n})_{RC}$ is provided by the [RhieChowMassFlux.md] object which uses pressure
+gradients and the discrete momentum equation to compute face velocities and mass fluxes.
+For more information on the expression that is used, see [SIMPLE.md].
+
+!syntax parameters /LinearFVKernels/LinearFVEnergyAdvection
+
+!syntax inputs /LinearFVKernels/LinearFVEnergyAdvection
+
+!syntax children /LinearFVKernels/LinearFVEnergyAdvection

modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVMomentumBoussinesq.md

@@ -0,0 +1,14 @@
+# LinearFVMomentumBoussinesq
+
+This kernel adds the contributions of the Boussinesq buoyancy treatment for density through a force/source term to the right hand side of the momentum equation system for the finite volume SIMPLE segregated solver [SIMPLE.md]. The Boussinesq buoyancy treatment is applicable for low changes in density, and assumes constant density value in all other equation terms.
+
+This term is described by $-\rho_{ref}\alpha\vec{g}(T - T_{ref})$ present in the momentum equation conservation when describing an incompressible fluid, where $\rho_{ref}$ is the reference density, $\alpha$ is the thermal expansion coefficient, $\vec{g}$ is the gravity vector, $T$ is the temperature, and $T_{ref}$ is a reference temperature. The Boussinesq buoyancy model assumes the changes in density as a function of temperature are linear and relevant only in the buoyant force term of the equation system. The Boussinesq kernel allows for modeling natural convection.
+
+This term deals only with the force due to the variation in density $\Delta \rho \vec{g}$, with the fluid density being $\rho = \rho_{ref}+\Delta\rho$. Thus, with no extra added terms to the conventional incompressible Navier Stokes equations, the system will solve for the total pressure minus the hydrostatic pressure.
+For natural convection simulations, it is advisable to compute relevant dimensionless numbers such as the Rayleigh number or the Richardson number to decide on the need for turbulence models, mesh refinement and stability considerations.
+
+!syntax parameters /LinearFVKernels/LinearFVMomentumBoussinesq
+
+!syntax inputs /LinearFVKernels/LinearFVMomentumBoussinesq
+
+!syntax children /LinearFVKernels/LinearFVMomentumBoussinesq

modules/navier_stokes/include/executioners/SIMPLE.h

@@ -52,18 +52,38 @@ class SIMPLE : public SegregatedSolverBase
   /// the pressure equation.
   std::pair<unsigned int, Real> solvePressureCorrector();
 
+  /// Solve an equation which contains an advection term that depends
+  /// on the solution of the segregated Navier-Stokes equations.
+  /// @param system_num The number of the system which is solved
+  /// @param system Reference to the system which is solved
+  /// @param relaxation_factor The relaxation factor for matrix relaxation
+  /// @param solver_config The solver configuration object for the linear solve
+  /// @param abs_tol The scaled absolute tolerance for the linear solve
+  /// @return The normalized residual norm of the equation.
+  std::pair<unsigned int, Real> solveAdvectedSystem(const unsigned int system_num,
+                                                    LinearSystem & system,
+                                                    const Real relaxation_factor,
+                                                    SolverConfiguration & solver_config,
+                                                    const Real abs_tol);
+
   /// The number(s) of the system(s) corresponding to the momentum equation(s)
   std::vector<unsigned int> _momentum_system_numbers;
 
   /// The number of the system corresponding to the pressure equation
   const unsigned int _pressure_sys_number;
 
+  /// The number of the system corresponding to the energy equation
+  const unsigned int _energy_sys_number;
+
   /// Pointer(s) to the system(s) corresponding to the momentum equation(s)
   std::vector<LinearSystem *> _momentum_systems;
 
   /// Reference to the nonlinear system corresponding to the pressure equation
   LinearSystem & _pressure_system;
 
+  /// Pointer to the nonlinear system corresponding to the fluid energy equation
+  LinearSystem * _energy_system;
+
   /// Pointer to the segregated RhieChow interpolation object
   RhieChowMassFlux * _rc_uo;
 };

modules/navier_stokes/include/linearfvkernels/LinearFVEnergyAdvection.h

@@ -0,0 +1,58 @@
+//* 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 "LinearFVFluxKernel.h"
+#include "RhieChowMassFlux.h"
+#include "LinearFVAdvectionDiffusionBC.h"
+
+/**
+ * An advection kernel that implements the advection term for the enthalpy in the
+ * energy equation.
+ */
+class LinearFVEnergyAdvection : public LinearFVFluxKernel
+{
+public:
+  static InputParameters validParams();
+  LinearFVEnergyAdvection(const InputParameters & params);
+
+  virtual Real computeElemMatrixContribution() override;
+
+  virtual Real computeNeighborMatrixContribution() override;
+
+  virtual Real computeElemRightHandSideContribution() override;
+
+  virtual Real computeNeighborRightHandSideContribution() override;
+
+  virtual Real computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc) override;
+
+  virtual Real computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc) override;
+
+  virtual void setupFaceData(const FaceInfo * face_info) override;
+
+protected:
+  /// The Rhie-Chow user object that provides us with the face velocity
+  const RhieChowMassFlux & _mass_flux_provider;
+
+  /// The constant specific heat value
+  const Real _cp;
+
+private:
+  /// Container for the current advected interpolation coefficients on the face to make sure
+  /// we don't compute it multiple times for different terms.
+  std::pair<Real, Real> _advected_interp_coeffs;
+
+  /// Container for the mass flux on the face which will be reused in the advection term's
+  /// matrix and right hand side contribution
+  Real _face_mass_flux;
+
+  /// The interpolation method to use for the advected quantity
+  Moose::FV::InterpMethod _advected_interp_method;
+};

modules/navier_stokes/include/linearfvkernels/LinearFVMomentumBoussinesq.h

@@ -0,0 +1,49 @@
+//* 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 "LinearFVElementalKernel.h"
+
+/**
+ * Kernel that adds the component of the Boussinesq term in the momentum
+ * equations to the right hand side.
+ */
+class LinearFVMomentumBoussinesq : public LinearFVElementalKernel
+{
+public:
+  static InputParameters validParams();
+
+  /**
+   * Class constructor.
+   * @param params The InputParameters for the kernel.
+   */
+  LinearFVMomentumBoussinesq(const InputParameters & params);
+
+  virtual Real computeMatrixContribution() override;
+
+  virtual Real computeRightHandSideContribution() override;
+
+protected:
+  /// Fluid Temperature
+  const MooseLinearVariableFV<Real> & getTemperatureVariable(const std::string & vname);
+
+  /// Index x|y|z of the momentum equation component
+  const unsigned int _index;
+  /// Pointer to the linear finite volume temperature variable
+  const MooseLinearVariableFV<Real> & _temperature_var;
+  /// The gravity vector
+  const RealVectorValue _gravity;
+  /// The thermal expansion coefficient
+  const Moose::Functor<Real> & _alpha;
+  /// Reference temperature at which the reference value of the density (_rho) was measured
+  const Real _ref_temperature;
+  /// the density
+  const Moose::Functor<Real> & _rho;
+};

modules/navier_stokes/src/executioners/SIMPLE.C

@@ -35,16 +35,11 @@ SIMPLE::validParams()
   /*
    * We suppress parameters which are not supported yet
    */
-  params.suppressParameter<SolverSystemName>("energy_system");
   params.suppressParameter<SolverSystemName>("solid_energy_system");
   params.suppressParameter<std::vector<SolverSystemName>>("passive_scalar_systems");
   params.suppressParameter<std::vector<SolverSystemName>>("turbulence_systems");
-  params.suppressParameter<Real>("energy_equation_relaxation");
   params.suppressParameter<std::vector<Real>>("passive_scalar_equation_relaxation");
   params.suppressParameter<std::vector<Real>>("turbulence_equation_relaxation");
-  params.suppressParameter<MultiMooseEnum>("energy_petsc_options");
-  params.suppressParameter<MultiMooseEnum>("energy_petsc_options_iname");
-  params.suppressParameter<std::vector<std::string>>("energy_petsc_options_value");
   params.suppressParameter<MultiMooseEnum>("solid_energy_petsc_options");
   params.suppressParameter<MultiMooseEnum>("solid_energy_petsc_options_iname");
   params.suppressParameter<std::vector<std::string>>("solid_energy_petsc_options_value");
@@ -54,13 +49,9 @@ SIMPLE::validParams()
   params.suppressParameter<MultiMooseEnum>("turbulence_petsc_options");
   params.suppressParameter<MultiMooseEnum>("turbulence_petsc_options_iname");
   params.suppressParameter<std::vector<std::string>>("turbulence_petsc_options_value");
-  params.suppressParameter<Real>("energy_absolute_tolerance");
   params.suppressParameter<Real>("solid_energy_absolute_tolerance");
   params.suppressParameter<std::vector<Real>>("passive_scalar_absolute_tolerance");
   params.suppressParameter<std::vector<Real>>("turbulence_absolute_tolerance");
-  params.suppressParameter<Real>("energy_l_tol");
-  params.suppressParameter<Real>("energy_l_abs_tol");
-  params.suppressParameter<unsigned int>("energy_l_max_its");
   params.suppressParameter<Real>("solid_energy_l_tol");
   params.suppressParameter<Real>("solid_energy_l_abs_tol");
   params.suppressParameter<unsigned int>("solid_energy_l_max_its");
@@ -77,7 +68,11 @@ SIMPLE::validParams()
 SIMPLE::SIMPLE(const InputParameters & parameters)
   : SegregatedSolverBase(parameters),
     _pressure_sys_number(_problem.linearSysNum(getParam<SolverSystemName>("pressure_system"))),
-    _pressure_system(_problem.getLinearSystem(_pressure_sys_number))
+    _energy_sys_number(_has_energy_system
+                           ? _problem.linearSysNum(getParam<SolverSystemName>("energy_system"))
+                           : libMesh::invalid_uint),
+    _pressure_system(_problem.getLinearSystem(_pressure_sys_number)),
+    _energy_system(_has_energy_system ? &_problem.getLinearSystem(_energy_sys_number) : nullptr)
 {
   // We fetch the system numbers for the momentum components plus add vectors
   // for removing the contribution from the pressure gradient terms.
@@ -253,6 +248,74 @@ SIMPLE::solvePressureCorrector()
   return std::make_pair(its_res_pair.first, pressure_solver.get_initial_residual() / norm_factor);
 }
 
+std::pair<unsigned int, Real>
+SIMPLE::solveAdvectedSystem(const unsigned int system_num,
+                            LinearSystem & system,
+                            const Real relaxation_factor,
+                            SolverConfiguration & solver_config,
+                            const Real absolute_tol)
+{
+  _problem.setCurrentLinearSystem(system_num);
+
+  // We will need some members from the implicit linear system
+  LinearImplicitSystem & li_system = libMesh::cast_ref<LinearImplicitSystem &>(system.system());
+
+  // We will need the solution, the right hand side and the matrix
+  NumericVector<Number> & current_local_solution = *(li_system.current_local_solution);
+  NumericVector<Number> & solution = *(li_system.solution);
+  SparseMatrix<Number> & mmat = *(li_system.matrix);
+  NumericVector<Number> & rhs = *(li_system.rhs);
+
+  mmat.zero();
+  rhs.zero();
+
+  // We need a vector that stores the (diagonal_relaxed-original_diagonal) vector
+  auto diff_diagonal = solution.zero_clone();
+
+  // Fetch the linear solver from the system
+  PetscLinearSolver<Real> & linear_solver =
+      libMesh::cast_ref<PetscLinearSolver<Real> &>(*li_system.get_linear_solver());
+
+  _problem.computeLinearSystemSys(li_system, mmat, rhs, true);
+
+  // Go and relax the system matrix and the right hand side
+  relaxMatrix(mmat, relaxation_factor, *diff_diagonal);
+  relaxRightHandSide(rhs, solution, *diff_diagonal);
+
+  if (_print_fields)
+  {
+    _console << system.name() << " system matrix" << std::endl;
+    mmat.print();
+  }
+
+  // We compute the normalization factors based on the fluxes
+  Real norm_factor = computeNormalizationFactor(solution, mmat, rhs);
+
+  // We need the non-preconditioned norm to be consistent with the norm factor
+  LIBMESH_CHKERR(KSPSetNormType(linear_solver.ksp(), KSP_NORM_UNPRECONDITIONED));
+
+  // Setting the linear tolerances and maximum iteration counts
+  solver_config.real_valued_data["abs_tol"] = absolute_tol * norm_factor;
+  linear_solver.set_solver_configuration(solver_config);
+
+  // Solve the system and update current local solution
+  auto its_res_pair = linear_solver.solve(mmat, mmat, solution, rhs);
+  li_system.update();
+
+  if (_print_fields)
+  {
+    _console << " rhs when we solve " << system.name() << std::endl;
+    rhs.print();
+    _console << system.name() << " solution " << std::endl;
+    solution.print();
+    _console << " Norm factor " << norm_factor << std::endl;
+  }
+
+  system.setSolution(current_local_solution);
+
+  return std::make_pair(its_res_pair.first, linear_solver.get_initial_residual() / norm_factor);
+}
+
 void
 SIMPLE::execute()
 {
@@ -307,15 +370,18 @@ SIMPLE::execute()
     unsigned int iteration_counter = 0;
 
     // Assign residuals to general residual vector
-    unsigned int no_systems = _momentum_systems.size() + 1;
+    unsigned int no_systems = _momentum_systems.size() + 1 + _has_energy_system;
     std::vector<std::pair<unsigned int, Real>> ns_residuals(no_systems, std::make_pair(0, 1.0));
     std::vector<Real> ns_abs_tols(_momentum_systems.size(), _momentum_absolute_tolerance);
     ns_abs_tols.push_back(_pressure_absolute_tolerance);
+    if (_has_energy_system)
+      ns_abs_tols.push_back(_energy_absolute_tolerance);
 
     // Loop until converged or hit the maximum allowed iteration number
     while (iteration_counter < _num_iterations && !converged(ns_residuals, ns_abs_tols))
     {
       iteration_counter++;
+      size_t residual_index = 0;
 
       // We set the preconditioner/controllable parameters through petsc options. Linear
       // tolerances will be overridden within the solver. In case of a segregated momentum
@@ -360,6 +426,25 @@ SIMPLE::execute()
       // Reconstruct the cell velocity as well to accelerate convergence
       _rc_uo->computeCellVelocity();
 
+      // Update residual index
+      residual_index = momentum_residual.size();
+
+      // If we have an energy equation, solve it here. We assume the material properties in the
+      // Navier-Stokes equations depend on temperature, therefore we can not solve for temperature
+      // outside of the velocity-pressure loop
+      if (_has_energy_system)
+      {
+        // We set the preconditioner/controllable parameters through petsc options. Linear
+        // tolerances will be overridden within the solver.
+        Moose::PetscSupport::petscSetOptions(_energy_petsc_options, solver_params);
+        residual_index += 1;
+        ns_residuals[residual_index] = solveAdvectedSystem(_energy_sys_number,
+                                                           *_energy_system,
+                                                           _energy_equation_relaxation,
+                                                           _energy_linear_control,
+                                                           _energy_l_abs_tol);
+      }
+      _problem.execute(EXEC_NONLINEAR);
       // Printing residuals
       _console << "Iteration " << iteration_counter << " Initial residual norms:" << std::endl;
       for (auto system_i : index_range(_momentum_systems))
@@ -373,6 +458,15 @@ SIMPLE::execute()
       _console << " Pressure equation: " << COLOR_GREEN
                << ns_residuals[momentum_residual.size()].second << COLOR_DEFAULT
                << " Linear its: " << ns_residuals[momentum_residual.size()].first << std::endl;
+      residual_index = momentum_residual.size();
+
+      if (_has_energy_system)
+      {
+        residual_index += 1;
+        _console << " Energy equation: " << COLOR_GREEN << ns_residuals[residual_index].second
+                 << COLOR_DEFAULT << " Linear its: " << ns_residuals[residual_index].first
+                 << std::endl;
+      }
     }
   }