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 "NSSUPGMomentum.h"
      11             : 
      12             : // FluidProperties includes
      13             : #include "IdealGasFluidProperties.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", NSSUPGMomentum);
      16             : 
      17             : InputParameters
      18          82 : NSSUPGMomentum::validParams()
      19             : {
      20          82 :   InputParameters params = NSSUPGBase::validParams();
      21          82 :   params.addClassDescription(
      22             :       "Compute residual and Jacobian terms form the SUPG terms in the momentum equation.");
      23         164 :   params.addRequiredParam<unsigned int>("component", "");
      24          82 :   return params;
      25           0 : }
      26             : 
      27          44 : NSSUPGMomentum::NSSUPGMomentum(const InputParameters & parameters)
      28          88 :   : NSSUPGBase(parameters), _component(getParam<unsigned>("component"))
      29             : {
      30          44 : }
      31             : 
      32             : Real
      33    10967040 : NSSUPGMomentum::computeQpResidual()
      34             : {
      35             :   // See "Component SUPG contributions" section of notes for details.
      36             : 
      37             :   // Values to be summed up and returned
      38             :   Real mass_term = 0.0;
      39             :   Real mom_term = 0.0;
      40             :   Real energy_term = 0.0;
      41             : 
      42             :   // Velocity vector
      43    10967040 :   RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      44             : 
      45             :   // Velocity vector magnitude squared
      46             :   Real velmag2 = vel.norm_sq();
      47             : 
      48             :   // Ratio of specific heats
      49    10967040 :   const Real gam = _fp.gamma();
      50             : 
      51             :   // Velocity vector, dotted with the test function gradient
      52    10967040 :   Real U_grad_phi = vel * _grad_test[_i][_qp];
      53             : 
      54             :   // _component'th entry of test function gradient
      55    10967040 :   Real dphi_dxk = _grad_test[_i][_qp](_component);
      56             : 
      57             :   // Vector object of momentum equation strong residuals
      58             :   RealVectorValue Ru(
      59    10967040 :       _strong_residuals[_qp][1], _strong_residuals[_qp][2], _strong_residuals[_qp][3]);
      60             : 
      61             :   // 1.) The mass-residual term:
      62    10967040 :   Real mass_coeff = 0.5 * (gam - 1.0) * velmag2 * dphi_dxk - vel(_component) * U_grad_phi;
      63             : 
      64    10967040 :   mass_term = _tauc[_qp] * mass_coeff * _strong_residuals[_qp][0];
      65             :   // Moose::out << "mass_term[" << _qp << "]=" << mass_term << std::endl;
      66             : 
      67             :   // 2.) The momentum-residual term:
      68             :   Real mom_term1 =
      69    10967040 :       U_grad_phi *
      70    10967040 :       _strong_residuals[_qp][_component + 1]; // <- momentum indices are 1,2,3, _component is 0,1,2
      71    10967040 :   Real mom_term2 = (1. - gam) * dphi_dxk * (vel * Ru);
      72    10967040 :   Real mom_term3 = vel(_component) * (_grad_test[_i][_qp] * Ru);
      73             : 
      74    10967040 :   mom_term = _taum[_qp] * (mom_term1 + mom_term2 + mom_term3);
      75             :   // Moose::out << "mom_term[" << _qp << "]=" << mom_term << std::endl;
      76             : 
      77             :   // 3.) The energy-residual term:
      78    10967040 :   energy_term = _taue[_qp] * (gam - 1.0) * dphi_dxk * _strong_residuals[_qp][4];
      79             : 
      80             :   // For printing purposes only
      81    10967040 :   Real result = mass_term + mom_term + energy_term;
      82             : 
      83    10967040 :   return result;
      84             : }
      85             : 
      86             : Real
      87     6709248 : NSSUPGMomentum::computeQpJacobian()
      88             : {
      89             :   // Set variable number for the computeJacobianHelper() function based on the
      90             :   // _component and the knowledge that this is the on-diagonal entry.
      91     6709248 :   unsigned int var_number[3] = {_rhou_var_number, _rhov_var_number, _rhow_var_number};
      92             : 
      93             :   // Call the common computeJacobianHelper() function
      94     6709248 :   return computeJacobianHelper(var_number[_component]);
      95             : }
      96             : 
      97             : Real
      98    20127744 : NSSUPGMomentum::computeQpOffDiagJacobian(unsigned int jvar)
      99             : {
     100    20127744 :   return computeJacobianHelper(jvar);
     101             : }
     102             : 
     103             : Real
     104    26836992 : NSSUPGMomentum::computeJacobianHelper(unsigned int var)
     105             : {
     106    26836992 :   if (isNSVariable(var))
     107             :   {
     108             :     // Convert the Moose numbering to canonical NS variable numbering.
     109    26836992 :     unsigned mapped_var_number = mapVarNumber(var);
     110             : 
     111             :     // Convenience vars
     112             : 
     113             :     // Velocity vector
     114    26836992 :     RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     115             : 
     116             :     // Velocity vector magnitude squared
     117             :     Real velmag2 = vel.norm_sq();
     118             : 
     119             :     // Ratio of specific heats
     120    26836992 :     const Real gam = _fp.gamma();
     121             : 
     122             :     // Shortcuts for shape function gradients at current qp.
     123    26836992 :     RealVectorValue grad_test_i = _grad_test[_i][_qp];
     124    26836992 :     RealVectorValue grad_phi_j = _grad_phi[_j][_qp];
     125             : 
     126             :     // ...
     127             : 
     128             :     // 1.) taum- and taue-proportional terms present for any variable:
     129             : 
     130             :     //
     131             :     // Art. Diffusion matrix for taum-proportional term = ( C_k + (1-gam)*C_k^T + diag(u_k) ) *
     132             :     // calA_{ell}
     133             :     //
     134             :     RealTensorValue mom_mat;
     135    26836992 :     mom_mat(0, 0) = mom_mat(1, 1) = mom_mat(2, 2) = vel(_component); // (diag(u_k)
     136    26836992 :     mom_mat += _calC[_qp][_component];                               //  + C_k
     137    26836992 :     mom_mat += (1.0 - gam) * _calC[_qp][_component].transpose();     //  + (1-gam)*C_k^T)
     138    26836992 :     mom_mat = mom_mat * _calA[_qp][mapped_var_number];               // * calA_{ell}
     139             :     Real mom_term =
     140    26836992 :         _taum[_qp] * grad_test_i * (mom_mat * grad_phi_j); // taum * grad(phi_i) * (M*grad(phi_j))
     141             : 
     142             :     //
     143             :     // Art. Diffusion matrix for taue-proportional term = (gam-1) * calE_km
     144             :     //
     145    26836992 :     RealTensorValue ene_mat = (gam - 1) * _calE[_qp][_component][mapped_var_number];
     146             :     Real ene_term =
     147    26836992 :         _taue[_qp] * grad_test_i * (ene_mat * grad_phi_j); // taue * grad(phi_i) * (M*grad(phi_j))
     148             : 
     149             :     // 2.) Terms only present if the variable is one of the momentums
     150             :     Real mass_term = 0.0;
     151             : 
     152    26836992 :     switch (mapped_var_number)
     153             :     {
     154    13418496 :       case 1:
     155             :       case 2:
     156             :       case 3:
     157             :       {
     158             :         // Variable for zero-based indexing into local matrices and vectors.
     159             :         unsigned m_local = mapped_var_number - 1;
     160             : 
     161             :         //
     162             :         // Art. Diffusion matrix for tauc-proportional term = 0.5*(gam - 1.0)*velmag2*D_km -
     163             :         // vel(_component)*C_m
     164             :         //
     165             :         RealTensorValue mass_mat;
     166    13418496 :         mass_mat(_component, m_local) = 0.5 * (gam - 1.0) * velmag2; // 0.5*(gam - 1.0)*velmag2*D_km
     167    13418496 :         mass_mat -= vel(_component) * _calC[_qp][m_local];           // vel(_component)*C_m
     168    13418496 :         mass_term = _tauc[_qp] * grad_test_i *
     169    13418496 :                     (mass_mat * grad_phi_j); // tauc * grad(phi_i) * (M*grad(phi_j))
     170             :       }
     171             :         // Nothing else to do if we are not a momentum...
     172             :     }
     173             : 
     174             :     // Sum up values and return
     175    26836992 :     return mass_term + mom_term + ene_term;
     176             :   }
     177             :   else
     178             :     return 0.0;
     179             : }