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 "RANSYPlusAux.h"
      11             : #include "NonlinearSystemBase.h"
      12             : #include "NavierStokesMethods.h"
      13             : #include "libmesh/nonlinear_solver.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", RANSYPlusAux);
      16             : 
      17             : InputParameters
      18         129 : RANSYPlusAux::validParams()
      19             : {
      20         129 :   InputParameters params = AuxKernel::validParams();
      21         129 :   params.addClassDescription("Calculates non-dimensional wall distance (y+) value.");
      22         258 :   params.addRequiredCoupledVar("u", "The velocity in the x direction.");
      23         258 :   params.addCoupledVar("v", "The velocity in the y direction.");
      24         258 :   params.addCoupledVar("w", "The velocity in the z direction.");
      25         258 :   params.addParam<MooseFunctorName>("k", "Turbulent kinetic energy functor.");
      26         258 :   params.deprecateParam("k", NS::TKE, "01/01/2025");
      27         129 :   params.addRequiredParam<MooseFunctorName>(NS::density, "Fluid density.");
      28         129 :   params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
      29         258 :   params.addParam<std::vector<BoundaryName>>(
      30             :       "walls", {}, "Boundaries that correspond to solid walls.");
      31         258 :   MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
      32         258 :   params.addParam<MooseEnum>("wall_treatment",
      33             :                              wall_treatment,
      34             :                              "The method used for computing the y_plus in the wall functions "
      35             :                              "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'");
      36         258 :   params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure coefficient.");
      37             : 
      38         129 :   return params;
      39         129 : }
      40             : 
      41          66 : RANSYPlusAux::RANSYPlusAux(const InputParameters & params)
      42             :   : AuxKernel(params),
      43          66 :     _dim(_subproblem.mesh().dimension()),
      44         132 :     _u_var(getFunctor<Real>("u")),
      45         198 :     _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
      46          66 :     _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
      47         132 :     _k(params.isParamValid(NS::TKE) ? &(getFunctor<Real>(NS::TKE)) : nullptr),
      48          66 :     _rho(getFunctor<Real>(NS::density)),
      49          66 :     _mu(getFunctor<Real>(NS::mu)),
      50         132 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
      51         132 :     _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
      52         198 :     _C_mu(getParam<Real>("C_mu"))
      53             : {
      54          66 :   if (_dim >= 2 && !_v_var)
      55           0 :     paramError("v", "In two or more dimensions, the v velocity must be supplied!");
      56             : 
      57          66 :   if (_dim >= 3 && !_w_var)
      58           0 :     paramError("w", "In three or more dimensions, the w velocity must be supplied!");
      59             : 
      60          66 :   if (_wall_treatment == NS::WallTreatmentEnum::NEQ && !_k)
      61           0 :     paramError(NS::TKE, "In the non-equilibrium wall treatment the TKE must be supplied!");
      62          66 : }
      63             : 
      64             : void
      65          66 : RANSYPlusAux::initialSetup()
      66             : {
      67             :   // Getting wall treatment maps
      68         132 :   NS::getWallBoundedElements(
      69          66 :       _wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _wall_bounded);
      70          66 :   NS::getWallDistance(_wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _dist);
      71          66 :   NS::getElementFaceArgs(_wall_boundary_names, _c_fe_problem, _subproblem, blockIDs(), _face_infos);
      72          66 : }
      73             : 
      74             : Real
      75     3851752 : RANSYPlusAux::computeValue()
      76             : {
      77     3851752 :   if (_wall_bounded.find(_current_elem) != _wall_bounded.end())
      78             :   {
      79     1004912 :     const auto state = determineState();
      80     1004912 :     const auto elem_arg = makeElemArg(_current_elem);
      81     1004912 :     const auto rho = _rho(elem_arg, state);
      82     1004912 :     const auto mu = _mu(elem_arg, state);
      83             :     Real y_plus;
      84             :     std::vector<Real> y_plus_vec;
      85             : 
      86     1004912 :     RealVectorValue velocity(_u_var(elem_arg, state));
      87     1004912 :     if (_v_var)
      88     1004912 :       velocity(1) = (*_v_var)(elem_arg, state);
      89     1004912 :     if (_w_var)
      90           0 :       velocity(2) = (*_w_var)(elem_arg, state);
      91             : 
      92     1004912 :     const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem);
      93     1004912 :     const auto & distance_vec = libmesh_map_find(_dist, _current_elem);
      94             : 
      95     2018184 :     for (unsigned int i = 0; i < distance_vec.size(); i++)
      96             :     {
      97     1013272 :       const auto parallel_speed = NS::computeSpeed<Real>(
      98     1013272 :           velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
      99     1013272 :       const auto distance = distance_vec[i];
     100             : 
     101     1013272 :       if (_wall_treatment == NS::WallTreatmentEnum::NEQ)
     102             :         // Non-equilibrium / Non-iterative
     103      400000 :         y_plus = std::pow(_C_mu, 0.25) * distance * std::sqrt((*_k)(elem_arg, state)) * rho / mu;
     104             :       else
     105             :         // Equilibrium / Iterative
     106      613296 :         y_plus = NS::findyPlus<Real>(mu, rho, std::max(parallel_speed, 1e-10), distance);
     107             : 
     108     1013272 :       y_plus_vec.push_back(y_plus);
     109             :     }
     110             :     // Return average of y+ for cells with multiple wall faces
     111     1004912 :     return std::accumulate(y_plus_vec.begin(), y_plus_vec.end(), 0.0) / y_plus_vec.size();
     112     1004912 :   }
     113             :   else
     114             :     return 0.;
     115             : }