Heat conduction tests descriptions

Simplifying PorousFlow's heat equation

PorousFlow heat conduction is governed by the equation (1) Here is the energy per unit volume in the rock-fluid system, and is the heat flux. In the PorousFlow module when there is no adsorbed species. Here is the temperature, is the rock porosity, and is the saturation of phase . The remainder of the notation is described in the next paragraph. When studying problems involving heat conduction (with no fluid convection) where is the tensorial thermal conductivity of the rock-fluid system.

The tests described in this page use the following simple forms for each term

  • , which is the rock-grain density (that is, the density of rock with zero porespace), measured in kg.m, is assumed constant in the current implementation of the PorousFlow module.

  • , which is the rock-grain specific heat capacity, measured in J.kg.K, is assumed constant in the current implementation of the PorousFlow module.

  • , which is the density of fluid phase , is assumed in here to be a function of the fluid pressure only. This is so Eq. (1) may be easily solved, but more general forms are allowed in the PorousFlow module.

  • , which is the specific internal energy of the fluid phase , and is measured in J.kg, is assumed here to be , where is the fluid's specific heat capacity at constant volume. This specific heat capacity is assumed constant, so that Eq. (1) may be easily solved — more general forms are allowed in the PorousFlow module.

  • is assumed to vary between and , depending on the aqueous saturation: , where is the aqueous saturation, and is a positive user-defined exponent. More general forms may be easily accommodated in the PorousFlow module, but to date none have been coded.

Under these conditions, Eq. (1) becomes (2) The tensor is For constant saturation and porepressure, is also constant.

Testing heat conduction in 1D

Consider the one-dimensional case where the spatial dimension is the semi-infinite line . Suppose that initially the temperature is constant, so that Then apply a fixed-pressure Dirichlet boundary condition at so that The solution of the above differential equation is well known to be (3) where Erf is the error function.

This is verified by using the following tests on a line of 10 elements.

  • A transient analysis with no fluids. The parameters chosen are , , and is chosen, so that .

# 0phase heat conduction.
# apply a boundary condition of T=300 to a bar that
# is initially at T=200, and observe the expected
# error-function response
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = 0
  xmax = 100
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [temp]
    initial_condition = 200
  []
[]

[Kernels]
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [heat_conduction]
    type = PorousFlowHeatConduction
    variable = temp
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp'
    number_fluid_phases = 0
    number_fluid_components = 0
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2.2 0 0  0 0 0  0 0 0'
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 2.2
    density = 0.5
  []
[]

[BCs]
  [left]
    type = DirichletBC
    boundary = left
    value = 300
    variable = temp
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E1
  end_time = 1E2
[]

[Postprocessors]
  [t000]
    type = PointValue
    variable = temp
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [t010]
    type = PointValue
    variable = temp
    point = '10 0 0'
    execute_on = 'initial timestep_end'
  []
  [t020]
    type = PointValue
    variable = temp
    point = '20 0 0'
    execute_on = 'initial timestep_end'
  []
  [t030]
    type = PointValue
    variable = temp
    point = '30 0 0'
    execute_on = 'initial timestep_end'
  []
  [t040]
    type = PointValue
    variable = temp
    point = '40 0 0'
    execute_on = 'initial timestep_end'
  []
  [t050]
    type = PointValue
    variable = temp
    point = '50 0 0'
    execute_on = 'initial timestep_end'
  []
  [t060]
    type = PointValue
    variable = temp
    point = '60 0 0'
    execute_on = 'initial timestep_end'
  []
  [t070]
    type = PointValue
    variable = temp
    point = '70 0 0'
    execute_on = 'initial timestep_end'
  []
  [t080]
    type = PointValue
    variable = temp
    point = '80 0 0'
    execute_on = 'initial timestep_end'
  []
  [t090]
    type = PointValue
    variable = temp
    point = '90 0 0'
    execute_on = 'initial timestep_end'
  []
  [t100]
    type = PointValue
    variable = temp
    point = '100 0 0'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = no_fluid
  [csv]
    type = CSV
  []
  exodus = false
