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