ver-1ie

Species Equilibration Problem in Lawdep Condition with Equal 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 similar to ver-1ia, except that, in the current case, the reaction is in lawdep condition. The case is simulated in (test/tests/ver-1ie/ver-1ie.i).

The problem considers the reaction between two isotopic species, A $var element = document.getElementById("moose-equation-6f234960-0314-4955-993d-117cc273eb18");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-29cca4a1-889b-477d-9bd1-6e9863436d0d");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , on a surface in lawdep condition. The reaction between AB, A $var element = document.getElementById("moose-equation-ff3c4cb9-2043-4ed6-996f-0b69099983f2");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , and B $var element = document.getElementById("moose-equation-1d6c977e-8796-4c49-9c77-71117782ca42");katex.render("_2", element, {displayMode:false,throwOnError:false});$ is the same as in ver-1ia. Therefore, the partial pressure of AB in equilibrium depends on the initial partial pressures of A $var element = document.getElementById("moose-equation-ddd578dc-6289-4742-836f-8c0c062314da");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-da2b4295-4cae-4e9f-a917-516800db4f5a");katex.render("_2", element, {displayMode:false,throwOnError:false});$ :

$(1)var element = document.getElementById("moose-equation-98b3476b-96b3-4a3a-883e-f954dd116ddd");katex.render(" P_{AB}^{eq} = \\frac{2 P_{A_2}^0 P_{B_2}^0}{P_{A_2}^0 + P_{B_2}^0}.", element, {displayMode:true,throwOnError:false});$

Just as in ver-1ia, we solve the net current of AB molecules from surface to the enclosure with

$(2)var element = document.getElementById("moose-equation-5c69efd6-0e55-4e03-85cd-ddbfa0f4a373");katex.render(" \\frac{d P_{AB}}{dt} = \\frac{S k_b T}{V} (2 K_r C_A C_B - K_d P_{AB}),", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-3f36d91a-453f-4262-b885-ba540270c2a6");katex.render("t", element, {displayMode:false,throwOnError:false});$ is the time, $var element = document.getElementById("moose-equation-9313219b-bb84-4d34-afaa-8ff540d40184");katex.render("S", element, {displayMode:false,throwOnError:false});$ is the surface area, $var element = document.getElementById("moose-equation-5a668415-42be-4286-aa4a-0f0f5764c3a1");katex.render("k_b", element, {displayMode:false,throwOnError:false});$ is the Boltzmann’s constant, $var element = document.getElementById("moose-equation-411ea8b6-1a54-4e5f-886d-8eb8a8b31bfe");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature, $var element = document.getElementById("moose-equation-ef31bc15-7312-4c31-910b-3fb50fff50c6");katex.render("V", element, {displayMode:false,throwOnError:false});$ is the volume in the enclosure, $var element = document.getElementById("moose-equation-416d3a50-5d83-4c7b-b1a8-b42071540aa5");katex.render("K_r", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-4b3bdde5-ed43-4c8a-b594-b1fb3de1de81");katex.render("K_d", element, {displayMode:false,throwOnError:false});$ are the recombination and dissociation coefficients, and $var element = document.getElementById("moose-equation-364dc276-b034-495c-af58-15cb1387b929");katex.render("C_A", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-7ca48d60-adc4-4ab3-a06a-6cdd27904d29");katex.render("C_B", element, {displayMode:false,throwOnError:false});$ are the concentration of atoms from A $var element = document.getElementById("moose-equation-d1e60200-ae3b-44bc-bccb-c1252d2c52db");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-261f38e3-4772-4d6c-929b-db739a76200d");katex.render("_2", element, {displayMode:false,throwOnError:false});$ on the reactive surface, respectively. In lawdep diffusion boundary condition, the concentration of A $var element = document.getElementById("moose-equation-4b0dedd0-0076-48d5-9b97-a54b156078ca");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-2d978125-6308-4393-ba71-f9699c6fe956");katex.render("_2", element, {displayMode:false,throwOnError:false});$ are always fixed relative to the partial pressures in the gas over the surface. When heteronuclear species formation is involved, TMAP8 uses logic similar to that used in the ratedep and surfdep condition for the arrival rate of gas atoms to the surface. However, there are no barriers to adsorption or release, and conversion is assumed to take place instantaneously. Any gas that does not diffuse away is immediately released from the surface. Therefore, the concentration of A $var element = document.getElementById("moose-equation-af240b22-3d67-4d6b-979d-8577e1dcf2c5");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-782fe619-8bed-48c1-95c1-d3889e2c2e77");katex.render("_2", element, {displayMode:false,throwOnError:false});$ from Sieverts' law are given by

