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