ver-1ib

Species Equilibration Problem in Ratedep Condition with Unequal Starting Pressures

General Case Description

This verification problem is taken from Ambrosek and Longhurst (2008) and builds on ver-1ia. The configuration and modeling parameters are the same as in ver-1ia, except that, in the current case, the starting pressures for A $var element = document.getElementById("moose-equation-1930b87f-67b6-4a24-8076-c5124de4ddf4");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-cc9e4a56-ce9b-4bd6-b978-c98dd2bb2287");katex.render("_2", element, {displayMode:false,throwOnError:false});$ are not equal. The case is simulated in (test/tests/ver-1ia/ver-1ia.i), but the starting pressures of A $var element = document.getElementById("moose-equation-53be190a-da93-4cc3-a8ad-fa7d9120634c");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-da6a3c17-0105-4d70-b320-52c8d6da7894");katex.render("_2", element, {displayMode:false,throwOnError:false});$ are $var element = document.getElementById("moose-equation-7cf3f7ee-c398-4a27-88bc-2df7b6609020");katex.render("1 \\times 10^4", element, {displayMode:false,throwOnError:false});$ Pa and $var element = document.getElementById("moose-equation-e8d930d4-9bec-4556-902c-65b787cf4a0e");katex.render("1 \\times 10^5", element, {displayMode:false,throwOnError:false});$ , respectively, and there is no AB initially present.

Analytical solution

Ambrosek and Longhurst (2008) provides the analytical equation for the partial pressure of AB when the conversion rate at the surface is high as

$(1)var element = document.getElementById("moose-equation-7eb2767e-7caa-44f4-ae93-0cf83e16f786");katex.render(" P_{AB} = \\frac{2 P_{A_2}^0 P_{B_2}^0}{P_{A_2}^0 + P_{B_2}^0} \\left(1 - \\exp \\left( -\\frac{S K_d k_b T}{V} t \\right)\\right).", element, {displayMode:true,throwOnError:false});$

Results

A comparison of the concentration of AB as a function of time is plotted in Figure 1. The TMAP8 calculations are found to be in good agreement with the analytical solution, with a root mean square percentage error (RMSPE) of RMSPE = 1.32%. The concentrations of A $var element = document.getElementById("moose-equation-57d8b87f-b928-4631-91fc-9f643a27ad4f");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-1a46a80b-3340-43a1-8474-dd7799a9965d");katex.render("_2", element, {displayMode:false,throwOnError:false});$ as a function of time are also plotted in Figure 1.

Figure 1: Comparison of concentration of AB as a function of time calculated through TMAP8 and analytically for the solution in ratedep condition when A $var element = document.getElementById("moose-equation-fb692268-a046-47e0-91b0-c16735cdc62d");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-6c2714d4-394b-4608-9cb3-1bf5be933b48");katex.render("_2", element, {displayMode:false,throwOnError:false});$ have unequal pressures (Ambrosek and Longhurst, 2008).

Input files

The input file for this case can be obtained by updating the initial pressure for B $var element = document.getElementById("moose-equation-2cbf1060-0a16-46de-9daf-8995349c0898");katex.render("_2", element, {displayMode:false,throwOnError:false});$ in (test/tests/ver-1ia/ver-1ia.i), which is the approach used to create tests in TMAP8 at (test/tests/ver-1ib/tests).

References

James Ambrosek and GR Longhurst. Verification and Validation of TMAP7. Technical Report INEEL/EXT-04-01657, Idaho National Engineering and Environmental Laboratory, December 2008.

@techreport{ambrosek2008verification,
    author = "Ambrosek, James and Longhurst, GR",
    title = "{Verification and Validation of TMAP7}",
    year = "2008",
    number = "INEEL/EXT-04-01657",
    month = "December",
    institution = "Idaho National Engineering and Environmental Laboratory"
}

(test/tests/ver-1ia/ver-1ia.i)

k_b = '${units 1.380649e-23 J/K}' # Boltzmann constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
T = '${units 1000 K}' # Temperature
V = '${units 1 m^3}' # Volume
S = '${units 25 cm^2 -> m^2}' # Area
p0_A2 = '${units 1e4 Pa}' # Initial pressure for A2
p0_B2 = '${units 1e4 Pa}' # Initial pressure for B2
peq_AB = '${units ${fparse 2 * ${p0_A2} * ${p0_B2} / ( ${p0_A2} + ${p0_B2} )} Pa}' # pressure in equilibration for AB
simulation_time = '${units 6 s}'
time_interval = '${units 0.01 s}'
K_r = '${units 5.88e-26 m^4/at/s}' # recombination rate for A2 or B2
K_d = '${units ${fparse 1.858e24 / sqrt( ${T} )} at/m^2/s/Pa}' # dissociation rate for AB

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[Variables]
  [p_AB]
    initial_condition = 0
  []
[]

[AuxVariables]
  [p_A2]
    initial_condition = ${p0_A2}
  []
  [p_B2]
    initial_condition = ${p0_B2}
  []
  [c_A_dot_c_B]
  []
[]

[AuxKernels]
  [p_A2_kernel]
    type = ParsedAux
    variable = p_A2
    coupled_variables = 'p_AB'
    expression = '${p0_A2} - p_AB / 2'
  []
  [p_B2_kernel]
    type = ParsedAux
    variable = p_B2
    coupled_variables = 'p_AB'
    expression = '${p0_B2} - p_AB / 2'
  []
  [c_A_dot_c_B_kernel]
    type = ParsedAux
    variable = c_A_dot_c_B
    expression = '${K_d} * ${peq_AB} / 2 / ${K_r}'
  []
