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 "INSFVTKEDWallFunctionBC.h"
      11             : #include "Function.h"
      12             : #include "NavierStokesMethods.h"
      13             : #include "NS.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", INSFVTKEDWallFunctionBC);
      16             : 
      17             : InputParameters
      18          82 : INSFVTKEDWallFunctionBC::validParams()
      19             : {
      20          82 :   InputParameters params = FVDirichletBCBase::validParams();
      21          82 :   params.addClassDescription("Adds Reichardt extrapolated wall values to set up directly the"
      22             :                              "Dirichlet BC for the turbulent kinetic energy dissipation rate.");
      23         164 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      24         164 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      25         164 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      26          82 :   params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
      27          82 :   params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
      28          82 :   params.addRequiredParam<MooseFunctorName>(NS::mu_t, "The turbulent viscosity.");
      29         164 :   params.addRequiredParam<MooseFunctorName>("k", "The turbulent kinetic energy.");
      30         164 :   params.deprecateParam("k", NS::TKE, "01/01/2025");
      31         164 :   params.addParam<MooseFunctorName>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
      32         164 :   params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
      33         164 :   params.addParamNamesToGroup("newton_solve", "Advanced");
      34          82 :   return params;
      35           0 : }
      36             : 
      37          44 : INSFVTKEDWallFunctionBC::INSFVTKEDWallFunctionBC(const InputParameters & params)
      38             :   : FVDirichletBCBase(params),
      39          44 :     _dim(_subproblem.mesh().dimension()),
      40          88 :     _u_var(getFunctor<ADReal>("u")),
      41          44 :     _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
      42          44 :     _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
      43          44 :     _rho(getFunctor<ADReal>(NS::density)),
      44          44 :     _mu(getFunctor<ADReal>(NS::mu)),
      45          44 :     _mu_t(getFunctor<ADReal>(NS::mu_t)),
      46          44 :     _k(getFunctor<ADReal>(NS::TKE)),
      47          88 :     _C_mu(getFunctor<ADReal>("C_mu")),
      48         132 :     _newton_solve(getParam<bool>("newton_solve"))
      49             : {
      50          44 : }
      51             : 
      52             : ADReal
      53       22600 : INSFVTKEDWallFunctionBC::boundaryValue(const FaceInfo & fi, const Moose::StateArg & state) const
      54             : {
      55       22600 :   const bool use_elem = (fi.faceType(std::make_pair(_var.number(), _var.sys().number())) ==
      56             :                          FaceInfo::VarFaceNeighbors::ELEM);
      57             :   const Real dist = std::abs(
      58       22600 :       ((use_elem ? fi.elemCentroid() : fi.neighborCentroid()) - fi.faceCentroid()) * fi.normal());
      59       22600 :   const Elem * const elem_ptr = use_elem ? fi.elemPtr() : fi.neighborPtr();
      60       22600 :   const auto elem_arg = makeElemArg(elem_ptr);
      61       22600 :   const auto mu = _mu(elem_arg, state);
      62       22600 :   const auto rho = _rho(elem_arg, state);
      63             : 
      64             :   // Assign boundary weights to element
      65             :   // This is based on the theory of linear TKE development for each boundary
      66             :   // This is, it assumes no interaction across turbulence production from boundaries
      67       22600 :   Real weight = 0.0;
      68      113000 :   for (unsigned int i_side = 0; i_side < elem_ptr->n_sides(); ++i_side)
      69       90400 :     weight += static_cast<Real>(_subproblem.mesh().getBoundaryIDs(elem_ptr, i_side).size());
      70             : 
      71             :   // Get the velocity vector
      72       45200 :   ADRealVectorValue velocity(_u_var(elem_arg, state));
      73       22600 :   if (_v_var)
      74           0 :     velocity(1) = (*_v_var)(elem_arg, state);
      75       22600 :   if (_w_var)
      76           0 :     velocity(2) = (*_w_var)(elem_arg, state);
      77             : 
      78             :   // Compute the velocity and direction of the velocity component that is parallel to the wall
      79             :   const ADReal parallel_speed =
      80       67800 :       NS::computeSpeed<ADReal>(velocity - velocity * (fi.normal()) * (fi.normal()));
      81             : 
      82             :   // Get friction velocity
      83       22600 :   const ADReal u_star = NS::findUStar<ADReal>(mu, rho, parallel_speed, dist);
      84             : 
      85             :   // Get associated non-dimensional wall distance
      86       22600 :   const ADReal y_plus = dist * u_star * rho / mu;
      87             : 
      88       22600 :   const auto TKE = _k(elem_arg, state);
      89             : 
      90       22600 :   if (y_plus <= 5.0) // sub-laminar layer
      91             :   {
      92       19222 :     const auto laminar_value = 2.0 * weight * TKE * mu / std::pow(dist, 2);
      93       19222 :     if (!_newton_solve)
      94       19222 :       return laminar_value;
      95             :     else
      96             :       // Additional zero term to make sure new derivatives are not introduced as y_plus
      97             :       // changes
      98           0 :       return laminar_value + 0 * _mu_t(elem_arg, state);
      99             :   }
     100        3378 :   else if (y_plus >= 30.0) // log-layer
     101             :   {
     102         248 :     const auto turbulent_value = weight * _C_mu(elem_arg, state) * std::pow(std::abs(TKE), 1.5) /
     103         248 :                                  (_mu_t(elem_arg, state) * dist);
     104         248 :     if (!_newton_solve)
     105         248 :       return turbulent_value;
     106             :     else
     107             :       // Additional zero term to make sure new derivatives are not introduced as y_plus changes
     108           0 :       return turbulent_value + 0 * mu;
     109             :   }
     110             :   else // blending function
     111             :   {
     112        3130 :     const auto laminar_value = 2.0 * weight * TKE * mu / std::pow(dist, 2);
     113        3130 :     const auto turbulent_value = weight * _C_mu(elem_arg, state) * std::pow(std::abs(TKE), 1.5) /
     114        3130 :                                  (_mu_t(elem_arg, state) * dist);
     115        6260 :     const auto interpolation_coef = (y_plus - 5.0) / 25.0;
     116        3130 :     return (interpolation_coef * (turbulent_value - laminar_value) + laminar_value);
     117             :   }
     118             : }