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 "NSSUPGEnergy.h"
      11             : 
      12             : // FluidProperties includes
      13             : #include "IdealGasFluidProperties.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", NSSUPGEnergy);
      16             : 
      17             : InputParameters
      18          19 : NSSUPGEnergy::validParams()
      19             : {
      20          19 :   InputParameters params = NSSUPGBase::validParams();
      21          19 :   params.addClassDescription(
      22             :       "Compute residual and Jacobian terms form the SUPG terms in the energy equation.");
      23          19 :   return params;
      24           0 : }
      25             : 
      26          10 : NSSUPGEnergy::NSSUPGEnergy(const InputParameters & parameters) : NSSUPGBase(parameters) {}
      27             : 
      28             : Real
      29     3698688 : NSSUPGEnergy::computeQpResidual()
      30             : {
      31             :   // See "Component SUPG contributions" section of notes for details.
      32             : 
      33             :   // Values to be summed up and returned
      34             :   Real mass_term = 0.0;
      35             :   Real mom_term = 0.0;
      36             :   Real energy_term = 0.0;
      37             : 
      38             :   // Ratio of specific heats
      39     3698688 :   const Real gam = _fp.gamma();
      40             : 
      41             :   // Velocity vector
      42     3698688 :   RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      43             : 
      44             :   // Velocity vector magnitude squared
      45             :   Real velmag2 = vel.norm_sq();
      46             : 
      47             :   // Velocity vector, dotted with the test function gradient
      48     3698688 :   Real U_grad_phi = vel * _grad_test[_i][_qp];
      49             : 
      50             :   // Vector object of momentum equation strong residuals
      51             :   RealVectorValue Ru(
      52     3698688 :       _strong_residuals[_qp][1], _strong_residuals[_qp][2], _strong_residuals[_qp][3]);
      53             : 
      54             :   // 1.) The mass-residual term:
      55     3698688 :   Real mass_coeff = (0.5 * (gam - 1.0) * velmag2 - _specific_total_enthalpy[_qp]) * U_grad_phi;
      56             : 
      57     3698688 :   mass_term = _tauc[_qp] * mass_coeff * _strong_residuals[_qp][0];
      58             : 
      59             :   // 2.) The momentum-residual term:
      60     3698688 :   Real mom_term1 = _specific_total_enthalpy[_qp] * (_grad_test[_i][_qp] * Ru);
      61     3698688 :   Real mom_term2 = (1.0 - gam) * U_grad_phi * (vel * Ru);
      62             : 
      63     3698688 :   mom_term = _taum[_qp] * (mom_term1 + mom_term2);
      64             : 
      65             :   // 3.) The energy-residual term:
      66     3698688 :   energy_term = _taue[_qp] * gam * U_grad_phi * _strong_residuals[_qp][4];
      67             : 
      68             :   // For printing purposes only
      69     3698688 :   Real result = mass_term + mom_term + energy_term;
      70             :   // Moose::out << "result[" << _qp << "]=" << result << std::endl;
      71             : 
      72     3698688 :   return result;
      73             : }
      74             : 
      75             : Real
      76     2248704 : NSSUPGEnergy::computeQpJacobian()
      77             : {
      78             :   // This is the energy equation, so pass the on-diagonal variable number.
      79     2248704 :   return computeJacobianHelper(_rho_et_var_number);
      80             : }
      81             : 
      82             : Real
      83     6746112 : NSSUPGEnergy::computeQpOffDiagJacobian(unsigned int jvar)
      84             : {
      85     6746112 :   return computeJacobianHelper(jvar);
      86             : }
      87             : 
      88             : Real
      89     8994816 : NSSUPGEnergy::computeJacobianHelper(unsigned var)
      90             : {
      91     8994816 :   if (isNSVariable(var))
      92             :   {
      93             :     // Convert the Moose numbering to canonical NS variable numbering.
      94     8994816 :     unsigned mapped_var_number = mapVarNumber(var);
      95             : 
      96             :     // Convenience vars
      97             : 
      98             :     // Velocity vector
      99     8994816 :     RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     100             : 
     101             :     // Velocity vector magnitude squared
     102             :     Real velmag2 = vel.norm_sq();
     103             : 
     104             :     // Ratio of specific heats
     105     8994816 :     const Real gam = _fp.gamma();
     106             : 
     107             :     // Shortcuts for shape function gradients at current qp.
     108     8994816 :     RealVectorValue grad_test_i = _grad_test[_i][_qp];
     109     8994816 :     RealVectorValue grad_phi_j = _grad_phi[_j][_qp];
     110             : 
     111             :     // ...
     112             : 
     113             :     // 1.) taum- and taue-proportional terms present for any variable:
     114             : 
     115             :     //
     116             :     // Art. Diffusion matrix for taum-proportional term = (diag(H) + (1-gam)*S) * A_{ell}
     117             :     //
     118             :     RealTensorValue mom_mat;
     119     8994816 :     mom_mat(0, 0) = mom_mat(1, 1) = mom_mat(2, 2) = _specific_total_enthalpy[_qp]; // (diag(H)
     120     8994816 :     mom_mat += (1. - gam) * _calC[_qp][0] * _calC[_qp][0].transpose();             //  + (1-gam)*S)
     121     8994816 :     mom_mat = mom_mat * _calA[_qp][mapped_var_number];                             // * A_{ell}
     122     8994816 :     Real mom_term = _taum[_qp] * grad_test_i * (mom_mat * grad_phi_j);
     123             : 
     124             :     //
     125             :     // Art. Diffusion matrix for taue-proportinal term = gam * E_{ell},
     126             :     // where E_{ell} = C_k * E_{k ell} for any k, summation over k *not* implied.
     127             :     //
     128     8994816 :     RealTensorValue ene_mat = gam * _calC[_qp][0] * _calE[_qp][0][mapped_var_number];
     129     8994816 :     Real ene_term = _taue[_qp] * grad_test_i * (ene_mat * grad_phi_j);
     130             : 
     131             :     // 2.) Terms only present if the variable is one of the momentums
     132             :     Real mass_term = 0.;
     133             : 
     134     8994816 :     switch (mapped_var_number)
     135             :     {
     136     4497408 :       case 1:
     137             :       case 2:
     138             :       case 3:
     139             :       {
     140             :         // Variable for zero-based indexing into local matrices and vectors.
     141             :         unsigned m_local = mapped_var_number - 1;
     142             : 
     143             :         //
     144             :         // Art. Diffusion matrix for tauc-proportional term = (0.5*(gam-1.)*velmag2 - H)*C_m
     145             :         //
     146             :         RealTensorValue mass_mat =
     147     4497408 :             (0.5 * (gam - 1.) * velmag2 - _specific_total_enthalpy[_qp]) * _calC[_qp][m_local];
     148     4497408 :         mass_term = _tauc[_qp] * grad_test_i * (mass_mat * grad_phi_j);
     149             : 
     150             :         // Don't even need to break, no other cases to fall through to...
     151             :         break;
     152             :       }
     153             :     }
     154             : 
     155             :     // Sum up values and return
     156     8994816 :     return mass_term + mom_term + ene_term;
     157             :   }
     158             :   else
     159             :     return 0.0;
     160             : }