ver-1kc-2

Sieverts’ Law Boundaries with Chemical Reactions and No Volumetric Source

General Case Description

Two enclosures are separated by a membrane that allows diffusion according to Sieverts' law and chemical reactions, with no volumetric source present. Enclosure 2 has twice the volume of Enclosure 1.

Case Set Up

This verification problem is taken from Ambrosek and Longhurst (2008).

Unlike the ver-1kc-1 case, which only considers tritium T $var element = document.getElementById("moose-equation-1f0633f3-2d18-48a2-b4cf-c220a3239944");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , this setup describes a diffusion system in which tritium T $var element = document.getElementById("moose-equation-d213a830-9d98-4340-aa1c-eab94f396340");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , dihydrogen H $var element = document.getElementById("moose-equation-90e96695-6cca-4344-ab83-4552f61e3a1f");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and HT are modeled across a one-dimensional domain split into two enclosures. The total system length is $var element = document.getElementById("moose-equation-a5fb7480-745e-4e31-b78a-962349ec96b7");katex.render("2.5 \\times 10^{-4}", element, {displayMode:false,throwOnError:false});$ m, divided into 100 segments. The system operates at a constant temperature of 500 Kelvin. Initial tritium T $var element = document.getElementById("moose-equation-bca7756d-20c1-43c6-b7b0-4df5cb90c0c8");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and dihydrogen H $var element = document.getElementById("moose-equation-c4c8220f-2f79-4e9e-b487-6e1895348367");katex.render("_2", element, {displayMode:false,throwOnError:false});$ pressures are specified as $var element = document.getElementById("moose-equation-a938b57a-9f34-4bc5-a757-3afc86fbc67b");katex.render("10^{5}", element, {displayMode:false,throwOnError:false});$ Pa for Enclosure 1 and $var element = document.getElementById("moose-equation-5da1a3a8-e783-42e8-93f3-463031883203");katex.render("10^{-10}", element, {displayMode:false,throwOnError:false});$ Pa for Enclosure 2. Initially, there is no HT in either enclosure.

The reaction between the species is described as follows

$var element = document.getElementById("moose-equation-af446c01-9fab-4d5f-a7d5-4a48cab0c97f");katex.render("\\text{H}_2 + \\text{T}_2 \\leftrightarrow 2\\text{HT}", element, {displayMode:true,throwOnError:false});$

The kinematic evolutions of the species are given by the following equations

