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 "INSFVTKEDSourceSink.h"
      11             : #include "NonlinearSystemBase.h"
      12             : #include "NavierStokesMethods.h"
      13             : #include "libmesh/nonlinear_solver.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", INSFVTKEDSourceSink);
      16             : 
      17             : InputParameters
      18         487 : INSFVTKEDSourceSink::validParams()
      19             : {
      20         487 :   InputParameters params = FVElementalKernel::validParams();
      21         487 :   params.addClassDescription("Elemental kernel to compute the production and destruction "
      22             :                              " terms of turbulent kinetic energy dissipation (TKED).");
      23         974 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      24         974 :   params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
      25         974 :   params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
      26         487 :   params.addRequiredParam<MooseFunctorName>(NS::TKE, "Coupled turbulent kinetic energy.");
      27         487 :   params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
      28         487 :   params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
      29         487 :   params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
      30         974 :   params.addParam<std::vector<BoundaryName>>(
      31             :       "walls", {}, "Boundaries that correspond to solid walls.");
      32         974 :   params.addParam<bool>(
      33             :       "linearized_model",
      34         974 :       true,
      35             :       "Boolean to determine if the problem should be used in a linear or nonlinear solve");
      36         974 :   MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
      37         974 :   params.addParam<MooseEnum>("wall_treatment",
      38             :                              wall_treatment,
      39             :                              "The method used for computing the wall functions "
      40             :                              "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'");
      41         974 :   params.addParam<MooseFunctorName>("C1_eps", 1.44, "First epsilon coefficient");
      42         974 :   params.addParam<MooseFunctorName>("C2_eps", 1.92, "Second epsilon coefficient");
      43         974 :   params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
      44         974 :   params.addParam<Real>("C_pl", 10.0, "Production limiter constant multiplier.");
      45         487 :   params.set<unsigned short>("ghost_layers") = 2;
      46         974 :   params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
      47         974 :   params.addParamNamesToGroup("newton_solve", "Advanced");
      48         487 :   return params;
      49         487 : }
      50             : 
      51         257 : INSFVTKEDSourceSink::INSFVTKEDSourceSink(const InputParameters & params)
      52             :   : FVElementalKernel(params),
      53         257 :     _dim(_subproblem.mesh().dimension()),
      54         514 :     _u_var(getFunctor<ADReal>("u")),
      55         771 :     _v_var(params.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
      56         257 :     _w_var(params.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
      57         257 :     _k(getFunctor<ADReal>(NS::TKE)),
      58         257 :     _rho(getFunctor<ADReal>(NS::density)),
      59         257 :     _mu(getFunctor<ADReal>(NS::mu)),
      60         257 :     _mu_t(getFunctor<ADReal>(NS::mu_t)),
      61         514 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
      62         514 :     _linearized_model(getParam<bool>("linearized_model")),
      63         514 :     _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
      64         514 :     _C1_eps(getFunctor<ADReal>("C1_eps")),
      65         514 :     _C2_eps(getFunctor<ADReal>("C2_eps")),
      66         514 :     _C_mu(getParam<Real>("C_mu")),
      67         514 :     _C_pl(getParam<Real>("C_pl")),
      68         771 :     _newton_solve(getParam<bool>("newton_solve"))
      69             : {
      70         257 :   if (_dim >= 2 && !_v_var)
      71           0 :     paramError("v", "In two or more dimensions, the v velocity must be supplied!");
      72             : 
      73         257 :   if (_dim >= 3 && !_w_var)
      74           0 :     paramError("w", "In three or more dimensions, the w velocity must be supplied!");
      75         257 : }
      76             : 
      77             : void
      78         929 : INSFVTKEDSourceSink::initialSetup()
      79             : {
      80        1858 :   NS::getWallBoundedElements(
      81         929 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      82         929 :   NS::getWallDistance(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _dist);
      83         929 :   NS::getElementFaceArgs(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _face_infos);
      84         929 : }
      85             : 
      86             : ADReal
      87      351348 : INSFVTKEDSourceSink::computeQpResidual()
      88             : {
      89             :   using std::max, std::sqrt, std::pow, std::min;
      90             : 
      91      351348 :   ADReal residual = 0.0;
      92      351348 :   ADReal production = 0.0;
      93      351348 :   ADReal destruction = 0.0;
      94      351348 :   const auto elem_arg = makeElemArg(_current_elem);
      95      351348 :   const auto state = determineState();
      96             :   const auto old_state =
      97      351348 :       _linearized_model ? Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear) : state;
      98      351348 :   const auto mu = _mu(elem_arg, state);
      99      351348 :   const auto rho = _rho(elem_arg, state);
     100             :   const auto TKE_old =
     101      351348 :       _newton_solve ? max(_k(elem_arg, old_state), 1e-10) : _k(elem_arg, old_state);
     102             :   ADReal y_plus;
     103             : 
     104      351348 :   if (_wall_bounded.find(_current_elem) != _wall_bounded.end())
     105             :   {
     106             :     std::vector<ADReal> y_plus_vec;
     107             : 
     108       83194 :     Real tot_weight = 0.0;
     109             : 
     110      166388 :     ADRealVectorValue velocity(_u_var(elem_arg, state));
     111       83194 :     if (_v_var)
     112       83194 :       velocity(1) = (*_v_var)(elem_arg, state);
     113       83194 :     if (_w_var)
     114           0 :       velocity(2) = (*_w_var)(elem_arg, state);
     115             : 
     116       83194 :     const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem);
     117       83194 :     const auto & distance_vec = libmesh_map_find(_dist, _current_elem);
     118             :     mooseAssert(distance_vec.size(), "Should have found a distance vector");
     119             :     mooseAssert(distance_vec.size() == face_info_vec.size(),
     120             :                 "Should be as many distance vectors as face info vectors");
     121             : 
     122      172110 :     for (unsigned int i = 0; i < distance_vec.size(); i++)
     123             :     {
     124       88916 :       const auto distance = distance_vec[i];
     125             :       mooseAssert(distance > 0, "Should be at a non-zero distance");
     126             : 
     127       88916 :       if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium / Non-iterative
     128        8640 :         y_plus = distance * sqrt(sqrt(_C_mu) * TKE_old) * rho / mu;
     129             :       else
     130             :       {
     131             :         // Equilibrium / Iterative
     132             :         const auto parallel_speed = NS::computeSpeed<ADReal>(
     133      258108 :             velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
     134             : 
     135       86036 :         y_plus = NS::findyPlus<ADReal>(mu, rho, max(parallel_speed, 1e-10), distance);
     136             :       }
     137             : 
     138       88916 :       y_plus_vec.push_back(y_plus);
     139             : 
     140       88916 :       tot_weight += 1.0;
     141             :     }
     142             : 
     143      172110 :     for (const auto i : index_range(y_plus_vec))
     144             :     {
     145       88916 :       const auto y_plus = y_plus_vec[i];
     146             : 
     147       88916 :       if (y_plus < 11.25)
     148             :       {
     149       61452 :         const auto fi = face_info_vec[i];
     150       61452 :         const bool defined_on_elem_side = _var.hasFaceSide(*fi, true);
     151       61452 :         const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr();
     152       61452 :         const Moose::FaceArg facearg = {
     153       61452 :             fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr};
     154       61452 :         destruction += 2.0 * TKE_old * _mu(facearg, state) / rho /
     155      122904 :                        Utility::pow<2>(distance_vec[i]) / tot_weight;
     156             :       }
     157             :       else
     158       27464 :         destruction += pow(_C_mu, 0.75) * pow(TKE_old, 1.5) /
     159       54928 :                        (NS::von_karman_constant * distance_vec[i]) / tot_weight;
     160             :     }
     161             : 
     162      166388 :     residual = _var(makeElemArg(_current_elem), state) - destruction;
     163       83194 :   }
     164             :   else
     165             :   {
     166      268154 :     const auto subdomain_id = _current_elem->subdomain_id();
     167      268154 :     const auto coord_sys = _subproblem.getCoordSystem(subdomain_id);
     168             :     const auto rz_radial_coord =
     169      268154 :         coord_sys == Moose::COORD_RZ ? _subproblem.getAxisymmetricRadialCoord() : 0;
     170             :     const auto symmetric_strain_tensor_sq_norm = NS::computeShearStrainRateNormSquared<ADReal>(
     171      268154 :         _u_var, _v_var, _w_var, elem_arg, state, coord_sys, rz_radial_coord);
     172             : 
     173      268154 :     ADReal production_k = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
     174             :     // Compute production limiter (needed for flows with stagnation zones)
     175             :     const auto eps_old =
     176      268154 :         _newton_solve ? max(_var(elem_arg, old_state), 1e-10) : _var(elem_arg, old_state);
     177      268154 :     const ADReal production_limit = _C_pl * rho * eps_old;
     178             :     // Apply production limiter
     179      268154 :     production_k = min(production_k, production_limit);
     180             : 
     181      268154 :     const auto time_scale = raw_value(TKE_old) / raw_value(eps_old);
     182      536308 :     production = _C1_eps(elem_arg, state) * production_k / time_scale;
     183      536308 :     destruction = _C2_eps(elem_arg, state) * rho * _var(elem_arg, state) / time_scale;
     184             : 
     185      268154 :     residual = destruction - production;
     186             :   }
     187             : 
     188      351348 :   return residual;
     189             : }