[]
(modules/porous_flow/test/tests/heat_conduction/no_fluid.i)

  • A transient analysis with 2 fluid phases. The parameters chosen are , and , so that . , , , , , , and . With these parameters, .

# 2phase heat conduction, with saturation fixed at 0.5
# apply a boundary condition of T=300 to a bar that
# is initially at T=200, and observe the expected
# error-function response
[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 10
  xmin = 0
  xmax = 100
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [phase0_porepressure]
    initial_condition = 0
  []
  [phase1_saturation]
    initial_condition = 0.5
  []
  [temp]
    initial_condition = 200
  []
[]

[Kernels]
  [dummy_p0]
    type = TimeDerivative
    variable = phase0_porepressure
  []
  [dummy_s1]
    type = TimeDerivative
    variable = phase1_saturation
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = temp
  []
  [heat_conduction]
    type = PorousFlowHeatConduction
    variable = temp
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'temp phase0_porepressure phase1_saturation'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
[]

[FluidProperties]
  [simple_fluid0]
    type = SimpleFluidProperties
    bulk_modulus = 1.5
    density0 = 0.4
    thermal_expansion = 0
    cv = 1
  []
  [simple_fluid1]
    type = SimpleFluidProperties
    bulk_modulus = 0.5
    density0 = 0.3
    thermal_expansion = 0
    cv = 2
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.3 0 0  0 0 0  0 0 0'
    wet_thermal_conductivity = '1.7 0 0  0 0 0  0 0 0'
    exponent = 1.0
    aqueous_phase_number = 1
  []
  [ppss]
    type = PorousFlow2PhasePS
    phase0_porepressure = phase0_porepressure
    phase1_saturation = phase1_saturation
    capillary_pressure = pc
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.8
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1.0
    density = 0.25
  []
  [simple_fluid0]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid0
    phase = 0
  []
  [simple_fluid1]
    type = PorousFlowSingleComponentFluid
    fp = simple_fluid1
    phase = 1
  []
[]

[BCs]
  [left]
    type = DirichletBC
    boundary = left
    value = 300
    variable = temp
  []
[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1E1
  end_time = 1E2
[]

[Postprocessors]
  [t000]
    type = PointValue
    variable = temp
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [t010]
    type = PointValue
    variable = temp
    point = '10 0 0'
    execute_on = 'initial timestep_end'
  []
  [t020]
    type = PointValue
    variable = temp
    point = '20 0 0'
    execute_on = 'initial timestep_end'
  []
  [t030]
    type = PointValue
    variable = temp
    point = '30 0 0'
    execute_on = 'initial timestep_end'
  []
  [t040]
    type = PointValue
    variable = temp
    point = '40 0 0'
    execute_on = 'initial timestep_end'
  []
  [t050]
    type = PointValue
    variable = temp
    point = '50 0 0'
    execute_on = 'initial timestep_end'
  []
  [t060]
    type = PointValue
    variable = temp
    point = '60 0 0'
    execute_on = 'initial timestep_end'
  []
  [t070]
    type = PointValue
    variable = temp
    point = '70 0 0'
    execute_on = 'initial timestep_end'
  []
  [t080]
    type = PointValue
    variable = temp
    point = '80 0 0'
    execute_on = 'initial timestep_end'
  []
  [t090]
    type = PointValue
    variable = temp
    point = '90 0 0'
    execute_on = 'initial timestep_end'
  []
  [t100]
    type = PointValue
    variable = temp
    point = '100 0 0'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = two_phase
  [csv]
    type = CSV
  []
  exodus = false
[]
(modules/porous_flow/test/tests/heat_conduction/two_phase.i)

PorousFlow yields the expected answer, as shown in Figure 1.

Figure 1: Comparison between the MOOSE result and the exact analytic expression given by Eq. (3).