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             : // Navier-Stokes includes
      11             : #include "NS.h"
      12             : #include "NSMomentumInviscidFlux.h"
      13             : 
      14             : // FluidProperties includes
      15             : #include "IdealGasFluidProperties.h"
      16             : 
      17             : registerMooseObject("NavierStokesApp", NSMomentumInviscidFlux);
      18             : 
      19             : InputParameters
      20          82 : NSMomentumInviscidFlux::validParams()
      21             : {
      22          82 :   InputParameters params = NSKernel::validParams();
      23          82 :   params.addClassDescription(
      24             :       "The inviscid flux (convective + pressure terms) for the momentum conservation equations.");
      25          82 :   params.addRequiredCoupledVar(NS::pressure, "pressure");
      26         164 :   params.addRequiredParam<unsigned int>(
      27             :       "component",
      28             :       "0,1,2 depending on if we are solving the x,y,z component of the momentum equation");
      29          82 :   return params;
      30           0 : }
      31             : 
      32          44 : NSMomentumInviscidFlux::NSMomentumInviscidFlux(const InputParameters & parameters)
      33             :   : NSKernel(parameters),
      34          88 :     _pressure(coupledValue(NS::pressure)),
      35          88 :     _component(getParam<unsigned int>("component"))
      36             : {
      37          44 : }
      38             : 
      39             : Real
      40    10967040 : NSMomentumInviscidFlux::computeQpResidual()
      41             : {
      42             :   // For _component = k,
      43             : 
      44             :   // (rho*u) * u_k = (rho*u_k) * u <- we write it this way
      45    10967040 :   RealVectorValue vec(_u[_qp] * _u_vel[_qp],  // (U_k) * u_1
      46    10967040 :                       _u[_qp] * _v_vel[_qp],  // (U_k) * u_2
      47    10967040 :                       _u[_qp] * _w_vel[_qp]); // (U_k) * u_3
      48             : 
      49             :   // (rho*u_k) * u + e_k * P [ e_k = unit vector in k-direction ]
      50    10967040 :   vec(_component) += _pressure[_qp];
      51             : 
      52             :   // -((rho*u_k) * u + e_k * P) * grad(phi)
      53    10967040 :   return -(vec * _grad_test[_i][_qp]);
      54             : }
      55             : 
      56             : Real
      57     6709248 : NSMomentumInviscidFlux::computeQpJacobian()
      58             : {
      59             :   // The on-diagonal entry corresponds to variable number _component+1.
      60     6709248 :   return computeJacobianHelper(_component + 1);
      61             : }
      62             : 
      63             : Real
      64    20127744 : NSMomentumInviscidFlux::computeQpOffDiagJacobian(unsigned int jvar)
      65             : {
      66    20127744 :   if (isNSVariable(jvar))
      67    20127744 :     return computeJacobianHelper(mapVarNumber(jvar));
      68             :   else
      69             :     return 0.0;
      70             : }
      71             : 
      72             : Real
      73    26836992 : NSMomentumInviscidFlux::computeJacobianHelper(unsigned int m)
      74             : {
      75    26836992 :   const RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      76             : 
      77             :   // Ratio of specific heats
      78    26836992 :   const Real gam = _fp.gamma();
      79             : 
      80    26836992 :   switch (m)
      81             :   {
      82             :     case 0: // density
      83             :     {
      84             :       const Real V2 = vel.norm_sq();
      85     6709248 :       return vel(_component) * (vel * _grad_test[_i][_qp]) -
      86     6709248 :              0.5 * (gam - 1.0) * V2 * _grad_test[_i][_qp](_component);
      87             :     }
      88             : 
      89    13418496 :     case 1:
      90             :     case 2:
      91             :     case 3: // momentums
      92             :     {
      93             :       // Map m into m_local = {0,1,2}
      94    13418496 :       unsigned int m_local = m - 1;
      95             : 
      96             :       // Kronecker delta
      97    13418496 :       const Real delta_kl = (_component == m_local ? 1. : 0.);
      98             : 
      99    13418496 :       return -1.0 *
     100    13418496 :              (vel(_component) * _grad_test[_i][_qp](m_local) +
     101    13418496 :               delta_kl * (vel * _grad_test[_i][_qp]) +
     102    13418496 :               (1. - gam) * vel(m_local) * _grad_test[_i][_qp](_component)) *
     103    13418496 :              _phi[_j][_qp];
     104             :     }
     105             : 
     106     6709248 :     case 4: // energy
     107     6709248 :       return -1.0 * (gam - 1.0) * _phi[_j][_qp] * _grad_test[_i][_qp](_component);
     108             :   }
     109             : 
     110           0 :   mooseError("Shouldn't get here!");
     111             : }