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 "INSFEFluidEnergyKernel.h"
      11             : 
      12             : registerMooseObject("NavierStokesApp", INSFEFluidEnergyKernel);
      13             : registerMooseObjectRenamed("NavierStokesApp",
      14             :                            MDFluidEnergyKernel,
      15             :                            "02/01/2024 00:00",
      16             :                            INSFEFluidEnergyKernel);
      17             : 
      18             : InputParameters
      19         123 : INSFEFluidEnergyKernel::validParams()
      20             : {
      21         123 :   InputParameters params = INSFEFluidKernelStabilization::validParams();
      22         123 :   params.addClassDescription("Adds advection, diffusion, and heat source terms to energy equation, "
      23             :                              "potentially with stabilization");
      24         246 :   params.addParam<bool>("conservative_form", false, "if conservative form is used");
      25         246 :   params.addCoupledVar("porosity_elem", "Element averaged porosity");
      26         246 :   params.addCoupledVar("power_density", "volumetric heat source");
      27         246 :   params.addCoupledVar("pke_power_var", "Normalized power from PKE");
      28         246 :   params.addParam<FunctionName>("power_shape_function", "Function that defines power profile");
      29             : 
      30         123 :   return params;
      31           0 : }
      32             : 
      33          66 : INSFEFluidEnergyKernel::INSFEFluidEnergyKernel(const InputParameters & parameters)
      34             :   : INSFEFluidKernelStabilization(parameters),
      35          66 :     _conservative_form(getParam<bool>("conservative_form")),
      36         132 :     _k_elem(getMaterialProperty<Real>("k_fluid_elem")),
      37         132 :     _cp(getMaterialProperty<Real>("cp_fluid")),
      38         132 :     _has_porosity_elem(isParamValid("porosity_elem")),
      39          66 :     _porosity_elem(_has_porosity_elem ? coupledValue("porosity_elem")
      40         132 :                                       : (_has_porosity ? coupledValue("porosity") : _zero)),
      41         132 :     _has_qv(isParamValid("power_density")),
      42          66 :     _qv(_has_qv ? coupledValue("power_density") : _zero),
      43         132 :     _has_pke(isParamValid("pke_power_var")),
      44          66 :     _pke_qv(_has_pke ? coupledScalarValue("pke_power_var") : _zero),
      45          66 :     _power_shape_function(
      46         198 :         isParamValid("power_shape_function") ? &getFunction("power_shape_function") : NULL)
      47             : {
      48             :   // Input sanity check: cannot have both 'power_density' and 'pke_power_var'
      49          66 :   if (_has_qv && _has_pke)
      50           0 :     mooseError(
      51             :         "'power_density' and 'pke_power_var' cannot be both provided in 'INSFEFluidEnergyKernel'.");
      52             : 
      53          66 :   if (_has_pke && (_power_shape_function == NULL))
      54           0 :     mooseError("'power_shape_function' is required if 'pke_power_var' is provided in "
      55             :                "'INSFEFluidEnergyKernel'.");
      56          66 : }
      57             : 
      58             : Real
      59    29772160 : INSFEFluidEnergyKernel::computeQpResidual()
      60             : {
      61    29772160 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
      62    29772160 :   Real porosity_elem = (_has_porosity_elem || _has_porosity) ? _porosity_elem[_qp] : 1.0;
      63    29772160 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      64    29772160 :   Real enthalpy = _eos.h_from_p_T(_pressure[_qp], _u[_qp]);
      65             : 
      66             :   // Residual weak form terms: convection + diffusion + volumetric heat
      67    29772160 :   Real convective_part = _conservative_form
      68    29772160 :                              ? -_rho[_qp] * enthalpy * vec_vel * _grad_test[_i][_qp]
      69     6700160 :                              : _rho[_qp] * _cp[_qp] * vec_vel * _grad_u[_qp] * _test[_i][_qp];
      70             : 
      71             :   // heat source, one of the two terms, or both two terms will be zero
      72    29772160 :   Real heat_source_part = -_qv[_qp] * _test[_i][_qp];
      73    29772160 :   if (_has_pke)
      74           0 :     heat_source_part +=
      75           0 :         -_pke_qv[0] * _power_shape_function->value(_t, _q_point[_qp]) * _test[_i][_qp];
      76             : 
      77    29772160 :   Real diffusion_part = porosity_elem * _k_elem[_qp] * _grad_u[_qp] * _grad_test[_i][_qp];
      78             : 
      79             :   // Residual strong form terms for SUPG contributions
      80    29772160 :   Real transient_supg = _bTransient ? porosity * _rho[_qp] * _cp[_qp] * _u_dot[_qp] : 0;
      81    29772160 :   Real convection_supg = _rho[_qp] * _cp[_qp] * vec_vel * _grad_u[_qp];
      82    29772160 :   Real heat_source_supg = -_qv[_qp];
      83    29772160 :   if (_has_pke)
      84           0 :     heat_source_supg += -_pke_qv[0] * _power_shape_function->value(_t, _q_point[_qp]);
      85             : 
      86    29772160 :   Real diffusion_supg = -porosity_elem * _k_elem[_qp] * _second_u[_qp].tr();
      87             : 
      88             :   // assemble residuals
      89    29772160 :   Real normal_part = convective_part + heat_source_part + diffusion_part;
      90    29772160 :   Real supg_part = (transient_supg + convection_supg + heat_source_supg + diffusion_supg) *
      91    29772160 :                    _taue[_qp] * _vel_elem * _grad_test[_i][_qp];
      92             : 
      93    29772160 :   return normal_part + supg_part;
      94             : }
      95             : 
      96             : Real
      97     2513920 : INSFEFluidEnergyKernel::computeQpJacobian()
      98             : {
      99     2513920 :   Real porosity = _has_porosity ? _porosity[_qp] : 1.0;
     100     2513920 :   Real porosity_elem = (_has_porosity_elem || _has_porosity) ? _porosity_elem[_qp] : 1.0;
     101             :   Real rho, drho_dp, drho_dT;
     102     2513920 :   _eos.rho_from_p_T(_pressure[_qp], _u[_qp], rho, drho_dp, drho_dT);
     103     2513920 :   Real enthalpy = _eos.h_from_p_T(_pressure[_qp], _u[_qp]);
     104     2513920 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     105             : 
     106             :   // convection part
     107     2513920 :   Real convection_part = _conservative_form
     108     2513920 :                              ? -(_rho[_qp] * _cp[_qp] + enthalpy * drho_dT) * vec_vel *
     109      486400 :                                    _phi[_j][_qp] * _grad_test[_i][_qp]
     110     2027520 :                              : _cp[_qp] * vec_vel *
     111     2027520 :                                    (drho_dT * _grad_u[_qp] + _rho[_qp] * _grad_phi[_j][_qp]) *
     112     2027520 :                                    _test[_i][_qp];
     113             : 
     114             :   // diffusion part
     115     2513920 :   Real diffusion_part = porosity_elem * _k_elem[_qp] * _grad_phi[_j][_qp] * _grad_test[_i][_qp];
     116             : 
     117     2513920 :   Real normal_part = convection_part + diffusion_part;
     118             : 
     119             :   // SUPG contributions (only transient and convection part, ignore diffusion part)
     120             :   Real transient_supg =
     121     2513920 :       _bTransient ? porosity * _rho[_qp] * _cp[_qp] * _du_dot_du[_qp] * _phi[_j][_qp] : 0;
     122     2513920 :   Real convection_supg = _rho[_qp] * _cp[_qp] * vec_vel * _grad_phi[_j][_qp];
     123             :   Real supg_part =
     124     2513920 :       _taue[_qp] * _vel_elem * _grad_test[_i][_qp] * (transient_supg + convection_supg);
     125             : 
     126             :   // Assembly total residuals
     127     2513920 :   return normal_part + supg_part;
     128             : }
     129             : 
     130             : Real
     131     7541760 : INSFEFluidEnergyKernel::computeQpOffDiagJacobian(unsigned int jvar)
     132             : {
     133     7541760 :   unsigned m = this->mapVarNumber(jvar);
     134             : 
     135     7541760 :   RealVectorValue vec_vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     136     7541760 :   Real enthalpy = _eos.h_from_p_T(_pressure[_qp], _u[_qp]);
     137             : 
     138             :   Real jac = 0.;
     139     7541760 :   switch (m)
     140             :   {
     141     5027840 :     case 1:
     142             :     case 2:
     143             :     case 3:
     144             :       // convection weak form part
     145     5027840 :       if (_conservative_form)
     146      972800 :         jac -= _rho[_qp] * _phi[_j][_qp] * enthalpy * _grad_test[_i][_qp](m - 1);
     147             :       else
     148     4055040 :         jac += _rho[_qp] * _cp[_qp] * _phi[_j][_qp] * _grad_u[_qp](m - 1) * _test[_i][_qp];
     149             :       // convection supg part
     150     5027840 :       jac += _rho[_qp] * _cp[_qp] * _phi[_j][_qp] * _grad_u[_qp](m - 1) * _taue[_qp] * _vel_elem *
     151     5027840 :              _grad_test[_i][_qp];
     152     5027840 :       break;
     153             : 
     154             :     default:
     155             :       jac = 0.0;
     156             :   }
     157             : 
     158     7541760 :   return jac;
     159             : }