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 "NSEnergyViscousBC.h"
      13             : 
      14             : registerMooseObject("NavierStokesApp", NSEnergyViscousBC);
      15             : 
      16             : InputParameters
      17           0 : NSEnergyViscousBC::validParams()
      18             : {
      19           0 :   InputParameters params = NSIntegratedBC::validParams();
      20           0 :   params.addRequiredCoupledVar(NS::temperature, "temperature");
      21           0 :   return params;
      22           0 : }
      23             : 
      24           0 : NSEnergyViscousBC::NSEnergyViscousBC(const InputParameters & parameters)
      25             :   : NSIntegratedBC(parameters),
      26           0 :     _grad_temperature(coupledGradient(NS::temperature)),
      27           0 :     _thermal_conductivity(getMaterialProperty<Real>("thermal_conductivity")),
      28             :     // Viscous stress tensor derivative computing object
      29             :     _vst_derivs(*this),
      30             :     // Temperature derivative computing object
      31           0 :     _temp_derivs(*this)
      32             : {
      33             :   // Store pointers to all variable gradients in a single vector.
      34           0 :   _gradU.resize(5);
      35           0 :   _gradU[0] = &_grad_rho;
      36           0 :   _gradU[1] = &_grad_rho_u;
      37           0 :   _gradU[2] = &_grad_rho_v;
      38           0 :   _gradU[3] = &_grad_rho_w;
      39           0 :   _gradU[4] = &_grad_rho_et;
      40           0 : }
      41             : 
      42             : Real
      43           0 : NSEnergyViscousBC::computeQpResidual()
      44             : {
      45             :   // n . (- k*grad(T) - tau*u) v
      46             : 
      47             :   // Velocity vector object
      48           0 :   RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      49             : 
      50             :   // k*grad(T)
      51           0 :   RealVectorValue thermal_vec = _thermal_conductivity[_qp] * _grad_temperature[_qp];
      52             : 
      53             :   // tau*u
      54           0 :   RealVectorValue visc_vec = _viscous_stress_tensor[_qp] * vel;
      55             : 
      56             :   // Add everything up, dot with normal, hit with test function.
      57           0 :   return ((-thermal_vec - visc_vec) * _normals[_qp]) * _test[_i][_qp];
      58             : }
      59             : 
      60             : Real
      61           0 : NSEnergyViscousBC::computeQpJacobian()
      62             : {
      63             :   // See notes for this term, involves temperature Hessian
      64             :   Real thermal_term = 0.0;
      65             : 
      66           0 :   for (unsigned int ell = 0; ell < LIBMESH_DIM; ++ell)
      67             :   {
      68             :     Real intermediate_result = 0.;
      69             : 
      70             :     // The temperature Hessian contribution
      71           0 :     for (unsigned n = 0; n < 5; ++n)
      72           0 :       intermediate_result += _temp_derivs.get_hess(/*m=*/4, n) * (*_gradU[n])[_qp](ell);
      73             : 
      74             :     // Hit Hessian contribution with test function
      75           0 :     intermediate_result *= _phi[_j][_qp];
      76             : 
      77             :     // Add in the temperature gradient contribution
      78           0 :     intermediate_result += _temp_derivs.get_grad(/*rhoE=*/4) * _grad_phi[_j][_qp](ell);
      79             : 
      80             :     // Hit the result with the normal component, accumulate in thermal_term
      81           0 :     thermal_term += intermediate_result * _normals[_qp](ell);
      82             :   }
      83             : 
      84             :   // Hit thermal_term with thermal conductivity
      85           0 :   thermal_term *= _thermal_conductivity[_qp];
      86             : 
      87           0 :   return (-thermal_term) * _test[_i][_qp];
      88             : }
      89             : 
      90             : Real
      91           0 : NSEnergyViscousBC::computeQpOffDiagJacobian(unsigned jvar)
      92             : {
      93           0 :   if (isNSVariable(jvar))
      94             :   {
      95             :     // Note: This function requires both _vst_derivs *and* _temp_derivs
      96             : 
      97             :     // Convenience variables
      98           0 :     const RealTensorValue & tau = _viscous_stress_tensor[_qp];
      99             : 
     100           0 :     Real rho = _rho[_qp];
     101           0 :     Real phij = _phi[_j][_qp];
     102           0 :     RealVectorValue U(_rho_u[_qp], _rho_v[_qp], _rho_w[_qp]);
     103             : 
     104             :     // Map jvar into the variable m for our problem, regardless of
     105             :     // how Moose has numbered things.
     106           0 :     unsigned m = mapVarNumber(jvar);
     107             : 
     108             :     //
     109             :     // 1.) Thermal term derivatives
     110             :     //
     111             : 
     112             :     // See notes for this term, involves temperature Hessian
     113             :     Real thermal_term = 0.;
     114             : 
     115           0 :     for (unsigned ell = 0; ell < LIBMESH_DIM; ++ell)
     116             :     {
     117             :       Real intermediate_result = 0.;
     118             : 
     119             :       // The temperature Hessian contribution
     120           0 :       for (unsigned n = 0; n < 5; ++n)
     121           0 :         intermediate_result += _temp_derivs.get_hess(m, n) * (*_gradU[n])[_qp](ell);
     122             : 
     123             :       // Hit Hessian contribution with test function
     124           0 :       intermediate_result *= _phi[_j][_qp];
     125             : 
     126             :       // Add in the temperature gradient contribution
     127           0 :       intermediate_result += _temp_derivs.get_grad(m) * _grad_phi[_j][_qp](ell);
     128             : 
     129             :       // Hit the result with the normal component, accumulate in thermal_term
     130           0 :       thermal_term += intermediate_result * _normals[_qp](ell);
     131             :     }
     132             : 
     133             :     // Hit thermal_term with thermal conductivity
     134           0 :     thermal_term *= _thermal_conductivity[_qp];
     135             : 
     136             :     //
     137             :     // 2.) Viscous term derivatives
     138             :     //
     139             : 
     140             :     // Compute viscous term derivatives
     141             :     Real visc_term = 0.;
     142             : 
     143           0 :     switch (m)
     144             :     {
     145             :       case 0: // density
     146             :       {
     147             :         // Loop over k and ell as in the notes...
     148           0 :         for (const auto k : make_range(Moose::dim))
     149             :         {
     150             :           Real intermediate_value = 0.0;
     151           0 :           for (unsigned int ell = 0; ell < LIBMESH_DIM; ++ell)
     152           0 :             intermediate_value +=
     153           0 :                 (U(ell) / rho * (-tau(k, ell) * phij / rho + _vst_derivs.dtau(k, ell, m)));
     154             : 
     155             :           // Hit accumulated value with normal component k.  We will multiply by test function at
     156             :           // the end of this routine...
     157           0 :           visc_term += intermediate_value * _normals[_qp](k);
     158             :         } // end for k
     159             : 
     160             :         break;
     161             :       } // end case 0
     162             : 
     163           0 :       case 1:
     164             :       case 2:
     165             :       case 3: // momentums
     166             :       {
     167             :         // Map m -> 0,1,2 as usual...
     168           0 :         unsigned int m_local = m - 1;
     169             : 
     170             :         // Loop over k and ell as in the notes...
     171           0 :         for (const auto k : make_range(Moose::dim))
     172             :         {
     173           0 :           Real intermediate_value = tau(k, m_local) * phij / rho;
     174             : 
     175           0 :           for (unsigned int ell = 0; ell < LIBMESH_DIM; ++ell)
     176           0 :             intermediate_value += _vst_derivs.dtau(k, ell, m) * U(ell) /
     177             :                                   rho; // Note: pass 'm' to dtau, it will convert it internally
     178             : 
     179             :           // Hit accumulated value with normal component k.
     180           0 :           visc_term += intermediate_value * _normals[_qp](k);
     181             :         } // end for k
     182             : 
     183             :         break;
     184             :       } // end case 1,2,3
     185             : 
     186           0 :       case 4: // energy
     187           0 :         mooseError("Shouldn't get here, this is the on-diagonal entry!");
     188             :         break;
     189             : 
     190           0 :       default:
     191           0 :         mooseError("Invalid m value.");
     192             :         break;
     193             :     }
     194             : 
     195             :     // Finally, sum up the different contributions (with appropriate
     196             :     // sign) multiply by the test function, and return.
     197           0 :     return (-thermal_term - visc_term) * _test[_i][_qp];
     198             :   }
     199             :   else
     200             :     return 0.0;
     201             : }