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 "INSFEFluidMomentumBC.h"
      11             : 
      12             : registerMooseObject("NavierStokesApp", INSFEFluidMomentumBC);
      13             : registerMooseObjectRenamed("NavierStokesApp",
      14             :                            MDFluidMomentumBC,
      15             :                            "02/01/2024 00:00",
      16             :                            INSFEFluidMomentumBC);
      17             : 
      18             : InputParameters
      19         246 : INSFEFluidMomentumBC::validParams()
      20             : {
      21         246 :   InputParameters params = INSFEFluidIntegratedBCBase::validParams();
      22         246 :   params.addClassDescription("Specifies flow of momentum through a boundary");
      23         492 :   params.addParam<bool>("conservative_form", false, "Whether the conservative form is used");
      24         492 :   params.addParam<bool>(
      25         492 :       "p_int_by_parts", false, "Whether integration by parts is applied to the pressure term");
      26         492 :   params.addRequiredParam<unsigned>("component", "0,1,or 2 for x-, y-, or z- direction");
      27         492 :   params.addParam<FunctionName>("p_fn", "Pressure function with time at the boundary");
      28         492 :   params.addParam<FunctionName>("v_fn", "Velocity function with time at the boundary");
      29             : 
      30             :   // coupled with branch pressure and density
      31             :   // The 'branch_center' is a little bit tricky, because SAM 1D and multi-D could be in
      32             :   // different mesh system.
      33             :   //   * The volume branch center is always defined in physical 3D XYZ coordinate system,
      34             :   //   * but multi-D flow could be simulated in 2D XY coordinate system,
      35             :   //   * the volume brance center needs be mapped to the 2D/3D flow mesh system
      36             :   // The pressure at the multi-D boundary and the branch pressure is related by:
      37             :   //   p_boundary = p_branch + rho_branch * (point_boundary - branch_center) * gravity
      38         492 :   params.addCoupledVar("p_branch", "Coupled scalar branch pressure");
      39         492 :   params.addCoupledVar("rho_branch", "Coupled scalar branch density for gravity head calculation");
      40         492 :   params.addParam<VectorValue<Real>>("gravity", "Gravity vector in 2D/3D flow mesh system");
      41         492 :   params.addParam<Point>("branch_center", "Position of branch center in 2D/3D flow mesh system");
      42         246 :   return params;
      43           0 : }
      44             : 
      45         132 : INSFEFluidMomentumBC::INSFEFluidMomentumBC(const InputParameters & parameters)
      46             :   : INSFEFluidIntegratedBCBase(parameters),
      47         132 :     _conservative_form(getParam<bool>("conservative_form")),
      48         264 :     _p_int_by_parts(getParam<bool>("p_int_by_parts")),
      49         264 :     _component(getParam<unsigned>("component")),
      50         264 :     _mu(getMaterialProperty<Real>("dynamic_viscosity")),
      51         264 :     _mu_t(getMaterialProperty<Real>("turbulence_viscosity")),
      52         132 :     _has_pbc(parameters.isParamValid("p_fn")),
      53         132 :     _has_vbc(parameters.isParamValid("v_fn")),
      54         132 :     _p_fn(_has_pbc ? &getFunction("p_fn") : nullptr),
      55         132 :     _v_fn(_has_vbc ? &getFunction("v_fn") : nullptr),
      56         132 :     _has_pbranch(parameters.isParamValid("p_branch")),
      57         132 :     _p_branch(_has_pbranch ? coupledScalarValue("p_branch") : _zero),
      58         132 :     _p_branch_var_number(_has_pbranch ? coupledScalar("p_branch") : libMesh::invalid_uint),
      59         264 :     _rho_branch(_has_pbranch ? coupledScalarValue("rho_branch") : _zero)
      60             : {
      61         132 :   if ((_has_pbc || _has_pbranch) && _has_vbc)
      62           0 :     mooseError("Pressure and velocity cannot be BOTH specified in INSFEFluidMomentumBC.");
      63             :   //
      64         132 :   if (_has_pbc && _has_pbranch)
      65           0 :     mooseError(
      66             :         "Pressure function and branch pressure cannot be BOTH specified in INSFEFluidMomentumBC.");
      67             :   //
      68         132 :   if (_has_pbranch)
      69             :   {
      70           0 :     if (!(isParamValid("gravity") && isParamValid("branch_center")))
      71             :     {
      72           0 :       mooseError(
      73             :           name(),
      74             :           ": this boundary is coupled to a volume branch, ",
      75             :           "please provide 'gravity' vector and 'branch_center' for gravity head calculation.");
      76             :     }
      77           0 :     _vec_g = getParam<VectorValue<Real>>("gravity");
      78           0 :     _branch_center = getParam<Point>("branch_center");
      79             :   }
      80         132 : }
      81             : 
      82             : Real
      83     1257216 : INSFEFluidMomentumBC::computeQpResidual()
      84             : {
      85     1257216 :   Real porosity = _has_porosity ? _porosity[_qp] : 1;
      86     1257216 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      87             : 
      88             :   // pressure bc value
      89             :   Real p_bc = 0.0;
      90     1257216 :   if (_has_pbc)
      91      680416 :     p_bc = _p_fn->value(_t, _q_point[_qp]);
      92      576800 :   else if (_has_pbranch)
      93             :   {
      94           0 :     Real dH = (_q_point[_qp] - _branch_center) * _vec_g;
      95           0 :     p_bc = _p_branch[0] + _rho_branch[0] * dH;
      96             :   }
      97             :   else
      98      576800 :     p_bc = _pressure[_qp];
      99             : 
     100             :   // velocity bc value
     101     1257216 :   Real v_bc = _has_vbc ? -_v_fn->value(_t, _q_point[_qp]) : vec_vel * _normals[_qp];
     102             : 
     103             :   Real viscous_part = (porosity > 0.99)
     104     1257216 :                           ? -(_mu[_qp] + _mu_t[_qp]) * _grad_u[_qp] * _normals[_qp] * _test[_i][_qp]
     105             :                           : 0;
     106     1257216 :   Real p_part = _p_int_by_parts ? porosity * p_bc * _normals[_qp](_component) * _test[_i][_qp] : 0;
     107     1257216 :   Real conv_part = _conservative_form ? _rho[_qp] * _u[_qp] * v_bc / porosity * _test[_i][_qp] : 0;
     108             : 
     109     1257216 :   return p_part + conv_part + viscous_part;
     110             : }
     111             : 
     112             : Real
     113       50176 : INSFEFluidMomentumBC::computeQpJacobian()
     114             : {
     115       50176 :   Real porosity = _has_porosity ? _porosity[_qp] : 1;
     116       50176 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     117       50176 :   Real v_bc = _has_vbc ? -_v_fn->value(_t, _q_point[_qp]) : vec_vel * _normals[_qp];
     118             : 
     119             :   Real conv_part = 0;
     120       50176 :   if (_conservative_form)
     121             :   {
     122       24320 :     if (_has_vbc)
     123       12160 :       conv_part = _rho[_qp] * _phi[_j][_qp] * v_bc / porosity * _test[_i][_qp];
     124             :     else
     125       12160 :       conv_part = _rho[_qp] * _phi[_j][_qp] * v_bc / porosity * _test[_i][_qp] +
     126       12160 :                   _rho[_qp] * _u[_qp] * _phi[_j][_qp] * _normals[_qp](_component) / porosity *
     127             :                       _test[_i][_qp];
     128             :   }
     129             :   Real viscous_part = (porosity > 0.99)
     130       50176 :                           ? -(_mu[_qp] + _mu_t[_qp]) * _grad_phi[_j][_qp](_component) *
     131           0 :                                 _normals[_qp](_component) * _test[_i][_qp]
     132             :                           : 0;
     133             : 
     134       50176 :   return conv_part + viscous_part;
     135             : }
     136             : 
     137             : Real
     138      124672 : INSFEFluidMomentumBC::computeQpOffDiagJacobian(unsigned int jvar)
     139             : {
     140      124672 :   unsigned m = this->mapVarNumber(jvar);
     141      124672 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     142      124672 :   Real porosity = _has_porosity ? _porosity[_qp] : 1;
     143             : 
     144             :   // this is the jacobian w.r.t branch pressure
     145      124672 :   if (jvar == _p_branch_var_number)
     146           0 :     return _p_int_by_parts ? porosity * _normals[_qp](_component) * _test[_i][_qp] : 0;
     147             : 
     148             :   Real jac = 0;
     149      124672 :   switch (m)
     150             :   {
     151       50176 :     case 0:
     152       50176 :       if (!(_has_pbc || _has_pbranch))
     153       12160 :         jac = _p_int_by_parts
     154       12160 :                   ? porosity * _phi[_j][_qp] * _normals[_qp](_component) * _test[_i][_qp]
     155             :                   : 0;
     156             :       break;
     157             : 
     158       50176 :     case 1:
     159             :     case 2:
     160             :     case 3:
     161       50176 :       if (m != (_component + 1))
     162       50176 :         jac = _conservative_form ? _rho[_qp] / porosity * _phi[_j][_qp] * _test[_i][_qp] * _u[_qp] *
     163       24320 :                                        _normals[_qp](m - 1)
     164             :                                  : 0;
     165             :       break;
     166             : 
     167             :     case 4:
     168             :     default:
     169             :       jac = 0;
     170             :   }
     171             : 
     172             :   return jac;
     173             : }