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