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 "INSFVSymmetryVelocityBC.h"
      11             : 
      12             : registerMooseObject("NavierStokesApp", INSFVSymmetryVelocityBC);
      13             : 
      14             : InputParameters
      15        6948 : INSFVSymmetryVelocityBC::validParams()
      16             : {
      17        6948 :   InputParameters params = INSFVFluxBC::validParams();
      18        6948 :   params += INSFVSymmetryBC::validParams();
      19        6948 :   params.addClassDescription("Implements a symmetry boundary condition for the velocity.");
      20       13896 :   params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
      21       13896 :   params.addParam<MooseFunctorName>("v", 0, "The velocity in the y direction.");
      22       13896 :   params.addParam<MooseFunctorName>("w", 0, "The velocity in the z direction.");
      23       13896 :   params.addRequiredParam<MooseFunctorName>("mu", "The viscosity");
      24        6948 :   return params;
      25           0 : }
      26             : 
      27        3904 : INSFVSymmetryVelocityBC::INSFVSymmetryVelocityBC(const InputParameters & params)
      28             :   : INSFVFluxBC(params),
      29             :     INSFVSymmetryBC(params),
      30        3904 :     _u_functor(getFunctor<ADReal>("u")),
      31        7808 :     _v_functor(getFunctor<ADReal>("v")),
      32        7808 :     _w_functor(getFunctor<ADReal>("w")),
      33        7808 :     _mu(getFunctor<ADReal>("mu")),
      34        7808 :     _dim(_subproblem.mesh().dimension())
      35             : {
      36        3904 : }
      37             : 
      38             : ADReal
      39      370192 : INSFVSymmetryVelocityBC::computeSegregatedContribution()
      40             : {
      41      370192 :   const bool use_elem = (_face_type == FaceInfo::VarFaceNeighbors::ELEM);
      42             :   const auto elem_arg =
      43      370192 :       use_elem ? makeElemArg(&_face_info->elem()) : makeElemArg(_face_info->neighborPtr());
      44             : 
      45      370192 :   const auto state = determineState();
      46             :   const auto state_old = Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear);
      47             : 
      48      370192 :   const auto normal = use_elem ? _face_info->normal() : Point(-_face_info->normal());
      49             :   const auto & cell_centroid =
      50      370192 :       use_elem ? _face_info->elemCentroid() : _face_info->neighborCentroid();
      51             : 
      52             :   ADRealVectorValue vel_C;
      53      370192 :   vel_C(0) = _u_functor(elem_arg, state);
      54      370192 :   if (_dim > 1)
      55             :   {
      56      370192 :     vel_C(1) = _v_functor(elem_arg, state);
      57      370192 :     if (_dim > 2)
      58           0 :       vel_C(2) = _w_functor(elem_arg, state);
      59             :   }
      60             : 
      61             :   ADRealVectorValue vel_C_old;
      62      370192 :   vel_C_old(0) = _u_functor(elem_arg, state_old);
      63      370192 :   if (_dim > 1)
      64             :   {
      65      370192 :     vel_C_old(1) = _v_functor(elem_arg, state_old);
      66      370192 :     if (_dim > 2)
      67           0 :       vel_C_old(2) = _w_functor(elem_arg, state_old);
      68             :   }
      69             : 
      70      370192 :   const auto d_perpendicular = std::abs((_face_info->faceCentroid() - cell_centroid) * normal);
      71             : 
      72             :   const auto face = singleSidedFaceArg();
      73      370192 :   const auto mu_b = _mu(face, state);
      74             : 
      75      370192 :   auto normal_x_normal = outer_product(normal, normal);
      76     1110576 :   for (const auto dim_i : make_range(_dim))
      77      740384 :     normal_x_normal(dim_i, dim_i) = 0.0;
      78             : 
      79      370192 :   const auto normal_squared = normal(_index) * normal(_index);
      80             : 
      81      370192 :   const auto face_flux = _rc_uo.getVelocity(Moose::FV::InterpMethod::RhieChow,
      82      370192 :                                             *_face_info,
      83             :                                             state,
      84             :                                             _tid,
      85             :                                             /*subtract_mesh_velocity=*/true) *
      86      370192 :                          normal;
      87             : 
      88             :   ADReal matrix_contribution;
      89             :   ADReal rhs_contribution;
      90             : 
      91             :   // First, add the advective flux
      92      740384 :   matrix_contribution = face_flux * (1 - normal_squared) * vel_C(_index);
      93      740384 :   rhs_contribution = face_flux * (normal_x_normal * vel_C_old)(_index);
      94             : 
      95             :   // Second, add the diffusive flux
      96      740384 :   matrix_contribution += mu_b / d_perpendicular * normal_squared * vel_C(_index);
      97      740384 :   rhs_contribution += mu_b / d_perpendicular * (normal_x_normal * vel_C_old)(_index);
      98             : 
      99      370192 :   return matrix_contribution - rhs_contribution;
     100             : }
     101             : 
     102             : void
     103      826204 : INSFVSymmetryVelocityBC::gatherRCData(const FaceInfo & fi)
     104             : {
     105      826204 :   _face_info = &fi;
     106      826204 :   _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
     107             : 
     108             :   const bool use_elem = _face_type == FaceInfo::VarFaceNeighbors::ELEM;
     109             :   const auto elem_arg =
     110      826204 :       use_elem ? makeElemArg(&_face_info->elem()) : makeElemArg(_face_info->neighborPtr());
     111      826204 :   const auto state = determineState();
     112      826204 :   const auto normal = use_elem ? _face_info->normal() : Point(-_face_info->normal());
     113             :   const Point & cell_centroid =
     114      826204 :       use_elem ? _face_info->elemCentroid() : _face_info->neighborCentroid();
     115      826204 :   const auto u_C = _u_functor(elem_arg, state);
     116      826204 :   const auto v_C = _v_functor(elem_arg, state);
     117      826204 :   const auto w_C = _w_functor(elem_arg, state);
     118             : 
     119      826204 :   const auto d_perpendicular = std::abs((_face_info->faceCentroid() - cell_centroid) * normal);
     120             : 
     121             :   // See Moukalled 15.150. Recall that we multiply by the area in the base class, so S_b ->
     122             :   // normal.norm() -> 1 here.
     123             : 
     124      826204 :   const auto face = singleSidedFaceArg();
     125      826204 :   const auto mu_b = _mu(face, state);
     126             : 
     127             :   ADReal v_dot_n = u_C * normal(0);
     128      826204 :   if (_dim > 1)
     129      826204 :     v_dot_n += v_C * normal(1);
     130      826204 :   if (_dim > 2)
     131           0 :     v_dot_n += w_C * normal(2);
     132             : 
     133      826204 :   const auto strong_resid = mu_b / d_perpendicular * normal(_index) * v_dot_n;
     134             : 
     135             :   // The strong residual for this object is a superposition of all the velocity components, however,
     136             :   // for the on-diagonal 'a' coefficient, we only care about the coefficient multiplying the
     137             :   // velocity component corresponding to _index, hence v_dot_n -> normal(_index) when moving from
     138             :   // strong_resid -> a
     139      826204 :   const auto a = mu_b / d_perpendicular * normal(_index) * normal(_index);
     140             : 
     141     2478612 :   _rc_uo.addToA((_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? &fi.elem() : fi.neighborPtr(),
     142             :                 _index,
     143      826204 :                 a * (fi.faceArea() * fi.faceCoord()));
     144             : 
     145     1652408 :   addResidualAndJacobian(strong_resid * (fi.faceArea() * fi.faceCoord()));
     146      826204 : }