Temperature-dependent density for heat structures violates energy conservation

[moosebuild]

Imported from idaholab/thm#164 by @andrsd at Nov 30, 2020 10:33AM MST

---

Bug Description

If the user specifies a temperature-dependent density for a heat structure, energy conservation is violated. This can be shown theoretically as well as numerically. The reason stems from the fixed control volume that does not expand with density - you don't even conserve mass when changing density.

Bug Impact

Energy balances are violated in this case, which can severely affect the accuracy of the solution.

Steps to Reproduce Bug

Run this input file:

T_hs = 300
heat_flux = 1000
t = 0.001

L = 2
D_i = 0.2
thickness = 0.5

# SS 316
density = 8.0272e3
density2 = ${fparse 0.9 * density}
specific_heat_capacity = 502.1
conductivity = 16.26

R_i = ${fparse 0.5 * D_i}
D_o = ${fparse D_i + 2 * thickness}
A = ${fparse pi * D_o * L}
scale = 0.8
E_change = ${fparse scale * heat_flux * A * t}

[Functions]
  [./q_fn]
    type = ConstantFunction
    value = ${heat_flux}
  [../]
  [rho_fn]
    type = PiecewiseLinear
    x = '300 301'
    y = '${density} ${density2}'
  []
[]

[HeatStructureMaterials]
  [./hs_mat]
    type = SolidMaterialProperties
    rho = rho_fn
    cp = ${specific_heat_capacity}
    k = ${conductivity}
  [../]
[]

[Components]
  [./hs]
    type = HeatStructureCylindrical
    orientation = '0 0 1'
    position = '0 0 0'
    length = ${L}
    n_elems = 10

    inner_radius = ${R_i}
    widths = '${thickness}'
    n_part_elems = '10'
    materials = 'hs_mat'
    names = 'region'

    initial_T = ${T_hs}
  [../]

  [./heat_flux_boundary]
    type = HSBoundaryHeatFlux
    boundary = 'hs:outer'
    hs = hs
    q = q_fn
    scale_pp = bc_scale_pp
  [../]
[]

[Postprocessors]
  [./bc_scale_pp]
    type = FunctionValuePostprocessor
    function = ${scale}
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./E_hs]
    type = HeatStructureEnergyRZ
    block = 'hs:region'
    axis_dir = '0 0 1'
    axis_point = '0 0 0'
    offset = ${R_i}
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./E_hs_change]
    type = ChangeOverTimePostprocessor
    postprocessor = E_hs
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
  [./E_change_relerr]
    type = RelativeDifferencePostprocessor
    value1 = E_hs_change
    value2 = ${E_change}
    execute_on = 'INITIAL TIMESTEP_END'
  [../]
[]

[Executioner]
  type = Transient

  [./TimeIntegrator]
    type = ActuallyExplicitEuler
    solve_type = lumped
  [../]
  dt = ${t}
  num_steps = 1
  abort_on_solve_fail = true
[]

[Outputs]
  [./out]
    type = CSV
    show = 'E_change_relerr'
    execute_on = 'FINAL'
  [../]
[]

Note E_change_rel_err is O(1) instead of O(-8).

[moosebuild]

Imported from idaholab/thm#164 by @andrsd at Nov 30, 2020 10:36AM MST

---

Note that this is a property of the heat conduction equation on a fixed-volume domain; thus, it doesn't have any solution in the case of THM (fixed-volume). The recommended response is to disallow temperature-dependent density in heat structures.

By Hansel, Joshua E on 2020-11-30T17:36:01 (imported from GitLab project)