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 "NSEnergyInviscidFlux.h"
      13             : 
      14             : // FluidProperties includes
      15             : #include "IdealGasFluidProperties.h"
      16             : 
      17             : registerMooseObject("NavierStokesApp", NSEnergyInviscidFlux);
      18             : 
      19             : InputParameters
      20          41 : NSEnergyInviscidFlux::validParams()
      21             : {
      22          41 :   InputParameters params = NSKernel::validParams();
      23          41 :   params.addClassDescription("This class computes the inviscid part of the energy flux.");
      24          41 :   params.addRequiredCoupledVar(NS::specific_total_enthalpy, "specific total enthalpy");
      25          41 :   return params;
      26           0 : }
      27             : 
      28          22 : NSEnergyInviscidFlux::NSEnergyInviscidFlux(const InputParameters & parameters)
      29          22 :   : NSKernel(parameters), _specific_total_enthalpy(coupledValue(NS::specific_total_enthalpy))
      30             : {
      31          22 : }
      32             : 
      33             : Real
      34     5483520 : NSEnergyInviscidFlux::computeQpResidual()
      35             : {
      36             :   // ht = specific total enthalpy = et + p/rho
      37             :   // => rho * u * ht = rho * u ( et + p/rho)
      38             :   //                =       u ( rho*et + p)
      39             : 
      40             :   // velocity vector
      41     5483520 :   RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      42             : 
      43             :   // Multiply vector U by the scalar value (rho*et + P) to get rho * U * ht
      44             :   // vel *= (_u[_qp] + _pressure[_qp]);
      45             : 
      46             :   // Multiply velocity vector by the scalar (rho * ht)
      47     5483520 :   vel *= (_rho[_qp] * _specific_total_enthalpy[_qp]);
      48             : 
      49             :   // Return -1 * vel * grad(phi_i)
      50     5483520 :   return -(vel * _grad_test[_i][_qp]);
      51             : }
      52             : 
      53             : Real
      54     3354624 : NSEnergyInviscidFlux::computeQpJacobian()
      55             : {
      56             :   // Derivative of this kernel wrt rho*et
      57     3354624 :   const RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      58             : 
      59             :   // Ratio of specific heats
      60     3354624 :   const Real gam = _fp.gamma();
      61             : 
      62             :   // -gamma * phi_j * (U*grad(phi_i))
      63     3354624 :   return -gam * _phi[_j][_qp] * (vel * _grad_test[_i][_qp]);
      64             : }
      65             : 
      66             : Real
      67    10063872 : NSEnergyInviscidFlux::computeQpOffDiagJacobian(unsigned int jvar)
      68             : {
      69    10063872 :   if (isNSVariable(jvar))
      70             :   {
      71    10063872 :     RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      72             :     Real V2 = vel.norm_sq();
      73             : 
      74             :     // Ratio of specific heats
      75    10063872 :     const Real gam = _fp.gamma();
      76             : 
      77             :     // Derivative wrt density
      78    10063872 :     if (jvar == _rho_var_number)
      79     3354624 :       return -((0.5 * (gam - 1) * V2 - _specific_total_enthalpy[_qp]) * _phi[_j][_qp] *
      80     3354624 :                (vel * _grad_test[_i][_qp]));
      81             : 
      82             :     // Derivatives wrt momentums
      83     6709248 :     else if ((jvar == _rhou_var_number) || (jvar == _rhov_var_number) || (jvar == _rhow_var_number))
      84             :     {
      85             :       // Map jvar into jlocal = {0,1,2}, regardless of how Moose has numbered things.
      86             :       unsigned jlocal = 0;
      87             : 
      88     6709248 :       if (jvar == _rhov_var_number)
      89             :         jlocal = 1;
      90     3354624 :       else if (jvar == _rhow_var_number)
      91             :         jlocal = 2;
      92             : 
      93             :       // Scale the velocity vector by the scalar (1-gamma)*vel(jlocal)
      94     6709248 :       vel *= (1.0 - gam) * vel(jlocal);
      95             : 
      96             :       // Add in the specific_total_enthalpy in the jlocal'th entry
      97     6709248 :       vel(jlocal) += _specific_total_enthalpy[_qp];
      98             : 
      99             :       // Return -1 * (vel * grad(phi_i)) * phi_j
     100     6709248 :       return -(vel * _grad_test[_i][_qp]) * _phi[_j][_qp];
     101             :     }
     102             : 
     103             :     else
     104             :     {
     105           0 :       std::ostringstream oss;
     106           0 :       oss << "Invalid jvar=" << jvar << " requested!\n"
     107             :           << "Did not match:\n"
     108           0 :           << " _rho_var_number =" << _rho_var_number << "\n"
     109           0 :           << " _rhou_var_number=" << _rhou_var_number << "\n"
     110           0 :           << " _rhov_var_number=" << _rhov_var_number << "\n"
     111           0 :           << " _rhow_var_number=" << _rhow_var_number << std::endl;
     112           0 :       mooseError(oss.str());
     113           0 :     }
     114             :   }
     115             :   else
     116             :     return 0.0;
     117             : }