How to set temperature-dependent specific heat capacity

[Ali1990dashti]

Dear community,
I want to make specific heat capacity in one of the blocks (the aquifer block) dependent on temperature. I have a .csv containing the temperature and specific heat capacity values. I have read this discussion but still, I could not solve my issue. Implementing the specific heat capacity values as a function also works for me because I can extract a rough function out of values stored in the .csv file. The following is the input file I want to implement the temperature-dependent specific heat capacity on it:

[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
  []
  [temperature]
    initial_condition = 293
    scaling = 1E-8
  []
[]

[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 373
    boundary = injection_area
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E9
    viscosity = 1.0E-3
    density0 = 1000.0
    thermal_expansion = 0.0002
    cp = 4194
    cv = 4186
    porepressure_coefficient = 0
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy_non_res]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
    block = 'caps'
  []
  [rock_internal_energy_res]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0 # this is what I want to change
    block = 'aquifer'
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 864000
  dt = 43200
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]

The .csv file is also attached.
Thanks for your support.

values.csv

[cpgr]

I think you will need to either make a new material that takes a function instead of a single value, or perhaps modify PorousFlowMatrixInternalEnergy to take a function

[Ali1990dashti]

Thanks for dedicating so much time to my problem. Then it complains about the multiplication at the end of .C file. This is the error:

Using HIT from /home/alid/projects/moose/framework/contrib/hit
Checking if header needs updating: /home/alid/projects/moose/framework/include/base/MooseRevision.h...
Creating Unity /home/alid/projects/moose/modules/porous_flow/build/unity_src/materials_Unity.C
Compiling C++ (in opt mode) /home/alid/projects/moose/modules/porous_flow/build/unity_src/materials_Unity.C...
In file included from /home/alid/projects/moose/modules/porous_flow/build/unity_src/materials_Unity.C:37:
/home/alid/projects/moose/modules/porous_flow/src/materials/PorousFlowMatrixInternalEnergy.C: In member function 'virtual void PorousFlowMatrixInternalEnergy::initQpStatefulProperties()':
/home/alid/projects/moose/modules/porous_flow/src/materials/PorousFlowMatrixInternalEnergy.C:48:20: error: invalid use of incomplete type 'const class Function'
   48 |   _en_nodal[_qp] = _cp.value(_t) * _density * _temperature_nodal[_qp];
...
...

I tried to find out how should I write the _cp correctly in equations but I could not find a similar example. I tried _cp.value(_t) and _cp.value(_t, _q_point[_qp]) but both led to the same error (error: invalid use of incomplete type 'const class Function').

[GiudGiud]

add the top of the file, you need to add

#include "Function.h"

because right now it knows that there exists a Function class, but it doesnt know what it s in it. Only the header will explain that to the compiler

[Ali1990dashti]

Many thanks for your help and time.
I finally got rid of the errors related to the compiler.
In the input file, I used your suggestion for extracting the function out of the .csv file, but how does MOOSE know that the relation is between specific heat capacity and temperature?
I have seen examples of MOOSE based on coordinates (x, y, z) and time (t), but how can say here that specific heat capacity should be a function of temperature?

[GiudGiud]

_cp.value(_t, _q_point[_qp]) isnt right if you want to use the function to describe the temperature dependence

you need to do something like

_cp.value(temperature[_qp], Point(0,0,0))

the point should not matter basically, there s only a 1D dependence

[Ali1990dashti]

Thanks a lot @GiudGiud , it worked with _cp.value(_temperature_nodal[_qp], Point(0,0,0)). This is the final veryion of the .C file:

#include "PorousFlowMatrixInternalEnergy.h"
#include "Function.h"

registerMooseObject("PorousFlowApp", PorousFlowMatrixInternalEnergy);

InputParameters
PorousFlowMatrixInternalEnergy::validParams()
{
  InputParameters params = PorousFlowMaterialVectorBase::validParams();
  params.addRequiredParam<FunctionName>("specific_heat_capacity",
                                "Specific heat capacity of the rock grains (J/kg/K).");
  params.addRequiredParam<Real>("density", "Density of the rock grains");
  params.set<bool>("at_nodes") = true;
  params.addPrivateParam<std::string>("pf_material_type", "matrix_internal_energy");
  params.addClassDescription("This Material calculates the internal energy of solid rock grains, "
                             "which is specific_heat_capacity * density * temperature.  Kernels "
                             "multiply this by (1 - porosity) to find the energy density of the "
                             "porous rock in a rock-fluid system");
  return params;
}

PorousFlowMatrixInternalEnergy::PorousFlowMatrixInternalEnergy(const InputParameters & parameters)
  : PorousFlowMaterialVectorBase(parameters),
    _cp(getFunction("specific_heat_capacity")),
    _density(getParam<Real>("density")),
    _temperature_nodal(getMaterialProperty<Real>("PorousFlow_temperature_nodal")),
    _dtemperature_nodal_dvar(
        getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")),
    _en_nodal(declareProperty<Real>("PorousFlow_matrix_internal_energy_nodal")),
    _den_nodal_dvar(
        declareProperty<std::vector<Real>>("dPorousFlow_matrix_internal_energy_nodal_dvar"))
{
  if (_nodal_material != true)
    mooseError("PorousFlowMatrixInternalEnergy classes are only defined for at_nodes = true");
}

void
PorousFlowMatrixInternalEnergy::initQpStatefulProperties()
{
  _en_nodal[_qp] = _cp.value(_temperature_nodal[_qp], Point(0,0,0)) * _density * _temperature_nodal[_qp];
}

void
PorousFlowMatrixInternalEnergy::computeQpProperties()
{
  _en_nodal[_qp] = _cp.value(_temperature_nodal[_qp], Point(0,0,0)) * _density * _temperature_nodal[_qp];

  _den_nodal_dvar[_qp].assign(_num_var, 0.0);
  for (unsigned v = 0; v < _num_var; ++v)
    _den_nodal_dvar[_qp][v] = _cp.value(_temperature_nodal[_qp], Point(0,0,0)) * _density * _dtemperature_nodal_dvar[_qp][v];
}