Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "LinearFVVelocitySymmetryBC.h"
      11             : #include "FEProblemBase.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", LinearFVVelocitySymmetryBC);
      14             : 
      15             : InputParameters
      16         474 : LinearFVVelocitySymmetryBC::validParams()
      17             : {
      18         474 :   InputParameters params = LinearFVAdvectionDiffusionBC::validParams();
      19         474 :   params.addClassDescription("Adds a symmetry boundary condition for the velocity.");
      20         948 :   params.addRequiredParam<SolverVariableName>("u", "The velocity in the x direction.");
      21         948 :   params.addParam<SolverVariableName>("v", "The velocity in the y direction.");
      22         948 :   params.addParam<SolverVariableName>("w", "The velocity in the z direction.");
      23         948 :   MooseEnum momentum_component("x=0 y=1 z=2");
      24         948 :   params.addRequiredParam<MooseEnum>(
      25             :       "momentum_component",
      26             :       momentum_component,
      27             :       "The component of the momentum equation that this kernel applies to.");
      28             : 
      29         474 :   return params;
      30         474 : }
      31             : 
      32         240 : LinearFVVelocitySymmetryBC::LinearFVVelocitySymmetryBC(const InputParameters & parameters)
      33             :   : LinearFVAdvectionDiffusionBC(parameters),
      34         240 :     _dim(_subproblem.mesh().dimension()),
      35         240 :     _u_var(dynamic_cast<const MooseLinearVariableFVReal *>(
      36         480 :         &_fv_problem.getVariable(_tid, getParam<SolverVariableName>("u")))),
      37         240 :     _v_var(parameters.isParamValid("v")
      38         480 :                ? dynamic_cast<const MooseLinearVariableFVReal *>(
      39         960 :                      &_fv_problem.getVariable(_tid, getParam<SolverVariableName>("v")))
      40             :                : nullptr),
      41         480 :     _w_var(parameters.isParamValid("w")
      42         240 :                ? dynamic_cast<const MooseLinearVariableFVReal *>(
      43         240 :                      &_fv_problem.getVariable(_tid, getParam<SolverVariableName>("w")))
      44             :                : nullptr),
      45         720 :     _index(getParam<MooseEnum>("momentum_component"))
      46             : {
      47         240 :   if (!_u_var)
      48           0 :     paramError("u", "the u velocity must be a MooseLinearVariableFVReal.");
      49             : 
      50         240 :   _vel_vars.push_back(_u_var);
      51             : 
      52         240 :   if (_dim >= 2 && !_v_var)
      53           0 :     paramError("v",
      54             :                "In two or more dimensions, the v velocity must be supplied and it must be a "
      55             :                "MooseLinearVariableFVReal.");
      56             : 
      57         240 :   _vel_vars.push_back(_v_var);
      58             : 
      59         240 :   if (_dim >= 3 && !_w_var)
      60           0 :     paramError("w",
      61             :                "In three-dimensions, the w velocity must be supplied and it must be a "
      62             :                "MooseLinearVariableFVReal.");
      63             : 
      64         240 :   _vel_vars.push_back(_w_var);
      65         240 : }
      66             : 
      67             : Real
      68      528660 : LinearFVVelocitySymmetryBC::computeBoundaryValue() const
      69             : {
      70             :   // We allow internal boundaries too so we need to check which side we are on
      71      528660 :   const auto elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
      72      528660 :                              ? _current_face_info->elemInfo()
      73           0 :                              : _current_face_info->neighborInfo();
      74             : 
      75             :   // By default we approximate the boundary value with the neighboring cell value
      76      528660 :   auto boundary_value = _var.getElemValue(*elem_info, determineState());
      77             :   auto reflected_boundary_value = boundary_value;
      78             : 
      79             :   // We don't have to flip the sign of the normal because we are subtacting normal*normal.
      80      528660 :   auto scaled_normal = _current_face_info->normal();
      81      528660 :   scaled_normal *= 2 * scaled_normal(_index);
      82             : 
      83     1585980 :   for (const auto dim_i : make_range(_dim))
      84     1057320 :     reflected_boundary_value -=
      85     1057320 :         scaled_normal(dim_i) * _vel_vars[dim_i]->getElemValue(*elem_info, determineState());
      86             : 
      87      528660 :   return 0.5 * (boundary_value + reflected_boundary_value);
      88             : }
      89             : 
      90             : Real
      91      382948 : LinearFVVelocitySymmetryBC::computeBoundaryNormalGradient() const
      92             : {
      93             :   // We allow internal boundaries too so we need to check which side we are on
      94      382948 :   const auto elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
      95      382948 :                              ? _current_face_info->elemInfo()
      96           0 :                              : _current_face_info->neighborInfo();
      97             : 
      98             :   Real boundary_normal_grad = 0.0;
      99             : 
     100             :   // We don't have to flip the sign of the normal because we are subtacting normal*normal.
     101      382948 :   auto scaled_normal = _current_face_info->normal();
     102      382948 :   scaled_normal *= scaled_normal(_index);
     103             : 
     104     1148844 :   for (const auto dim_i : make_range(_dim))
     105      765896 :     boundary_normal_grad +=
     106      765896 :         scaled_normal(dim_i) * _vel_vars[dim_i]->getElemValue(*elem_info, determineState());
     107             : 
     108      382948 :   return boundary_normal_grad / computeCellToFaceDistance();
     109             : }
     110             : 
     111             : Real
     112      373124 : LinearFVVelocitySymmetryBC::computeBoundaryValueMatrixContribution() const
     113             : {
     114             :   // No matter if we have a one-term or two-term expansion we will always
     115             :   // have a contribution to the matrix
     116      373124 :   const auto normal_component = _current_face_info->normal()(_index);
     117      373124 :   const auto normal_component_sq = normal_component * normal_component;
     118      373124 :   return 1.0 - normal_component_sq;
     119             : }
     120             : 
     121             : Real
     122      373124 : LinearFVVelocitySymmetryBC::computeBoundaryValueRHSContribution() const
     123             : {
     124             :   // We allow internal boundaries too so we need to check which side we are on
     125      373124 :   const auto elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     126      373124 :                              ? _current_face_info->elemInfo()
     127           0 :                              : _current_face_info->neighborInfo();
     128             : 
     129             :   // We don't have to flip the sign of the normal because we are subtacting normal*normal.
     130      373124 :   auto scaled_normal = _current_face_info->normal();
     131      373124 :   scaled_normal *= scaled_normal(_index);
     132             : 
     133      373124 :   auto current_bd_value = computeBoundaryValue();
     134      373124 :   const auto normal_component = _current_face_info->normal()(_index);
     135      373124 :   const auto normal_component_sq = normal_component * normal_component;
     136             : 
     137             :   return current_bd_value -
     138      373124 :          (1.0 - normal_component_sq) * _var.getElemValue(*elem_info, determineState());
     139             : }
     140             : 
     141             : Real
     142      382948 : LinearFVVelocitySymmetryBC::computeBoundaryGradientMatrixContribution() const
     143             : {
     144             :   // We don't have to flip the sign of the normal because we are subtacting normal*normal.
     145      382948 :   const auto normal_component = _current_face_info->normal()(_index);
     146      382948 :   const auto normal_component_sq = normal_component * normal_component;
     147             : 
     148      382948 :   return normal_component_sq / computeCellToFaceDistance();
     149             : }
     150             : 
     151             : Real
     152      382948 : LinearFVVelocitySymmetryBC::computeBoundaryGradientRHSContribution() const
     153             : {
     154             :   // We allow internal boundaries too so we need to check which side we are on
     155      382948 :   const auto elem_info = _current_face_type == FaceInfo::VarFaceNeighbors::ELEM
     156      382948 :                              ? _current_face_info->elemInfo()
     157           0 :                              : _current_face_info->neighborInfo();
     158             : 
     159      382948 :   auto boundary_value = _var.getElemValue(*elem_info, determineState());
     160             : 
     161             :   // We don't have to flip the sign of the normal because we are subtacting normal*normal.
     162      382948 :   const auto normal_component = _current_face_info->normal()(_index);
     163      382948 :   const auto normal_component_sq = normal_component * normal_component;
     164      382948 :   return computeBoundaryNormalGradient() -
     165      382948 :          normal_component_sq / computeCellToFaceDistance() * boundary_value;
     166             : }