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 "INSFVMixingLengthReynoldsStress.h"
      11             : #include "INSFVVelocityVariable.h"
      12             : #include "NS.h"
      13             : #include "SystemBase.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", INSFVMixingLengthReynoldsStress);
      16             : 
      17             : InputParameters
      18       17057 : INSFVMixingLengthReynoldsStress::validParams()
      19             : {
      20       17057 :   InputParameters params = INSFVFluxKernel::validParams();
      21       17057 :   params.addClassDescription(
      22             :       "Computes the force due to the Reynolds stress term in the incompressible"
      23             :       " Reynolds-averaged Navier-Stokes equations.");
      24       34114 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      25       34114 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      26       34114 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      27       17057 :   params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
      28       34114 :   params.addRequiredParam<MooseFunctorName>("mixing_length", "Turbulent eddy mixing length.");
      29       34114 :   MooseEnum momentum_component("x=0 y=1 z=2");
      30       34114 :   params.addRequiredParam<MooseEnum>(
      31             :       "momentum_component",
      32             :       momentum_component,
      33             :       "The component of the momentum equation that this kernel applies to.");
      34             :   // We assume the worst, e.g. we are doing Rhie-Chow. In that case we need three layers. An 'a'
      35             :   // coefficient evaluation at a face will necessitate evaluation of *this* object at every face of
      36             :   // the adjoining element, necessitating a face gradient evaluation at those faces, necessitating a
      37             :   // cell gradient evaluation in neighboring elements, necessitating cell value evaluations in
      38             :   // neighbors of those neighbor elements
      39       17057 :   params.set<unsigned short>("ghost_layers") = 3;
      40       17057 :   return params;
      41       17057 : }
      42             : 
      43         212 : INSFVMixingLengthReynoldsStress::INSFVMixingLengthReynoldsStress(const InputParameters & params)
      44             :   : INSFVFluxKernel(params),
      45         424 :     _dim(blocksMaxDimension()),
      46         424 :     _axis_index(getParam<MooseEnum>("momentum_component")),
      47         424 :     _u(getFunctor<ADReal>("u")),
      48         636 :     _v(params.isParamValid("v") ? &getFunctor<ADReal>("v") : nullptr),
      49         212 :     _w(params.isParamValid("w") ? &getFunctor<ADReal>("w") : nullptr),
      50         212 :     _rho(getFunctor<ADReal>(NS::density)),
      51         636 :     _mixing_len(getFunctor<ADReal>("mixing_length"))
      52             : {
      53         212 :   if (_dim >= 2 && !_v)
      54           0 :     mooseError(
      55             :         "In two or more dimensions, the v velocity must be supplied using the 'v' parameter");
      56         212 :   if (_dim >= 3 && !_w)
      57           0 :     mooseError("In three dimensions, the w velocity must be supplied using the 'w' parameter");
      58         212 : }
      59             : 
      60             : ADReal
      61     5175768 : INSFVMixingLengthReynoldsStress::computeStrongResidual(const bool populate_a_coeffs)
      62             : {
      63             :   using std::sqrt;
      64             : 
      65             :   constexpr Real offset = 1e-15; // prevents explosion of sqrt(x) derivative to infinity
      66             : 
      67     5175768 :   const auto face = makeCDFace(*_face_info);
      68     5175768 :   const auto state = determineState();
      69             : 
      70     5175768 :   const auto grad_u = _u.gradient(face, state);
      71             :   // Compute the dot product of the strain rate tensor and the normal vector
      72             :   // aka (grad_v + grad_v^T) * n_hat
      73     5175768 :   ADReal norm_strain_rate = grad_u(_axis_index) * _normal(0);
      74             :   ADRealVectorValue grad_v;
      75             :   ADRealVectorValue grad_w;
      76     5175768 :   if (_dim >= 2)
      77             :   {
      78    10351536 :     grad_v = _v->gradient(face, state);
      79    10351536 :     norm_strain_rate += grad_v(_axis_index) * _normal(1);
      80     5175768 :     if (_dim >= 3)
      81             :     {
      82           0 :       grad_w = _w->gradient(face, state);
      83           0 :       norm_strain_rate += grad_w(_axis_index) * _normal(2);
      84             :     }
      85             :   }
      86     5175768 :   const ADRealVectorValue & var_grad = _index == 0 ? grad_u : (_index == 1 ? grad_v : grad_w);
      87     5175768 :   norm_strain_rate += var_grad * _normal;
      88             : 
      89    10351536 :   ADReal symmetric_strain_tensor_norm = 2.0 * Utility::pow<2>(grad_u(0));
      90     5175768 :   if (_dim >= 2)
      91             :   {
      92             :     symmetric_strain_tensor_norm +=
      93    10351536 :         2.0 * Utility::pow<2>(grad_v(1)) + Utility::pow<2>(grad_v(0) + grad_u(1));
      94     5175768 :     if (_dim >= 3)
      95           0 :       symmetric_strain_tensor_norm += 2.0 * Utility::pow<2>(grad_w(2)) +
      96           0 :                                       Utility::pow<2>(grad_u(2) + grad_w(0)) +
      97           0 :                                       Utility::pow<2>(grad_v(2) + grad_w(1));
      98             :   }
      99             : 
     100     5175768 :   symmetric_strain_tensor_norm = sqrt(symmetric_strain_tensor_norm + offset);
     101             : 
     102             :   // Interpolate the mixing length to the face
     103     5175768 :   const ADReal mixing_len = _mixing_len(face, state);
     104             : 
     105             :   // Compute the eddy diffusivity
     106     5175768 :   ADReal eddy_diff = symmetric_strain_tensor_norm * mixing_len * mixing_len;
     107             : 
     108     5175768 :   const ADReal rho = _rho(face, state);
     109             : 
     110     5175768 :   if (populate_a_coeffs)
     111             :   {
     112     5175768 :     if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     113             :         _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     114             :     {
     115     5175768 :       const auto dof_number = _face_info->elem().dof_number(_sys.number(), _var.number(), 0);
     116             :       // norm_strain_rate is a linear combination of degrees of freedom so it's safe to straight-up
     117             :       // index into the derivatives vector at the dof we care about
     118     5175768 :       _ae = norm_strain_rate.derivatives()[dof_number];
     119    10351536 :       _ae *= -rho * eddy_diff;
     120             :     }
     121     5175768 :     if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
     122             :         _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     123             :     {
     124     5016520 :       const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
     125     5016520 :       _an = norm_strain_rate.derivatives()[dof_number];
     126     5016520 :       _an *= rho * eddy_diff;
     127             :     }
     128             :   }
     129             : 
     130             :   // Return the turbulent stress contribution to the momentum equation
     131    10351536 :   return -1 * rho * eddy_diff * norm_strain_rate;
     132             : }
     133             : 
     134             : ADReal
     135           0 : INSFVMixingLengthReynoldsStress::computeSegregatedContribution()
     136             : {
     137           0 :   return computeStrongResidual(false);
     138             : }
     139             : 
     140             : void
     141     5623128 : INSFVMixingLengthReynoldsStress::gatherRCData(const FaceInfo & fi)
     142             : {
     143     5623128 :   if (skipForBoundary(fi))
     144             :     return;
     145             : 
     146     5175768 :   _face_info = &fi;
     147     5175768 :   _normal = fi.normal();
     148     5175768 :   _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
     149             : 
     150    10351536 :   addResidualAndJacobian(computeStrongResidual(true) * (fi.faceArea() * fi.faceCoord()));
     151             : 
     152     5175768 :   if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     153             :       _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     154    10351536 :     _rc_uo.addToA(&fi.elem(), _index, _ae * (fi.faceArea() * fi.faceCoord()));
     155     5175768 :   if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
     156             :       _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     157    15049560 :     _rc_uo.addToA(fi.neighborPtr(), _index, _an * (fi.faceArea() * fi.faceCoord()));
     158             : }