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 "INSFEFluidMassKernel.h"
      11             : #include "MooseMesh.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", INSFEFluidMassKernel);
      14             : registerMooseObjectRenamed("NavierStokesApp",
      15             :                            MDFluidMassKernel,
      16             :                            "02/01/2024 00:00",
      17             :                            INSFEFluidMassKernel);
      18             : 
      19             : InputParameters
      20         246 : INSFEFluidMassKernel::validParams()
      21             : {
      22         246 :   InputParameters params = INSFEFluidKernelStabilization::validParams();
      23         246 :   params.addClassDescription("Adds advective term of mass conservation equation along with "
      24             :                              "pressure-stabilized Petrov-Galerkin terms");
      25         246 :   return params;
      26           0 : }
      27             : 
      28         132 : INSFEFluidMassKernel::INSFEFluidMassKernel(const InputParameters & parameters)
      29             :   : INSFEFluidKernelStabilization(parameters),
      30         264 :     _u_vel_second(_u_var.secondSln()),
      31         132 :     _v_vel_second(_v_var.secondSln()),
      32         264 :     _w_vel_second(_mesh.dimension() == 3 ? _w_var.secondSln() : _second_zero)
      33             : {
      34         132 : }
      35             : 
      36             : Real
      37    59166720 : INSFEFluidMassKernel::computeQpResidual()
      38             : {
      39    59166720 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      40    59166720 :   Real masseq_part = -_rho[_qp] * vec_vel * _grad_test[_i][_qp];
      41             : 
      42             :   // Compute PSPG term due to momentum equation residual vector
      43    59166720 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
      44             :   RealVectorValue convection_vec(
      45    59166720 :       vec_vel * _grad_u_vel[_qp], vec_vel * _grad_v_vel[_qp], vec_vel * _grad_w_vel[_qp]);
      46    59166720 :   RealVectorValue psi_pspg = _tauc[_qp] * _grad_test[_i][_qp];
      47             : 
      48    59166720 :   Real transient_pspg = _bTransient ? _rho[_qp] * velocityDot() * psi_pspg : 0;
      49    59166720 :   Real convection_pspg = _rho[_qp] / porosity * convection_vec * psi_pspg;
      50    59166720 :   Real pressure_pspg = porosity * _grad_u[_qp] * psi_pspg;
      51    59166720 :   Real gravity_pspg = -porosity * _rho[_qp] * _vec_g * psi_pspg;
      52             : 
      53             :   Real viscous_pspg = 0;
      54             :   Real pm_friction_pspg = 0;
      55    59166720 :   if (porosity > 0.99)
      56             :   {
      57             :     RealVectorValue vec_vel_second(
      58    27923200 :         _u_vel_second[_qp].tr(), _v_vel_second[_qp].tr(), _w_vel_second[_qp].tr());
      59             :     viscous_pspg =
      60    27923200 :         -(_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) * vec_vel_second * psi_pspg;
      61             :   }
      62             :   else
      63             :   {
      64    31243520 :     Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
      65    31243520 :                             _w_vel[_qp] * _w_vel[_qp]);
      66             : 
      67    31243520 :     pm_friction_pspg += _inertia_resistance_coeff[_qp] * vec_vel * velmag * psi_pspg;
      68    31243520 :     pm_friction_pspg += _viscous_resistance_coeff[_qp] * vec_vel * psi_pspg;
      69             :   }
      70             : 
      71             :   // Assemble PSPG terms
      72    59166720 :   Real momeq_part = transient_pspg + convection_pspg + pressure_pspg + viscous_pspg + gravity_pspg +
      73             :                     pm_friction_pspg;
      74             : 
      75             :   // Assemble final residual
      76    59166720 :   return masseq_part + momeq_part;
      77             : }
      78             : 
      79             : Real
      80     6423040 : INSFEFluidMassKernel::computeQpJacobian()
      81             : {
      82     6423040 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
      83     6423040 :   return porosity * _tauc[_qp] * _grad_test[_i][_qp] * _grad_phi[_j][_qp];
      84             : }
      85             : 
      86             : Real
      87    15360000 : INSFEFluidMassKernel::computeQpOffDiagJacobian(unsigned int jvar)
      88             : {
      89             :   // Convert the Moose numbering to internal porous medium model variable numbering.
      90    15360000 :   unsigned m = this->mapVarNumber(jvar);
      91             : 
      92    15360000 :   switch (m)
      93             :   {
      94    12846080 :     case 1:
      95             :     case 2:
      96             :     case 3:
      97             :     {
      98    12846080 :       RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      99    12846080 :       Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
     100    12846080 :                               _w_vel[_qp] * _w_vel[_qp]);
     101             : 
     102    12846080 :       Real mass_eqn_part = -_rho[_qp] * _phi[_j][_qp] * _grad_test[_i][_qp](m - 1);
     103             : 
     104             :       Real pm_inertial_pspg = 0;
     105             :       Real pm_viscous_pspg = 0;
     106    12846080 :       if (velmag < 1e-3)
     107             :         pm_inertial_pspg = 0.;
     108             :       else
     109    12154336 :         pm_inertial_pspg = _inertia_resistance_coeff[_qp](m - 1, m - 1) *
     110    12154336 :                            (velmag + vec_vel(m - 1) * vec_vel(m - 1) / velmag) * _phi[_j][_qp] *
     111    12154336 :                            _tauc[_qp] * _grad_test[_i][_qp](m - 1);
     112    12846080 :       pm_viscous_pspg = _viscous_resistance_coeff[_qp](m - 1, m - 1) * _phi[_j][_qp] * _tauc[_qp] *
     113             :                         _grad_test[_i][_qp](m - 1);
     114             : 
     115    12846080 :       return mass_eqn_part + pm_inertial_pspg + pm_viscous_pspg;
     116             :     }
     117             : 
     118             :     default:
     119             :       return 0;
     120             :   }
     121             : }