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 "INSFVMixingLengthScalarDiffusion.h"
      11             : #include "MathFVUtils.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", INSFVMixingLengthScalarDiffusion);
      14             : 
      15             : InputParameters
      16          54 : INSFVMixingLengthScalarDiffusion::validParams()
      17             : {
      18          54 :   InputParameters params = FVFluxKernel::validParams();
      19          54 :   params += FVDiffusionInterpolationInterface::validParams();
      20          54 :   params.addClassDescription("Computes the turbulent diffusive flux that appears in "
      21             :                              "Reynolds-averaged fluid conservation equations.");
      22         108 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      23         108 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      24         108 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      25         108 :   params.addRequiredParam<MooseFunctorName>("mixing_length", "The turbulent mixing length.");
      26         108 :   params.addRequiredParam<Real>(
      27             :       "schmidt_number",
      28             :       "The turbulent Schmidt number (or turbulent Prandtl number if the passive scalar is energy) "
      29             :       "that relates the turbulent scalar diffusivity to the turbulent momentum diffusivity.");
      30          54 :   params.set<unsigned short>("ghost_layers") = 2;
      31             : 
      32             :   // We add the relationship manager here, this will select the right number of
      33             :   // ghosting layers depending on the chosen interpolation method
      34         108 :   params.addRelationshipManager(
      35             :       "ElementSideNeighborLayers",
      36             :       Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
      37             :           Moose::RelationshipManagerType::COUPLING,
      38             :       [](const InputParameters & obj_params, InputParameters & rm_params)
      39          33 :       { FVRelationshipManagerInterface::setRMParamsDiffusion(obj_params, rm_params, 3); });
      40             : 
      41          54 :   return params;
      42           0 : }
      43             : 
      44          29 : INSFVMixingLengthScalarDiffusion::INSFVMixingLengthScalarDiffusion(const InputParameters & params)
      45             :   : FVFluxKernel(params),
      46             :     FVDiffusionInterpolationInterface(params),
      47          29 :     _dim(_subproblem.mesh().dimension()),
      48          58 :     _u(getFunctor<ADReal>("u")),
      49         116 :     _v(isParamValid("v") ? &getFunctor<ADReal>("v") : nullptr),
      50          58 :     _w(isParamValid("w") ? &getFunctor<ADReal>("w") : nullptr),
      51          58 :     _mixing_len(getFunctor<ADReal>("mixing_length")),
      52          87 :     _schmidt_number(getParam<Real>("schmidt_number"))
      53             : {
      54          29 :   if (_dim >= 2 && !_v)
      55           0 :     mooseError(
      56             :         "In two or more dimensions, the v velocity must be supplied using the 'v' parameter");
      57          29 :   if (_dim >= 3 && !_w)
      58           0 :     mooseError("In three dimensions, the w velocity must be supplied using the 'w' parameter");
      59          29 : }
      60             : 
      61             : ADReal
      62      481338 : INSFVMixingLengthScalarDiffusion::computeQpResidual()
      63             : {
      64             :   using std::sqrt;
      65             : 
      66             :   constexpr Real offset = 1e-15; // prevents explosion of sqrt(x) derivative to infinity
      67             : 
      68      481338 :   auto face = makeCDFace(*_face_info);
      69      481338 :   const auto state = determineState();
      70             : 
      71      481338 :   const auto grad_u = _u.gradient(face, state);
      72      962676 :   ADReal symmetric_strain_tensor_norm = 2.0 * Utility::pow<2>(grad_u(0));
      73      481338 :   if (_dim >= 2)
      74             :   {
      75      481338 :     const auto grad_v = _v->gradient(face, state);
      76             :     symmetric_strain_tensor_norm +=
      77      962676 :         2.0 * Utility::pow<2>(grad_v(1)) + Utility::pow<2>(grad_v(0) + grad_u(1));
      78      481338 :     if (_dim >= 3)
      79             :     {
      80           0 :       const auto grad_w = _w->gradient(face, state);
      81           0 :       symmetric_strain_tensor_norm += 2.0 * Utility::pow<2>(grad_w(2)) +
      82           0 :                                       Utility::pow<2>(grad_u(2) + grad_w(0)) +
      83           0 :                                       Utility::pow<2>(grad_v(2) + grad_w(1));
      84             :     }
      85             :   }
      86             : 
      87      481338 :   symmetric_strain_tensor_norm = sqrt(symmetric_strain_tensor_norm + offset);
      88             : 
      89             :   // Interpolate the mixing length to the face
      90      481338 :   ADReal mixing_len = _mixing_len(face, state);
      91             : 
      92             :   // Compute the eddy diffusivity for momentum
      93      481338 :   ADReal eddy_diff = symmetric_strain_tensor_norm * mixing_len * mixing_len;
      94             : 
      95             :   // Use the turbulent Schmidt/Prandtl number to get the eddy diffusivity for
      96             :   // the scalar variable
      97      481338 :   eddy_diff /= _schmidt_number;
      98             : 
      99             :   // Compute the diffusive flux of the scalar variable, skewness correction is not used here
     100      481338 :   auto dudn = gradUDotNormal(state, _correct_skewness);
     101      481338 :   return -1 * eddy_diff * dudn;
     102             : }