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 "NavierStokesMethods.h"
      11             : #include "MooseError.h"
      12             : #include "libmesh/vector_value.h"
      13             : #include "NS.h"
      14             : 
      15             : namespace NS
      16             : {
      17             : int
      18           0 : delta(unsigned int i, unsigned int j)
      19             : {
      20           0 :   if (i == j)
      21             :     return 1;
      22             :   else
      23           0 :     return 0;
      24             : }
      25             : 
      26             : int
      27           0 : computeSign(const Real & a)
      28             : {
      29           0 :   return a > 0 ? 1 : (a < 0 ? -1 : 0);
      30             : }
      31             : 
      32             : unsigned int
      33           0 : getIndex(const Real & p, const std::vector<Real> & bounds)
      34             : {
      35           0 :   if (p < bounds.front() || p > bounds.back())
      36           0 :     mooseError("Point exceeds bounds of domain!");
      37             : 
      38           0 :   for (unsigned int i = 1; i < bounds.size(); ++i)
      39           0 :     if (p <= bounds[i])
      40           0 :       return i - 1;
      41             : 
      42           0 :   return bounds.size() - 2;
      43             : }
      44             : 
      45             : Real
      46      356838 : reynoldsPropertyDerivative(
      47             :     const Real & Re, const Real & rho, const Real & mu, const Real & drho, const Real & dmu)
      48             : {
      49      356838 :   return Re * (drho / std::max(rho, 1e-6) - dmu / std::max(mu, 1e-8));
      50             : }
      51             : 
      52             : Real
      53      356838 : prandtlPropertyDerivative(const Real & mu,
      54             :                           const Real & cp,
      55             :                           const Real & k,
      56             :                           const Real & dmu,
      57             :                           const Real & dcp,
      58             :                           const Real & dk)
      59             : {
      60      356838 :   return (k * (mu * dcp + cp * dmu) - mu * cp * dk) / std::max(k * k, 1e-8);
      61             : }
      62             : 
      63             : template <typename T>
      64             : T
      65     1208304 : findUStar(const T & mu, const T & rho, const T & u, const Real dist)
      66             : {
      67             :   // usually takes about 3-4 iterations
      68             :   constexpr int MAX_ITERS{50};
      69             :   constexpr Real REL_TOLERANCE{1e-6};
      70             : 
      71             :   // Check inputs
      72             :   mooseAssert(mu > 0, "Need a strictly positive viscosity");
      73             :   mooseAssert(rho > 0, "Need a strictly positive density");
      74             :   mooseAssert(u > 0, "Need a strictly positive velocity");
      75             :   mooseAssert(dist > 0, "Need a strictly positive wall distance");
      76             : 
      77      678144 :   const T nu = mu / rho;
      78             : 
      79             :   // Wall-function linearized guess
      80             :   const Real a_c = 1 / NS::von_karman_constant;
      81     2798784 :   const T b_c = 1.0 / NS::von_karman_constant * (std::log(NS::E_turb_constant * dist / mu) + 1.0);
      82             :   const T & c_c = u;
      83             : 
      84             :   /// This is important to reduce the number of nonlinear iterations
      85     2798784 :   T u_star = std::max(1e-20, (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c));
      86             : 
      87             :   // Newton-Raphson method to solve for u_star (friction velocity).
      88     4774622 :   for (int i = 0; i < MAX_ITERS; ++i)
      89             :   {
      90     6843162 :     T residual =
      91     2706082 :         u_star / NS::von_karman_constant * std::log(NS::E_turb_constant * u_star * dist / nu) - u;
      92     6843162 :     T deriv = (1.0 + std::log(NS::E_turb_constant * u_star * dist / nu)) / NS::von_karman_constant;
      93     4774622 :     T new_u_star = std::max(1e-20, u_star - residual / deriv);
      94             : 
      95             :     Real rel_err =
      96     4774622 :         std::abs(MetaPhysicL::raw_value(new_u_star - u_star) / MetaPhysicL::raw_value(new_u_star));
      97             : 
      98     2068540 :     u_star = new_u_star;
      99     4774622 :     if (rel_err < REL_TOLERANCE)
     100     1208304 :       return u_star;
     101             :   }
     102             : 
     103           0 :   mooseException("Could not find the wall friction velocity (mu: ",
     104             :                  mu,
     105             :                  " rho: ",
     106             :                  rho,
     107             :                  " velocity: ",
     108             :                  u,
     109             :                  " wall distance: ",
     110             :                  dist,
     111             :                  ")");
     112             : }
     113             : template Real findUStar<Real>(const Real & mu, const Real & rho, const Real & u, const Real dist);
     114             : template ADReal
     115             : findUStar<ADReal>(const ADReal & mu, const ADReal & rho, const ADReal & u, const Real dist);
     116             : 
     117             : template <typename T>
     118             : T
     119     1078536 : findyPlus(const T & mu, const T & rho, const T & u, const Real dist)
     120             : {
     121             :   // Fixed point iteration method to find y_plus
     122             :   // It should take 3 or 4 iterations
     123             :   constexpr int MAX_ITERS{10};
     124             :   constexpr Real REL_TOLERANCE{1e-2};
     125             : 
     126             :   // Check inputs
     127             :   mooseAssert(mu > 0, "Need a strictly positive viscosity");
     128             :   mooseAssert(u > 0, "Need a strictly positive velocity");
     129             :   mooseAssert(rho > 0, "Need a strictly positive density");
     130             :   mooseAssert(dist > 0, "Need a strictly positive wall distance");
     131             : 
     132             :   Real yPlusLast = 0.0;
     133     1078536 :   T yPlus = dist * u * rho / mu; // Assign initial value to laminar
     134     1078536 :   const Real rev_yPlusLam = 1.0 / MetaPhysicL::raw_value(yPlus);
     135     1078536 :   const T kappa_time_Re = NS::von_karman_constant * u * dist / (mu / rho);
     136             :   unsigned int iters = 0;
     137             : 
     138             :   do
     139             :   {
     140             :     yPlusLast = MetaPhysicL::raw_value(yPlus);
     141             :     // Negative y plus does not make sense
     142     3238375 :     yPlus = std::max(NS::min_y_plus, yPlus);
     143     3997262 :     yPlus = (kappa_time_Re + yPlus) / (1.0 + std::log(NS::E_turb_constant * yPlus));
     144     3238375 :   } while (std::abs(rev_yPlusLam * (MetaPhysicL::raw_value(yPlus) - yPlusLast)) > REL_TOLERANCE &&
     145             :            ++iters < MAX_ITERS);
     146             : 
     147     1078536 :   return std::max(NS::min_y_plus, yPlus);
     148             : }
     149             : template Real findyPlus<Real>(const Real & mu, const Real & rho, const Real & u, Real dist);
     150             : template ADReal
     151             : findyPlus<ADReal>(const ADReal & mu, const ADReal & rho, const ADReal & u, Real dist);
     152             : 
     153             : template <typename T>
     154             : T
     155    81772228 : computeSpeed(const libMesh::VectorValue<T> & velocity)
     156             : {
     157             :   // if the velocity is zero, then the norm function call fails because AD tries to calculate the
     158             :   // derivatives which causes a divide by zero - because d/dx(sqrt(f(x))) = 1/2/sqrt(f(x))*df/dx.
     159             :   // So add a bit of noise (based on hitchhiker's guide to the galaxy's meaning of life number) to
     160             :   // avoid this failure mode.
     161    81772228 :   return isZero(velocity) ? 1e-42 : velocity.norm();
     162             : }
     163             : template Real computeSpeed<Real>(const libMesh::VectorValue<Real> & velocity);
     164             : template ADReal computeSpeed<ADReal>(const libMesh::VectorValue<ADReal> & velocity);
     165             : 
     166             : template <typename T>
     167             : T
     168      807880 : computeShearStrainRateNormSquared(const Moose::Functor<T> & u,
     169             :                                   const Moose::Functor<T> * v,
     170             :                                   const Moose::Functor<T> * w,
     171             :                                   const Moose::ElemArg & elem_arg,
     172             :                                   const Moose::StateArg & state)
     173             : {
     174      807880 :   const auto & grad_u = u.gradient(elem_arg, state);
     175      807880 :   const T Sij_xx = 2.0 * grad_u(0);
     176           0 :   T Sij_xy = 0.0;
     177           0 :   T Sij_xz = 0.0;
     178           0 :   T Sij_yy = 0.0;
     179           0 :   T Sij_yz = 0.0;
     180           0 :   T Sij_zz = 0.0;
     181             : 
     182           0 :   const T grad_xx = grad_u(0);
     183           0 :   T grad_xy = 0.0;
     184           0 :   T grad_xz = 0.0;
     185           0 :   T grad_yx = 0.0;
     186           0 :   T grad_yy = 0.0;
     187           0 :   T grad_yz = 0.0;
     188           0 :   T grad_zx = 0.0;
     189           0 :   T grad_zy = 0.0;
     190           0 :   T grad_zz = 0.0;
     191             : 
     192      807880 :   T trace = Sij_xx / 3.0;
     193             : 
     194      807880 :   if (v) // dim >= 2
     195             :   {
     196      807880 :     const auto & grad_v = (*v).gradient(elem_arg, state);
     197      807880 :     Sij_xy = grad_u(1) + grad_v(0);
     198      807880 :     Sij_yy = 2.0 * grad_v(1);
     199             : 
     200           0 :     grad_xy = grad_u(1);
     201           0 :     grad_yx = grad_v(0);
     202           0 :     grad_yy = grad_v(1);
     203             : 
     204      807880 :     trace += Sij_yy / 3.0;
     205             : 
     206      807880 :     if (w) // dim >= 3
     207             :     {
     208           0 :       const auto & grad_w = (*w).gradient(elem_arg, state);
     209             : 
     210           0 :       Sij_xz = grad_u(2) + grad_w(0);
     211           0 :       Sij_yz = grad_v(2) + grad_w(1);
     212           0 :       Sij_zz = 2.0 * grad_w(2);
     213             : 
     214           0 :       grad_xz = grad_u(2);
     215           0 :       grad_yz = grad_v(2);
     216           0 :       grad_zx = grad_w(0);
     217           0 :       grad_zy = grad_w(1);
     218           0 :       grad_zz = grad_w(2);
     219             : 
     220           0 :       trace += Sij_zz / 3.0;
     221             :     }
     222             :   }
     223             : 
     224      807880 :   return (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
     225      807880 :          (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
     226      807880 :          (Sij_zz - trace) * grad_zz;
     227             : }
     228             : template Real computeShearStrainRateNormSquared<Real>(const Moose::Functor<Real> & u,
     229             :                                                       const Moose::Functor<Real> * v,
     230             :                                                       const Moose::Functor<Real> * w,
     231             :                                                       const Moose::ElemArg & elem_arg,
     232             :                                                       const Moose::StateArg & state);
     233             : template ADReal computeShearStrainRateNormSquared<ADReal>(const Moose::Functor<ADReal> & u,
     234             :                                                           const Moose::Functor<ADReal> * v,
     235             :                                                           const Moose::Functor<ADReal> * w,
     236             :                                                           const Moose::ElemArg & elem_arg,
     237             :                                                           const Moose::StateArg & state);
     238             : 
     239             : /// Bounded element maps for wall treatment
     240             : void
     241       18243 : getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_names,
     242             :                        const FEProblemBase & fe_problem,
     243             :                        const SubProblem & subproblem,
     244             :                        const std::set<SubdomainID> & block_ids,
     245             :                        std::unordered_set<const Elem *> & wall_bounded)
     246             : {
     247             : 
     248             :   wall_bounded.clear();
     249       18243 :   const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_names);
     250             : 
     251     6089940 :   for (const auto & elem : fe_problem.mesh().getMesh().active_element_ptr_range())
     252             :   {
     253     3026727 :     if (block_ids.find(elem->subdomain_id()) != block_ids.end())
     254    15098355 :       for (const auto i_side : elem->side_index_range())
     255             :       {
     256    12078684 :         const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
     257    34131940 :         for (const auto & wall_id : wall_boundary_ids)
     258             :         {
     259    23553162 :           for (const auto side_id : side_bnds)
     260     1499906 :             if (side_id == wall_id)
     261             :               wall_bounded.insert(elem);
     262             :         }
     263    12078684 :       }
     264       18243 :   }
     265       18243 : }
     266             : 
     267             : /// Bounded element face distances for wall treatment
     268             : void
     269        4941 : getWallDistance(const std::vector<BoundaryName> & wall_boundary_name,
     270             :                 const FEProblemBase & fe_problem,
     271             :                 const SubProblem & subproblem,
     272             :                 const std::set<SubdomainID> & block_ids,
     273             :                 std::map<const Elem *, std::vector<Real>> & dist_map)
     274             : {
     275             :   dist_map.clear();
     276             : 
     277     1697412 :   for (const auto & elem : fe_problem.mesh().getMesh().active_element_ptr_range())
     278      843765 :     if (block_ids.find(elem->subdomain_id()) != block_ids.end())
     279     4208745 :       for (const auto i_side : elem->side_index_range())
     280             :       {
     281     3366996 :         const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
     282    12936736 :         for (const auto & name : wall_boundary_name)
     283             :         {
     284     9569740 :           const auto wall_id = subproblem.mesh().getBoundaryID(name);
     285    10214457 :           for (const auto side_id : side_bnds)
     286      644717 :             if (side_id == wall_id)
     287             :             {
     288             :               // The list below stores the face infos with respect to their owning elements,
     289             :               // depending on the block restriction we might encounter situations where the
     290             :               // element outside of the block owns the face info.
     291      173957 :               const auto & neighbor = elem->neighbor_ptr(i_side);
     292      173957 :               const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
     293      173957 :               const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
     294         336 :               const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
     295             : 
     296      173957 :               const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
     297             :               const auto & elem_centroid =
     298      173957 :                   elem_has_fi ? fi->elemCentroid() : fi->neighborCentroid();
     299      173957 :               const Real dist = std::abs((elem_centroid - fi->faceCentroid()) * fi->normal());
     300      173957 :               dist_map[elem].push_back(dist);
     301             :             }
     302             :         }
     303     3371937 :       }
     304        4941 : }
     305             : 
     306             : /// Face arguments to wall-bounded faces for wall treatment
     307             : void
     308        4941 : getElementFaceArgs(const std::vector<BoundaryName> & wall_boundary_name,
     309             :                    const FEProblemBase & fe_problem,
     310             :                    const SubProblem & subproblem,
     311             :                    const std::set<SubdomainID> & block_ids,
     312             :                    std::map<const Elem *, std::vector<const FaceInfo *>> & face_info_map)
     313             : {
     314             :   face_info_map.clear();
     315             : 
     316     1697412 :   for (const auto & elem : fe_problem.mesh().getMesh().active_element_ptr_range())
     317      843765 :     if (block_ids.find(elem->subdomain_id()) != block_ids.end())
     318     4208745 :       for (const auto i_side : elem->side_index_range())
     319             :       {
     320     3366996 :         const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
     321    12936736 :         for (const auto & name : wall_boundary_name)
     322             :         {
     323     9569740 :           const auto wall_id = subproblem.mesh().getBoundaryID(name);
     324    10214457 :           for (const auto side_id : side_bnds)
     325      644717 :             if (side_id == wall_id)
     326             :             {
     327             :               // The list below stores the face infos with respect to their owning elements,
     328             :               // depending on the block restriction we might encounter situations where the
     329             :               // element outside of the block owns the face info.
     330      173957 :               const auto & neighbor = elem->neighbor_ptr(i_side);
     331      173957 :               const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
     332      173957 :               const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
     333         336 :               const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
     334             : 
     335      173957 :               const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
     336      173957 :               face_info_map[elem].push_back(fi);
     337             :             }
     338             :         }
     339     3371937 :       }
     340        4941 : }
     341             : }