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 "INSFVMixingLengthTurbulentViscosityAux.h"
      11             : #include "INSFVVelocityVariable.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", INSFVMixingLengthTurbulentViscosityAux);
      14             : 
      15             : InputParameters
      16          11 : INSFVMixingLengthTurbulentViscosityAux::validParams()
      17             : {
      18          11 :   InputParameters params = AuxKernel::validParams();
      19          11 :   params.addClassDescription("Computes the turbulent viscosity for the mixing length model.");
      20          22 :   params.addRequiredCoupledVar("u", "The velocity in the x direction.");
      21          22 :   params.addCoupledVar("v", "The velocity in the y direction.");
      22          22 :   params.addCoupledVar("w", "The velocity in the z direction.");
      23          22 :   params.addRequiredCoupledVar("mixing_length", "Turbulent eddy mixing length.");
      24          11 :   return params;
      25           0 : }
      26             : 
      27           6 : INSFVMixingLengthTurbulentViscosityAux::INSFVMixingLengthTurbulentViscosityAux(
      28           6 :     const InputParameters & params)
      29             :   : AuxKernel(params),
      30           6 :     _dim(_subproblem.mesh().dimension()),
      31           6 :     _u_var(dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("u", 0))),
      32           6 :     _v_var(params.isParamValid("v")
      33          12 :                ? dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("v", 0))
      34             :                : nullptr),
      35           6 :     _w_var(params.isParamValid("w")
      36           6 :                ? dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("w", 0))
      37             :                : nullptr),
      38          12 :     _mixing_len(coupledValue("mixing_length"))
      39             : {
      40           6 :   if (!_u_var)
      41           0 :     paramError("u", "the u velocity must be an INSFVVelocityVariable.");
      42             : 
      43           6 :   if (_dim >= 2 && !_v_var)
      44           0 :     paramError("v",
      45             :                "In two or more dimensions, the v velocity must be supplied and it must be an "
      46             :                "INSFVVelocityVariable.");
      47             : 
      48           6 :   if (_dim >= 3 && !_w_var)
      49           0 :     paramError("w",
      50             :                "In three-dimensions, the w velocity must be supplied and it must be an "
      51             :                "INSFVVelocityVariable.");
      52           6 : }
      53             : 
      54             : Real
      55       96800 : INSFVMixingLengthTurbulentViscosityAux::computeValue()
      56             : {
      57             :   constexpr Real offset = 1e-15; // prevents explosion of sqrt(x) derivative to infinity
      58       96800 :   const Elem & elem = *_current_elem;
      59       96800 :   const auto state = determineState();
      60             : 
      61       96800 :   const auto grad_u = MetaPhysicL::raw_value(_u_var->adGradSln(&elem, state));
      62       96800 :   Real symmetric_strain_tensor_norm = 2.0 * Utility::pow<2>(grad_u(0));
      63       96800 :   if (_dim >= 2)
      64             :   {
      65       96800 :     const auto grad_v = MetaPhysicL::raw_value(_v_var->adGradSln(&elem, state));
      66       96800 :     symmetric_strain_tensor_norm +=
      67       96800 :         2.0 * Utility::pow<2>(grad_v(1)) + Utility::pow<2>(grad_v(0) + grad_u(1));
      68       96800 :     if (_dim >= 3)
      69             :     {
      70           0 :       const auto grad_w = MetaPhysicL::raw_value(_w_var->adGradSln(&elem, state));
      71           0 :       symmetric_strain_tensor_norm += 2.0 * Utility::pow<2>(grad_w(2)) +
      72           0 :                                       Utility::pow<2>(grad_u(2) + grad_w(0)) +
      73           0 :                                       Utility::pow<2>(grad_v(2) + grad_w(1));
      74             :     }
      75             :   }
      76             : 
      77       96800 :   symmetric_strain_tensor_norm = std::sqrt(symmetric_strain_tensor_norm + offset);
      78             : 
      79             :   // Compute the eddy diffusivitiy
      80       96800 :   const Real eddy_diff = symmetric_strain_tensor_norm * _mixing_len[_qp] * _mixing_len[_qp];
      81             : 
      82             :   // Return the turbulent stress contribution to the momentum equation
      83       96800 :   return eddy_diff;
      84             : }