diff --git a/framework/include/linearfvkernels/LinearFVSource.h b/framework/include/linearfvkernels/LinearFVSource.h index 500dcf015a27..60ccdebfe28b 100644 --- a/framework/include/linearfvkernels/LinearFVSource.h +++ b/framework/include/linearfvkernels/LinearFVSource.h @@ -33,4 +33,7 @@ class LinearFVSource : public LinearFVElementalKernel protected: /// The functor for the source density const Moose::Functor & _source_density; + + /// Scale factor + const Real _scale; }; diff --git a/framework/src/linearfvkernels/LinearFVSource.C b/framework/src/linearfvkernels/LinearFVSource.C index b32fda7e382e..080282cfddb7 100644 --- a/framework/src/linearfvkernels/LinearFVSource.C +++ b/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("source_density", 1.0, "The source density."); + params.addParam("scaling_factor", 1.0, "Coefficient to multiply the body force term with."); return params; } LinearFVSource::LinearFVSource(const InputParameters & params) - : LinearFVElementalKernel(params), _source_density(getFunctor("source_density")) + : LinearFVElementalKernel(params), + _source_density(getFunctor("source_density")), + _scale(getParam("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; } diff --git a/modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVEnergyAdvection.md b/modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVEnergyAdvection.md new file mode 100644 index 000000000000..cf8564ea4405 --- /dev/null +++ b/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 diff --git a/modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVMomentumBoussinesq.md b/modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVMomentumBoussinesq.md new file mode 100644 index 000000000000..f74feca039f6 --- /dev/null +++ b/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 diff --git a/modules/navier_stokes/include/executioners/SIMPLE.h b/modules/navier_stokes/include/executioners/SIMPLE.h index 51bfad4decc7..7ffc5ee7d345 100644 --- a/modules/navier_stokes/include/executioners/SIMPLE.h +++ b/modules/navier_stokes/include/executioners/SIMPLE.h @@ -52,18 +52,38 @@ class SIMPLE : public SegregatedSolverBase /// the pressure equation. std::pair 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 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 _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 _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; }; diff --git a/modules/navier_stokes/include/linearfvkernels/LinearFVEnergyAdvection.h b/modules/navier_stokes/include/linearfvkernels/LinearFVEnergyAdvection.h new file mode 100644 index 000000000000..b315e78e878f --- /dev/null +++ b/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 _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; +}; diff --git a/modules/navier_stokes/include/linearfvkernels/LinearFVMomentumBoussinesq.h b/modules/navier_stokes/include/linearfvkernels/LinearFVMomentumBoussinesq.h new file mode 100644 index 000000000000..8722cae3b40c --- /dev/null +++ b/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 & 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 & _temperature_var; + /// The gravity vector + const RealVectorValue _gravity; + /// The thermal expansion coefficient + const Moose::Functor & _alpha; + /// Reference temperature at which the reference value of the density (_rho) was measured + const Real _ref_temperature; + /// the density + const Moose::Functor & _rho; +}; diff --git a/modules/navier_stokes/src/executioners/SIMPLE.C b/modules/navier_stokes/src/executioners/SIMPLE.C index f751e198ec12..214c7e5baa8c 100644 --- a/modules/navier_stokes/src/executioners/SIMPLE.C +++ b/modules/navier_stokes/src/executioners/SIMPLE.C @@ -35,16 +35,11 @@ SIMPLE::validParams() /* * We suppress parameters which are not supported yet */ - params.suppressParameter("energy_system"); params.suppressParameter("solid_energy_system"); params.suppressParameter>("passive_scalar_systems"); params.suppressParameter>("turbulence_systems"); - params.suppressParameter("energy_equation_relaxation"); params.suppressParameter>("passive_scalar_equation_relaxation"); params.suppressParameter>("turbulence_equation_relaxation"); - params.suppressParameter("energy_petsc_options"); - params.suppressParameter("energy_petsc_options_iname"); - params.suppressParameter>("energy_petsc_options_value"); params.suppressParameter("solid_energy_petsc_options"); params.suppressParameter("solid_energy_petsc_options_iname"); params.suppressParameter>("solid_energy_petsc_options_value"); @@ -54,13 +49,9 @@ SIMPLE::validParams() params.suppressParameter("turbulence_petsc_options"); params.suppressParameter("turbulence_petsc_options_iname"); params.suppressParameter>("turbulence_petsc_options_value"); - params.suppressParameter("energy_absolute_tolerance"); params.suppressParameter("solid_energy_absolute_tolerance"); params.suppressParameter>("passive_scalar_absolute_tolerance"); params.suppressParameter>("turbulence_absolute_tolerance"); - params.suppressParameter("energy_l_tol"); - params.suppressParameter("energy_l_abs_tol"); - params.suppressParameter("energy_l_max_its"); params.suppressParameter("solid_energy_l_tol"); params.suppressParameter("solid_energy_l_abs_tol"); params.suppressParameter("solid_energy_l_max_its"); @@ -77,7 +68,11 @@ SIMPLE::validParams() SIMPLE::SIMPLE(const InputParameters & parameters) : SegregatedSolverBase(parameters), _pressure_sys_number(_problem.linearSysNum(getParam("pressure_system"))), - _pressure_system(_problem.getLinearSystem(_pressure_sys_number)) + _energy_sys_number(_has_energy_system + ? _problem.linearSysNum(getParam("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 +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(system.system()); + + // We will need the solution, the right hand side and the matrix + NumericVector & current_local_solution = *(li_system.current_local_solution); + NumericVector & solution = *(li_system.solution); + SparseMatrix & mmat = *(li_system.matrix); + NumericVector & 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 & linear_solver = + libMesh::cast_ref &>(*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> ns_residuals(no_systems, std::make_pair(0, 1.0)); std::vector 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; + } } } diff --git a/modules/navier_stokes/src/linearfvkernels/LinearFVEnergyAdvection.C b/modules/navier_stokes/src/linearfvkernels/LinearFVEnergyAdvection.C new file mode 100644 index 000000000000..a8fb51986640 --- /dev/null +++ b/modules/navier_stokes/src/linearfvkernels/LinearFVEnergyAdvection.C @@ -0,0 +1,106 @@ +//* 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 "LinearFVEnergyAdvection.h" +#include "MooseLinearVariableFV.h" +#include "NSFVUtils.h" +#include "NS.h" + +registerMooseObject("NavierStokesApp", LinearFVEnergyAdvection); + +InputParameters +LinearFVEnergyAdvection::validParams() +{ + InputParameters params = LinearFVFluxKernel::validParams(); + params.addClassDescription("Represents the matrix and right hand side contributions of an " + "advection term for the energy e.g. h=int(cp dT). A user may still " + "override what quantity is advected, but the default is temperature."); + params.addParam("cp", 1, "Specific heat value"); + params.addRequiredParam( + "rhie_chow_user_object", + "The rhie-chow user-object which is used to determine the face velocity."); + params += Moose::FV::advectedInterpolationParameter(); + return params; +} + +LinearFVEnergyAdvection::LinearFVEnergyAdvection(const InputParameters & params) + : LinearFVFluxKernel(params), + _mass_flux_provider(getUserObject("rhie_chow_user_object")), + _cp(getParam("cp")), + _advected_interp_coeffs(std::make_pair(0, 0)), + _face_mass_flux(0.0) +{ + Moose::FV::setInterpolationMethod(*this, _advected_interp_method, "advected_interp_method"); +} + +Real +LinearFVEnergyAdvection::computeElemMatrixContribution() +{ + return _cp * _advected_interp_coeffs.first * _face_mass_flux * _current_face_area; +} + +Real +LinearFVEnergyAdvection::computeNeighborMatrixContribution() +{ + return _cp * _advected_interp_coeffs.second * _face_mass_flux * _current_face_area; +} + +Real +LinearFVEnergyAdvection::computeElemRightHandSideContribution() +{ + return 0.0; +} + +Real +LinearFVEnergyAdvection::computeNeighborRightHandSideContribution() +{ + return 0.0; +} + +Real +LinearFVEnergyAdvection::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc) +{ + const auto * const adv_bc = static_cast(&bc); + mooseAssert(adv_bc, "This should be a valid BC!"); + + const auto boundary_value_matrix_contrib = adv_bc->computeBoundaryValueMatrixContribution(); + + // We support internal boundaries too so we have to make sure the normal points always outward + const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0; + + return _cp * boundary_value_matrix_contrib * factor * _face_mass_flux * _current_face_area; +} + +Real +LinearFVEnergyAdvection::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc) +{ + const auto * const adv_bc = static_cast(&bc); + mooseAssert(adv_bc, "This should be a valid BC!"); + + // We support internal boundaries too so we have to make sure the normal points always outward + const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0); + + const auto boundary_value_rhs_contrib = adv_bc->computeBoundaryValueRHSContribution(); + return -_cp * boundary_value_rhs_contrib * factor * _face_mass_flux * _current_face_area; +} + +void +LinearFVEnergyAdvection::setupFaceData(const FaceInfo * face_info) +{ + LinearFVFluxKernel::setupFaceData(face_info); + + // Caching the mass flux on the face which will be reused in the advection term's matrix and right + // hand side contributions + _face_mass_flux = _mass_flux_provider.getMassFlux(*face_info); + + // Caching the interpolation coefficients so they will be reused for the matrix and right hand + // side terms + _advected_interp_coeffs = + interpCoeffs(_advected_interp_method, *_current_face_info, true, _face_mass_flux); +} diff --git a/modules/navier_stokes/src/linearfvkernels/LinearFVMomentumBoussinesq.C b/modules/navier_stokes/src/linearfvkernels/LinearFVMomentumBoussinesq.C new file mode 100644 index 000000000000..2bfd4533a2f3 --- /dev/null +++ b/modules/navier_stokes/src/linearfvkernels/LinearFVMomentumBoussinesq.C @@ -0,0 +1,79 @@ +//* 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 "LinearFVMomentumBoussinesq.h" +#include "Assembly.h" +#include "SubProblem.h" +#include "NS.h" + +registerMooseObject("NavierStokesApp", LinearFVMomentumBoussinesq); + +InputParameters +LinearFVMomentumBoussinesq::validParams() +{ + InputParameters params = LinearFVElementalKernel::validParams(); + params.addClassDescription("Represents the Boussinesq term in the Navier Stokes momentum " + "equations, added to the right hand side."); + params.addParam(NS::T_fluid, "The fluid temperature variable."); + params.addRequiredParam("gravity", "Gravitational acceleration vector."); + params.addParam("alpha_name", + NS::alpha_boussinesq, + "The name of the thermal expansion coefficient" + "this is of the form rho = rho_ref*(1-alpha (T-T_ref))"); + params.addRequiredParam("ref_temperature", "The value for the reference temperature."); + params.addRequiredParam(NS::density, "The value for the density"); + MooseEnum momentum_component("x=0 y=1 z=2"); + params.addRequiredParam( + "momentum_component", + momentum_component, + "The component of the momentum equation that this kernel applies to."); + return params; +} + +LinearFVMomentumBoussinesq::LinearFVMomentumBoussinesq(const InputParameters & params) + : LinearFVElementalKernel(params), + _index(getParam("momentum_component")), + _temperature_var(getTemperatureVariable(NS::T_fluid)), + _gravity(getParam("gravity")), + _alpha(getFunctor("alpha_name")), + _ref_temperature(getParam("ref_temperature")), + _rho(getFunctor(NS::density)) +{ + if (!_rho.isConstant()) + paramError(NS::density, "The density in the boussinesq term is not constant!"); +} + +const MooseLinearVariableFV & +LinearFVMomentumBoussinesq::getTemperatureVariable(const std::string & vname) +{ + auto * ptr = dynamic_cast *>( + &_fe_problem.getVariable(_tid, getParam(vname))); + + if (!ptr) + paramError(NS::T_fluid, + "The fluid temperature variable should be of type MooseLinearVariableFVReal!"); + + return *ptr; +} + +Real +LinearFVMomentumBoussinesq::computeMatrixContribution() +{ + return 0.0; +} + +Real +LinearFVMomentumBoussinesq::computeRightHandSideContribution() +{ + const auto elem = makeElemArg(_current_elem_info->elem()); + const auto state = determineState(); + return -_alpha(elem, state) * _gravity(_index) * _rho(elem, state) * + (_temperature_var.getElemValue(*_current_elem_info, state) - _ref_temperature) * + _current_elem_volume; +} diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq.i b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq.i new file mode 100644 index 000000000000..d4ed4dfd09b0 --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq.i @@ -0,0 +1,243 @@ +mu = 2.6 +rho = 1.0 +advected_interp_method = 'average' +cp = 300 +k = 20 +alpha_b = 1e-4 + +[Mesh] + [mesh] + type = CartesianMeshGenerator + dim = 2 + dx = '1.' + dy = '0.2' + ix = '10' + iy = '5' + [] +[] + +[Problem] + linear_sys_names = 'u_system v_system pressure_system energy_system' + previous_nl_solution_required = true +[] + +[UserObjects] + [rc] + type = RhieChowMassFlux + u = vel_x + v = vel_y + pressure = pressure + rho = ${rho} + p_diffusion_kernel = p_diffusion + [] +[] + +[Variables] + [vel_x] + type = MooseLinearVariableFVReal + initial_condition = 0.5 + solver_sys = u_system + [] + [vel_y] + type = MooseLinearVariableFVReal + solver_sys = v_system + initial_condition = 0.0 + [] + [pressure] + type = MooseLinearVariableFVReal + solver_sys = pressure_system + initial_condition = 0.2 + [] + [T] + type = MooseLinearVariableFVReal + solver_sys = energy_system + initial_condition = 300 + [] +[] + +[LinearFVKernels] + [u_advection_stress] + type = LinearWCNSFVMomentumFlux + variable = vel_x + advected_interp_method = ${advected_interp_method} + mu = ${mu} + u = vel_x + v = vel_y + momentum_component = 'x' + rhie_chow_user_object = 'rc' + use_nonorthogonal_correction = false + [] + [v_advection_stress] + type = LinearWCNSFVMomentumFlux + variable = vel_y + advected_interp_method = ${advected_interp_method} + mu = ${mu} + u = vel_x + v = vel_y + momentum_component = 'y' + rhie_chow_user_object = 'rc' + use_nonorthogonal_correction = false + [] + [u_pressure] + type = LinearFVMomentumPressure + variable = vel_x + pressure = pressure + momentum_component = 'x' + [] + [v_pressure] + type = LinearFVMomentumPressure + variable = vel_y + pressure = pressure + momentum_component = 'y' + [] + [u_boussinesq] + type = LinearFVMomentumBoussinesq + variable = vel_x + rho = ${rho} + gravity = '0 -9.81 0' + alpha_name = ${alpha_b} + ref_temperature = 300.0 + T_fluid = T + momentum_component = 'x' + [] + [v_boussinesq] + type = LinearFVMomentumBoussinesq + variable = vel_y + rho = ${rho} + gravity = '0 -9.81 0' + alpha_name = ${alpha_b} + ref_temperature = 300.0 + T_fluid = T + momentum_component = 'y' + [] + + [p_diffusion] + type = LinearFVAnisotropicDiffusion + variable = pressure + diffusion_tensor = Ainv + use_nonorthogonal_correction = false + [] + [HbyA_divergence] + type = LinearFVDivergence + variable = pressure + face_flux = HbyA + force_boundary_execution = true + [] + + [h_advection] + type = LinearFVEnergyAdvection + variable = T + cp = ${cp} + advected_interp_method = ${advected_interp_method} + rhie_chow_user_object = 'rc' + [] + [conduction] + type = LinearFVDiffusion + variable = T + diffusion_coeff = ${k} + use_nonorthogonal_correction = false + [] +[] + +[FunctorMaterials] + [constant_functors] + type = GenericFunctorMaterial + prop_names = 'cp alpha_b' + prop_values = '${cp} ${alpha_b}' + [] +[] + +[LinearFVBCs] + [inlet-u] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'left' + variable = vel_x + functor = '1.1' + [] + [inlet-v] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'left' + variable = vel_y + functor = '0.0' + [] + [walls-u] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'top bottom' + variable = vel_x + functor = 0.0 + [] + [walls-v] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'top bottom' + variable = vel_y + functor = 0.0 + [] + [outlet_p] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'right' + variable = pressure + functor = 1.4 + [] + [outlet_u] + type = LinearFVAdvectionDiffusionOutflowBC + variable = vel_x + use_two_term_expansion = false + boundary = right + [] + [outlet_v] + type = LinearFVAdvectionDiffusionOutflowBC + variable = vel_y + use_two_term_expansion = false + boundary = right + [] + [inlet_top_T] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = T + functor = ${fparse 300.0} + boundary = 'left top' + [] + [bottom_T] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = T + functor = ${fparse 350.0} + boundary = bottom + [] + [outlet_T] + type = LinearFVAdvectionDiffusionOutflowBC + variable = T + use_two_term_expansion = false + boundary = right + [] +[] + +[Executioner] + type = SIMPLE + momentum_l_abs_tol = 1e-10 + pressure_l_abs_tol = 1e-10 + energy_l_abs_tol = 1e-10 + momentum_l_tol = 0 + pressure_l_tol = 0 + energy_l_tol = 0 + rhie_chow_user_object = 'rc' + momentum_systems = 'u_system v_system' + pressure_system = 'pressure_system' + energy_system = 'energy_system' + momentum_equation_relaxation = 0.8 + energy_equation_relaxation = 0.9 + pressure_variable_relaxation = 0.3 + num_iterations = 200 + pressure_absolute_tolerance = 1e-10 + momentum_absolute_tolerance = 1e-10 + energy_absolute_tolerance = 1e-10 + momentum_petsc_options_iname = '-pc_type -pc_hypre_type' + momentum_petsc_options_value = 'hypre boomeramg' + pressure_petsc_options_iname = '-pc_type -pc_hypre_type' + pressure_petsc_options_value = 'hypre boomeramg' + energy_petsc_options_iname = '-pc_type -pc_hypre_type' + energy_petsc_options_value = 'hypre boomeramg' + print_fields = false +[] + +[Outputs] + exodus = true +[] diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/gold/2d-boussinesq_out.e b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/gold/2d-boussinesq_out.e new file mode 100644 index 000000000000..f12a320c14b4 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/gold/2d-boussinesq_out.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/tests b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/tests index 912c2b751882..261269687237 100644 --- a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/tests +++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/tests @@ -1,6 +1,6 @@ [Tests] - issues = '#27280' - design = 'SIMPLE.md LinearFVDivergence.md LinearWCNSFVMomentumFlux.md LinearFVMomentumPressure.md' + issues = '#27280 #28951' + design = 'SIMPLE.md LinearFVDivergence.md LinearWCNSFVMomentumFlux.md LinearFVMomentumPressure.md LinearFVEnergyAdvection.md LinearFVMomentumBoussinesq.md' [momentum-pressure] type = 'Exodiff' input = 2d-velocity-pressure.i @@ -10,4 +10,14 @@ recover = false # we don't support recovery for SIMPLE yet max_threads = 1 # see libmesh issue #3808 [] + [energy_boussinesq] + type = 'Exodiff' + input = 2d-boussinesq.i + exodiff = 2d-boussinesq_out.e + requirement = "The system shall be able to solve the steady-state Navier-Stokes equations with coupled fluid energy equations in a 2D channel using the Boussinesq approximation for buoyancy with the SIMPLE algorithm using the linear finite volume system." + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + abs_zero = 1e-5 + rel_err = 1e-5 + [] [] diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_segregated.i b/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_segregated.i new file mode 100644 index 000000000000..cfc2104c0b8c --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_segregated.i @@ -0,0 +1,250 @@ +################################################################################ +# MATERIAL PROPERTIES +################################################################################ +rho = 3279. +T_0 = 875.0 +mu = 1. +k_cond = 38.0 +cp = 640. +alpha = 3.26e-4 + +walls = 'right left top bottom' + +[GlobalParams] + rhie_chow_user_object = 'ins_rhie_chow_interpolator' + advected_interp_method = 'upwind' + u = superficial_vel_x + v = superficial_vel_y +[] + +[Problem] + linear_sys_names = 'u_system v_system pressure_system energy_system' + previous_nl_solution_required = true +[] + +################################################################################ +# GEOMETRY +################################################################################ + +[Mesh] + [gen] + type = GeneratedMeshGenerator + dim = 2 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + nx = 30 + ny = 30 + [] +[] + +################################################################################ +# EQUATIONS: VARIABLES, KERNELS & BCS +################################################################################ + +[UserObjects] + [ins_rhie_chow_interpolator] + type = RhieChowMassFlux + u = superficial_vel_x + v = superficial_vel_y + pressure = pressure + rho = ${rho} + p_diffusion_kernel = p_diffusion + [] +[] + +[Variables] + [superficial_vel_x] + type = MooseLinearVariableFVReal + solver_sys = u_system + [] + [superficial_vel_y] + type = MooseLinearVariableFVReal + solver_sys = v_system + [] + [pressure] + type = MooseLinearVariableFVReal + initial_condition = 0 + solver_sys = pressure_system + [] + [T_fluid] + type = MooseLinearVariableFVReal + solver_sys = energy_system + initial_condition = 875 + [] +[] + +[LinearFVKernels] + [u_advection_stress] + type = LinearWCNSFVMomentumFlux + variable = superficial_vel_x + mu = ${mu} + momentum_component = 'x' + use_nonorthogonal_correction = true + [] + [u_pressure] + type = LinearFVMomentumPressure + variable = superficial_vel_x + pressure = pressure + momentum_component = 'x' + [] + [u_buoyancy] + type = LinearFVMomentumBoussinesq + variable = superficial_vel_x + T_fluid = T_fluid + gravity = '0 -9.81 0' + rho = ${rho} + ref_temperature = ${T_0} + alpha_name = ${alpha} + momentum_component = 'x' + [] + + [v_advection_stress] + type = LinearWCNSFVMomentumFlux + variable = superficial_vel_y + mu = ${mu} + momentum_component = 'y' + use_nonorthogonal_correction = true + [] + [v_pressure] + type = LinearFVMomentumPressure + variable = superficial_vel_y + pressure = pressure + momentum_component = 'y' + [] + [v_buoyancy] + type = LinearFVMomentumBoussinesq + variable = superficial_vel_y + T_fluid = T_fluid + gravity = '0 -9.81 0' + rho = ${rho} + ref_temperature = ${T_0} + alpha_name = ${alpha} + momentum_component = 'y' + [] + + [p_diffusion] + type = LinearFVAnisotropicDiffusion + variable = pressure + diffusion_tensor = Ainv + use_nonorthogonal_correction = true + [] + [HbyA_divergence] + type = LinearFVDivergence + variable = pressure + face_flux = HbyA + force_boundary_execution = true + [] + + ####### FUEL ENERGY EQUATION ####### + + [heat_advection] + type = LinearFVEnergyAdvection + variable = T_fluid + cp = ${cp} + [] + [conduction] + type = LinearFVDiffusion + variable = T_fluid + diffusion_coeff = ${fparse k_cond} + [] +[] + + +[LinearFVBCs] + [no-slip-u] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = superficial_vel_x + boundary = ${walls} + functor = 0 + [] + [no-slip-v] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = superficial_vel_y + boundary = ${walls} + functor = 0 + [] + [T_cold] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = T_fluid + boundary = 'right' + functor = 870.0 + [] + [T_hot] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = T_fluid + boundary = 'left' + functor = 880.0 + [] + [T_all] + type = LinearFVAdvectionDiffusionExtrapolatedBC + variable = T_fluid + boundary = 'top bottom' + use_two_term_expansion = false + [] + [pressure-extrapolation] + type = LinearFVExtrapolatedPressureBC + boundary = ${walls} + variable = pressure + use_two_term_expansion = false + [] +[] + +[FunctorMaterials] + [constant_functors] + type = GenericFunctorMaterial + prop_names = 'cp alpha_b' + prop_values = '${cp} ${alpha}' + [] +[] + +################################################################################ +# EXECUTION / SOLVE +################################################################################ + +[Executioner] + type = SIMPLE + momentum_l_abs_tol = 1e-11 + pressure_l_abs_tol = 1e-11 + energy_l_abs_tol = 1e-11 + momentum_l_tol = 0 + pressure_l_tol = 0 + energy_l_tol = 0 + rhie_chow_user_object = 'ins_rhie_chow_interpolator' + momentum_systems = 'u_system v_system' + pressure_system = 'pressure_system' + energy_system = 'energy_system' + momentum_equation_relaxation = 0.3 + pressure_variable_relaxation = 0.7 + energy_equation_relaxation = 0.95 + num_iterations = 3000 + pressure_absolute_tolerance = 1e-8 + momentum_absolute_tolerance = 1e-8 + energy_absolute_tolerance = 1e-8 + print_fields = false + momentum_l_max_its = 300 + + pin_pressure = true + pressure_pin_value = 0.0 + pressure_pin_point = '0.5 0.0 0.0' + + # momentum_petsc_options = '-ksp_monitor' + momentum_petsc_options_iname = '-pc_type -pc_hypre_type' + momentum_petsc_options_value = 'hypre boomeramg' + + pressure_petsc_options_iname = '-pc_type -pc_hypre_type' + pressure_petsc_options_value = 'hypre boomeramg' + + energy_petsc_options_iname = '-pc_type -pc_hypre_type' + energy_petsc_options_value = 'hypre boomeramg' +[] + +################################################################################ +# SIMULATION OUTPUTS +################################################################################ + +[Outputs] + exodus = true +[] + diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/gold/diff_heated_cavity_linear_segregated_out.e b/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/gold/diff_heated_cavity_linear_segregated_out.e new file mode 100644 index 000000000000..f4a0ac45798f Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/gold/diff_heated_cavity_linear_segregated_out.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/tests b/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/tests new file mode 100644 index 000000000000..21c43ef37d68 --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/tests @@ -0,0 +1,13 @@ +[Tests] + design = 'LinearFVEnergyAdvection.md LinearFVMomentumBoussinesq.md' + issues = '#28951' + [diff_heated_cavity_linear_segregated] + type = 'Exodiff' + input = 'diff_heated_cavity_linear_segregated.i' + exodiff = 'diff_heated_cavity_linear_segregated_out.e' + abs_zero = 1e-5 + rel_err = 1e-3 + max_threads = 1 + requirement = 'The system shall be able to use the density Boussinesq approximation to solve for a differentially heated 2D cavity.' + [] +[]