- variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on
ADMatReactionFlexible
Kernel to add -coeff*K*vs, where coeff=coefficient, K=reaction rate, vs=variables product
Overview
This class adds to the residual a contribution from
where are species concentrations provided as nonlinear or coupled variables, is a constant coefficient, and is a material property. can be used to include the stoichiometry of a reaction for specific species. For example, for a reaction such as would be the same for both species, but would be defined as -1 for species X, and as 2 for species Y.
Input Parameters
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- coeff1A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:A coefficient for multiplying the reaction term. It can be used to include the stoichiometry of a reaction for specific species.
- displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
- matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
- reaction_rate_nameKThe reaction rate used with the kernel
Default:K
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:The reaction rate used with the kernel
- vsSet this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Set this to make vs a list of coupled variables, otherwise it will use the kernel's nonlinear variable for v
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill
Contribution To Tagged Field Data Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (test/tests/ver-1kc-2/ver-1kc-2.i)
- (test/tests/pore_scale_transport/2D_absorption_base.i)
- (test/tests/ver-1gc/ver-1gc.i)
- (test/tests/val-2c/val-2c_base.i)
- (test/tests/kernels/ad_mat_reaction_flexible/ad_mat_reaction_flexible_test.i)
- (test/tests/ver-1ia/ver-1ia.i)
- (test/tests/ver-1ic/ver-1ic.i)
- (test/tests/ver-1kd/ver-1kd.i)
- (test/tests/ver-1g/ver-1g.i)
- (test/tests/ver-1ie/ver-1ie.i)
- (test/tests/val-2e/val-2ee.i)
(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'
[]
(test/tests/pore_scale_transport/2D_absorption_base.i)
# TMAP8 input file
# Written by Pierre-Clément Simon - Idaho National Laboratory
#
# Published with:
# P.-C. A. Simon, P. W. Humrickhouse, A. D. Lindsay,
# "Tritium Transport Modeling at the Pore Scale in Ceramic Breeder Materials Using TMAP8,"
# in IEEE Transactions on Plasma Science, 2022, doi: 10.1109/TPS.2022.3183525.
#
# Phase Field Model: Isotropic bulk diffusion, pore surface reactions, isotropic pore diffusion
# Type: Transient
# Phase structure: Contains two phases: the bulk of a material, and the pores. An interface separates the two
# Physics: tritium trapping and detrapping, surface reactions, diffusion
# Species variable: tritium_s (Ts), tritium_trap (trapped T), tritium_2g (T2(g))
# Species AuxVariable: --
# Chemistry: at the pore surface: (s) means in the solid, (g) means in gas form
# Diatomic molecule dissolving at surface
# T2(g) -> 2 T(s) rate K [T2] (1-theta)^2
# Diatomic molecule combining at surface
# 2 T(s) -> T2(g) rate K [T]^2
#
# BCs: No flux for tritium and tritium_trap, zero at pore center for gaseous tritium in every form
#
#
# Info:
# - This input file uses the microstructure generated by 2D_microstructure_reader_smoothing_base.i
# and performs tritium transport simulation on this microstructure.
# - Basic 2D non-dimensional input file for the transport of a solute in a microstructure with pores.
# - The goal of this input base and associated parameter files is to have the same simulation with different pore configurations
# to show the effect of microstructure on tritium release.
#
#
# Once TMAP8 is installed and built (see instructions on the TMAP8 website), run with
# cd ~/projects/TMAP8/test/tests/pore_scale_transport/
# mpirun -np 8 ~/projects/TMAP8/tmap8-opt -i 2D_absorption_base.i pore_structure_open_absorption.params
# Scaling factors
scale_quantity = ${units 1e-18 moles}
scale_time = ${units 1 s}
scale_length = ${units 1 mum}
tritium_trap_scaling = 1e3 # (-)
tritium_s_scaling = 1e2 # (-)
# Conditions
tritium_2g_concentration = ${units 0.1121 mol/mum^3}
# Material properties
tritium_trap_0 = ${units 3.472805925e-18 mol/mum^3}
tritium_trap_0_scaled = ${fparse tritium_trap_0*scale_length^3/scale_quantity} # non-dimensional
available_sites_mult = 4.313540691 # (-) should be >1, otherwise there are fewer 'available sites' than 'trapping sites'
Diffusion_min = ${units 0 mum^2/s} # used as diffusion coefficients where I don't want the species to move (more stable than 0)
Diffusion_coefficient_Db_D0 = ${units 1649969.728 mum^2/s}
Diffusion_coefficient_D_pore_D0 = ${units 224767738.6 mum^2/s}
detrapping_rate_0 = ${units -4.065104516 1/s} # should be negative
trapping_rate_0 = ${units 8.333146645e+16 1/s/mol} # should be positive
reactionRateSurface_sXYg_ref_0 = ${units 1549103878000000.0 1/s/mol}
reactionRateSurface_gXYs_ref_0 = ${units 34.69205021 1/s}
[Mesh]
file = ${input_name}
[]
[UserObjects]
[initial_microstructure]
type = SolutionUserObject
mesh = ${input_name}
timestep = LATEST
[]
[]
[Variables]
[tritium_s]
scaling = ${tritium_s_scaling}
[]
[tritium_trap]
scaling = ${tritium_trap_scaling}
[]
[tritium_2g]
[]
[]
[ICs]
[tritium_2g_ic]
type = FunctionIC
variable = tritium_2g
function = 'gaseous_concentration_1_function_ic'
[]
[]
[Functions]
[pore_position_func]
type = SolutionFunction
solution = initial_microstructure
from_variable = gr0
[]
[gaseous_concentration_1_function_ic]
type = ParsedFunction
symbol_names = 'len_pore position_pore_int'
symbol_values = '25 4521'
expression = '${tritium_2g_concentration}*0.5*(1.0+tanh((sqrt(x*x+y*y)-position_pore_int)/(len_pore/2.2)))'
[]
[]
[BCs]
[tritium_2g_sides_d] # Fix concentration on the sides
type = DirichletBC
variable = tritium_2g
boundary = 'right top'
value = ${tritium_2g_concentration}
[]
[]
[AuxVariables]
[pore]
initial_from_file_var = gr0
initial_from_file_timestep = LATEST
[]
# Used to prevent negative concentrations
[bounds_dummy_tritium_s]
order = FIRST
family = LAGRANGE
[]
[bounds_dummy_tritium_2g]
order = FIRST
family = LAGRANGE
[]
[]
[Bounds]
# To prevent negative concentrations
[tritium_s_lower_bound]
type = ConstantBounds
variable = bounds_dummy_tritium_s
bounded_variable = tritium_s
bound_type = lower
bound_value = 0.
[]
[tritium_g_lower_bound]
type = ConstantBounds
variable = bounds_dummy_tritium_2g
bounded_variable = tritium_2g
bound_type = lower
bound_value = 0.
[]
[]
[Kernels]
# ========================================== Time derivative for each variable
[dtritium_sdt]
type = TimeDerivative
variable = 'tritium_s'
[]
[dtritium_2gdt]
type = TimeDerivative
variable = 'tritium_2g'
[]
[dtritium_trapdt]
type = TimeDerivative
variable = 'tritium_trap'
[]
# ============================================================= Bulk diffusion
[Diff_tritium]
type = ADMatDiffusion
variable = tritium_s
diffusivity = D_tritium
[]
# ============================================================= Pore diffusion
[Diff_tritium2g]
type = ADMatDiffusion
variable = tritium_2g
diffusivity = D_pore
[]
# ================================================ Tritium trapping/detrapping
[Trapping]
type = CoupledTimeDerivative
variable = tritium_s
v = 'tritium_trap'
[]
[Detrapping_trap]
type = ADMatReaction
variable = 'tritium_trap'
reaction_rate = ${fparse scale_time * detrapping_rate_0} # from 1/s to non-dimensional
[]
[Trapping_trap]
type = ADMatCoupledDefectAnnihilation
variable = 'tritium_trap'
eq_concentration = ${tritium_trap_0_scaled}
v = tritium_s
annihilation_rate = ${fparse scale_time*scale_quantity/scale_length^3*trapping_rate_0} # from um^3/s/mol to non-dimensional
[]
# ========================= Surface Reactions ================================
# ======================================= Diatomic molecule dissolve at surface
# T2(g) -> 2 T(s) rate K [T2] (1-theta) ======================================
[Surface_Reaction_dissolution_tritium_2g] # T2(g) -> 2 T(s) rate -K*tritium_2g
type = ADMatReactionFlexible
variable = tritium_2g
vs = 'tritium_2g'
coeff = -1
reaction_rate_name = ReactionRateSurface_gXYs
[]
[Surface_Reaction_dissolution_tritium_2g_tritium_s] # T2(g) -> 2 T(s) rate 2*K*tritium_2g
type = ADMatReactionFlexible
variable = tritium_s
vs = 'tritium_2g'
coeff = 2
reaction_rate_name = ReactionRateSurface_gXYs
[]
# ======================================= Diatomic molecule combine at surface
# 2 T(s) -> T2(g) rate K [T]^2 ===============================================
[Surface_Reaction_release_tritium_s_tritium_2g] # 2 T(s) -> T2(g) rate -2*K*tritium_s^2
type = ADMatReactionFlexible
variable = tritium_s
vs = 'tritium_s tritium_s'
coeff = -2
reaction_rate_name = ReactionRateSurface_sXYg
[]
[Surface_Reaction_formation_tritium_2g] # 2 T(s) -> T2(g) rate K*tritium_s^2
type = ADMatReactionFlexible
variable = tritium_2g
vs = 'tritium_s tritium_s'
coeff = 1
reaction_rate_name = ReactionRateSurface_sXYg
[]
[]
[AuxKernels]
[pore_aux]
type = FunctionAux
variable = pore
function = 'pore_position_func'
execute_on = 'INITIAL'
[]
[]
[Materials]
#===================================================== Interpolation functions
[hsurfaceAD] # equal to 1 at pore surface, 0 elsewhere. The values of hsurface should be smaller than the ones for hmat or hpore to allow the diffusion of the species out of the surface
type = ADDerivativeParsedMaterial
coupled_variables = 'pore'
expression = '16*(1-pore)^2*pore^2'
property_name = 'hsurfaceAD'
[]
[hmatAD] # equal to 1 in material, and 0 in pore.
type = ADDerivativeParsedMaterial
coupled_variables = 'pore'
expression = '(1-pore)*(1-pore)'
property_name = 'hmatAD'
[]
[hporeAD] # equal to 1 in pore, and 1 in material. # notice that hpore overlaps a bit with hmat to prevent some agglomeration of chemicals on the edges of the surface.
type = ADDerivativeParsedMaterial
coupled_variables = 'pore'
expression = 'pore*pore'
property_name = 'hporeAD'
[]
#================================================= Local species concentration
[Tritium_ceramic]
type = ADParsedMaterial
property_name = 'tritium_ceramic'
coupled_variables = 'tritium_s tritium_trap'
expression = 'tritium_s+tritium_trap'
[]
[Tritium_pore]
type = ADParsedMaterial
property_name = 'tritium_pore'
coupled_variables = 'tritium_2g'
expression = '2*tritium_2g'
[]
[Tritium_total]
type = ADParsedMaterial
property_name = 'tritium_tot'
coupled_variables = 'tritium_s tritium_trap tritium_2g'
expression = 'tritium_s+tritium_trap+2*tritium_2g'
[]
#=============================================== Occupied sites on the surface
[theta]
type = ADParsedMaterial
property_name = 'theta'
coupled_variables = 'tritium_s tritium_trap'
expression = '(tritium_s + tritium_trap)/(${available_sites_mult}*${tritium_trap_0_scaled})'
[]
#====================================================== Diffusion coefficients
[Diffusion_coefficient_D] # from um^2/s to non-dimensional
type = ADDerivativeParsedMaterial
property_name = 'D_tritium'
coupled_variables = 'pore'
material_property_names = 'hmatAD(pore)'
expression = '${scale_time}/${scale_length}^2*(hmatAD*${Diffusion_coefficient_Db_D0} + (1-hmatAD)*${Diffusion_min})'
derivative_order = 2
[]
[Diffusion_coefficient_D_pore] # from um^2/s to non-dimensional
type = ADDerivativeParsedMaterial
property_name = 'D_pore'
material_property_names = 'hporeAD(pore)'
expression = '${scale_time}/${scale_length}^2*(hporeAD*${Diffusion_coefficient_D_pore_D0} + (1-hporeAD)*${Diffusion_min})'
[]
#============================================ Surface reaction rate of tritium
[ReactionRateSurface_sXYg_ref] # from surface to gaseous for diatomic molecules # from um^3/s/mol to non-dimensional
type = ADDerivativeParsedMaterial
property_name = ReactionRateSurface_sXYg
material_property_names = 'hsurfaceAD(pore)'
expression = '${scale_time}*${scale_quantity}/${scale_length}^3*${reactionRateSurface_sXYg_ref_0}*hsurfaceAD'
[]
[ReactionRateSurface_gXYs_ref] # from gaseous to surface for diatomic molecules # from 1/s to non-dimensional
type = ADDerivativeParsedMaterial
property_name = ReactionRateSurface_gXYs
coupled_variables = 'pore tritium_s tritium_trap'
material_property_names = 'hsurfaceAD(pore) theta(tritium_s,tritium_trap)'
expression = '${scale_time}*${reactionRateSurface_gXYs_ref_0}*hsurfaceAD*(1-theta)^2'
[]
[]
[Postprocessors]
[Tritium_amount_free]
type = ElementIntegralVariablePostprocessor
variable = tritium_s
execute_on = 'INITIAL TIMESTEP_END'
[]
[Tritium_amount_trapped]
type = ElementIntegralVariablePostprocessor
variable = tritium_trap
execute_on = 'INITIAL TIMESTEP_END'
[]
[Tritium_amount_ceramic]
type = ADElementIntegralMaterialProperty
mat_prop = tritium_ceramic
execute_on = 'INITIAL TIMESTEP_END'
[]
[Tritium2g_amount]
type = ElementIntegralVariablePostprocessor
variable = tritium_2g
execute_on = 'INITIAL TIMESTEP_END'
[]
[Tritium_amount_pore]
type = ADElementIntegralMaterialProperty
mat_prop = tritium_pore
execute_on = 'INITIAL TIMESTEP_END'
[]
[Tritium_amount_total]
type = ADElementIntegralMaterialProperty
mat_prop = tritium_tot
execute_on = 'INITIAL TIMESTEP_END'
[]
[Surface_tot]
type = ElementIntegralMaterialProperty
mat_prop = 1
[]
[dt]
type = TimestepSize
[]
[]
[Preconditioning]
[Newtonlu]
type = SMP
full = true
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -sub_pc_type -snes_type'
petsc_options_value = 'asm lu vinewtonrsls' # This petsc option helps prevent negative concentrations'
[]
[]
[Executioner]
type = Transient
scheme = bdf2
nl_rel_tol = 1e-10
end_time = 2000
dtmax = 50
nl_max_its = 14
[TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 12
iteration_window = 1
growth_factor = 1.2
dt = 2.37e-7
cutback_factor = 0.75
cutback_factor_at_failure = 0.75
[]
[]
[Outputs]
[exodus]
type = Exodus
time_step_interval = 5
[]
csv = true
file_base = ${output_name}
[]
(test/tests/ver-1gc/ver-1gc.i)
concentration_A_0 = ${units 2.415e14 at/m^3 -> at/mum^3} # atoms/microns^3 initial concentration of species A
k_1 = 0.0125 # 1/s reaction rate for first reaction
k_2 = 0.0025 # 1/s reaction rate for second reaction
end_time = 1500 # s
[Mesh]
type = GeneratedMesh
dim = 1
[]
[Variables]
[c_A]
initial_condition = ${concentration_A_0}
scaling = 1e4
[]
[c_B]
scaling = 1e4
[]
[c_C]
scaling = 1e4
[]
[]
[Kernels]
# Equation for species A
[timeDerivative_c_A]
type = ADTimeDerivative
variable = c_A
[]
[MatReaction_A]
type = ADMatReactionFlexible
variable = c_A
vs = 'c_A'
coeff = -1
reaction_rate_name = K1
[]
# Equation for species B
[timeDerivative_c_B]
type = TimeDerivative
variable = c_B
[]
[MatReaction_B_1]
type = ADMatReactionFlexible
variable = c_B
vs = 'c_A'
coeff = 1
reaction_rate_name = K1
[]
[MatReaction_B_2]
type = ADMatReactionFlexible
variable = c_B
vs = 'c_B'
coeff = -1
reaction_rate_name = K2
[]
# Equation for species C
[timeDerivative_c_C]
type = TimeDerivative
variable = c_C
[]
[MatReaction_C]
type = ADMatReactionFlexible
variable = c_C
vs = 'c_B'
coeff = 1
reaction_rate_name = K2
[]
[]
[Materials]
[K1]
type = ADParsedMaterial
property_name = 'K1'
expression = ${k_1}
[]
[K2]
type = ADParsedMaterial
property_name = 'K2'
expression = ${k_2}
[]
[]
[BCs]
[c_A_neumann] # No flux on the sides
type = NeumannBC
variable = c_A
boundary = 'left right'
value = 0
[]
[c_B_neumann] # No flux on the sides
type = NeumannBC
variable = c_B
boundary = 'left right'
value = 0
[]
[c_C_neumann] # No flux on the sides
type = NeumannBC
variable = c_C
boundary = 'left right'
value = 0
[]
[]
[Postprocessors]
[concentration_A]
type = ElementAverageValue
variable = c_A
execute_on = 'INITIAL TIMESTEP_END'
[]
[concentration_B]
type = ElementAverageValue
variable = c_B
execute_on = 'INITIAL TIMESTEP_END'
[]
[concentration_C]
type = ElementAverageValue
variable = c_C
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Executioner]
type = Transient
scheme = bdf2
nl_rel_tol = 1e-11
nl_abs_tol = 1e-50
l_tol = 1e-10
solve_type = 'NEWTON'
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
end_time = ${end_time}
dtmax = 50
[TimeStepper]
type = IterationAdaptiveDT
dt = 1.0
optimal_iterations = 15
growth_factor = 1.25
cutback_factor = 0.8
[]
[]
[Outputs]
exodus = true
csv = true
[]
(test/tests/val-2c/val-2c_base.i)
# Validation problem #2c from TMAP4/TMAP7 V&V document
# Test Cell Release Experiment based on
# D. F. Holland and R. A. Jalbert, "A Model for Tritium Concentration Following Tritium
# Release into a Test Cell and Subsequent Operation of an Atmospheric Cleanup Systen,"
# Proceedings, Eleventh Symposium of Fusion Engineering, November 18-22, 1985,. Austin,
# TX, Vol I, pp. 638-43, IEEE Cat. No. CH2251-7.
# Note that the approach to model this validation case is different in TMAP4 and TMAP7.
[Mesh]
[base_mesh]
type = GeneratedMeshGenerator
dim = 1
nx = ${fparse mesh_num_nodes_paint + 2}
xmax = ${length_domain}
[]
[subdomain_id]
input = base_mesh
type = SubdomainPerElementGenerator
subdomain_ids = '0 0 0 0 0 0 0 0 0 0 0 0 1 1'
[]
[interface]
type = SideSetsBetweenSubdomainsGenerator
input = subdomain_id
primary_block = '0' # paint
paired_block = '1' # enclosure
new_boundary = 'interface'
[]
[interface_other_side]
type = SideSetsBetweenSubdomainsGenerator
input = interface
primary_block = '1' # enclosure
paired_block = '0' # paint
new_boundary = 'interface_other'
[]
[]
[Variables]
# T2 concentration in the enclosure in molecules/microns^3
[t2_enclosure_concentration]
block = 1
[]
# HT concentration in the enclosure in molecules/microns^3
[ht_enclosure_concentration]
block = 1
[]
# HTO concentration in the enclosure in molecules/microns^3
[hto_enclosure_concentration]
block = 1
[]
# H2O concentration in the enclosure in molecules/microns^3
[h2o_enclosure_concentration]
block = 1
initial_condition = ${initial_H2O_concentration}
[]
# concentration of T2 in the paint in molecules/microns^3
[t2_paint_concentration]
block = 0
[]
# concentration of HT in the paint in molecules/microns^3
[ht_paint_concentration]
block = 0
[]
# concentration of HTO in the paint in molecules/microns^3
[hto_paint_concentration]
block = 0
[]
# concentration of H2O in the paint in molecules/microns^3
[h2o_paint_concentration]
block = 0
[]
[]
[AuxVariables]
# Used to prevent negative concentrations
[bounds_dummy_t2_paint_concentration]
order = FIRST
family = LAGRANGE
[]
[bounds_dummy_t2_enclosure_concentration]
order = FIRST
family = LAGRANGE
[]
[bounds_dummy_hto_paint_concentration]
order = FIRST
family = LAGRANGE
[]
[bounds_dummy_hto_enclosure_concentration]
order = FIRST
family = LAGRANGE
[]
[bounds_dummy_h2o_paint_concentration]
order = FIRST
family = LAGRANGE
[]
[bounds_dummy_h2o_enclosure_concentration]
order = FIRST
family = LAGRANGE
[]
[]
[Bounds]
# To prevent negative concentrations
[t2_paint_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy_t2_paint_concentration
bounded_variable = t2_paint_concentration
bound_type = lower
bound_value = ${lower_value_threshold}
[]
[t2_enclosure_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy_t2_enclosure_concentration
bounded_variable = t2_enclosure_concentration
bound_type = lower
bound_value = ${lower_value_threshold}
[]
[hto_paint_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy_hto_paint_concentration
bounded_variable = hto_paint_concentration
bound_type = lower
bound_value = ${lower_value_threshold}
[]
[hto_enclosure_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy_hto_enclosure_concentration
bounded_variable = hto_enclosure_concentration
bound_type = lower
bound_value = ${lower_value_threshold}
[]
[h2o_paint_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy_h2o_paint_concentration
bounded_variable = h2o_paint_concentration
bound_type = lower
bound_value = -1e-4
[]
[h2o_enclosure_concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy_h2o_enclosure_concentration
bounded_variable = h2o_enclosure_concentration
bound_type = lower
bound_value = -1e-4
[]
[]
[Kernels]
# In the enclosure
[t2_time_derivative]
type = TimeDerivative
variable = t2_enclosure_concentration
block = 1
[]
[t2_outflow]
type = MaskedBodyForce
variable = t2_enclosure_concentration
value = '-1'
mask = t2_enclosure_concentration_outflow
block = 1
[]
[t2_diffusion]
type = MatDiffusion
variable = t2_enclosure_concentration
block = 1
diffusivity = ${diffusivity_artificial_enclosure}
[]
[ht_time_derivative]
type = TimeDerivative
variable = ht_enclosure_concentration
block = 1
[]
[ht_outflow]
type = MaskedBodyForce
variable = ht_enclosure_concentration
value = '-1'
mask = ht_enclosure_concentration_outflow
block = 1
[]
[ht_diffusion]
type = MatDiffusion
variable = ht_enclosure_concentration
block = 1
diffusivity = ${diffusivity_artificial_enclosure}
[]
[hto_time_derivative]
type = TimeDerivative
variable = hto_enclosure_concentration
block = 1
[]
[hto_outflow]
type = MaskedBodyForce
variable = hto_enclosure_concentration
value = '-1'
mask = hto_enclosure_concentration_outflow
block = 1
[]
[hto_diffusion]
type = MatDiffusion
variable = hto_enclosure_concentration
block = 1
diffusivity = ${diffusivity_artificial_enclosure}
[]
[h2o_time_derivative]
type = TimeDerivative
variable = h2o_enclosure_concentration
block = 1
[]
[h2o_inflow]
type = MaskedBodyForce
variable = h2o_enclosure_concentration
mask = ${inflow_concentration}
block = 1
[]
[h2o_outflow]
type = MaskedBodyForce
variable = h2o_enclosure_concentration
value = '-1'
mask = h2o_enclosure_concentration_outflow
block = 1
[]
[h2o_diffusion]
type = MatDiffusion
variable = h2o_enclosure_concentration
block = 1
diffusivity = ${diffusivity_artificial_enclosure}
[]
# reaction T2+H2O->HTO+HT
[reaction_1_t2]
type = ADMatReactionFlexible
variable = t2_enclosure_concentration
block = 1
coeff = -1
reaction_rate_name = reaction_rate_t2
[]
[reaction_1_h2o]
type = ADMatReactionFlexible
variable = h2o_enclosure_concentration
block = 1
coeff = -1
reaction_rate_name = reaction_rate_t2
[]
[reaction_1_hto]
type = ADMatReactionFlexible
variable = hto_enclosure_concentration
block = 1
coeff = 1
reaction_rate_name = reaction_rate_t2
[]
[reaction_1_ht]
type = ADMatReactionFlexible
variable = ht_enclosure_concentration
block = 1
coeff = 1
reaction_rate_name = reaction_rate_t2
[]
# reaction HT+H2O->HTO+H2
[reaction_2_HT]
type = ADMatReactionFlexible
variable = ht_enclosure_concentration
block = 1
coeff = -1
reaction_rate_name = reaction_rate_ht
[]
[reaction_2_h2o]
type = ADMatReactionFlexible
variable = h2o_enclosure_concentration
block = 1
coeff = -1
reaction_rate_name = reaction_rate_ht
[]
[reaction_2_hto]
type = ADMatReactionFlexible
variable = hto_enclosure_concentration
block = 1
coeff = 1
reaction_rate_name = reaction_rate_ht
[]
# In the paint
[t2_paint_time]
type = TimeDerivative
variable = t2_paint_concentration
block = 0
[]
[t2_paint_diffusion]
type = MatDiffusion
variable = t2_paint_concentration
block = 0
diffusivity = '${diffusivity_elemental_tritium}'
[]
[ht_paint_time]
type = TimeDerivative
variable = ht_paint_concentration
block = 0
[]
[ht_paint_diffusion]
type = MatDiffusion
variable = ht_paint_concentration
block = 0
diffusivity = '${diffusivity_elemental_tritium}'
[]
[hto_paint_time]
type = TimeDerivative
variable = hto_paint_concentration
block = 0
[]
[hto_paint_diffusion]
type = MatDiffusion
variable = hto_paint_concentration
block = 0
diffusivity = '${diffusivity_tritiated_water}'
[]
[h2o_paint_time]
type = TimeDerivative
variable = h2o_paint_concentration
block = 0
[]
[h2o_paint_diffusion]
type = MatDiffusion
variable = h2o_paint_concentration
block = 0
diffusivity = '${diffusivity_tritiated_water}'
[]
[]
[InterfaceKernels]
# solubility at the surface of the paint
[t2_solubility]
type = ADInterfaceSorption
variable = t2_paint_concentration
neighbor_var = t2_enclosure_concentration # molecules/microns^3
unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
K0 = ${solubility_elemental_tritium}
Ea = 0
n_sorption = 1 # Henry's law
temperature = ${temperature}
diffusivity = diffusivity_elemental_tritium
boundary = 'interface'
[]
[ht_solubility]
type = ADInterfaceSorption
variable = ht_paint_concentration
neighbor_var = ht_enclosure_concentration # molecules/microns^3
unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
K0 = ${solubility_elemental_tritium}
Ea = 0
n_sorption = 1 # Henry's law
temperature = ${temperature}
diffusivity = diffusivity_elemental_tritium
boundary = 'interface'
[]
[hto_solubility]
type = ADInterfaceSorption
variable = hto_paint_concentration
neighbor_var = hto_enclosure_concentration # molecules/microns^3
unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
K0 = ${solubility_tritiated_water}
Ea = 0
n_sorption = 1 # Henry's law
temperature = ${temperature}
diffusivity = diffusivity_tritiated_water
boundary = 'interface'
[]
[h2o_solubility]
type = ADInterfaceSorption
variable = h2o_paint_concentration
neighbor_var = h2o_enclosure_concentration # molecules/microns^3
unit_scale_neighbor = ${fparse 1e18/NA} # to convert neighbor concentration to mol/m^3
K0 = ${solubility_tritiated_water}
Ea = 0
n_sorption = 1 # Henry's law
temperature = ${temperature}
diffusivity = diffusivity_tritiated_water
boundary = 'interface'
[]
[]
[Materials]
[reaction_rate_t2]
type = ADDerivativeParsedMaterial
coupled_variables = 't2_enclosure_concentration ht_enclosure_concentration hto_enclosure_concentration'
expression = '${reaction_rate} * 2 * t2_enclosure_concentration * (2*t2_enclosure_concentration + ht_enclosure_concentration + hto_enclosure_concentration)'
property_name = reaction_rate_t2
block = 1
[]
[reaction_rate_ht]
type = ADDerivativeParsedMaterial
coupled_variables = 't2_enclosure_concentration ht_enclosure_concentration hto_enclosure_concentration'
expression = '${reaction_rate} * ht_enclosure_concentration * (2*t2_enclosure_concentration + ht_enclosure_concentration + hto_enclosure_concentration)'
property_name = reaction_rate_ht
block = 1
[]
[diffusivity_elemental_tritium_paint]
type = ADConstantMaterial
value = ${diffusivity_elemental_tritium}
property_name = 'diffusivity_elemental_tritium'
block = 0
[]
[diffusivity_elemental_tritium_enclosure]
type = ADConstantMaterial
value = ${diffusivity_artificial_enclosure}
property_name = 'diffusivity_elemental_tritium'
block = 1
[]
[diffusivity_tritiated_water_paint]
type = ADConstantMaterial
value = ${diffusivity_tritiated_water}
property_name = 'diffusivity_tritiated_water'
block = 0
[]
[diffusivity_tritiated_water_enclosure]
type = ADConstantMaterial
value = ${diffusivity_artificial_enclosure}
property_name = 'diffusivity_tritiated_water'
block = 1
[]
[t2_enclosure_concentration_outflow]
type = DerivativeParsedMaterial
coupled_variables = 't2_enclosure_concentration'
expression = 't2_enclosure_concentration * ${outflow} / ${volume_enclosure}'
property_name = 't2_enclosure_concentration_outflow'
block = 1
[]
[ht_enclosure_concentration_outflow]
type = DerivativeParsedMaterial
coupled_variables = 'ht_enclosure_concentration'
expression = 'ht_enclosure_concentration * ${outflow} / ${volume_enclosure}'
property_name = 'ht_enclosure_concentration_outflow'
block = 1
[]
[hto_enclosure_concentration_outflow]
type = DerivativeParsedMaterial
coupled_variables = 'hto_enclosure_concentration'
expression = 'hto_enclosure_concentration * ${outflow} / ${volume_enclosure}'
property_name = 'hto_enclosure_concentration_outflow'
block = 1
[]
[h2o_enclosure_concentration_outflow]
type = DerivativeParsedMaterial
coupled_variables = 'h2o_enclosure_concentration'
expression = 'h2o_enclosure_concentration * ${outflow} / ${volume_enclosure}'
property_name = 'h2o_enclosure_concentration_outflow'
block = 1
[]
[]
[Postprocessors]
# Pressures in enclosure
[t2_enclosure_edge_concentration] # (molecules/microns^3)
type = PointValue
point = '${length_domain} 0 0' # on the far side of the enclosure
variable = t2_enclosure_concentration
execute_on = 'initial timestep_end'
[]
[ht_enclosure_edge_concentration] # (molecules/microns^3)
type = PointValue
point = '${length_domain} 0 0' # on the far side of the enclosure
variable = ht_enclosure_concentration
execute_on = 'initial timestep_end'
[]
[hto_enclosure_edge_concentration] # (molecules/microns^3)
type = PointValue
point = '${length_domain} 0 0' # on the far side of the enclosure
variable = hto_enclosure_concentration
execute_on = 'initial timestep_end'
[]
[h2o_enclosure_edge_concentration] # (molecules/microns^3)
type = PointValue
point = '${length_domain} 0 0' # on the far side of the enclosure
variable = h2o_enclosure_concentration
execute_on = 'initial timestep_end'
[]
# Inventory in enclosure
[t2_enclosure_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = t2_enclosure_concentration
execute_on = 'initial timestep_end'
block = 1
[]
[ht_enclosure_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = ht_enclosure_concentration
execute_on = 'initial timestep_end'
block = 1
[]
[hto_enclosure_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = hto_enclosure_concentration
execute_on = 'initial timestep_end'
block = 1
[]
[h2o_enclosure_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = h2o_enclosure_concentration
execute_on = 'initial timestep_end'
block = 1
[]
# Inventory in paint
[t2_paint_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = t2_paint_concentration
execute_on = 'initial timestep_end'
block = 0
[]
[ht_paint_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = ht_paint_concentration
execute_on = 'initial timestep_end'
block = 0
[]
[hto_paint_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = hto_paint_concentration
execute_on = 'initial timestep_end'
block = 0
[]
[h2o_paint_inventory] # (molecules/microns^2)
type = ElementIntegralVariablePostprocessor
variable = h2o_paint_concentration
execute_on = 'initial timestep_end'
block = 0
[]
[tritium_total_inventory]
type = LinearCombinationPostprocessor
pp_names = 't2_enclosure_inventory ht_enclosure_inventory hto_enclosure_inventory t2_paint_inventory ht_paint_inventory hto_paint_inventory'
pp_coefs = '2 1 1 2 1 1'
execute_on = 'initial timestep_end'
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
scheme = 'bdf2'
petsc_options_iname = '-pc_type -sub_pc_type -snes_type'
petsc_options_value = 'asm lu vinewtonrsls' # This petsc option helps prevent negative concentrations with bounds'
nl_rel_tol = 1.e-10
automatic_scaling = true
compute_scaling_once = false
end_time = ${time_end}
dtmax = ${dtmax}
nl_max_its = 16
[TimeStepper]
type = IterationAdaptiveDT
dt = ${time_step}
optimal_iterations = 12
iteration_window = 1
growth_factor = 1.1
cutback_factor = 0.9
[]
[]
[Outputs]
exodus = true
perf_graph = true
[csv]
type = CSV
execute_on = 'initial timestep_end'
[]
[]
(test/tests/kernels/ad_mat_reaction_flexible/ad_mat_reaction_flexible_test.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 10
ny = 10
[]
[Variables]
[c_a]
[]
[c_b]
[]
[]
[ICs]
[c_a_IC]
type = ConstantIC
variable = c_a
value = 1
[]
[c_b_IC]
type = ConstantIC
variable = c_b
value = 1
[]
[]
[Kernels]
[timeDerivative_c_a]
type = ADTimeDerivative
variable = c_a
[]
[timeDerivative_c_b]
type = TimeDerivative
variable = c_b
[]
[MatReaction]
type = ADMatReactionFlexible
variable = c_b
vs = 'c_a'
coeff = 0.5
reaction_rate_name = K
[]
[]
[Materials]
[K]
type = ADParsedMaterial
property_name = 'K'
expression = '10'
[]
[]
[BCs]
[c_a_neumann] # No flux on the sides
type = NeumannBC
variable = c_a
boundary = 'left right bottom top'
value = 0
[]
[c_b_neumann] # No flux on the sides
type = NeumannBC
variable = c_b
boundary = 'left right bottom top'
value = 0
[]
[]
[Executioner]
type = Transient
scheme = bdf2
nl_rel_tol = 1e-10
solve_type = 'NEWTON'
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
start_time = 0.0
end_time = 1.
num_steps = 60000
dt = .2
n_startup_steps = 0
[]
[Outputs]
exodus = 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-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-1kd/ver-1kd.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 2 s}'
N_a = '${units 6.02214076e23 1/mol}' # Avogadro's number from PhysicalConstants.h
source = '${units ${fparse 1e23 / N_a} mol/m^3/s}' # Source term for T2 in enclosure 1
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]
# Source term for T2 in enclosure 1
[T2_source_term_enclosure_1]
type = BodyForce
variable = concentration_T2_enclosure_1
block = '1'
value = ${source}
[]
# 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-1kd_out_k10'
csv = true
exodus = true
execute_on = 'initial timestep_end'
[]
(test/tests/ver-1g/ver-1g.i)
## Run this file with one of three input files to simulate various verification cases
## a. equal_conc.i -> TMAP4 and TMAP7 equal concentration case
## b. diff_conc_TMAP4.i -> TMAP4 different concentration case
## c. diff_conc_TMAP7.i -> TMAP7 different concentration case
## Example: ~/projects/TMAP8/tmap8-opt -i ver-1g.i equal_conc.i
R = 8.31446261815324 # Gas constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
T = '${units 25 degC -> K}' # Temperature
Na = 6.02214076E23 # Avogadro's constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?na)
[Mesh]
type = GeneratedMesh
dim = 2
[]
[Variables]
[c_a]
[]
[c_b]
[]
[c_ab]
initial_condition = 0
[]
[]
[Kernels]
[timeDerivative_c_a]
type = ADTimeDerivative
variable = c_a
[]
[timeDerivative_c_b]
type = TimeDerivative
variable = c_b
[]
[timeDerivative_c_ab]
type = TimeDerivative
variable = c_ab
[]
[MatReaction_b]
type = ADMatReactionFlexible
variable = c_b
vs = 'c_a c_b'
coeff = -1
reaction_rate_name = K
[]
[MatReaction_a]
type = ADMatReactionFlexible
variable = c_a
vs = 'c_b c_a'
coeff = -1
reaction_rate_name = K
[]
[MatReaction_ab]
type = ADMatReactionFlexible
variable = c_ab
vs = 'c_b c_a'
coeff = 1
reaction_rate_name = K
[]
[]
[Materials]
[K]
type = ADParsedMaterial
property_name = 'K'
expression = '4.14e3' # units: molecule.micrometer^3/atom/second
[]
[]
[BCs]
[c_a_neumann] # No flux on the sides
type = NeumannBC
variable = c_a
boundary = 'left right bottom top'
value = 0
[]
[c_b_neumann] # No flux on the sides
type = NeumannBC
variable = c_b
boundary = 'left right bottom top'
value = 0
[]
[c_ab_neumann] # No flux on the sides
type = NeumannBC
variable = c_ab
boundary = 'left right bottom top'
value = 0
[]
[]
[Postprocessors]
[conc_a]
type = ElementAverageValue
variable = c_a
[]
[conc_b]
type = ElementAverageValue
variable = c_b
[]
[conc_ab]
type = ElementAverageValue
variable = c_ab
[]
[]
[Executioner]
type = Transient
scheme = bdf2
nl_rel_tol = 1e-10
nl_abs_tol = 1e-15
solve_type = 'NEWTON'
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
start_time = 0.0
end_time = 40
num_steps = 60000
dt = .2
n_startup_steps = 0
[]
[Outputs]
exodus = true
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/val-2e/val-2ee.i)
# Validation Problem #2ee from TMAP7's V&V document
# Deuterium permeation through 0.05-mm Pd at 825 K.
# No Soret effect, or trapping included.
# Necessary physical and model parameters (kb, R, temperature)
!include parameters_three_gases.params
# Surface reaction data used in TMAP7 case
K_r_pre_D2 = '${units ${fparse 2.502e-24 / sqrt(4 * temperature)} m^4/s -> mum^4/s}'
K_r_pre_H2 = '${units ${fparse 2.502e-24 / sqrt(2 * temperature)} m^4/s -> mum^4/s}'
K_r_pre_HD = '${units ${fparse 2.502e-24 / sqrt(3 * temperature)} m^4/s -> mum^4/s}'
K_r_energy = '${units ${fparse -11836 * R} J/mol}'
K_d_D2 = '${units ${fparse 2.1897e22 / sqrt(4 * temperature)} at/m^2/Pa -> at/mum^2/Pa}'
K_d_H2 = '${units ${fparse 2.1897e22 / sqrt(2 * temperature)} at/m^2/Pa -> at/mum^2/Pa}'
K_d_HD = '${units ${fparse 2.1897e22 / sqrt(3 * temperature)} at/m^2/Pa -> at/mum^2/Pa}'
# Modeling data used in current case
file_name = 'val-2ee_out'
!include val-2e_base_three_gases.i
[Kernels]
# Gas flow kernels
# Equation for D2 in enclosure upstream
[MatReaction_upstream_D2_reaction]
type = ADMatReactionFlexible
variable = D2_pressure_upstream
vs = 'D_concentration H_concentration'
reaction_rate_name = 'K_r_HD'
coeff = '${fparse -0.5 * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
[MatReaction_upstream_D2_re_reaction]
type = ADMatReaction
variable = D2_pressure_upstream
v = 'HD_pressure_upstream'
reaction_rate = '${fparse 0.5 * K_d_HD * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
# Equation for H2 in enclosure upstream
[MatReaction_upstream_H2_reaction]
type = ADMatReactionFlexible
variable = H2_pressure_upstream
vs = 'D_concentration H_concentration'
reaction_rate_name = 'K_r_HD'
coeff = '${fparse -0.5 * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
[MatReaction_upstream_H2_re_reaction]
type = ADMatReaction
variable = H2_pressure_upstream
v = 'HD_pressure_upstream'
reaction_rate = '${fparse 0.5 * K_d_HD * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
# Equation for HD in enclosure upstream
[MatReaction_upstream_HD_reaction]
type = ADMatReactionFlexible
variable = HD_pressure_upstream
vs = 'D_concentration H_concentration'
reaction_rate_name = 'K_r_HD'
coeff = '${fparse kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
[MatReaction_upstream_HD_re_reaction]
type = ADMatReaction
variable = HD_pressure_upstream
v = 'HD_pressure_upstream'
reaction_rate = '${fparse - K_d_HD * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
# Membrane upstream
[MatReaction_upstream_outflux_membrane_D2]
type = ADMatBodyForce
variable = D2_pressure_upstream
material_property = 'membrane_reaction_rate_right_D2'
extra_vector_tags = 'ref'
[]
[MatReaction_upstream_outflux_membrane_H2]
type = ADMatBodyForce
variable = H2_pressure_upstream
material_property = 'membrane_reaction_rate_right_H2'
extra_vector_tags = 'ref'
[]
[MatReaction_upstream_outflux_membrane_HD]
type = ADMatBodyForce
variable = HD_pressure_upstream
material_property = 'membrane_reaction_rate_right_HD'
extra_vector_tags = 'ref'
[]
# Equation for D2 enclosure downstream
[MatReaction_downstream_D2_reaction]
type = ADMatReactionFlexible
variable = D2_pressure_downstream
vs = 'D_concentration H_concentration'
reaction_rate_name = 'K_r_HD'
coeff = '${fparse -0.5 * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
[MatReaction_downstream_D2_re_reaction]
type = ADMatReaction
variable = D2_pressure_downstream
v = 'HD_pressure_downstream'
reaction_rate = '${fparse 0.5 * K_d_HD * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
# Equation for H2 enclosure downstream
[MatReaction_downstream_H2_reaction]
type = ADMatReactionFlexible
variable = H2_pressure_downstream
vs = 'D_concentration H_concentration'
reaction_rate_name = 'K_r_HD'
coeff = '${fparse -0.5 * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
[MatReaction_downstream_H2_re_reaction]
type = ADMatReaction
variable = H2_pressure_downstream
v = 'HD_pressure_downstream'
reaction_rate = '${fparse 0.5 * K_d_HD * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
# Equation for HD enclosure downstream
[MatReaction_downstream_HD_reaction]
type = ADMatReactionFlexible
variable = HD_pressure_downstream
vs = 'D_concentration H_concentration'
reaction_rate_name = 'K_r_HD'
coeff = '${fparse kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
[MatReaction_downstream_HD_re_reaction]
type = ADMatReaction
variable = HD_pressure_downstream
v = 'HD_pressure_downstream'
reaction_rate = '${fparse - K_d_HD * kb * temperature * surface_area / volume_enclosure}'
extra_vector_tags = 'ref'
[]
# Membrane downstream
[MatReaction_downstream_influx_membrane_D2]
type = ADMatBodyForce
variable = D2_pressure_downstream
material_property = 'membrane_reaction_rate_left_D2'
extra_vector_tags = 'ref'
[]
[MatReaction_downstream_influx_membrane_H2]
type = ADMatBodyForce
variable = H2_pressure_downstream
material_property = 'membrane_reaction_rate_left_H2'
extra_vector_tags = 'ref'
[]
[MatReaction_downstream_influx_membrane_HD]
type = ADMatBodyForce
variable = HD_pressure_downstream
material_property = 'membrane_reaction_rate_left_HD'
extra_vector_tags = 'ref'
[]
[]
[Materials]
[K_r_D2]
type = ADParsedMaterial
property_name = 'K_r_D2'
expression = '${K_r_pre_D2} * exp( - ${K_r_energy} / ${R} / ${temperature})'
[]
[K_r_H2]
type = ADParsedMaterial
property_name = 'K_r_H2'
expression = '${K_r_pre_H2} * exp( - ${K_r_energy} / ${R} / ${temperature})'
[]
[K_r_HD]
type = ADParsedMaterial
property_name = 'K_r_HD'
expression = '${K_r_pre_HD} * exp( - ${K_r_energy} / ${R} / ${temperature})'
[]
[membrane_reaction_rate_right_D2]
type = ADParsedMaterial
property_name = 'membrane_reaction_rate_right_D2'
material_property_names = flux_on_right_D2
expression = '-flux_on_right_D2 * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}'
[]
[membrane_reaction_rate_left_D2]
type = ADParsedMaterial
property_name = 'membrane_reaction_rate_left_D2'
material_property_names = flux_on_left_D2
expression = '-flux_on_left_D2 * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}'
[]
[membrane_reaction_rate_right_H2]
type = ADParsedMaterial
property_name = 'membrane_reaction_rate_right_H2'
material_property_names = flux_on_right_H2
expression = '-flux_on_right_H2 * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}'
[]
[membrane_reaction_rate_left_H2]
type = ADParsedMaterial
property_name = 'membrane_reaction_rate_left_H2'
material_property_names = flux_on_left_H2
expression = '-flux_on_left_H2 * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}'
[]
[membrane_reaction_rate_right_HD]
type = ADParsedMaterial
property_name = 'membrane_reaction_rate_right_HD'
material_property_names = flux_on_right_HD
expression = '-flux_on_right_HD * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}'
[]
[membrane_reaction_rate_left_HD]
type = ADParsedMaterial
property_name = 'membrane_reaction_rate_left_HD'
material_property_names = flux_on_left_HD
expression = '-flux_on_left_HD * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}'
[]
[flux_on_left_D]
type = ADParsedMaterial
coupled_variables = 'D_concentration H_concentration D2_pressure_downstream HD_pressure_downstream'
property_name = 'flux_on_left_D'
material_property_names = 'K_r_D2 K_r_HD'
expression = '- 2 * K_r_D2 * D_concentration ^ 2 + 2 * ${K_d_D2} * D2_pressure_downstream - K_r_HD * D_concentration * H_concentration + ${K_d_HD} * HD_pressure_downstream'
outputs = 'exodus'
[]
[flux_on_right_D]
type = ADParsedMaterial
coupled_variables = 'D_concentration H_concentration D2_pressure_upstream HD_pressure_upstream'
property_name = 'flux_on_right_D'
material_property_names = 'K_r_D2 K_r_HD'
expression = '- 2 * K_r_D2 * D_concentration ^ 2 + 2 * ${K_d_D2} * D2_pressure_upstream - K_r_HD * D_concentration * H_concentration + ${K_d_HD} * HD_pressure_upstream'
outputs = 'exodus'
[]
[flux_on_left_H]
type = ADParsedMaterial
coupled_variables = 'D_concentration H_concentration H2_pressure_downstream HD_pressure_downstream'
property_name = 'flux_on_left_H'
material_property_names = 'K_r_H2 K_r_HD'
expression = '- 2 * K_r_H2 * H_concentration ^ 2 + 2 * ${K_d_H2} * H2_pressure_downstream - K_r_HD * D_concentration * H_concentration + ${K_d_HD} * HD_pressure_downstream'
outputs = 'exodus'
[]
[flux_on_right_H]
type = ADParsedMaterial
coupled_variables = 'D_concentration H_concentration H2_pressure_upstream HD_pressure_upstream'
property_name = 'flux_on_right_H'
material_property_names = 'K_r_H2 K_r_HD'
expression = '- 2 * K_r_H2 * H_concentration ^ 2 + 2 * ${K_d_H2} * H2_pressure_upstream - K_r_HD * D_concentration * H_concentration + ${K_d_HD} * HD_pressure_upstream'
outputs = 'exodus'
[]
[flux_on_left_D2]
type = ADParsedMaterial
coupled_variables = 'D_concentration D2_pressure_downstream'
property_name = 'flux_on_left_D2'
material_property_names = 'K_r_D2'
expression = '- K_r_D2 * D_concentration ^ 2 + ${K_d_D2} * D2_pressure_downstream'
outputs = 'exodus'
[]
[flux_on_left_HD]
type = ADParsedMaterial
coupled_variables = 'D_concentration H_concentration HD_pressure_downstream'
property_name = 'flux_on_left_HD'
material_property_names = 'K_r_HD'
expression = '- K_r_HD * D_concentration * H_concentration + ${K_d_HD} * HD_pressure_downstream'
outputs = 'exodus'
[]
[flux_on_left_H2]
type = ADParsedMaterial
coupled_variables = 'H_concentration H2_pressure_downstream'
property_name = 'flux_on_left_H2'
material_property_names = 'K_r_H2'
expression = '- K_r_H2 * H_concentration ^ 2 + ${K_d_H2} * H2_pressure_downstream'
outputs = 'exodus'
[]
[flux_on_right_D2]
type = ADParsedMaterial
coupled_variables = 'D_concentration D2_pressure_upstream'
property_name = 'flux_on_right_D2'
material_property_names = 'K_r_D2'
expression = '- K_r_D2 * D_concentration ^ 2 + ${K_d_D2} * D2_pressure_upstream'
outputs = 'exodus'
[]
[flux_on_right_HD]
type = ADParsedMaterial
coupled_variables = 'D_concentration H_concentration HD_pressure_upstream'
property_name = 'flux_on_right_HD'
material_property_names = 'K_r_HD'
expression = '- K_r_HD * D_concentration * H_concentration + ${K_d_HD} * HD_pressure_upstream'
outputs = 'exodus'
[]
[flux_on_right_H2]
type = ADParsedMaterial
coupled_variables = 'H_concentration H2_pressure_upstream'
property_name = 'flux_on_right_H2'
material_property_names = 'K_r_H2'
expression = '- K_r_H2 * H_concentration ^ 2 + ${K_d_H2} * H2_pressure_upstream'
outputs = 'exodus'
[]
[]
[BCs]
# recombination BCs
[left_diffusion_D]
type = ADMatNeumannBC
variable = D_concentration
boundary = left
value = 1
boundary_material = flux_on_left_D
[]
[right_diffusion_D]
type = ADMatNeumannBC
variable = D_concentration
boundary = right
value = 1
boundary_material = flux_on_right_D
[]
[left_diffusion_H]
type = ADMatNeumannBC
variable = H_concentration
boundary = left
value = 1
boundary_material = flux_on_left_H
[]
[right_diffusion_H]
type = ADMatNeumannBC
variable = H_concentration
boundary = right
value = 1
boundary_material = flux_on_right_H
[]
[]
[Postprocessors]
# Flux
[dcdx_right_D]
type = ADSideAverageMaterialProperty
boundary = right
property = flux_on_right_D
outputs = 'console csv exodus'
[]
[dcdx_left_D]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_left_D
outputs = 'console csv exodus'
[]
[dcdx_right_H]
type = ADSideAverageMaterialProperty
boundary = right
property = flux_on_right_H
outputs = 'console csv exodus'
[]
[dcdx_left_H]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_left_H
outputs = 'console csv exodus'
[]
[flux_on_left_H2]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_left_H2
outputs = 'console csv exodus'
[]
[flux_on_left_D2]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_left_D2
outputs = 'console csv exodus'
[]
[flux_on_left_HD]
type = ADSideAverageMaterialProperty
boundary = left
property = flux_on_left_HD
outputs = 'console csv exodus'
[]
[flux_on_right_H2]
type = ADSideAverageMaterialProperty
boundary = right
property = flux_on_right_H2
outputs = 'console csv exodus'
[]
[flux_on_right_D2]
type = ADSideAverageMaterialProperty
boundary = right
property = flux_on_right_D2
outputs = 'console csv exodus'
[]
[flux_on_right_HD]
type = ADSideAverageMaterialProperty
boundary = right
property = flux_on_right_HD
outputs = 'console csv exodus'
[]
[]