ver-1if

Species Equilibration Problem in Lawdep Condition with Unequal Starting Pressures

General Case Description

This verification problem is taken from Ambrosek and Longhurst (2008) and builds on ver-1ie. The configuration and modeling parameters are the same as in ver-1ie, except that, in the current case, the starting pressures for A $var element = document.getElementById("moose-equation-aba841d7-dd29-4aa9-aa97-a3645c152031");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-2287d82b-e1b7-4353-8150-d5b67712cbbd");katex.render("_2", element, {displayMode:false,throwOnError:false});$ are not equal. The case is simulated in (test/tests/ver-1ie/ver-1ie.i), but the starting pressures of A $var element = document.getElementById("moose-equation-f386da3b-f80c-45cf-93cb-8be049a1a777");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-feb8c08a-3d60-4564-942b-fdf55a484c1f");katex.render("_2", element, {displayMode:false,throwOnError:false});$ are $var element = document.getElementById("moose-equation-74c1b865-647c-4f7e-bc7e-7cfe180851b1");katex.render("1 \\times 10^4", element, {displayMode:false,throwOnError:false});$ Pa and $var element = document.getElementById("moose-equation-d5384f34-790e-4c8d-887a-9cd211ae3369");katex.render("1 \\times 10^5", element, {displayMode:false,throwOnError:false});$ Pa, respectively, and there is no AB initially present.

Analytical solution

Similar with ver-1ie, the governing equation becomes

$(1)var element = document.getElementById("moose-equation-72be28fe-00d6-468e-bc98-54a8a66aef07");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, the analytical solution is hard to calculate. Ambrosek and Longhurst (2008) uses the solution from ver-1ia with a different saturation time constant to compare with the numerical solution. The equation is defined by

$(2)var element = document.getElementById("moose-equation-0fb3477e-b405-4497-965b-0229f846a5f4");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-88dd57a8-8397-40b8-9227-89655579826e");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ is the saturation time constant. The value of $var element = document.getElementById("moose-equation-61c2b197-1e75-4a0f-8506-c85b9fd17074");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ is selected as 0.123 to match the numerical solution in Eq. (2).

note:No exact analytical solution

Because Eq. (1) is a complex non-linear function, we have no exact analytical solution to compare against. This case is therefore not truly a verification case. However, we compare the numerical solution provided by TMAP8 against an adapted analytical solution from ver-1ia and we benchmark the TMAP8 predictions against the TMAP7 results from Ambrosek and Longhurst (2008).

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 reasonable agreement with the solution from TMAP7, with a root mean square percentage error (RMSPE) of RMSPE = 13.56%. It makes sense for us because the solution from TMAP7 is not the real analytical solution. Also, the results from TMAP8 represent an improvement from those from TMAP7, which showed some numerical instability around $var element = document.getElementById("moose-equation-3b5d26db-abe4-4299-8f60-6f4c3d4586bb");katex.render("t=2.5", element, {displayMode:false,throwOnError:false});$ s. The concentrations of A $var element = document.getElementById("moose-equation-204c08b4-ef9d-4aed-a252-73894f3d1bb8");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-ed7c126b-e3bc-4ff4-8929-16a077b349cb");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-ff9376f7-414e-40fb-8665-98c453c8530c");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and B $var element = document.getElementById("moose-equation-a201c702-3656-4c52-b26e-ee1958c4a551");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-52b27f35-f2ee-499c-a3c5-f6196dfc4fb7");katex.render("_2", element, {displayMode:false,throwOnError:false});$ in (test/tests/ver-1ie/ver-1ie.i), which is the approach used to create tests in TMAP8 at (test/tests/ver-1if/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-1if/tests)

[Tests]
  design = 'ADMatReactionFlexible.md'
  issues = '#12'
  verification = 'ver-1if.md'
  [ver-1if_csv]
    type = CSVDiff
    input = '../ver-1ie/ver-1ie.i'
    cli_args = 'p0_B2="${units 1e5 Pa}" Outputs/file_base="ver-1if_out"' # these are the only things that need to be updated from ver-1ie to ver-1if.
    csvdiff = ver-1if_out.csv
    requirement = 'The system shall be able to model a equilibration problem on a reactive surface in lawdep condition with unequal starting pressures and to generate CSV data for use in comparisons with the analytic solution over time.'
  []
  [ver-1if_comparison]
    type = RunCommand
    command = 'python3 comparison_ver-1if.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 unequal starting pressures.'
    required_python_packages = 'matplotlib numpy pandas scipy os'
  []
[]

General Case Description
Analytical solution
Results
Input files
References