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 "INSFEFluidMomentumKernel.h"
      11             : 
      12             : registerMooseObject("NavierStokesApp", INSFEFluidMomentumKernel);
      13             : registerMooseObjectRenamed("NavierStokesApp",
      14             :                            MDFluidMomentumKernel,
      15             :                            "02/01/2024 00:00",
      16             :                            INSFEFluidMomentumKernel);
      17             : 
      18             : InputParameters
      19         228 : INSFEFluidMomentumKernel::validParams()
      20             : {
      21         228 :   InputParameters params = INSFEFluidKernelStabilization::validParams();
      22         228 :   params.addClassDescription("Adds advection, viscous, pressure, friction, and gravity terms to "
      23             :                              "the Navier-Stokes momentum equation, potentially with stabilization");
      24         456 :   params.addParam<bool>("conservative_form", false, "Whether conservative form is used");
      25         456 :   params.addParam<bool>(
      26         456 :       "p_int_by_parts", false, "Whether integration by parts is applied to the pressure term");
      27         456 :   params.addRequiredParam<unsigned>("component", "0,1,or 2 for x-, y-, or z- direction");
      28         228 :   return params;
      29           0 : }
      30             : 
      31         120 : INSFEFluidMomentumKernel::INSFEFluidMomentumKernel(const InputParameters & parameters)
      32             :   : INSFEFluidKernelStabilization(parameters),
      33         120 :     _grad_eps(coupledGradient("porosity")),
      34         240 :     _conservative_form(getParam<bool>("conservative_form")),
      35         240 :     _p_int_by_parts(getParam<bool>("p_int_by_parts")),
      36         360 :     _component(getParam<unsigned>("component"))
      37             : {
      38         120 : }
      39             : 
      40             : Real
      41    78095360 : INSFEFluidMomentumKernel::computeQpResidual()
      42             : {
      43    78095360 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
      44    78095360 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      45             : 
      46             :   // Weak form residual of the momentum equation
      47    78095360 :   Real convection_part = _conservative_form
      48    78095360 :                              ? -_rho[_qp] / porosity * _u[_qp] * vec_vel * _grad_test[_i][_qp]
      49    48232960 :                              : _rho[_qp] / porosity * vec_vel * _grad_u[_qp] * _test[_i][_qp];
      50    78095360 :   Real gravity_part = -porosity * _rho[_qp] * _vec_g(_component) * _test[_i][_qp];
      51    78095360 :   Real pressure_part = _p_int_by_parts
      52    78095360 :                            ? -porosity * _pressure[_qp] * _grad_test[_i][_qp](_component) -
      53    31257600 :                                  _pressure[_qp] * _grad_eps[_qp](_component) * _test[_i][_qp]
      54    46837760 :                            : porosity * _grad_pressure[_qp](_component) * _test[_i][_qp];
      55             : 
      56             :   Real viscous_part = 0;
      57             :   Real friction = 0;
      58             :   Real friction_part = 0;
      59    78095360 :   if (porosity > 0.99)
      60             :   {
      61             :     RealVectorValue viscous_stress_vec(_viscous_stress_tensor[_qp](_component, 0),
      62             :                                        _viscous_stress_tensor[_qp](_component, 1),
      63    37260800 :                                        _viscous_stress_tensor[_qp](_component, 2));
      64    37260800 :     viscous_part = viscous_stress_vec * _grad_test[_i][_qp];
      65             :   }
      66             :   else
      67             :   {
      68    40834560 :     Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
      69    40834560 :                             _w_vel[_qp] * _w_vel[_qp]);
      70             :     Real pm_inertial_part =
      71    40834560 :         _inertia_resistance_coeff[_qp](_component, _component) * vec_vel(_component) * velmag;
      72             :     Real pm_viscous_part =
      73    40834560 :         _viscous_resistance_coeff[_qp](_component, _component) * vec_vel(_component);
      74    40834560 :     friction = pm_inertial_part + pm_viscous_part;
      75    40834560 :     friction_part = friction * _test[_i][_qp];
      76             :   }
      77             : 
      78             :   // Assemble weak form terms
      79    78095360 :   Real normal_part = convection_part + pressure_part + gravity_part + viscous_part + friction_part;
      80             : 
      81             :   // SUPG contributions
      82    78095360 :   Real psi_supg = _taum[_qp] * _vel_elem * _grad_test[_i][_qp];
      83             : 
      84    78095360 :   Real transient_supg = _bTransient ? _rho[_qp] * velocityDot()(_component) : 0;
      85    78095360 :   Real convection_supg = _rho[_qp] / porosity * vec_vel * _grad_u[_qp];
      86    78095360 :   Real pressure_supg = porosity * _grad_pressure[_qp](_component);
      87             :   Real gravity_supg = -porosity * _rho[_qp] * _vec_g(_component);
      88             :   Real viscous_supg = 0;
      89             :   Real friction_supg = friction;
      90    78095360 :   if (porosity > 0.99)
      91    37260800 :     viscous_supg = -(_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) * _second_u[_qp].tr();
      92             : 
      93             :   // Assemble SUPG terms
      94    78095360 :   Real supg_part = psi_supg * (transient_supg + convection_supg + pressure_supg + viscous_supg +
      95    78095360 :                                friction_supg + gravity_supg);
      96             : 
      97             :   // Assemble weak form and SUPG contributions
      98    78095360 :   Real res = normal_part + supg_part;
      99    93724160 :   if (_p_int_by_parts && (_component == 0) &&
     100    15628800 :       (_fe_problem.getCoordSystem(_current_elem->subdomain_id()) == Moose::COORD_RZ))
     101           0 :     res -= porosity * _pressure[_qp] / _q_point[_qp](0) * _test[_i][_qp];
     102             : 
     103    78095360 :   return res;
     104             : }
     105             : 
     106             : Real
     107     8657920 : INSFEFluidMomentumKernel::computeQpJacobian()
     108             : {
     109     8657920 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
     110     8657920 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     111             : 
     112             :   // supg parts (some terms can be reused in the weak form)
     113     8657920 :   Real psi_supg = _taum[_qp] * _vel_elem * _grad_test[_i][_qp];
     114     8657920 :   Real transient_supg = _bTransient ? _rho[_qp] * _du_dot_du[_qp] * _phi[_j][_qp] : 0;
     115     8657920 :   Real convection_supg = _rho[_qp] / porosity *
     116     8657920 :                          (_phi[_j][_qp] * _grad_u[_qp](_component) + _grad_phi[_j][_qp] * vec_vel);
     117             : 
     118             :   // Weak form parts
     119             :   Real convection_part = 0;
     120     8657920 :   if (_conservative_form)
     121      665600 :     convection_part = -_rho[_qp] / porosity * _phi[_j][_qp] *
     122      665600 :                       (vec_vel * _grad_test[_i][_qp] + _u[_qp] * _grad_test[_i][_qp](_component));
     123             :   else
     124     7992320 :     convection_part = convection_supg * _test[_i][_qp];
     125             : 
     126             :   Real viscous_part = 0;
     127             :   Real friction_part = 0;
     128             :   Real friction_supg = 0;
     129     8657920 :   if (porosity > 0.99)
     130     4761600 :     viscous_part = (_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) *
     131     4761600 :                    _grad_phi[_j][_qp](_component) * _grad_test[_i][_qp](_component);
     132             :   else
     133             :   {
     134     3896320 :     Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
     135     3896320 :                             _w_vel[_qp] * _w_vel[_qp]);
     136             : 
     137             :     Real pm_inertial_part = 0;
     138     3896320 :     if (velmag > 1e-3)
     139     3793408 :       pm_inertial_part = _inertia_resistance_coeff[_qp](_component, _component) * _phi[_j][_qp] *
     140     3793408 :                          (velmag + vec_vel(_component) * vec_vel(_component) / velmag);
     141     3896320 :     Real pm_viscous_part = _viscous_resistance_coeff[_qp](_component, _component) * _phi[_j][_qp];
     142             : 
     143     3896320 :     friction_part = (pm_inertial_part + pm_viscous_part) * _test[_i][_qp];
     144     3896320 :     friction_supg = (pm_inertial_part + pm_viscous_part) * psi_supg;
     145             :   }
     146             : 
     147             :   // Assemble weak form and SUPG contributions
     148     8657920 :   Real normal_part = convection_part + viscous_part + friction_part;
     149     8657920 :   Real supg_part = psi_supg * (transient_supg + convection_supg) + friction_supg;
     150             : 
     151     8657920 :   return normal_part + supg_part;
     152             : }
     153             : 
     154             : Real
     155    20705280 : INSFEFluidMomentumKernel::computeQpOffDiagJacobian(unsigned int jvar)
     156             : {
     157    20705280 :   unsigned m = this->mapVarNumber(jvar);
     158    20705280 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
     159    20705280 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     160             : 
     161    20705280 :   Real psi_supg = _taum[_qp] * _vel_elem * _grad_test[_i][_qp];
     162             : 
     163             :   Real jac = 0.;
     164    20705280 :   switch (m)
     165             :   {
     166     8657920 :     case 0:
     167     8657920 :       jac = _p_int_by_parts ? -porosity * _phi[_j][_qp] * _grad_test[_i][_qp](_component) -
     168     1013760 :                                   _phi[_j][_qp] * _grad_eps[_qp](_component) * _test[_i][_qp]
     169     7644160 :                             : porosity * _grad_phi[_j][_qp](_component) * _test[_i][_qp];
     170     9164800 :       if (_p_int_by_parts && (_component == 0) &&
     171      506880 :           (_fe_problem.getCoordSystem(_current_elem->subdomain_id()) == Moose::COORD_RZ))
     172           0 :         jac -= porosity * _phi[_j][_qp] / _q_point[_qp](0) * _test[_i][_qp];
     173             :       break;
     174     8657920 :     case 1:
     175             :     case 2:
     176             :     case 3:
     177     8657920 :       if (m != (_component + 1))
     178             :       {
     179             :         // convection term weak form
     180     8657920 :         if (_conservative_form)
     181      665600 :           jac = -_rho[_qp] / porosity * _phi[_j][_qp] * _u[_qp] * _grad_test[_i][_qp](m - 1);
     182             :         else
     183     7992320 :           jac = _rho[_qp] / porosity * _phi[_j][_qp] * _grad_u[_qp](m - 1) * _test[_i][_qp];
     184             :         // convection SUPG term
     185     8657920 :         jac += _rho[_qp] / porosity * _phi[_j][_qp] * _grad_u[_qp](m - 1) * psi_supg;
     186             :         // pm friction
     187     8657920 :         if (porosity <= 0.99)
     188             :         {
     189     3896320 :           Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
     190     3896320 :                                   _w_vel[_qp] * _w_vel[_qp]);
     191     3896320 :           if (velmag > 1e-3)
     192             :           {
     193             :             // inertia_resistance, both normal and supg parts
     194     3793408 :             jac += _inertia_resistance_coeff[_qp](_component, _component) * _phi[_j][_qp] *
     195     3793408 :                    vec_vel(_component) * vec_vel(m - 1) / velmag * (_test[_i][_qp] + psi_supg);
     196             :           }
     197             :         }
     198             :       }
     199             :       break;
     200             : 
     201             :     case 4:
     202             :     default:
     203             :       jac = 0.;
     204             :   }
     205    20705280 :   return jac;
     206             : }