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             :   using std::pow, std::abs;
      55             : 
      56       22600 :   const bool use_elem = (fi.faceType(std::make_pair(_var.number(), _var.sys().number())) ==
      57             :                          FaceInfo::VarFaceNeighbors::ELEM);
      58             :   const Real dist = abs(
      59       22600 :       ((use_elem ? fi.elemCentroid() : fi.neighborCentroid()) - fi.faceCentroid()) * fi.normal());
      60       22600 :   const Elem * const elem_ptr = use_elem ? fi.elemPtr() : fi.neighborPtr();
      61       22600 :   const auto elem_arg = makeElemArg(elem_ptr);
      62       22600 :   const auto mu = _mu(elem_arg, state);
      63       22600 :   const auto rho = _rho(elem_arg, state);
      64             : 
      65             :   // Assign boundary weights to element
      66             :   // This is based on the theory of linear TKE development for each boundary
      67             :   // This is, it assumes no interaction across turbulence production from boundaries
      68       22600 :   Real weight = 0.0;
      69      113000 :   for (unsigned int i_side = 0; i_side < elem_ptr->n_sides(); ++i_side)
      70       90400 :     weight += static_cast<Real>(_subproblem.mesh().getBoundaryIDs(elem_ptr, i_side).size());
      71             : 
      72             :   // Get the velocity vector
      73       45200 :   ADRealVectorValue velocity(_u_var(elem_arg, state));
      74       22600 :   if (_v_var)
      75           0 :     velocity(1) = (*_v_var)(elem_arg, state);
      76       22600 :   if (_w_var)
      77           0 :     velocity(2) = (*_w_var)(elem_arg, state);
      78             : 
      79             :   // Compute the velocity and direction of the velocity component that is parallel to the wall
      80             :   const ADReal parallel_speed =
      81       67800 :       NS::computeSpeed<ADReal>(velocity - velocity * (fi.normal()) * (fi.normal()));
      82             : 
      83             :   // Get friction velocity
      84       22600 :   const ADReal u_star = NS::findUStar<ADReal>(mu, rho, parallel_speed, dist);
      85             : 
      86             :   // Get associated non-dimensional wall distance
      87       22600 :   const ADReal y_plus = dist * u_star * rho / mu;
      88             : 
      89       22600 :   const auto TKE = _k(elem_arg, state);
      90             : 
      91       22600 :   if (y_plus <= 5.0) // sub-laminar layer
      92             :   {
      93       19222 :     const auto laminar_value = 2.0 * weight * TKE * mu / pow(dist, 2);
      94       19222 :     if (!_newton_solve)
      95       19222 :       return laminar_value;
      96             :     else
      97             :       // Additional zero term to make sure new derivatives are not introduced as y_plus
      98             :       // changes
      99           0 :       return laminar_value + 0 * _mu_t(elem_arg, state);
     100             :   }
     101        3378 :   else if (y_plus >= 30.0) // log-layer
     102             :   {
     103             :     const auto turbulent_value =
     104         496 :         weight * _C_mu(elem_arg, state) * pow(abs(TKE), 1.5) / (_mu_t(elem_arg, state) * dist);
     105         248 :     if (!_newton_solve)
     106         248 :       return turbulent_value;
     107             :     else
     108             :       // Additional zero term to make sure new derivatives are not introduced as y_plus changes
     109           0 :       return turbulent_value + 0 * mu;
     110             :   }
     111             :   else // blending function
     112             :   {
     113        3130 :     const auto laminar_value = 2.0 * weight * TKE * mu / pow(dist, 2);
     114             :     const auto turbulent_value =
     115        6260 :         weight * _C_mu(elem_arg, state) * pow(abs(TKE), 1.5) / (_mu_t(elem_arg, state) * dist);
     116        6260 :     const auto interpolation_coef = (y_plus - 5.0) / 25.0;
     117        3130 :     return (interpolation_coef * (turbulent_value - laminar_value) + laminar_value);
     118             :   }
     119             : }