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