$(3)var element = document.getElementById("moose-equation-6285c940-c245-499b-b0a7-b5d59a9f24e3");katex.render(" C_A = K_s \\sqrt{P_{A_2}},", element, {displayMode:true,throwOnError:false});$

and

$(4)var element = document.getElementById("moose-equation-ab4712a6-51fc-4279-a1b3-edf8932f8aee");katex.render(" C_B = K_s \\sqrt{P_{B_2}},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-9aa00312-f155-4684-8267-3cf47948fcd8");katex.render("K_s", element, {displayMode:false,throwOnError:false});$ is Sieverts’ solubility. Due to in the isotopic reaction, $var element = document.getElementById("moose-equation-e7e0705d-f0e8-4877-a1b2-461a817dc6b9");katex.render("K_s", element, {displayMode:false,throwOnError:false});$ is the same for each homonuclear species. The relationship between $var element = document.getElementById("moose-equation-46f9331f-d703-4f07-b61e-814388fa1d39");katex.render("K_s", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-d7af99b2-4f03-4f14-8222-682b560b0a0d");katex.render("K_r", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-08aaf536-6b47-4c7a-bbd8-68bf4ee31363");katex.render("K_d", element, {displayMode:false,throwOnError:false});$ is given by

$(5)var element = document.getElementById("moose-equation-cfcf1c7e-fc42-4062-bd1b-64a51a6fc986");katex.render(" K_d = {K_s}^2 K_r.", element, {displayMode:true,throwOnError:false});$

This case uses equal starting pressures of $var element = document.getElementById("moose-equation-c41e968a-2a1c-4abc-b1c4-1b2f776308f8");katex.render("1 \\times 10^{4}", element, {displayMode:false,throwOnError:false});$ Pa of A $var element = document.getElementById("moose-equation-6074d985-5166-43f5-a503-a1dcc9b4fb65");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-b9347c6e-7964-4cb8-ac86-a2ce39a6e034");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and no AB. $var element = document.getElementById("moose-equation-4b6886f5-b489-49cf-9b8f-4d90ce94159b");katex.render("K_d", element, {displayMode:false,throwOnError:false});$ is specified to be $var element = document.getElementById("moose-equation-6fdc40bc-d1a1-4212-828b-be1f4697ccb2");katex.render("1.858 \\times 10^{24}/\\sqrt{T}", element, {displayMode:false,throwOnError:false});$ atom/m $var element = document.getElementById("moose-equation-ac3f5e17-2aa3-4ba4-a79e-2ac333bdb0fd");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s/pa. $var element = document.getElementById("moose-equation-5e514035-1d44-4b54-878b-4721728ccd2f");katex.render("K_s", element, {displayMode:false,throwOnError:false});$ is specified to be $var element = document.getElementById("moose-equation-6c24c3cc-1e31-48f7-bb26-7e4508d59d85");katex.render("1 \\times 10^{24}", element, {displayMode:false,throwOnError:false});$ atom/m $var element = document.getElementById("moose-equation-b20affe9-f16c-4aaa-81ca-10b33db9f4a2");katex.render("^3", element, {displayMode:false,throwOnError:false});$ /Pa $var element = document.getElementById("moose-equation-0d90974c-f8db-4ac4-9b1c-63f89bb6c825");katex.render("^{0.5}", element, {displayMode:false,throwOnError:false});$ , the temperature is 1000 K, the surface area for reaction is 0.05 m $var element = document.getElementById("moose-equation-41504a41-7bec-49ef-bd64-3711e94c11cf");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 0.05 m square, and the enclosure volume is 1 m $var element = document.getElementById("moose-equation-8933d46c-52e0-4fcd-804f-38b4a4850529");katex.render("^3", element, {displayMode:false,throwOnError:false});$ .

Analytical solution

After combining Eq. (3) and Eq. (4), Eq. (2) becomes

$(6)var element = document.getElementById("moose-equation-96a65f94-0fa5-4572-8a1d-31f0cbccaba2");katex.render(" \\frac{d P_{AB}}{dt} = \\frac{S k_b T K_d}{V} \\left(2 \\sqrt{P^0_{A_2} - \\frac{P_{AB}}{2}} \\sqrt{P^0_{B_2} - \\frac{P_{AB}}{2}} - P_{AB} \\right).", element, {displayMode:true,throwOnError:false});$

This is a non-linear function, but it has a special solution when $var element = document.getElementById("moose-equation-bc8539c1-9720-44d2-aa53-3a67d469b1d1");katex.render("P^0_{A_2} = P^0_{B_2}", element, {displayMode:false,throwOnError:false});$ , which is true in the current case. Thus, the analytical solution for the partial pressure of AB is given by

