Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "LinearFVTurbulentViscosityWallFunctionBC.h"
      11             : #include "NavierStokesMethods.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", LinearFVTurbulentViscosityWallFunctionBC);
      14             : 
      15             : InputParameters
      16          88 : LinearFVTurbulentViscosityWallFunctionBC::validParams()
      17             : {
      18          88 :   InputParameters params = LinearFVAdvectionDiffusionBC::validParams();
      19          88 :   params.addClassDescription("Adds Dirichlet BC for wall values of the turbulent viscosity.");
      20         176 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      21         176 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      22         176 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      23          88 :   params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
      24          88 :   params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
      25         176 :   params.addParam<MooseFunctorName>("k", "The turbulent kinetic energy.");
      26         176 :   params.deprecateParam("k", NS::TKE, "01/01/2025");
      27         176 :   params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
      28             : 
      29         176 :   MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
      30         176 :   params.addParam<MooseEnum>(
      31             :       "wall_treatment", wall_treatment, "The method used for computing the wall functions");
      32          88 :   return params;
      33          88 : }
      34             : 
      35          44 : LinearFVTurbulentViscosityWallFunctionBC::LinearFVTurbulentViscosityWallFunctionBC(
      36          44 :     const InputParameters & params)
      37             :   : LinearFVAdvectionDiffusionBC(params),
      38          44 :     _dim(_subproblem.mesh().dimension()),
      39          88 :     _u_var(getFunctor<Real>("u")),
      40         132 :     _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
      41          44 :     _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
      42          44 :     _rho(getFunctor<Real>(NS::density)),
      43          44 :     _mu(getFunctor<Real>(NS::mu)),
      44          44 :     _k(getFunctor<Real>(NS::TKE)),
      45          88 :     _C_mu(getParam<Real>("C_mu")),
      46         132 :     _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>())
      47             : {
      48          44 : }
      49             : 
      50             : Real
      51     1022784 : LinearFVTurbulentViscosityWallFunctionBC::computeTurbulentViscosity() const
      52             : {
      53             :   // Utility functions
      54     1022784 :   const auto wall_dist = computeCellToFaceDistance();
      55     1022784 :   const auto & elem = _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR
      56     1022784 :                           ? _current_face_info->neighborPtr()
      57      966532 :                           : _current_face_info->elemPtr();
      58     1022784 :   const auto re = makeElemArg(elem);
      59             :   const auto old_state = Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear);
      60     1022784 :   const auto mu = _mu(re, old_state);
      61     1022784 :   const auto rho = _rho(re, old_state);
      62             : 
      63             :   // Get the velocity vector
      64     1022784 :   RealVectorValue velocity(_u_var(re, old_state));
      65     1022784 :   if (_v_var)
      66     1022784 :     velocity(1) = (*_v_var)(re, old_state);
      67     1022784 :   if (_w_var)
      68           0 :     velocity(2) = (*_w_var)(re, old_state);
      69             : 
      70             :   // Compute the velocity and direction of the velocity component that is parallel to the wall
      71     1022784 :   const auto parallel_speed = NS::computeSpeed<Real>(
      72     1022784 :       velocity - velocity * (_current_face_info->normal()) * (_current_face_info->normal()));
      73             : 
      74             :   // Switch for determining the near wall quantities
      75             :   // wall_treatment can be: "eq_newton eq_incremental eq_linearized neq"
      76             :   Real y_plus = 0.0;
      77             :   Real mu_wall = 0.0; // total wall viscosity to obtain the shear stress at the wall
      78             : 
      79     1022784 :   if (_wall_treatment == NS::WallTreatmentEnum::EQ_NEWTON)
      80             :   {
      81             :     // Full Newton-Raphson solve to find the wall quantities from the law of the wall
      82      651784 :     const auto u_tau = NS::findUStar<Real>(mu, rho, parallel_speed, wall_dist);
      83      651784 :     y_plus = wall_dist * u_tau * rho / mu;
      84      651784 :     mu_wall = rho * Utility::pow<2>(u_tau) * wall_dist / parallel_speed;
      85             :   }
      86      371000 :   else if (_wall_treatment == NS::WallTreatmentEnum::EQ_INCREMENTAL)
      87             :   {
      88             :     // Incremental solve on y_plus to get the near-wall quantities
      89       21072 :     y_plus = NS::findyPlus<Real>(mu, rho, std::max(parallel_speed, 1e-10), wall_dist);
      90       42000 :     mu_wall = mu * (NS::von_karman_constant * y_plus /
      91       21072 :                     std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4)));
      92             :   }
      93      350000 :   else if (_wall_treatment == NS::WallTreatmentEnum::EQ_LINEARIZED)
      94             :   {
      95             :     // Linearized approximation to the wall function to find the near-wall quantities faster
      96             :     const Real a_c = 1 / NS::von_karman_constant;
      97             :     const Real b_c = 1 / NS::von_karman_constant *
      98           0 :                      (std::log(NS::E_turb_constant * std::max(wall_dist, 1.0) / mu) + 1.0);
      99             :     const Real c_c = parallel_speed;
     100             : 
     101           0 :     const auto u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c);
     102           0 :     y_plus = wall_dist * u_tau * rho / mu;
     103           0 :     mu_wall = rho * Utility::pow<2>(u_tau) * wall_dist / parallel_speed;
     104             :   }
     105      350000 :   else if (_wall_treatment == NS::WallTreatmentEnum::NEQ)
     106             :   {
     107             :     // Assign non-equilibrium wall function value
     108      350000 :     y_plus = std::pow(_C_mu, 0.25) * wall_dist * std::sqrt(_k(re, old_state)) * rho / mu;
     109      700000 :     mu_wall = mu * (NS::von_karman_constant * y_plus /
     110      350000 :                     std::log(std::max(NS::E_turb_constant * y_plus, 1.0 + 1e-4)));
     111             :   }
     112             :   else
     113             :     mooseAssert(false,
     114             :                 "For `INSFVTurbulentViscosityWallFunction` , wall treatment should not reach here");
     115             : 
     116     1022784 :   const Real mut_log = mu_wall - mu; // turbulent log-layer viscosity
     117             : 
     118             :   Real mu_t = 0;
     119             : 
     120     1022784 :   if (y_plus <= 5.0)
     121             :     // sub-laminar layer
     122             :     mu_t += 0.0;
     123     1022586 :   else if (y_plus >= 30.0)
     124             :     // log-layer
     125      452534 :     mu_t += std::max(mut_log, NS::mu_t_low_limit);
     126             :   else
     127             :   {
     128             :     // buffer layer
     129      570052 :     const auto blending_function = (y_plus - 5.0) / 25.0;
     130             :     // the blending depends on the mut_log at y+=30
     131      570052 :     const auto mut_log = mu * _mut_30;
     132      570052 :     mu_t += std::max(blending_function * mut_log, NS::mu_t_low_limit);
     133             :   }
     134     1022784 :   return mu_t;
     135             : }
     136             : 
     137             : Real
     138     1022784 : LinearFVTurbulentViscosityWallFunctionBC::computeBoundaryValue() const
     139             : {
     140     1022784 :   return this->computeTurbulentViscosity();
     141             : }
     142             : 
     143             : Real
     144           0 : LinearFVTurbulentViscosityWallFunctionBC::computeBoundaryNormalGradient() const
     145             : {
     146           0 :   const auto elem_arg = makeElemArg(_current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     147           0 :                                         ? _current_face_info->elemPtr()
     148           0 :                                         : _current_face_info->neighborPtr());
     149           0 :   const Real distance = computeCellToFaceDistance();
     150           0 :   return (this->computeTurbulentViscosity() - raw_value(_var(elem_arg, determineState()))) /
     151           0 :          distance;
     152             : }
     153             : 
     154             : Real
     155           0 : LinearFVTurbulentViscosityWallFunctionBC::computeBoundaryValueMatrixContribution() const
     156             : {
     157           0 :   mooseError("We should not solve for the turbulent viscosity directly meaning that this should "
     158             :              "contribute to neither vector nor a right hand side.");
     159             : }
     160             : 
     161             : Real
     162           0 : LinearFVTurbulentViscosityWallFunctionBC::computeBoundaryValueRHSContribution() const
     163             : {
     164           0 :   mooseError("We should not solve for the turbulent viscosity directly meaning that this should "
     165             :              "contribute to neither vector nor a right hand side.");
     166             : }
     167             : 
     168             : Real
     169           0 : LinearFVTurbulentViscosityWallFunctionBC::computeBoundaryGradientMatrixContribution() const
     170             : {
     171           0 :   mooseError("We should not solve for the turbulent viscosity directly meaning that this should "
     172             :              "contribute to neither vector nor a right hand side.");
     173             : }
     174             : 
     175             : Real
     176           0 : LinearFVTurbulentViscosityWallFunctionBC::computeBoundaryGradientRHSContribution() const
     177             : {
     178           0 :   mooseError("We should not solve for the turbulent viscosity directly meaning that this should "
     179             :              "contribute to neither vector nor a right hand side.");
     180             : }