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         492 : INSFEFluidMomentumKernel::validParams()
      20             : {
      21         492 :   InputParameters params = INSFEFluidKernelStabilization::validParams();
      22         492 :   params.addClassDescription("Adds advection, viscous, pressure, friction, and gravity terms to "
      23             :                              "the Navier-Stokes momentum equation, potentially with stabilization");
      24         984 :   params.addParam<bool>("conservative_form", false, "Whether conservative form is used");
      25         984 :   params.addParam<bool>(
      26         984 :       "p_int_by_parts", false, "Whether integration by parts is applied to the pressure term");
      27         984 :   params.addRequiredParam<unsigned>("component", "0,1,or 2 for x-, y-, or z- direction");
      28         492 :   return params;
      29           0 : }
      30             : 
      31         264 : INSFEFluidMomentumKernel::INSFEFluidMomentumKernel(const InputParameters & parameters)
      32             :   : INSFEFluidKernelStabilization(parameters),
      33         264 :     _grad_eps(coupledGradient("porosity")),
      34         528 :     _conservative_form(getParam<bool>("conservative_form")),
      35         528 :     _p_int_by_parts(getParam<bool>("p_int_by_parts")),
      36         792 :     _component(getParam<unsigned>("component"))
      37             : {
      38         264 : }
      39             : 
      40             : Real
      41   118333440 : INSFEFluidMomentumKernel::computeQpResidual()
      42             : {
      43   118333440 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
      44   118333440 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      45             : 
      46             :   // Weak form residual of the momentum equation
      47   118333440 :   Real convection_part = _conservative_form
      48   118333440 :                              ? -_rho[_qp] / porosity * _u[_qp] * vec_vel * _grad_test[_i][_qp]
      49    72189440 :                              : _rho[_qp] / porosity * vec_vel * _grad_u[_qp] * _test[_i][_qp];
      50   118333440 :   Real gravity_part = -porosity * _rho[_qp] * _vec_g(_component) * _test[_i][_qp];
      51   118333440 :   Real pressure_part = _p_int_by_parts
      52   118333440 :                            ? -porosity * _pressure[_qp] * _grad_test[_i][_qp](_component) -
      53    48216320 :                                  _pressure[_qp] * _grad_eps[_qp](_component) * _test[_i][_qp]
      54    70117120 :                            : porosity * _grad_pressure[_qp](_component) * _test[_i][_qp];
      55             : 
      56             :   Real viscous_part = 0;
      57             :   Real friction = 0;
      58             :   Real friction_part = 0;
      59   118333440 :   if (porosity > 0.99)
      60             :   {
      61             :     RealVectorValue viscous_stress_vec(_viscous_stress_tensor[_qp](_component, 0),
      62             :                                        _viscous_stress_tensor[_qp](_component, 1),
      63    55846400 :                                        _viscous_stress_tensor[_qp](_component, 2));
      64    55846400 :     viscous_part = viscous_stress_vec * _grad_test[_i][_qp];
      65             :   }
      66             :   else
      67             :   {
      68    62487040 :     Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
      69    62487040 :                             _w_vel[_qp] * _w_vel[_qp]);
      70             :     Real pm_inertial_part =
      71    62487040 :         _inertia_resistance_coeff[_qp](_component, _component) * vec_vel(_component) * velmag;
      72             :     Real pm_viscous_part =
      73    62487040 :         _viscous_resistance_coeff[_qp](_component, _component) * vec_vel(_component);
      74    62487040 :     friction = pm_inertial_part + pm_viscous_part;
      75    62487040 :     friction_part = friction * _test[_i][_qp];
      76             :   }
      77             : 
      78             :   // Assemble weak form terms
      79   118333440 :   Real normal_part = convection_part + pressure_part + gravity_part + viscous_part + friction_part;
      80             : 
      81             :   // SUPG contributions
      82   118333440 :   Real psi_supg = _taum[_qp] * _vel_elem * _grad_test[_i][_qp];
      83             : 
      84   118333440 :   Real transient_supg = _bTransient ? _rho[_qp] * velocityDot()(_component) : 0;
      85   118333440 :   Real convection_supg = _rho[_qp] / porosity * vec_vel * _grad_u[_qp];
      86   118333440 :   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   118333440 :   if (porosity > 0.99)
      91    55846400 :     viscous_supg = -(_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) * _second_u[_qp].tr();
      92             : 
      93             :   // Assemble SUPG terms
      94   118333440 :   Real supg_part = psi_supg * (transient_supg + convection_supg + pressure_supg + viscous_supg +
      95   118333440 :                                friction_supg + gravity_supg);
      96             : 
      97             :   // Assemble weak form and SUPG contributions
      98   118333440 :   Real res = normal_part + supg_part;
      99   142441600 :   if (_p_int_by_parts && (_component == 0) &&
     100    24108160 :       (_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   118333440 :   return res;
     104             : }
     105             : 
     106             : Real
     107    12846080 : INSFEFluidMomentumKernel::computeQpJacobian()
     108             : {
     109    12846080 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
     110    12846080 :   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    12846080 :   Real psi_supg = _taum[_qp] * _vel_elem * _grad_test[_i][_qp];
     114    12846080 :   Real transient_supg = _bTransient ? _rho[_qp] * _du_dot_du[_qp] * _phi[_j][_qp] : 0;
     115    12846080 :   Real convection_supg = _rho[_qp] / porosity *
     116    12846080 :                          (_phi[_j][_qp] * _grad_u[_qp](_component) + _grad_phi[_j][_qp] * vec_vel);
     117             : 
     118             :   // Weak form parts
     119             :   Real convection_part = 0;
     120    12846080 :   if (_conservative_form)
     121      972800 :     convection_part = -_rho[_qp] / porosity * _phi[_j][_qp] *
     122      972800 :                       (vec_vel * _grad_test[_i][_qp] + _u[_qp] * _grad_test[_i][_qp](_component));
     123             :   else
     124    11873280 :     convection_part = convection_supg * _test[_i][_qp];
     125             : 
     126             :   Real viscous_part = 0;
     127             :   Real friction_part = 0;
     128             :   Real friction_supg = 0;
     129    12846080 :   if (porosity > 0.99)
     130     7065600 :     viscous_part = (_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) *
     131     7065600 :                    _grad_phi[_j][_qp](_component) * _grad_test[_i][_qp](_component);
     132             :   else
     133             :   {
     134     5780480 :     Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
     135     5780480 :                             _w_vel[_qp] * _w_vel[_qp]);
     136             : 
     137             :     Real pm_inertial_part = 0;
     138     5780480 :     if (velmag > 1e-3)
     139     5631488 :       pm_inertial_part = _inertia_resistance_coeff[_qp](_component, _component) * _phi[_j][_qp] *
     140     5631488 :                          (velmag + vec_vel(_component) * vec_vel(_component) / velmag);
     141     5780480 :     Real pm_viscous_part = _viscous_resistance_coeff[_qp](_component, _component) * _phi[_j][_qp];
     142             : 
     143     5780480 :     friction_part = (pm_inertial_part + pm_viscous_part) * _test[_i][_qp];
     144     5780480 :     friction_supg = (pm_inertial_part + pm_viscous_part) * psi_supg;
     145             :   }
     146             : 
     147             :   // Assemble weak form and SUPG contributions
     148    12846080 :   Real normal_part = convection_part + viscous_part + friction_part;
     149    12846080 :   Real supg_part = psi_supg * (transient_supg + convection_supg) + friction_supg;
     150             : 
     151    12846080 :   return normal_part + supg_part;
     152             : }
     153             : 
     154             : Real
     155    30720000 : INSFEFluidMomentumKernel::computeQpOffDiagJacobian(unsigned int jvar)
     156             : {
     157    30720000 :   unsigned m = this->mapVarNumber(jvar);
     158    30720000 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
     159    30720000 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     160             : 
     161    30720000 :   Real psi_supg = _taum[_qp] * _vel_elem * _grad_test[_i][_qp];
     162             : 
     163             :   Real jac = 0.;
     164    30720000 :   switch (m)
     165             :   {
     166    12846080 :     case 0:
     167    12846080 :       jac = _p_int_by_parts ? -porosity * _phi[_j][_qp] * _grad_test[_i][_qp](_component) -
     168     1489920 :                                   _phi[_j][_qp] * _grad_eps[_qp](_component) * _test[_i][_qp]
     169    11356160 :                             : porosity * _grad_phi[_j][_qp](_component) * _test[_i][_qp];
     170    13591040 :       if (_p_int_by_parts && (_component == 0) &&
     171      744960 :           (_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    12846080 :     case 1:
     175             :     case 2:
     176             :     case 3:
     177    12846080 :       if (m != (_component + 1))
     178             :       {
     179             :         // convection term weak form
     180    12846080 :         if (_conservative_form)
     181      972800 :           jac = -_rho[_qp] / porosity * _phi[_j][_qp] * _u[_qp] * _grad_test[_i][_qp](m - 1);
     182             :         else
     183    11873280 :           jac = _rho[_qp] / porosity * _phi[_j][_qp] * _grad_u[_qp](m - 1) * _test[_i][_qp];
     184             :         // convection SUPG term
     185    12846080 :         jac += _rho[_qp] / porosity * _phi[_j][_qp] * _grad_u[_qp](m - 1) * psi_supg;
     186             :         // pm friction
     187    12846080 :         if (porosity <= 0.99)
     188             :         {
     189     5780480 :           Real velmag = std::sqrt(_u_vel[_qp] * _u_vel[_qp] + _v_vel[_qp] * _v_vel[_qp] +
     190     5780480 :                                   _w_vel[_qp] * _w_vel[_qp]);
     191     5780480 :           if (velmag > 1e-3)
     192             :           {
     193             :             // inertia_resistance, both normal and supg parts
     194     5631488 :             jac += _inertia_resistance_coeff[_qp](_component, _component) * _phi[_j][_qp] *
     195     5631488 :                    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    30720000 :   return jac;
     206             : }