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 "INSFVTurbulentViscosityWallFunction.h"
      11             : #include "Function.h"
      12             : #include "NavierStokesMethods.h"
      13             : 
      14             : registerMooseObject("NavierStokesApp", INSFVTurbulentViscosityWallFunction);
      15             : 
      16             : InputParameters
      17        1916 : INSFVTurbulentViscosityWallFunction::validParams()
      18             : {
      19        1916 :   InputParameters params = FVDirichletBCBase::validParams();
      20        1916 :   params.addClassDescription("Adds Dirichlet BC for wall values of the turbulent viscosity.");
      21        3832 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      22        3832 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      23        3832 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      24        1916 :   params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
      25        1916 :   params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
      26        1916 :   params.addRequiredParam<MooseFunctorName>(NS::mu_t, "The turbulent viscosity.");
      27        3832 :   params.addParam<MooseFunctorName>("k", "The turbulent kinetic energy.");
      28        3832 :   params.deprecateParam("k", NS::TKE, "01/01/2025");
      29        3832 :   params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
      30             : 
      31        3832 :   MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
      32        3832 :   params.addParam<MooseEnum>(
      33             :       "wall_treatment", wall_treatment, "The method used for computing the wall functions");
      34        1916 :   return params;
      35        1916 : }
      36             : 
      37         441 : INSFVTurbulentViscosityWallFunction::INSFVTurbulentViscosityWallFunction(
      38         441 :     const InputParameters & params)
      39             :   : FVDirichletBCBase(params),
      40         441 :     _dim(_subproblem.mesh().dimension()),
      41         882 :     _u_var(getFunctor<ADReal>("u")),
      42        1323 :     _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
      43         441 :     _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
      44         441 :     _rho(getFunctor<ADReal>(NS::density)),
      45         441 :     _mu(getFunctor<ADReal>(NS::mu)),
      46         441 :     _mu_t(getFunctor<ADReal>(NS::mu_t)),
      47         441 :     _k(getFunctor<ADReal>(NS::TKE)),
      48         882 :     _C_mu(getParam<Real>("C_mu")),
      49         882 :     _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
      50         882 :     _preserve_sparsity_pattern(_fv_problem.preserveMatrixSparsityPattern())
      51             : {
      52         441 : }
      53             : 
      54             : ADReal
      55      389040 : INSFVTurbulentViscosityWallFunction::boundaryValue(const FaceInfo & fi,
      56             :                                                    const Moose::StateArg & /* state */) const
      57             : {
      58             :   // if on an internal face (internal to the mesh, but an external boundary of the flow area),
      59             :   // we have to select the element on which the flow variables / viscosity is defined
      60             :   const bool use_elem = (fi.faceType(_var_sys_numbers_pair) == FaceInfo::VarFaceNeighbors::ELEM);
      61             :   const Real wall_dist = std::abs(
      62      389040 :       ((use_elem ? fi.elemCentroid() : fi.neighborCentroid()) - fi.faceCentroid()) * fi.normal());
      63      389040 :   const Elem * const elem_ptr = use_elem ? fi.elemPtr() : fi.neighborPtr();
      64      389040 :   const auto elem_arg = makeElemArg(elem_ptr);
      65             : 
      66             :   // Get the previous non linear iteration values
      67             :   const auto old_state = Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear);
      68      389040 :   const auto mu = _mu(elem_arg, old_state);
      69      389040 :   const auto rho = _rho(elem_arg, old_state);
      70             : 
      71             :   // Get the velocity vector
      72      778080 :   ADRealVectorValue velocity(_u_var(elem_arg, old_state));
      73      389040 :   if (_v_var)
      74      389040 :     velocity(1) = (*_v_var)(elem_arg, old_state);
      75      389040 :   if (_w_var)
      76           0 :     velocity(2) = (*_w_var)(elem_arg, old_state);
      77             : 
      78             :   // Compute the velocity and direction of the velocity component that is parallel to the wall
      79             :   const auto parallel_speed =
      80     1167120 :       NS::computeSpeed<ADReal>(velocity - velocity * (fi.normal()) * (fi.normal()));
      81             : 
      82             :   // Switch for determining the near wall quantities
      83             :   // wall_treatment can be: "eq_newton eq_incremental eq_linearized neq"
      84             :   ADReal y_plus;
      85             :   ADReal mut_log; // turbulent log-layer viscosity
      86             :   ADReal mu_wall; // total wall viscosity to obtain the shear stress at the wall
      87             : 
      88      389040 :   if (_wall_treatment == NS::WallTreatmentEnum::EQ_NEWTON)
      89             :   {
      90             :     // Full Newton-Raphson solve to find the wall quantities from the law of the wall
      91      341160 :     const auto u_tau = NS::findUStar<ADReal>(mu, rho, parallel_speed, wall_dist);
      92      341160 :     y_plus = wall_dist * u_tau * rho / mu;
      93      341160 :     mu_wall = rho * Utility::pow<2>(u_tau) * wall_dist / parallel_speed;
      94      341160 :     mut_log = mu_wall - mu;
      95             :   }
      96       47880 :   else if (_wall_treatment == NS::WallTreatmentEnum::EQ_INCREMENTAL)
      97             :   {
      98             :     // Incremental solve on y_plus to get the near-wall quantities
      99       21960 :     y_plus = NS::findyPlus<ADReal>(mu, rho, std::max(parallel_speed, 1e-10), wall_dist);
     100           0 :     mu_wall = mu * (NS::von_karman_constant * y_plus /
     101       65880 :                     std::log(std::max(NS::E_turb_constant * y_plus, 1 + 1e-4)));
     102       21960 :     mut_log = mu_wall - mu;
     103             :   }
     104       25920 :   else if (_wall_treatment == NS::WallTreatmentEnum::EQ_LINEARIZED)
     105             :   {
     106             :     // Linearized approximation to the wall function to find the near-wall quantities faster
     107       12960 :     const ADReal a_c = 1 / NS::von_karman_constant;
     108       25920 :     const ADReal b_c = 1 / NS::von_karman_constant *
     109       38880 :                        (std::log(NS::E_turb_constant * std::max(wall_dist, 1.0) / mu) + 1.0);
     110       12960 :     const ADReal c_c = parallel_speed;
     111             : 
     112       64800 :     const auto u_tau = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c);
     113       12960 :     y_plus = wall_dist * u_tau * rho / mu;
     114       12960 :     mu_wall = rho * Utility::pow<2>(u_tau) * wall_dist / parallel_speed;
     115       12960 :     mut_log = mu_wall - mu;
     116             :   }
     117       12960 :   else if (_wall_treatment == NS::WallTreatmentEnum::NEQ)
     118             :   {
     119             :     // Assign non-equilibrium wall function value
     120       25920 :     y_plus = std::pow(_C_mu, 0.25) * wall_dist * std::sqrt(_k(elem_arg, old_state)) * rho / mu;
     121           0 :     mu_wall = mu * (NS::von_karman_constant * y_plus /
     122       38880 :                     std::log(std::max(NS::E_turb_constant * y_plus, 1.0 + 1e-4)));
     123       12960 :     mut_log = mu_wall - mu;
     124             :   }
     125             :   else
     126             :     mooseAssert(false,
     127             :                 "For `INSFVTurbulentViscosityWallFunction` , wall treatment should not reach here");
     128             : 
     129      389040 :   ADReal residual = 0;
     130             :   // To keep the same sparsity pattern for all y_plus
     131      389040 :   if (_preserve_sparsity_pattern)
     132      778080 :     residual = 0 * mut_log * y_plus;
     133             : 
     134      389040 :   if (y_plus <= 5.0)
     135             :     // sub-laminar layer
     136             :     residual += 0.0;
     137      135408 :   else if (y_plus >= 30.0)
     138             :     // log-layer
     139      119466 :     residual += std::max(mut_log, NS::mu_t_low_limit);
     140             :   else
     141             :   {
     142             :     // buffer layer
     143       15942 :     const auto blending_function = (y_plus - 5.0) / 25.0;
     144             :     // the blending depends on the mut_log at y+=30
     145       15942 :     const auto mut_log = mu * _mut_30;
     146       15942 :     residual += std::max(blending_function * mut_log, NS::mu_t_low_limit);
     147             :   }
     148      389040 :   return residual;
     149             : }