$var element = document.getElementById("moose-equation-699dc6ac-54ec-405c-9301-9d71a1b7bfb3");katex.render("\\frac{d C_{\\text{HT}}}{dt} = 2K_1 C_{\\text{H}_2} C_{\\text{T}_2} - K_2 C_{\\text{HT}}^2", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-e6542929-e826-4e6e-b38b-81eb5b00d598");katex.render("\\frac{d C_{\\text{H}_2}}{dt} = -K_1 C_{\\text{H}_2} C_{\\text{T}_2} + \\frac{1}{2} K_2 C_{\\text{HT}}^2", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-768e1bfd-01c4-4c7c-bbe1-eeeac29cc9d3");katex.render("\\frac{d C_{\\text{T}_2}}{dt} = -K_1 c_{\\text{H}_2} C_{\\text{T}_2} + \\frac{1}{2} K_2 C_{\\text{HT}}^2", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-07087b0c-1579-4533-b34c-6f591b57d68a");katex.render("K_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-f97d2b7a-be80-4a78-a9d6-6d17e8014a5f");katex.render("K_2", element, {displayMode:false,throwOnError:false});$ represent the reaction rates for the forward and reverse reactions, respectively.

At equilibrium, the time derivatives are zero

$var element = document.getElementById("moose-equation-d70ca0cb-a7ac-44c9-8b46-057b7688e569");katex.render("2K_1 C_{\\text{H}_2} C_{\\text{T}_2} - K_2 C_{\\text{HT}}^2 = 0", element, {displayMode:true,throwOnError:false});$

From this, we can derive the same equilibrium condition as used in TMAP7:

$var element = document.getElementById("moose-equation-eb2c445e-4f1a-4654-9e94-dd2c73a115dd");katex.render("P_{\\text{HT}} = \\eta \\sqrt{P_{\\text{H}_2} P_{\\text{T}_2}}", element, {displayMode:true,throwOnError:false});$

where the equilibrium constant $var element = document.getElementById("moose-equation-1c233c66-06d5-41be-bb85-3ecb709624d3");katex.render("\\eta", element, {displayMode:false,throwOnError:false});$ is defined as

$(1)var element = document.getElementById("moose-equation-92601546-b89e-49dd-bce7-5640bd04847e");katex.render(" \\eta = \\sqrt{\\frac{2K_1}{K_2}}", element, {displayMode:true,throwOnError:false});$

Similarly to TMAP7, the equilibrium constant $var element = document.getElementById("moose-equation-770dc47a-c0a4-44d5-8853-3c769470bb96");katex.render("\\eta", element, {displayMode:false,throwOnError:false});$ has been set to a fixed value of $var element = document.getElementById("moose-equation-b559fc6a-f2c3-4216-849a-ad9c9958b208");katex.render("\\eta = 2", element, {displayMode:false,throwOnError:false});$ .

The diffusion process for each species in the two enclosures can be expressed by

$var element = document.getElementById("moose-equation-fcbc838b-8821-4025-8c7e-94e46273dfdc");katex.render("\\frac{\\partial C_1}{\\partial t} = \\nabla D \\nabla C_1,", element, {displayMode:true,throwOnError:false});$ and $var element = document.getElementById("moose-equation-0cef4e96-940e-4636-a3c9-4b29149e6d27");katex.render("\\frac{\\partial C_2}{\\partial t} = \\nabla D \\nabla C_2,", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-415a4c6a-5a55-48f7-b48d-db0763f76c61");katex.render("C_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-6935e880-7dbf-404f-bc74-202d901c07c4");katex.render("C_2", element, {displayMode:false,throwOnError:false});$ represent the concentration fields in enclosures 1 and 2 respectively, $var element = document.getElementById("moose-equation-bab73020-895d-4f30-a440-fb72f6147c5f");katex.render("t", element, {displayMode:false,throwOnError:false});$ is the time, and $var element = document.getElementById("moose-equation-959afbf6-e343-4a37-92a0-2e148690908b");katex.render("D", element, {displayMode:false,throwOnError:false});$ denotes the diffusivity. Note that the diffusivity may vary across different species and enclosures. However, in this case, it is assumed to be identical for all.

The concentration in Enclosure 1 is related to the partial pressure and concentration in Enclosure 2 via the interface sorption law:

$var element = document.getElementById("moose-equation-f4957050-c0b7-485e-85bf-af266989b85b");katex.render("C_1 = K P_2^n = K \\left( C_2 RT \\right)^n", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-a82bb40a-6ce2-442c-a984-71f6ec4ea59e");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the ideal gas constant in J/mol/K, $var element = document.getElementById("moose-equation-cbfe25ad-67f1-44d0-ad58-80fb19187595");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K, $var element = document.getElementById("moose-equation-be77c47a-1067-46ee-9829-55a6b156d93c");katex.render("K", element, {displayMode:false,throwOnError:false});$ is the solubility, and $var element = document.getElementById("moose-equation-d893b699-44ff-4cc0-874b-80868e80f53b");katex.render("n", element, {displayMode:false,throwOnError:false});$ is the exponent of the sorption law. For Sieverts' law, $var element = document.getElementById("moose-equation-28119f77-6aa4-48c8-972c-9258e1b64de0");katex.render("n=0.5", element, {displayMode:false,throwOnError:false});$ .

Results

We assume that $var element = document.getElementById("moose-equation-afe7d218-aa8a-4dce-afaa-a4ebe16bacbb");katex.render("K = 10/\\sqrt{RT}", element, {displayMode:false,throwOnError:false});$ , which is expected to result in $var element = document.getElementById("moose-equation-31a7e063-4a46-4fdf-a108-bd1be8ea7e66");katex.render("C_1 = 10 \\sqrt{C_2}", element, {displayMode:false,throwOnError:false});$ at equilibrium. As illustrated in Figure 1, similarly to ver-1kc-1, T $var element = document.getElementById("moose-equation-58419337-11f8-45cb-8bad-da5b529caac7");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and H $var element = document.getElementById("moose-equation-4810cea5-cfc4-430c-b090-75df14dcc6b6");katex.render("_2", element, {displayMode:false,throwOnError:false});$ pressures reach equilibrium in both enclosures. What is new, however, is that HT is produced in both enclosures following the sorption law.

Thus, it is crucial to ensure that the chemical equilibrium between HT, T $var element = document.getElementById("moose-equation-de44ab4f-9b8d-46ec-bc53-637c23b9cf0e");katex.render("_2", element, {displayMode:false,throwOnError:false});$ and H $var element = document.getElementById("moose-equation-b000ffcc-7e85-4b21-9784-a567116981c2");katex.render("_2", element, {displayMode:false,throwOnError:false});$ is achieved. This can be verified in both enclosures by examining the ratio between $var element = document.getElementById("moose-equation-a3b88da0-dbb3-4d14-9b4c-69a668ea8096");katex.render("P_{\\text{HT}}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-270caf53-1fbc-4fe6-ab0f-8e05f47efa0f");katex.render("\\sqrt{P_{\\text{H}_2} P_{\\text{T}_2}}", element, {displayMode:false,throwOnError:false});$ , which must equal $var element = document.getElementById("moose-equation-f04404c5-5b22-4956-9994-74895b54131a");katex.render("\\eta=2", element, {displayMode:false,throwOnError:false});$ . As shown in Figure 2, this ratio approaches $var element = document.getElementById("moose-equation-31bbfca8-ed18-4d1d-8d0f-2d10dfbd62c0");katex.render("\\eta=2", element, {displayMode:false,throwOnError:false});$ for both enclosures, as observed in TMAP7. However, achieving this balance involves a compromise to reproduce TMAP7's results. The ratio of $var element = document.getElementById("moose-equation-d111bed7-9354-437d-b17c-fe8859eec37c");katex.render("K_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-cd518c98-46b6-4b89-8a16-e8b62c958128");katex.render("K_2", element, {displayMode:false,throwOnError:false});$ must respect Eq. (1). The values of $var element = document.getElementById("moose-equation-bbce412e-6acf-43d7-8c91-f6e13e505764");katex.render("K_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-d1e18b83-8d9a-46ad-968e-c5a9bffcc9af");katex.render("K_2", element, {displayMode:false,throwOnError:false});$ must also be large enough to ensure that the kinetics of chemical reactions are faster than diffusion or surface permeation to be closer to the equilibrium assumption used in TMAP7. Here, the equilibrium in enclosure 1 is achieved rapidly. Increasing $var element = document.getElementById("moose-equation-2d0f0e12-5a5c-48de-a703-c33e3a769eae");katex.render("K_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-1813785e-2805-404e-9ba3-ccfffbaf0a28");katex.render("K_2", element, {displayMode:false,throwOnError:false});$ would also enable a quicker attainment of equilibrium in enclosure 2. However, using very high values for $var element = document.getElementById("moose-equation-d8be0f7a-31e1-480f-9fc0-719784b40cdf");katex.render("K_1", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-06944009-b944-4fc7-98a3-da6ddaee61a5");katex.render("K_2", element, {displayMode:false,throwOnError:false});$ would lead to an unnecessary increase in computational costs.

The concentration ratios for T $var element = document.getElementById("moose-equation-9aa08e4e-51cd-4601-ab46-e9ee2fe1b1e5");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , H $var element = document.getElementById("moose-equation-1bd822bf-33b8-4925-90e9-4900fc55f204");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , and HT between enclosures 1 and 2, shown in Figure 3, Figure 4, and Figure 5, demonstrate that the results obtained with TMAP8 are consistent with the analytical results derived from the sorption law for $var element = document.getElementById("moose-equation-72ecf11c-d297-42e6-aae7-5290743f251e");katex.render("K \\sqrt{RT} = 10", element, {displayMode:false,throwOnError:false});$ .

As shown in Figure 6, mass is conserved between the two enclosures over time. The variation in mass is only $var element = document.getElementById("moose-equation-0e5c2135-98c3-4f59-851e-ee3e91373cdd");katex.render("0.6", element, {displayMode:false,throwOnError:false});$ %. This variation in mass can be further minimized by refining the mesh, i.e., increasing the number of segments in the domain.

Figure 1: Evolution of species concentration over time governed by Sieverts' law with $var element = document.getElementById("moose-equation-41a32c90-e2d9-431a-9b8a-1f3f263a695b");katex.render("K = 10/\\sqrt{RT}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-15648fff-11a3-4c95-ac08-57d6d4da9c4d");katex.render("\\eta = \\sqrt{2K_1/K_2}", element, {displayMode:false,throwOnError:false});$ .

Figure 2: Equilibrium constant as a function of time for $var element = document.getElementById("moose-equation-29d65d33-f417-4429-ba2c-741f9959d551");katex.render("\\eta = \\sqrt{2K_1/K_2}=2", element, {displayMode:false,throwOnError:false});$ .

Figure 3: T $var element = document.getElementById("moose-equation-fd4dd127-a085-4b8c-8e9e-387bf373d223");katex.render("_2", element, {displayMode:false,throwOnError:false});$ concentration ratio between enclosures 1 and 2 at the interface for $var element = document.getElementById("moose-equation-f0ea9961-bd07-4c76-afd9-3500390bb30d");katex.render("K = 10/\\sqrt{RT}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-d8f63366-90c0-446f-9868-f7c3d4cc07b2");katex.render("\\eta = \\sqrt{2K_1/K_2}", element, {displayMode:false,throwOnError:false});$ . This verifies TMAP8's ability to apply Sieverts' law across the interface.

Figure 4: H $var element = document.getElementById("moose-equation-d1fc9ee0-deda-4a01-9249-db771c9ac950");katex.render("_2", element, {displayMode:false,throwOnError:false});$ concentration ratio between enclosures 1 and 2 at the interface for $var element = document.getElementById("moose-equation-a814dbec-2d01-4b6d-929c-af35222290b9");katex.render("K = 10/\\sqrt{RT}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-0a8cd3b9-2f88-4982-aa11-4f2bdfdfa0f7");katex.render("\\eta = \\sqrt{2K_1/K_2}", element, {displayMode:false,throwOnError:false});$ . This verifies TMAP8's ability to apply Sieverts' law across the interface.

Figure 5: HT concentration ratio between enclosures 1 and 2 at the interface for $var element = document.getElementById("moose-equation-1d83f68d-3bb2-4751-bc20-48109fb9ccea");katex.render("K = 10/\\sqrt{RT}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-09798991-2361-4284-b89f-621f1013c837");katex.render("\\eta = \\sqrt{2K_1/K_2}", element, {displayMode:false,throwOnError:false});$ . This verifies TMAP8's ability to apply Sieverts' law across the interface.

Figure 6: Total mass conservation across both enclosures over time for $var element = document.getElementById("moose-equation-1c6d5613-fc95-4870-882f-80d906231c94");katex.render("K = 10/\\sqrt{RT}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-c1365a39-301a-489e-aa93-3d3473017143");katex.render("\\eta = \\sqrt{2K_1/K_2}", element, {displayMode:false,throwOnError:false});$ .

Input files

The input file for this case can be found at (test/tests/ver-1kc-2/ver-1kc-2.i).

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-1kc-2/ver-1kc-2.i)

nb_segments_TMAP7 = 20
node_size_TMAP7 = '${units 1.25e-5 m}'
long_total = '${units ${fparse nb_segments_TMAP7 * node_size_TMAP7} m}'
nb_segments_TMAP8 = 1e2
simulation_time = '${units 0.25 s}'
temperature = '${units 500 K}'
R = '${units 8.31446261815324 J/mol/K}' # ideal gas constant from PhysicalConstants.h
initial_pressure_1 = '${units 1e5 Pa}'
initial_pressure_2 = '${units 1e-10 Pa}'
initial_concentration_1 = '${units ${fparse initial_pressure_1 / (R*temperature)} mol/m^3}'
initial_concentration_2 = '${units ${fparse initial_pressure_2 / (R*temperature)} mol/m^3}'
solubility = '${units ${fparse 10/sqrt(R*temperature)} mol/m^3/Pa^(1/2)}' # Sieverts' law solubility
diffusivity = '${units ${fparse 4.31e-6 * exp(-2818/temperature)} m^2/s}'
n_sorption = 0.5 # Sieverts' Law
K1 = '${units 4000 mol/m^3/s}' # reaction rate for H2+T2->2HT
equilibrium_constant = 2.0
K2 = '${units ${fparse (2*K1) / (equilibrium_constant)^2} mol/m^3/s}' # reaction rate for 2HT->H2+T2
unit_scale = 1
unit_scale_neighbor = 1

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = ${nb_segments_TMAP8}
    xmax = ${long_total}
  []
  [enclosure_1]
    type = SubdomainBoundingBoxGenerator
    input = generated
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '${fparse 1/3 * long_total} 0 0'
  []
  [enclosure_2]
    type = SubdomainBoundingBoxGenerator
    input = enclosure_1
    block_id = 2
    bottom_left = '${fparse 1/3 * long_total} 0 0'
    top_right = '${fparse long_total} 0 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = enclosure_2
    primary_block = 1
    paired_block = 2
    new_boundary = interface
  []
  [interface2]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = 2
    paired_block = 1
    new_boundary = interface2
  []
[]

[Variables]

  # Variables for H2, T2, and HT in enclosure 1
  [concentration_H2_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_T2_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_1}'
  []
  [concentration_HT_enclosure_1]
    block = 1
    initial_condition = '${initial_concentration_2}'
  []

  # Variables for H2, T2, and HT in enclosure 2
  [concentration_H2_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
  [concentration_T2_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
  [concentration_HT_enclosure_2]
    block = 2
    initial_condition = '${initial_concentration_2}'
  []
[]

[AuxVariables]

  # Auxiliary variables for H2, T2, and HT in enclosure 1
  [pressure_H2_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_T2_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_1}'
  []
  [pressure_HT_enclosure_1]
    order = CONSTANT
    family = MONOMIAL
    block = 1
    initial_condition = '${initial_pressure_2}'
  []

  # Auxiliary variables for H2, T2, and HT in enclosure 2
  [pressure_H2_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
  [pressure_T2_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
  [pressure_HT_enclosure_2]
    order = CONSTANT
    family = MONOMIAL
    block = 2
    initial_condition = '${initial_pressure_2}'
  []
[]

[Kernels]
  # Diffusion equation for H2
  [H2_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_H2_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [H2_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_H2_enclosure_1
    block = '1'
  []
  [H2_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_H2_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [H2_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_H2_enclosure_2
    block = '2'
  []

  # Diffusion equation for T2
  [T2_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_T2_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [T2_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_T2_enclosure_1
    block = '1'
  []
  [T2_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_T2_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [T2_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_T2_enclosure_2
    block = '2'
  []

  # Diffusion equation for HT
  [HT_diffusion_enclosure_1]
    type = MatDiffusion
    variable = concentration_HT_enclosure_1
    diffusivity = ${diffusivity}
    block = '1'
  []
  [HT_time_derivative_enclosure_1]
    type = TimeDerivative
    variable = concentration_HT_enclosure_1
    block = '1'
  []
  [HT_diffusion_enclosure_2]
    type = MatDiffusion
    variable = concentration_HT_enclosure_2
    diffusivity = ${diffusivity}
    block = '2'
  []
  [HT_time_derivative_enclosure_2]
    type = TimeDerivative
    variable = concentration_HT_enclosure_2
    block = '2'
  []

  # Reaction H2+T2->2HT in enclosure 1
  [reaction_H2_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_T2_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_HT_encl_1_1]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_1
    vs = 'concentration_H2_enclosure_1 concentration_T2_enclosure_1'
    block = 1
    coeff = 2
    reaction_rate_name = ${K1}
  []

  # Reaction 2HT->H2+T2 in enclosure 1
  [reaction_H2_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_T2_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_HT_encl_1_2]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_1
    vs = 'concentration_HT_enclosure_1 concentration_HT_enclosure_1'
    block = 1
    coeff = -1
    reaction_rate_name = ${K2}
  []

  # Reaction H2+T2->2HT in enclosure 2
  [reaction_H2_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_T2_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K1}
  []
  [reaction_HT_encl_2_1]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_2
    vs = 'concentration_H2_enclosure_2 concentration_T2_enclosure_2'
    block = 2
    coeff = 2
    reaction_rate_name = ${K1}
  []

  # Reaction 2HT->H2+T2 in enclosure 2
  [reaction_H2_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_H2_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_T2_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_T2_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = 0.5
    reaction_rate_name = ${K2}
  []
  [reaction_HT_encl_2_2]
    type = ADMatReactionFlexible
    variable = concentration_HT_enclosure_2
    vs = 'concentration_HT_enclosure_2 concentration_HT_enclosure_2'
    block = 2
    coeff = -1
    reaction_rate_name = ${K2}
  []
[]

[AuxKernels]

  # Auxiliary kernels for H2, T2, and HT in enclosure 1
  [pressure_H2_enclosure_1]
    type = ParsedAux
    variable = pressure_H2_enclosure_1
    coupled_variables = 'concentration_H2_enclosure_1'
    expression = '${fparse R*temperature}*concentration_H2_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_1]
    type = ParsedAux
    variable = pressure_T2_enclosure_1
    coupled_variables = 'concentration_T2_enclosure_1'
    expression = '${fparse R*temperature}*concentration_T2_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_1]
    type = ParsedAux
    variable = pressure_HT_enclosure_1
    coupled_variables = 'concentration_HT_enclosure_1'
    expression = '${fparse R*temperature}*concentration_HT_enclosure_1'
    block = 1
    execute_on = 'initial timestep_end'
  []

  # Auxiliary kernels for H2, T2, and HT in enclosure 2
  [pressure_H2_enclosure_2]
    type = ParsedAux
    variable = pressure_H2_enclosure_2
    coupled_variables = 'concentration_H2_enclosure_2'
    expression = '${fparse R*temperature}*concentration_H2_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_2]
    type = ParsedAux
    variable = pressure_T2_enclosure_2
    coupled_variables = 'concentration_T2_enclosure_2'
    expression = '${fparse R*temperature}*concentration_T2_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_2]
    type = ParsedAux
    variable = pressure_HT_enclosure_2
    coupled_variables = 'concentration_HT_enclosure_2'
    expression = '${fparse R*temperature}*concentration_HT_enclosure_2'
    block = 2
    execute_on = 'initial timestep_end'
  []
[]

[InterfaceKernels]
  [interface_sorption_H2]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_H2_enclosure_1
    neighbor_var = concentration_H2_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
  [interface_sorption_T2]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_T2_enclosure_1
    neighbor_var = concentration_T2_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
  [interface_sorption_HT]
    type = InterfaceSorption
    K0 = ${solubility}
    Ea = 0
    n_sorption = ${n_sorption}
    diffusivity = ${diffusivity}
    unit_scale = ${unit_scale}
    unit_scale_neighbor = ${unit_scale_neighbor}
    temperature = ${temperature}
    variable = concentration_HT_enclosure_1
    neighbor_var = concentration_HT_enclosure_2
    sorption_penalty = 1e1
    boundary = interface
  []
[]

[Postprocessors]

  # postprocessors for H2
  [pressure_H2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_2]
    type = ElementAverageValue
    variable = pressure_H2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_H2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_H2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_H2]
    type = ParsedPostprocessor
    expression = 'concentration_H2_enclosure_1_at_interface / sqrt(concentration_H2_enclosure_2_at_interface)'
    pp_names = 'concentration_H2_enclosure_1_at_interface concentration_H2_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [pressure_H2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = pressure_H2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [pressure_H2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = pressure_H2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_H2_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_H2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for T2
  [pressure_T2_enclosure_1]
    type = ElementAverageValue
    variable = pressure_T2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_T2_enclosure_2]
    type = ElementAverageValue
    variable = pressure_T2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_T2_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_T2_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_T2]
    type = ParsedPostprocessor
    expression = 'concentration_T2_enclosure_1_at_interface / sqrt(concentration_T2_enclosure_2_at_interface)'
    pp_names = 'concentration_T2_enclosure_1_at_interface concentration_T2_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [concentration_T2_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_T2_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_T2_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_T2_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for HT
  [pressure_HT_enclosure_1]
    type = ElementAverageValue
    variable = pressure_HT_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_HT_enclosure_2]
    type = ElementAverageValue
    variable = pressure_HT_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_enclosure_1_at_interface]
    type = SideAverageValue
    boundary = interface
    variable = concentration_HT_enclosure_1
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_enclosure_2_at_interface]
    type = SideAverageValue
    boundary = interface2
    variable = concentration_HT_enclosure_2
    outputs = 'csv console'
    execute_on = 'initial timestep_end'
  []
  [concentration_ratio_HT]
    type = ParsedPostprocessor
    expression = 'concentration_HT_enclosure_1_at_interface / sqrt(concentration_HT_enclosure_2_at_interface)'
    pp_names = 'concentration_HT_enclosure_1_at_interface concentration_HT_enclosure_2_at_interface'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [concentration_HT_encl_1_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_HT_enclosure_1
    block = 1
    execute_on = 'initial timestep_end'
  []
  [concentration_HT_encl_2_inventory]
    type = ElementIntegralVariablePostprocessor
    variable = concentration_HT_enclosure_2
    block = 2
    execute_on = 'initial timestep_end'
  []

  # postprocessors for mass conservation
  [mass_conservation_sum_encl1_encl2]
    type = LinearCombinationPostprocessor
    pp_names = 'concentration_HT_encl_1_inventory concentration_HT_encl_2_inventory concentration_H2_encl_1_inventory concentration_H2_encl_2_inventory concentration_T2_encl_1_inventory concentration_T2_encl_2_inventory'
    pp_coefs = '1            1           1            1            1            1'
    execute_on = 'initial timestep_end'
  []

  # postprocessors for equilibrium constant
  [equilibrium_constant_encl_1]
    type = ParsedPostprocessor
    expression = 'pressure_HT_enclosure_1 / sqrt(pressure_H2_enclosure_1 * pressure_T2_enclosure_1)'
    pp_names = 'pressure_HT_enclosure_1 pressure_H2_enclosure_1 pressure_T2_enclosure_1'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [equilibrium_constant_encl_2]
    type = ParsedPostprocessor
    expression = 'pressure_HT_enclosure_2 / sqrt(pressure_H2_enclosure_2 * pressure_T2_enclosure_2)'
    pp_names = 'pressure_HT_enclosure_2 pressure_H2_enclosure_2 pressure_T2_enclosure_2'
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  end_time = ${simulation_time}
  dtmax = 1e-2
  nl_max_its = 9
  l_max_its = 30
  scheme = 'bdf2'
  solve_type = PJFNK
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-6
    optimal_iterations = 7
    iteration_window = 1
    growth_factor = 1.1
    cutback_factor = 0.9
    cutback_factor_at_failure = 0.9
  []
[]

[Outputs]
  file_base = 'ver-1kc-2_out_k10'
  csv = true
  exodus = true
  execute_on = 'initial timestep_end'
[]

General Case Description
Case Set Up
Results
Input files
References