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 "WallFunctionWallShearStressAux.h"
      11             : #include "NavierStokesMethods.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", WallFunctionWallShearStressAux);
      14             : 
      15             : InputParameters
      16          11 : WallFunctionWallShearStressAux::validParams()
      17             : {
      18          11 :   InputParameters params = AuxKernel::validParams();
      19          11 :   params.addClassDescription(
      20             :       "Calculates the wall shear stress based on algebraic standard velocity wall functions.");
      21          22 :   params.addRequiredCoupledVar("u", "The velocity in the x direction.");
      22          22 :   params.addCoupledVar("v", "The velocity in the y direction.");
      23          22 :   params.addCoupledVar("w", "The velocity in the z direction.");
      24          22 :   params.addRequiredParam<MooseFunctorName>("rho", "fluid density");
      25          22 :   params.addRequiredParam<MooseFunctorName>("mu", "Dynamic viscosity");
      26          22 :   params.addRequiredParam<std::vector<BoundaryName>>("walls",
      27             :                                                      "Boundaries that correspond to solid walls");
      28          11 :   return params;
      29           0 : }
      30             : 
      31           6 : WallFunctionWallShearStressAux::WallFunctionWallShearStressAux(const InputParameters & params)
      32             :   : AuxKernel(params),
      33           6 :     _dim(_subproblem.mesh().dimension()),
      34           6 :     _u_var(dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("u", 0))),
      35           6 :     _v_var(params.isParamValid("v")
      36          12 :                ? dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("v", 0))
      37             :                : nullptr),
      38           6 :     _w_var(params.isParamValid("w")
      39           6 :                ? dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("w", 0))
      40             :                : nullptr),
      41          12 :     _rho(getFunctor<Real>("rho")),
      42          12 :     _mu(getFunctor<Real>("mu")),
      43          18 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls"))
      44             : {
      45           6 :   if (!_u_var)
      46           0 :     paramError("u", "the u velocity must be an INSFVVelocityVariable.");
      47             : 
      48           6 :   if (_dim >= 2 && !_v_var)
      49           0 :     paramError("v",
      50             :                "In two or more dimensions, the v velocity must be supplied and it must be an "
      51             :                "INSFVVelocityVariable.");
      52             : 
      53           6 :   if (_dim >= 3 && !params.isParamValid("w"))
      54           0 :     paramError("w",
      55             :                "In three-dimensions, the w velocity must be supplied and it must be an "
      56             :                "INSFVVelocityVariable.");
      57           6 : }
      58             : 
      59             : Real
      60      123200 : WallFunctionWallShearStressAux::computeValue()
      61             : {
      62      123200 :   const Elem & elem = *_current_elem;
      63             : 
      64             :   bool wall_bounded = false;
      65             :   Real min_wall_dist = 1e10;
      66             :   Point wall_vec;
      67             :   Point normal;
      68      616000 :   for (unsigned int i_side = 0; i_side < elem.n_sides(); ++i_side)
      69             :   {
      70             :     const std::vector<BoundaryID> side_bnds =
      71      492800 :         _subproblem.mesh().getBoundaryIDs(_current_elem, i_side);
      72      985600 :     for (const BoundaryName & name : _wall_boundary_names)
      73             :     {
      74      492800 :       BoundaryID wall_id = _subproblem.mesh().getBoundaryID(name);
      75      516432 :       for (BoundaryID side_id : side_bnds)
      76             :       {
      77       23632 :         if (side_id == wall_id)
      78             :         {
      79       11200 :           const FaceInfo * const fi = _mesh.faceInfo(&elem, i_side);
      80             :           const Point & this_normal = fi->normal();
      81       11200 :           Point this_wall_vec = (elem.vertex_average() - fi->faceCentroid());
      82             :           Real dist = std::abs(this_wall_vec * normal);
      83       11200 :           if (dist < min_wall_dist)
      84             :           {
      85             :             min_wall_dist = dist;
      86             :             wall_vec = this_wall_vec;
      87       11200 :             normal = this_normal;
      88             :           }
      89             :           wall_bounded = true;
      90             :         }
      91             :       }
      92             :     }
      93      492800 :   }
      94             : 
      95      123200 :   if (!wall_bounded)
      96             :     return 0;
      97             : 
      98       11200 :   const auto state = determineState();
      99             : 
     100             :   // Get the velocity vector
     101       11200 :   RealVectorValue velocity(MetaPhysicL::raw_value(_u_var->getElemValue(&elem, state)));
     102       11200 :   if (_v_var)
     103       11200 :     velocity(1) = MetaPhysicL::raw_value(_v_var->getElemValue(&elem, state));
     104       11200 :   if (_w_var)
     105           0 :     velocity(2) = MetaPhysicL::raw_value(_w_var->getElemValue(&elem, state));
     106             : 
     107             :   // Compute the velocity and direction of the velocity component that is parallel to the wall
     108             :   Real dist = std::abs(wall_vec * normal);
     109             :   Real perpendicular_speed = velocity * normal;
     110             :   RealVectorValue parallel_velocity = velocity - perpendicular_speed * normal;
     111       11200 :   Real parallel_speed = parallel_velocity.norm();
     112             : 
     113       11200 :   if (parallel_speed < 1e-6)
     114             :     return 0;
     115             : 
     116       11200 :   if (!std::isfinite(parallel_speed))
     117             :     return parallel_speed;
     118             : 
     119             :   // Compute the friction velocity and the wall shear stress
     120       11200 :   const auto elem_arg = makeElemArg(_current_elem);
     121       11200 :   const auto rho = _rho(elem_arg, state);
     122       11200 :   const auto mu = _mu(elem_arg, state);
     123       11200 :   const Real u_star = NS::findUStar<Real>(mu, rho, parallel_speed, dist);
     124       11200 :   const Real tau = u_star * u_star * rho;
     125             : 
     126       11200 :   return tau;
     127             : }