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 "MomentumFreeSlipBC.h"
      11             : #include "MooseMesh.h"
      12             : #include "MooseVariable.h"
      13             : #include "SystemBase.h"
      14             : #include "FEProblemBase.h"
      15             : #include "libmesh/numeric_vector.h"
      16             : 
      17             : registerMooseObject("NavierStokesApp", MomentumFreeSlipBC);
      18             : 
      19             : InputParameters
      20          38 : MomentumFreeSlipBC::validParams()
      21             : {
      22          38 :   InputParameters params = NodalNormalBC::validParams();
      23          76 :   params.addRequiredCoupledVar("rho_u", "x-component of velocity");
      24          76 :   params.addRequiredCoupledVar("rho_v", "y-component of velocity");
      25          76 :   params.addCoupledVar("rho_w", "z-component of velocity");
      26          38 :   params.addClassDescription("Implements free slip boundary conditions for the Navier Stokes"
      27             :                              "momentum equation.");
      28          38 :   return params;
      29           0 : }
      30             : 
      31          19 : MomentumFreeSlipBC::MomentumFreeSlipBC(const InputParameters & parameters)
      32             :   : NodalNormalBC(parameters),
      33          38 :     _mesh_dimension(_mesh.dimension()),
      34          19 :     _rho_u(coupledValue("rho_u")),
      35          19 :     _rho_v(_mesh_dimension >= 2 ? coupledValue("rho_v") : _zero),
      36          19 :     _rho_w(_mesh_dimension >= 3 ? coupledValue("rho_w") : _zero),
      37          19 :     _rho_u_var(getVar("rho_u", 0)),
      38          19 :     _rho_v_var(getVar("rho_v", 0)),
      39          38 :     _rho_w_var(getVar("rho_w", 0))
      40             : {
      41          19 :   if (_mesh_dimension == 1)
      42           0 :     mooseError(type(), " is not applicable for one-dimensional mesh.");
      43          19 :   else if (_mesh_dimension == 3)
      44           0 :     mooseError(type(), " has not been implemented for three-dimensional mesh.");
      45          19 :   else if (_mesh_dimension != 2)
      46           0 :     mooseError("Mesh dimension ", std::to_string(_mesh_dimension), " not supported.");
      47             : 
      48          38 :   auto check_var = [this](const auto & var_name, const auto * const var_ptr)
      49             :   {
      50          38 :     if (var_ptr)
      51             :       return;
      52             : 
      53           0 :     if (isCoupledConstant(var_name))
      54           0 :       paramError(var_name, "A coupled constant for this variable is not supported in this class");
      55             :     else
      56           0 :       mooseError(var_name, "This variable must be supplied.");
      57             : 
      58             :     if (!var_ptr->isNodal())
      59             :       paramError(var_name, "Only nodal variables supported");
      60          19 :   };
      61             : 
      62          19 :   check_var("rho_u", _rho_u_var);
      63          19 :   check_var("rho_v", _rho_v_var);
      64          19 :   if (_mesh_dimension == 3)
      65           0 :     check_var("rho_w", _rho_w_var);
      66          19 : }
      67             : 
      68          19 : MomentumFreeSlipBC::~MomentumFreeSlipBC() {}
      69             : 
      70             : bool
      71      298726 : MomentumFreeSlipBC::shouldApply() const
      72             : {
      73             :   // this prevents zeroing out the row
      74      298726 :   return !_fe_problem.currentlyComputingJacobian();
      75             : }
      76             : 
      77             : void
      78      295610 : MomentumFreeSlipBC::computeResidual()
      79             : {
      80      295610 :   _normal = Point(_nx[_qp], _ny[_qp], _nz[_qp]);
      81             : 
      82      295610 :   auto set_residual = [this](auto & residual)
      83             :   {
      84      295610 :     const auto rho_u_dof_idx = _rho_u_var->nodalDofIndex();
      85      295610 :     const auto rho_v_dof_idx = _rho_v_var->nodalDofIndex();
      86             : 
      87      295610 :     const auto rho_un = _normal(0) * _rho_u[0] + _normal(1) * _rho_v[0];
      88      295610 :     const auto Re_u = residual(rho_u_dof_idx);
      89      295610 :     const auto Re_v = residual(rho_v_dof_idx);
      90             : 
      91             :     // We obtain these contributions in 3 steps:
      92             :     // 1.) Tranform the momentum residuals into (tangential, normal)
      93             :     //     components by left-multiplying the residual by:
      94             :     //     R = [tx ty] = [-ny nx]
      95             :     //         [nx ny]   [ nx ny]
      96             :     // 2.) Impose the no-normal-flow BC, rho_un = 0, in the normal momentum component's
      97             :     // residual.
      98             :     // 3.) Transform back to (x,y) momentum components by left-multiplying the residual by
      99             :     // R^{-1} = R^T.
     100      295610 :     const auto rho_u_val =
     101      295610 :         (Re_u * _normal(1) * _normal(1) - Re_v * _normal(0) * _normal(1)) + rho_un * _normal(0);
     102      295610 :     const auto rho_v_val =
     103      295610 :         (-Re_u * _normal(0) * _normal(1) + Re_v * _normal(0) * _normal(0)) + rho_un * _normal(1);
     104             : 
     105             :     // NOTE: we have to handle all components at the same time, otherwise we'd work with the
     106             :     // modified residual and we do not want that
     107      295610 :     residual.set(rho_u_dof_idx, rho_u_val);
     108      295610 :     residual.set(rho_v_dof_idx, rho_v_val);
     109      591220 :   };
     110             : 
     111      295610 :   setResidual(_sys, set_residual);
     112      295610 : }