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