Merge pull request idaholab#26648 from tophmatthews/spheat_dt_26647

Add spheat and density deriv to heat conduction

modules/heat_transfer/include/kernels/HeatConductionTimeDerivative.h

@@ -39,5 +39,7 @@ class HeatConductionTimeDerivative : public TimeDerivative
   virtual Real computeQpJacobian();
 
   const MaterialProperty<Real> & _specific_heat;
+  const MaterialProperty<Real> * const _specific_heat_dT;
   const MaterialProperty<Real> & _density;
+  const MaterialProperty<Real> * const _density_dT;
 };

modules/heat_transfer/include/materials/HeatConductionMaterial.h

@@ -41,6 +41,7 @@ class HeatConductionMaterialTempl : public Material
 
   GenericMaterialProperty<Real, is_ad> & _specific_heat;
   const Function * const _specific_heat_temperature_function;
+  GenericMaterialProperty<Real, is_ad> * const _specific_heat_dT;
 
   /// Minimum temperature, below which temperature is "clipped" before evaluating functions
   const Real * _min_T;

modules/heat_transfer/src/kernels/HeatConductionTimeDerivative.C

@@ -24,7 +24,11 @@ HeatConductionTimeDerivative::validParams()
   params.set<bool>("use_displaced_mesh") = true;
 
   params.addParam<MaterialPropertyName>(
-      "specific_heat", "specific_heat", "Property name of the specific heat material property");
+      "specific_heat", "specific_heat", "Name of the specific heat material property");
+  params.addParam<MaterialPropertyName>(
+      "specific_heat_dT",
+      "Name of the material property for the derivative of the specific heat with respect "
+      "to the variable.");
 
   /**
    * We would like to rename this input parameter to 'density' but gratuitous use of
@@ -33,13 +37,21 @@ HeatConductionTimeDerivative::validParams()
    */
   params.addParam<MaterialPropertyName>(
       "density_name", "density", "Property name of the density material property");
+  params.addParam<MaterialPropertyName>(
+      "density_name_dT",
+      "Name of material property for the derivative of the density with respect to the variable.");
   return params;
 }
 
 HeatConductionTimeDerivative::HeatConductionTimeDerivative(const InputParameters & parameters)
   : TimeDerivative(parameters),
     _specific_heat(getMaterialProperty<Real>("specific_heat")),
-    _density(getMaterialProperty<Real>("density_name"))
+    _specific_heat_dT(isParamValid("specific_heat_dT")
+                          ? &getMaterialProperty<Real>("specific_heat_dT")
+                          : nullptr),
+    _density(getMaterialProperty<Real>("density_name")),
+    _density_dT(isParamValid("density_name_dT") ? &getMaterialProperty<Real>("density_name_dT")
+                                                : nullptr)
 {
 }
 
@@ -52,5 +64,12 @@ HeatConductionTimeDerivative::computeQpResidual()
 Real
 HeatConductionTimeDerivative::computeQpJacobian()
 {
-  return _specific_heat[_qp] * _density[_qp] * TimeDerivative::computeQpJacobian();
+  auto jac = _specific_heat[_qp] * _density[_qp] * TimeDerivative::computeQpJacobian();
+  if (_specific_heat_dT)
+    jac += (*_specific_heat_dT)[_qp] * _density[_qp] * _phi[_j][_qp] *
+           TimeDerivative::computeQpResidual();
+  if (_density_dT)
+    jac += _specific_heat[_qp] * (*_density_dT)[_qp] * _phi[_j][_qp] *
+           TimeDerivative::computeQpResidual();
+  return jac;
 }

modules/heat_transfer/src/materials/HeatConductionMaterial.C

@@ -60,6 +60,9 @@ HeatConductionMaterialTempl<is_ad>::HeatConductionMaterialTempl(const InputParam
         getParam<FunctionName>("specific_heat_temperature_function") != ""
             ? &getFunction("specific_heat_temperature_function")
             : nullptr),
+    _specific_heat_dT(is_ad && !_specific_heat_temperature_function
+                          ? nullptr
+                          : &declareGenericProperty<Real, is_ad>("specific_heat_dT")),
     _min_T(isParamValid("min_T") ? &getParam<Real>("min_T") : nullptr)
 {
   if (_thermal_conductivity_temperature_function && !_has_temp)
@@ -86,16 +89,15 @@ HeatConductionMaterialTempl<is_ad>::computeQpProperties()
 
   if (_has_temp && _min_T)
   {
-    if (qp_temperature < *_min_T)
+    if (_temperature[_qp] < *_min_T)
     {
-      std::stringstream msg;
-      msg << "WARNING:  In HeatConductionMaterial:  negative temperature!\n"
-          << "\tResetting to 'min_T'.\n"
-          << "\t_qp: " << _qp << "\n"
-          << "\ttemp: " << qp_temperature << "\n"
-          << "\telem: " << _current_elem->id() << "\n"
-          << "\tproc: " << processor_id() << "\n";
-      mooseWarning(msg.str());
+      mooseWarning("Temperature (",
+                   _temperature[_qp],
+                   ") is below min_T (",
+                   *_min_T,
+                   ") at ",
+                   _q_point[_qp],
+                   ". \nUsing min_T instead of temperature.");
       qp_temperature = *_min_T;
     }
   }
@@ -113,7 +115,12 @@ HeatConductionMaterialTempl<is_ad>::computeQpProperties()
   }
 
   if (_specific_heat_temperature_function)
