ver-1id

Species Equilibration Model in Surfdep Conditions with High Barrier Energy

General Case Description

This verification problem is taken from Ambrosek and Longhurst (2008) and builds on ver-1ic. The configuration and modeling parameters are the same as in ver-1ic, except that, in the current case, the reaction has a high barrier energy. The case is simulated in (test/tests/ver-1ic/ver-1ic.i), but $var element = document.getElementById("moose-equation-59e7724c-7fa6-4146-8181-ab4bea36dc4a");katex.render("E_x", element, {displayMode:false,throwOnError:false});$ is set to 0.20 eV instead of 0.05 eV.

Analytical solution

Ambrosek and Longhurst (2008) provides the analytical equation for the partial pressure of AB as

$(1)var element = document.getElementById("moose-equation-c1558759-f2fa-4407-84df-76ce29b193cd");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{t}{\\tau}\\right)\\right),", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-77ba57fd-9145-4535-bc6a-2aed581e6aaa");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ is defined as

$(2)var element = document.getElementById("moose-equation-0ca2598f-c245-488f-bc67-64e4c7054363");katex.render(" \\tau = \\frac{V (\\hat{K_r} + \\hat{K_b})}{S k_b T \\hat{K_d} \\hat{K_b}}.", 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 = 0.19%. The concentrations of A $var element = document.getElementById("moose-equation-9f50229f-c54c-4731-a571-ea31e2abe290");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-ec4e9f38-1c28-4007-a62c-a178203d24fa");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 surfdep condition with high barrier energy (Ambrosek and Longhurst, 2008).

Input files

The input file for this case can be obtained by updating the adsorption barrier energy $var element = document.getElementById("moose-equation-d21cb63a-fc00-41ca-bfc0-8d1c98c191d1");katex.render("E_x", element, {displayMode:false,throwOnError:false});$ in (test/tests/ver-1ic/ver-1ic.i), which is the approach used to create tests in TMAP8 at (test/tests/ver-1id/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-1ic/ver-1ic.i)

k = '${units 1.380649e-23 J/K}' # Boltzmann constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
amu = '${units 1.6605390666e-27 kg}' # Atomic mass unit
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
end_time = '${units 10 s}'
time_interval = '${units 0.1 s}'
E_x = '${units 0.05 eV -> J}'
E_c = '${units -0.01 eV -> J}'
E_b = '${units 0.00 eV -> J}'
nu = '${units 8.4e12 m/s}' # Debye frequency
M = '${units ${fparse 2 * amu} kg}'
K_d = '${units ${fparse 1 / sqrt(2 * pi * ${M} * ${k} * ${T}) * exp( - ${E_x} / ( ${k} * ${T} ) )} s/kg/m}' # deposition rate
K_r = '${units ${fparse ${nu} * exp(( ${E_c} - ${E_x} ) / ${k} / ${T})} m/s}' # release rate
K_b = '${units ${fparse ${nu} * exp( - E_b / ${k} / ${T} ) } m/s}' # dissociation rate
D_s_lamda = '${units ${fparse 5.3167e-7 * exp( -4529 / ${T} )} m^4/at/s}'


[Mesh]
  type = GeneratedMesh
  dim = 2
[]

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

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

[AuxKernels]
  [c_A_dot_c_B_kernel]
    type = ParsedAux
    variable = c_A_dot_c_B
    # coupled_variables = 'p_AB c_A c_B'
    expression = ' ${peq_AB} * ${K_d} * ${K_b} / ${K_r} / 2 / ${D_s_lamda} '
  []
  [c_AB_kernel]
    type = ParsedAux
    variable = c_AB
    coupled_variables = 'p_AB c_A_dot_c_B'
    expression = ' (p_AB * ${K_d} + c_A_dot_c_B * 2 * ${D_s_lamda}) / ( ${K_r} + ${K_b} ) '
  []
  [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'
  []
[]

[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 * ${S} * ${k} * ${T} / ${V} }'
    reaction_rate_name = '${D_s_lamda}'
  []
  [MatReaction_p_AB_dissociation]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'c_AB'
    coeff = '${fparse - ${S} * ${k} * ${T} / ${V} }'
    reaction_rate_name = '${K_b}'
  []
[]

[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-15

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

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

[Outputs]
  csv = true
[]

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

k = '${units 1.380649e-23 J/K}' # Boltzmann constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
amu = '${units 1.6605390666e-27 kg}' # Atomic mass unit
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
end_time = '${units 10 s}'
time_interval = '${units 0.1 s}'
E_x = '${units 0.05 eV -> J}'
E_c = '${units -0.01 eV -> J}'
E_b = '${units 0.00 eV -> J}'
nu = '${units 8.4e12 m/s}' # Debye frequency
M = '${units ${fparse 2 * amu} kg}'
K_d = '${units ${fparse 1 / sqrt(2 * pi * ${M} * ${k} * ${T}) * exp( - ${E_x} / ( ${k} * ${T} ) )} s/kg/m}' # deposition rate
K_r = '${units ${fparse ${nu} * exp(( ${E_c} - ${E_x} ) / ${k} / ${T})} m/s}' # release rate
K_b = '${units ${fparse ${nu} * exp( - E_b / ${k} / ${T} ) } m/s}' # dissociation rate
D_s_lamda = '${units ${fparse 5.3167e-7 * exp( -4529 / ${T} )} m^4/at/s}'


[Mesh]
  type = GeneratedMesh
  dim = 2
[]

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

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

[AuxKernels]
  [c_A_dot_c_B_kernel]
    type = ParsedAux
    variable = c_A_dot_c_B
    # coupled_variables = 'p_AB c_A c_B'
    expression = ' ${peq_AB} * ${K_d} * ${K_b} / ${K_r} / 2 / ${D_s_lamda} '
  []
  [c_AB_kernel]
    type = ParsedAux
    variable = c_AB
    coupled_variables = 'p_AB c_A_dot_c_B'
    expression = ' (p_AB * ${K_d} + c_A_dot_c_B * 2 * ${D_s_lamda}) / ( ${K_r} + ${K_b} ) '
  []
  [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'
  []
[]

[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 * ${S} * ${k} * ${T} / ${V} }'
    reaction_rate_name = '${D_s_lamda}'
  []
  [MatReaction_p_AB_dissociation]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'c_AB'
    coeff = '${fparse - ${S} * ${k} * ${T} / ${V} }'
    reaction_rate_name = '${K_b}'
  []
[]

[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-15

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

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

[Outputs]
  csv = true
[]

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

[Tests]
  design = 'ADMatReactionFlexible.md'
  issues = '#12'
  verification = 'ver-1id.md'
  [ver-1id_csv]
    type = CSVDiff
    input = '../ver-1ic/ver-1ic.i'
    cli_args = "E_x='${units 0.2 eV -> J}' Outputs/file_base='ver-1id_out'" # these are the only things that need to be updated from ver-1ic to ver-1id.
    csvdiff = ver-1id_out.csv
    requirement = 'The system shall be able to model a equilibration problem on a reactive surface in surfdep conditions with high barrier energy and to generate CSV data for use in comparisons with the analytic solution over time.'
  []
  [ver-1id_comparison]
    type = RunCommand
    command = 'python3 comparison_ver-1id.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 surfdep condition with high barrier energy.'
    required_python_packages = 'matplotlib numpy pandas scipy os'
  []
[]

General Case Description
Analytical solution
Results
Input files
References