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       16544 : INSFVMixingLengthReynoldsStress::validParams()
      19             : {
      20       16544 :   InputParameters params = INSFVFluxKernel::validParams();
      21       16544 :   params.addClassDescription(
      22             :       "Computes the force due to the Reynolds stress term in the incompressible"
      23             :       " Reynolds-averaged Navier-Stokes equations.");
      24       33088 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      25       33088 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      26       33088 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      27       16544 :   params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
      28       33088 :   params.addRequiredParam<MooseFunctorName>("mixing_length", "Turbulent eddy mixing length.");
      29       33088 :   MooseEnum momentum_component("x=0 y=1 z=2");
      30       33088 :   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       16544 :   params.set<unsigned short>("ghost_layers") = 3;
      40       16544 :   return params;
      41       16544 : }
      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             :   constexpr Real offset = 1e-15; // prevents explosion of sqrt(x) derivative to infinity
      64             : 
      65     5175768 :   const auto face = makeCDFace(*_face_info);
      66     5175768 :   const auto state = determineState();
      67             : 
      68     5175768 :   const auto grad_u = _u.gradient(face, state);
      69             :   // Compute the dot product of the strain rate tensor and the normal vector
      70             :   // aka (grad_v + grad_v^T) * n_hat
      71     5175768 :   ADReal norm_strain_rate = grad_u(_axis_index) * _normal(0);
      72             :   ADRealVectorValue grad_v;
      73             :   ADRealVectorValue grad_w;
      74     5175768 :   if (_dim >= 2)
      75             :   {
      76    10351536 :     grad_v = _v->gradient(face, state);
      77    10351536 :     norm_strain_rate += grad_v(_axis_index) * _normal(1);
      78     5175768 :     if (_dim >= 3)
      79             :     {
      80           0 :       grad_w = _w->gradient(face, state);
      81           0 :       norm_strain_rate += grad_w(_axis_index) * _normal(2);
      82             :     }
      83             :   }
      84     5175768 :   const ADRealVectorValue & var_grad = _index == 0 ? grad_u : (_index == 1 ? grad_v : grad_w);
      85     5175768 :   norm_strain_rate += var_grad * _normal;
      86             : 
      87    10351536 :   ADReal symmetric_strain_tensor_norm = 2.0 * Utility::pow<2>(grad_u(0));
      88     5175768 :   if (_dim >= 2)
      89             :   {
      90             :     symmetric_strain_tensor_norm +=
      91    10351536 :         2.0 * Utility::pow<2>(grad_v(1)) + Utility::pow<2>(grad_v(0) + grad_u(1));
      92     5175768 :     if (_dim >= 3)
      93           0 :       symmetric_strain_tensor_norm += 2.0 * Utility::pow<2>(grad_w(2)) +
      94           0 :                                       Utility::pow<2>(grad_u(2) + grad_w(0)) +
      95           0 :                                       Utility::pow<2>(grad_v(2) + grad_w(1));
      96             :   }
      97             : 
      98     5175768 :   symmetric_strain_tensor_norm = std::sqrt(symmetric_strain_tensor_norm + offset);
      99             : 
     100             :   // Interpolate the mixing length to the face
     101     5175768 :   const ADReal mixing_len = _mixing_len(face, state);
     102             : 
     103             :   // Compute the eddy diffusivity
     104     5175768 :   ADReal eddy_diff = symmetric_strain_tensor_norm * mixing_len * mixing_len;
     105             : 
     106     5175768 :   const ADReal rho = _rho(face, state);
     107             : 
     108     5175768 :   if (populate_a_coeffs)
     109             :   {
     110     5175768 :     if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     111             :         _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     112             :     {
     113     5175768 :       const auto dof_number = _face_info->elem().dof_number(_sys.number(), _var.number(), 0);
     114             :       // norm_strain_rate is a linear combination of degrees of freedom so it's safe to straight-up
     115             :       // index into the derivatives vector at the dof we care about
     116     5175768 :       _ae = norm_strain_rate.derivatives()[dof_number];
     117    10351536 :       _ae *= -rho * eddy_diff;
     118             :     }
     119     5175768 :     if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
     120             :         _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     121             :     {
     122     5016520 :       const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
     123     5016520 :       _an = norm_strain_rate.derivatives()[dof_number];
     124     5016520 :       _an *= rho * eddy_diff;
     125             :     }
     126             :   }
     127             : 
     128             :   // Return the turbulent stress contribution to the momentum equation
     129    10351536 :   return -1 * rho * eddy_diff * norm_strain_rate;
     130             : }
     131             : 
     132             : ADReal
     133           0 : INSFVMixingLengthReynoldsStress::computeSegregatedContribution()
     134             : {
     135           0 :   return computeStrongResidual(false);
     136             : }
     137             : 
     138             : void
     139     5623128 : INSFVMixingLengthReynoldsStress::gatherRCData(const FaceInfo & fi)
     140             : {
     141     5623128 :   if (skipForBoundary(fi))
     142             :     return;
     143             : 
     144     5175768 :   _face_info = &fi;
     145     5175768 :   _normal = fi.normal();
     146     5175768 :   _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
     147             : 
     148    10351536 :   addResidualAndJacobian(computeStrongResidual(true) * (fi.faceArea() * fi.faceCoord()));
     149             : 
     150     5175768 :   if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     151             :       _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     152    10351536 :     _rc_uo.addToA(&fi.elem(), _index, _ae * (fi.faceArea() * fi.faceCoord()));
     153     5175768 :   if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
     154             :       _face_type == FaceInfo::VarFaceNeighbors::BOTH)
     155    15049560 :     _rc_uo.addToA(fi.neighborPtr(), _index, _an * (fi.faceArea() * fi.faceCoord()));
     156             : }