+  {
     _specific_heat[_qp] = _specific_heat_temperature_function->value(qp_temperature);
+    if (_specific_heat_dT)
+      (*_specific_heat_dT)[_qp] = _specific_heat_temperature_function->timeDerivative(
+          MetaPhysicL::raw_value(qp_temperature));
+  }
   else
     _specific_heat[_qp] = _my_specific_heat;
 }

modules/heat_transfer/test/tests/transient_heat/gold/transient_heat_derivatives_out.csv

@@ -0,0 +1,3 @@
+time,avg
+0,0
+0.1,1.810271319407

modules/heat_transfer/test/tests/transient_heat/tests

@@ -1,10 +1,29 @@
 [Tests]
-  issues = '#7759'
-  design = 'SpecificHeatConductionTimeDerivative.md'
-  [./test]
+  [transient_heat]
     type = 'Exodiff'
     input = 'transient_heat.i'
     exodiff = 'transient_heat_out.e'
     requirement = 'The system shall compute the time derivative term of the heat equation'
-  [../]
+    issues = '#7759'
+    design = 'SpecificHeatConductionTimeDerivative.md'
+  []
+  [transient_heat_derivatives]
+    type = 'CSVDiff'
+    input = 'transient_heat_derivatives.i'
+    csvdiff = 'transient_heat_derivatives_out.csv'
+    requirement = 'The system shall compute the time derivative term of the heat equation and '
+                  'utilize the derivative of the specific heat, thermal conductivity, and density '
+                  'to solve a heat conduction problem'
+    issues = '#26647'
+    design = 'HeatConductionTimeDerivative.md'
+  []
+  [transient_heat_t_limit]
+    type = 'RunException'
+    input = 'transient_heat_derivatives.i'
+    expect_err = 'is below min_T'
+    cli_args = 'Materials/constant/min_T=10'
+    requirement = 'The system shall use a lower temperature limit to compute the specific heat'
+    issues = '#26647'
+    design = 'HeatConductionTimeDerivative.md'
+  []
 []

modules/heat_transfer/test/tests/transient_heat/transient_heat_derivatives.i

@@ -0,0 +1,98 @@
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 10
+    ny = 10
+  []
+[]
+
+[Variables]
+  [temp]
+    order = FIRST
+    family = LAGRANGE
+    initial_condition = 2
+  []
+[]
+
+[Kernels]
+  [heat]
+    type = HeatConduction
+    variable = temp
+  []
+
+  [ie]
+    type = HeatConductionTimeDerivative
+    variable = temp
+    specific_heat_dT = specific_heat_dT
+    density_name_dT = density_dT
+  []
+[]
+
+[Functions]
+  [spheat]
+    type = ParsedFunction
+    expression = 't^4'
+  []
+  [thcond]
+    type = ParsedFunction
+    expression = 'exp(t)'
+  []
+[]
+
+[BCs]
+  [bottom]
+    type = DirichletBC
+    variable = temp
+    boundary = 1
+    value = 4
+  []
+
+  [top]
+    type = DirichletBC
+    variable = temp
+    boundary = 2
+    value = 1
+  []
+[]
+
+[Materials]
+  [constant]
+    type = HeatConductionMaterial
+    thermal_conductivity_temperature_function = thcond
+    specific_heat_temperature_function = spheat
+    temp = temp
+  []
+  [density]
+    type = ParsedMaterial
+    property_name = density
+    coupled_variables = temp
+    expression = 'temp^3 + 2/temp'
+  []
+  [density_dT]
+    type = ParsedMaterial
+    property_name = density_dT
+    coupled_variables = temp
+    expression = '3 * temp^2 - 2/temp/temp'
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = NEWTON
+  num_steps = 1
+  dt = .1
+  nl_max_its = 10
+  dtmin = .1
+[]
+
+[Postprocessors]
+  [avg]
+    type = ElementAverageValue
+    variable = temp
+  []
+[]
+
+[Outputs]
+  csv = true
+[]