$(7)var element = document.getElementById("moose-equation-67126c56-b2d1-44ba-9e39-0f6da0bf39d8");katex.render(" P_{AB} = P_{A_2}^0 \\left(1 - \\exp \\left( -\\frac{ 2 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 = 0.36%. The concentrations of A $var element = document.getElementById("moose-equation-e99db8d2-7299-4773-ad7f-cf6001bb29f6");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-897bf1d1-c43b-47c0-881a-626e31a47e68");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 lawdep condition when A $var element = document.getElementById("moose-equation-a3105692-0d3b-491f-97c8-f1b2908e9aab");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-11891ce9-f33a-4cbd-9db4-ab57af25c6a8");katex.render("_2", element, {displayMode:false,throwOnError:false});$ have equal pressures (Ambrosek and Longhurst, 2008).

Input files

The input file for this case can be found at (test/tests/ver-1ie/ver-1ie.i), which is also used as tests in TMAP8 at (test/tests/ver-1ie/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-1ie/ver-1ie.i)

k = '${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
simulation_time = '${units 3 s}'
time_interval = '${units 0.01 s}'
K_s = '${units 1.0e24 at/m^3/Pa^0.5}' # atom/m^3/pa^0.5 recombination rate for A2 or B2
K_d = '${units ${fparse 1.858e24 / sqrt( ${T} )} at/m^2/s/Pa}' # dissociation rate for AB
K_r = '${units ${fparse K_d / K_s / K_s} m^4/at/s}'

[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]
  []
  [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_kernel]
    type = ParsedAux
    variable = c_A
    coupled_variables = 'p_A2'
    expression = '${K_s} * sqrt(p_A2)'
  []
  [c_B_kernel]
    type = ParsedAux
    variable = c_B
    coupled_variables = 'p_B2'
    expression = '${K_s} * sqrt(p_B2)'
  []
[]

[Kernels]
  [timeDerivative_p_AB]
    type = ADTimeDerivative
    variable = p_AB
  []
  [MatReaction_p_AB_recombination]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'c_A c_B'
    coeff = '${fparse 2 * ${k} * ${T} * ${S} / ${V}}'
    reaction_rate_name = '${K_r}'
  []
  [MatReaction_p_AB_dissociation]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'p_AB'
    coeff = '${fparse -1 * ${k} * ${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-1ie/ver-1ie.i)

k = '${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
simulation_time = '${units 3 s}'
time_interval = '${units 0.01 s}'
K_s = '${units 1.0e24 at/m^3/Pa^0.5}' # atom/m^3/pa^0.5 recombination rate for A2 or B2
K_d = '${units ${fparse 1.858e24 / sqrt( ${T} )} at/m^2/s/Pa}' # dissociation rate for AB
K_r = '${units ${fparse K_d / K_s / K_s} m^4/at/s}'

[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]
  []
  [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_kernel]
    type = ParsedAux
    variable = c_A
    coupled_variables = 'p_A2'
    expression = '${K_s} * sqrt(p_A2)'
  []
  [c_B_kernel]
    type = ParsedAux
    variable = c_B
    coupled_variables = 'p_B2'
    expression = '${K_s} * sqrt(p_B2)'
  []
[]

[Kernels]
  [timeDerivative_p_AB]
    type = ADTimeDerivative
    variable = p_AB
  []
  [MatReaction_p_AB_recombination]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'c_A c_B'
    coeff = '${fparse 2 * ${k} * ${T} * ${S} / ${V}}'
    reaction_rate_name = '${K_r}'
  []
  [MatReaction_p_AB_dissociation]
    type = ADMatReactionFlexible
    variable = p_AB
    vs = 'p_AB'
    coeff = '${fparse -1 * ${k} * ${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-1ie/tests)

[Tests]
  design = 'ADMatReactionFlexible.md'
  issues = '#12'
  verification = 'ver-1ie.md'
  [ver-1ie_csv]
    type = CSVDiff
    input = 'ver-1ie.i'
    csvdiff = ver-1ie_out.csv
    requirement = 'The system shall be able to model a equilibration problem on a reactive surface in lawdep condition with equal starting pressures and to generate CSV data for use in comparisons with the analytic solution over time.'
  []
  [ver-1ie_comparison]
    type = RunCommand
    command = 'python3 comparison_ver-1ie.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 lawdep condition with equal starting pressures.'
    required_python_packages = 'matplotlib numpy pandas scipy os'
  []
[]

General Case Description
Analytical solution
Results
Input files
References