[]

[Kernels]
  [timeDerivative_p_AB]
    type = ADTimeDerivative
    variable = p_AB
  []
  [MatReaction_p_AB_recombination]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'c_A_dot_c_B'
    coeff = '${fparse 2 * ${k_b} * ${T} * ${S} / ${V}}'
    reaction_rate_name = '${K_r}'
  []
  [MatReaction_p_AB_dissociation]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'p_AB'
    coeff = '${fparse -1 * ${k_b} * ${T} * ${S} / ${V}}'
    reaction_rate_name = '${K_d}'
  []
[]

[BCs]
  [p_AB_neumann] # No flux on the sides
    type = NeumannBC
    variable = p_AB
    boundary = 'left right bottom top'
    value = 0
  []
[]

[Postprocessors]
  [pressure_A2]
    type = ElementAverageValue
    variable = p_A2
    execute_on = 'initial timestep_end'
  []
  [pressure_B2]
    type = ElementAverageValue
    variable = p_B2
    execute_on = 'initial timestep_end'
  []
  [pressure_AB]
    type = ElementAverageValue
    variable = p_AB
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-10

  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  end_time = ${simulation_time}
  dt = ${time_interval}
  automatic_scaling = true
[]

[Outputs]
  csv = true
[]

(test/tests/ver-1ia/ver-1ia.i)

k_b = '${units 1.380649e-23 J/K}' # Boltzmann constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
T = '${units 1000 K}' # Temperature
V = '${units 1 m^3}' # Volume
S = '${units 25 cm^2 -> m^2}' # Area
p0_A2 = '${units 1e4 Pa}' # Initial pressure for A2
p0_B2 = '${units 1e4 Pa}' # Initial pressure for B2
peq_AB = '${units ${fparse 2 * ${p0_A2} * ${p0_B2} / ( ${p0_A2} + ${p0_B2} )} Pa}' # pressure in equilibration for AB
simulation_time = '${units 6 s}'
time_interval = '${units 0.01 s}'
K_r = '${units 5.88e-26 m^4/at/s}' # recombination rate for A2 or B2
K_d = '${units ${fparse 1.858e24 / sqrt( ${T} )} at/m^2/s/Pa}' # dissociation rate for AB

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[Variables]
  [p_AB]
    initial_condition = 0
  []
[]

[AuxVariables]
  [p_A2]
    initial_condition = ${p0_A2}
  []
  [p_B2]
    initial_condition = ${p0_B2}
  []
  [c_A_dot_c_B]
  []
[]

[AuxKernels]
  [p_A2_kernel]
    type = ParsedAux
    variable = p_A2
    coupled_variables = 'p_AB'
    expression = '${p0_A2} - p_AB / 2'
  []
  [p_B2_kernel]
    type = ParsedAux
    variable = p_B2
    coupled_variables = 'p_AB'
    expression = '${p0_B2} - p_AB / 2'
  []
  [c_A_dot_c_B_kernel]
    type = ParsedAux
    variable = c_A_dot_c_B
    expression = '${K_d} * ${peq_AB} / 2 / ${K_r}'
  []
[]

[Kernels]
  [timeDerivative_p_AB]
    type = ADTimeDerivative
    variable = p_AB
  []
  [MatReaction_p_AB_recombination]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'c_A_dot_c_B'
    coeff = '${fparse 2 * ${k_b} * ${T} * ${S} / ${V}}'
    reaction_rate_name = '${K_r}'
  []
  [MatReaction_p_AB_dissociation]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'p_AB'
    coeff = '${fparse -1 * ${k_b} * ${T} * ${S} / ${V}}'
    reaction_rate_name = '${K_d}'
  []
[]

[BCs]
  [p_AB_neumann] # No flux on the sides
    type = NeumannBC
    variable = p_AB
    boundary = 'left right bottom top'
    value = 0
  []
[]

[Postprocessors]
  [pressure_A2]
    type = ElementAverageValue
    variable = p_A2
    execute_on = 'initial timestep_end'
  []
  [pressure_B2]
    type = ElementAverageValue
    variable = p_B2
    execute_on = 'initial timestep_end'
  []
  [pressure_AB]
    type = ElementAverageValue
    variable = p_AB
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-10

  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  end_time = ${simulation_time}
  dt = ${time_interval}
  automatic_scaling = true
[]

[Outputs]
  csv = true
[]

(test/tests/ver-1ib/tests)

[Tests]
  design = 'ADMatReactionFlexible.md'
  issues = '#12'
  verification = 'ver-1ib.md'
  [ver-1ib_csv]
    type = CSVDiff
    input = '../ver-1ia/ver-1ia.i'
    cli_args = 'p0_B2="${units 1e5 Pa}" time_interval="${units 0.05 s}" Outputs/file_base="ver-1ib_out"' # these are the only things that need to be updated from ver-1ia to ver-1ib.
    csvdiff = ver-1ib_out.csv
    requirement = 'The system shall be able to model a equilibration problem on a reactive surface with unequal starting pressures in ratedep condition and to generate CSV data for use in comparisons with the analytic solution over time.'
  []
  [ver-1ib_comparison]
    type = RunCommand
    command = 'python3 comparison_ver-1ib.py'
    requirement = 'The system shall be able to generate comparison plots between the analytical solution and simulated solution of a equilibration on a reactive surface in ratedep condition with unequal starting pressures.'
    required_python_packages = 'matplotlib numpy pandas scipy os'
  []
[]

General Case Description
Analytical solution
Results
Input files
References