EquilibriumBC

Enforces a species equilibrium condition between an enclosure and an adjacent diffusion structure.

Overview

This class strongly enforces the Dirichlet boundary condition

$var element = document.getElementById("moose-equation-c11dcab2-79ea-4ea6-a14a-2e3ed053f377");katex.render("C_i = K_o \\exp{\\left( \\frac{-E_a}{RT}\\right)} P_i^p", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-387bd5f8-8078-49aa-974f-febf67dc0132");katex.render("C_i", element, {displayMode:false,throwOnError:false});$ , represented by the variable parameter, is the concentration of specie $var element = document.getElementById("moose-equation-58ed74fe-5116-47ab-9155-33c8c8902389");katex.render("i", element, {displayMode:false,throwOnError:false});$ in a diffusion structure, $var element = document.getElementById("moose-equation-f3aae43e-e842-4d19-974a-227f07039d49");katex.render("P_i", element, {displayMode:false,throwOnError:false});$ is the partial pressure of specie $var element = document.getElementById("moose-equation-4946a486-b935-4e7d-9115-752f993459b3");katex.render("i", element, {displayMode:false,throwOnError:false});$ in the gas phase in the enclosure adjacent to the diffusion structure, $var element = document.getElementById("moose-equation-39d7725e-b9ed-400e-971b-46500cf7718f");katex.render("K_o", element, {displayMode:false,throwOnError:false});$ is a solubility constant, $var element = document.getElementById("moose-equation-c69b8b7d-9a34-4507-8895-f5df15a8445a");katex.render("E_a", element, {displayMode:false,throwOnError:false});$ is the activation energy, $var element = document.getElementById("moose-equation-507fed67-372f-4417-8d9d-0ff6ad258763");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the universal gas constant, $var element = document.getElementById("moose-equation-59a8f054-aa97-4af0-801a-748418a9c968");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature and $var element = document.getElementById("moose-equation-3c5cb6c8-8685-4dfd-9a02-186739baf811");katex.render("p", element, {displayMode:false,throwOnError:false});$ is the exponent of the solution law.

Input Parameters

KoThe solubility coefficient $Ko$ for the relationship $C_i = Ko exp{-Ea/RT} P_i^p$
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The solubility coefficient $Ko$ for the relationship $C_i = Ko exp{-Ea/RT} P_i^p$
boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
enclosure_varThe coupled enclosure variable $P_i$ in the relationship $C_i = Ko exp{-Ea/RT} P_i^p$. Can be a either a field or scalar variable.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled enclosure variable $P_i$ in the relationship $C_i = Ko exp{-Ea/RT} P_i^p$. Can be a either a field or scalar variable.
temperatureThe temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The temperature
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

Required Parameters

activation_energy0The activation energy $Ea$ for the relationship $C_i = Ko exp{-Ea/RT} P_i^p$ (J/mol)
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The activation energy $Ea$ for the relationship $C_i = Ko exp{-Ea/RT} P_i^p$ (J/mol)
diag_save_inThe name of auxiliary variables to save this BC'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 BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
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)
p1The exponent $p$ in the relationship $C_i = Ko exp{-Ea/RT} P_i^p$
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The exponent $p$ in the relationship $C_i = Ko exp{-Ea/RT} P_i^p$
save_inThe name of auxiliary variables to save this BC'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 BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
var_scaling_factor1The number of atoms that compose our arbitrary unit for quantity
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The number of atoms that compose our arbitrary unit for quantity

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

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.
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
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

matrix_tagssystem timeThe tag for the matrices this Kernel should fill
Default:system time
C++ Type:MultiMooseEnum
Options:nontime, system, time
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsresidualThe tag for the vectors this Kernel should fill
Default:residual
C++ Type:MultiMooseEnum
Options:nontime, time, residual
Controllable:No
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

Input Files

(test/tests/val-2e/val-2ea.i)
(test/tests/val-2b/val-2b.i)
(test/tests/bcs/EquilibriumBC/equilibriumBC.i)
(test/tests/ver-1a/ver-1a.i)
(test/tests/val-2e/val-2ed.i)

(test/tests/val-2e/val-2ea.i)

# Validation Problem #2ea 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_deuterium_gas.params

# Enclosure data used in TMAP7 case
pressure_enclosure1 = '${units 1e-6 Pa}'
pressure_initial = '${units 1e-6 Pa}'

# Modeling data used in current case
slab_thickness = '${units 5e-5 m -> mum}'
file_name = 'val-2ea_out'
simulation_time = '${units 1900 s}'

!include val-2e_base.i

# Diffusion data used in TMAP7 case
solubility_exponent = 0.9297 # -
solubility_pre = '${units ${fparse 1.511e23 / 1e18} at/mum^3/Pa^0.9297}'
solubility_energy = '${units ${fparse 5918 * R} J/mol}'

[Variables]
  [D2_pressure_upstream]
    initial_condition = '${pressure_initial}'
  []
  [D2_pressure_downstream]
    initial_condition = '${pressure_initial}'
  []
  [D_concentration]
    initial_condition = 0
  []
[]

[AuxVariables]
  [D2_pressure_enclosure1]
    initial_condition = '${pressure_enclosure1}'
  []
  [D2_pressure_enclosure4]
    initial_condition = '${pressure_enclosure4}'
  []
  [D2_pressure_enclosure5]
  []
[]

[AuxKernels]
  [D2_pressure_enclosure5_kernel]
    type = FunctionAux
    variable = D2_pressure_enclosure5
    function = D2_pressure_enclosure5_function
  []
[]

[Functions]
  [D2_pressure_enclosure5_function]
    type = ParsedFunction
    expression = 'if(t < 150, 1.20e-4,
                  if(t < 250, 2.41e-4,
                  if(t < 350, 6.06e-4,
                  if(t < 450, 1.30e-3,
                  if(t < 550, 2.53e-3,
                  if(t < 650, 7.08e-3,
                  if(t < 750, 1.45e-2,
                  if(t < 850, 2.63e-2,
                  if(t < 950, 6.51e-2,
                  if(t < 1050, 0.116,
                  if(t < 1150, 0.297,
                  if(t < 1250, 0.760,
                  if(t < 1350, 1.550,
                  if(t < 1900, 3.370, 3.370))))))))))))))'
  []
[]

[Kernels]
  # Gas flow kernels
  # Equation for enclosure upstream
  [MatReaction_upstream_D2_outflux_membrane]
    type = ADMatBodyForce
    variable = D2_pressure_upstream
    material_property = 'membrane_reaction_rate_right'
    extra_vector_tags = 'ref'
  []
  # Equation for enclosure downstream
  [MatReaction_downstream_D2_influx_membrane]
    type = ADMatBodyForce
    variable = D2_pressure_downstream
    material_property = 'membrane_reaction_rate_left'
    extra_vector_tags = 'ref'
  []
[]

[Materials]
  [membrane_reaction_rate_right]
    type = ADParsedMaterial
    property_name = 'membrane_reaction_rate_right'
    postprocessor_names = flux_surface_right_D
    expression = 'flux_surface_right_D * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor} / 2'
    outputs = 'exodus'
  []
  [membrane_reaction_rate_left]
    type = ADParsedMaterial
    property_name = 'membrane_reaction_rate_left'
    postprocessor_names = flux_surface_left_D
    expression = 'flux_surface_left_D * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor} / 2'
    outputs = 'exodus'
  []
  [converter_to_regular]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_D'
    reg_props_out = 'diffusivity_D_nonAD'
    outputs = none
  []
[]

[BCs]
  # The surface of the slab in contact with the source is assumed to be in equilibrium with the source enclosure
  [right_diffusion]
    type = EquilibriumBC
    variable = D_concentration
    enclosure_var = D2_pressure_upstream
    boundary = 'right'
    Ko = '${solubility_pre}'
    activation_energy = '${solubility_energy}'
    p = '${solubility_exponent}'
    temperature = ${temperature}
  []
  [left_diffusion]
    type = EquilibriumBC
    variable = D_concentration
    enclosure_var = D2_pressure_downstream
    boundary = 'left'
    Ko = '${solubility_pre}'
    activation_energy = '${solubility_energy}'
    p = '${solubility_exponent}'
    temperature = ${temperature}
  []
[]

(test/tests/val-2b/val-2b.i)


# Physical constants
R = ${units 8.31446261815324 J/mol/K} # ideal gas constant based on number used in include/utils/PhysicalConstants.h

# Pressure conditions
pressure_enclosure_init = ${units 13300 Pa}
pressure_enclosure_cooldown = ${units 1e-6 Pa} # vaccum
pressure_enclosure_desorption = ${units 1e-3 Pa} # vaccum, assumed

# Temperature conditions
temperature_initial = ${units 773 K}
temperature_desorption_min = ${units 300 K}
temperature_desorption_max = ${units 1073 K}
temperature_cooldown_min = ${temperature_desorption_min}
desorption_heating_rate = ${units ${fparse 3/60} K/s}

# Important times
charge_time = ${units 50 h -> s}
cooldown_time_constant = ${units ${fparse 45*60} s}
# TMAP4 and TMAP7 used 40 minutes for the cooldown duration,
# We use a 5 hour cooldown period to let the temperature decrease to around 300 K for the start of the desorption.
# R.G. Macaulay-Newcombe et al. (1991) is not very clear on how long samples cooled down.
cooldown_duration = ${units 5 h -> s}
desorption_duration = ${fparse (temperature_desorption_max-temperature_desorption_min)/desorption_heating_rate}
endtime = ${fparse charge_time + cooldown_duration + desorption_duration}

# Materials properties
concentration_scaling = 1e10 # (-)
diffusion_Be_preexponential = ${units 8.0e-9 m^2/s -> mum^2/s}
diffusion_Be_energy = ${units 4220 K}
diffusion_BeO_preexponential_charging = ${units 1.40e-4 m^2/s -> mum^2/s}
diffusion_BeO_energy_charging = ${units 24408 K}
diffusion_BeO_preexponential_desorption = ${units 7e-5 m^2/s -> mum^2/s}
diffusion_BeO_energy_desorption = ${units 27000 K}
solubility_order = .5 # order of the solubility law (Here, we use Sievert's law)
solubility_constant_BeO = ${fparse 5.00e20 / 1e18 / concentration_scaling} # at/m^3/Pa^0.5 -> at/mum^3/Pa^0.5}
solubility_energy_BeO = ${units -9377.7 K}
solubility_constant_Be = ${fparse 7.156e27 / 1e18 / concentration_scaling} # at/m^3/Pa^0.5 -> at/mum^3/Pa^0.5}
solubility_energy_Be = ${units 11606 K}
jump_penalty = 1e0 # (-)

# Numerical parameters
dt_max_large = ${units 100 s}
dt_max_small = ${units 10 s}
dt_start_charging = ${units 1 s}
dt_start_cooldown = ${units 10 s}
dt_start_desorption = ${units 1 s}

# Geometry and mesh
length_BeO = ${units 18 nm -> mum}
num_nodes_BeO = 18
node_length_BeO = ${fparse length_BeO / num_nodes_BeO}
length_Be = ${units 0.4 mm -> mum}
length_Be_modeled = ${fparse length_Be/2}
num_nodes_Be = 40
node_length_Be = ${fparse length_Be_modeled / num_nodes_Be}

[Mesh]
  [cmg]
    type = CartesianMeshGenerator
    dim = 1

    #     0                  1                  2                  3                  4                  5                  6                  7                  8                  9                  10                 11                 12                 13                 14                 15                 16                 17
    dx = '${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO} ${node_length_BeO}
          ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be}
          ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be} ${node_length_Be}'
    #     18                19                20                21                22                23                24                25                26                27                28                29                30                31                32                33                34                35                36                37

    #               0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
    subdomain_id = '0 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0
                    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
                    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1'
    #               18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37'

  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = cmg
    primary_block = '0' #BeO
    paired_block = '1' # Be
    new_boundary = 'interface'
  []
  [interface_other_side]
    type = SideSetsBetweenSubdomainsGenerator
    input = interface
    primary_block = '1' #BeO
    paired_block = '0' # Be
    new_boundary = 'interface_other'
  []
[]

[Variables]
  [deuterium_concentration_Be] # (atoms/microns^3) / concentration_scaling
    block = 1
  []
  [deuterium_concentration_BeO] # (atoms/microns^3) / concentration_scaling
    block = 0
  []
[]

[AuxVariables]
  [enclosure_pressure]
    family = SCALAR
    initial_condition = ${pressure_enclosure_init}
  []
  [temperature]
    initial_condition = ${temperature_initial}
  []
  [flux_x]
    order = FIRST
    family = MONOMIAL
  []
[]

[Kernels]
  [time_BeO]
    type = TimeDerivative
    variable = deuterium_concentration_BeO
    block = 0
  []
  [diffusion_BeO]
    type = ADMatDiffusion
    variable = deuterium_concentration_BeO
    diffusivity = diffusivity_BeO
    block = 0
  []
  [time_Be]
    type = TimeDerivative
    variable = deuterium_concentration_Be
    block = 1
  []
  [diffusion_Be]
    type = ADMatDiffusion
    variable = deuterium_concentration_Be
    diffusivity = diffusivity_Be
    block = 1
  []
[]

[InterfaceKernels]
  [tied]
    type = ADPenaltyInterfaceDiffusion
    variable = deuterium_concentration_BeO
    neighbor_var = deuterium_concentration_Be
    penalty = ${jump_penalty}
    jump_prop_name = solubility_ratio
    boundary = 'interface'
  []
[]

[AuxScalarKernels]
  [enclosure_pressure_aux]
    type = FunctionScalarAux
    variable = enclosure_pressure
    function = enclosure_pressure_func
  []
[]

[AuxKernels]
  [temperature_aux]
    type = FunctionAux
    variable = temperature
    function = temperature_bc_func
    execute_on = 'INITIAL LINEAR'
  []
  [flux_x_Be]
    type = DiffusionFluxAux
    diffusivity = diffusivity_Be
    variable = flux_x
    diffusion_variable = deuterium_concentration_Be
    component = x
    block = 1
  []
  [flux_x_BeO]
    type = DiffusionFluxAux
    diffusivity = diffusivity_BeO
    variable = flux_x
    diffusion_variable = deuterium_concentration_BeO
    component = x
    block = 0
  []
[]

[BCs]
  [left_flux]
    type = EquilibriumBC
    Ko = ${solubility_constant_BeO}
    activation_energy = '${fparse solubility_energy_BeO * R}'
    boundary = left
    enclosure_var = enclosure_pressure
    temperature = temperature
    variable = deuterium_concentration_BeO
    p = ${solubility_order}
  []
  [right_flux]
    type = ADNeumannBC
    boundary = right
    variable = deuterium_concentration_Be
    value = 0
  []
[]

[Functions]
  [temperature_bc_func]
    type = ParsedFunction
    expression = 'if(t<${charge_time}, ${temperature_initial}, if(t<${fparse charge_time + cooldown_duration}, ${temperature_initial}-((1-exp(-(t-${charge_time})/${cooldown_time_constant}))*${fparse temperature_initial - temperature_cooldown_min}), ${temperature_desorption_min}+${desorption_heating_rate}*(t-${fparse charge_time + cooldown_duration})))'
  []
  [diffusivity_BeO_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = 'if(t<${fparse charge_time + cooldown_duration}, ${diffusion_BeO_preexponential_charging}*exp(-${diffusion_BeO_energy_charging}/T), ${diffusion_BeO_preexponential_desorption}*exp(-${diffusion_BeO_energy_desorption}/T))'
  []
  [diffusivity_Be_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${diffusion_Be_preexponential}*exp(-${diffusion_Be_energy}/T)'
  []
  [enclosure_pressure_func]
    type = ParsedFunction
    expression = 'if(t<${charge_time}, ${pressure_enclosure_init}, if(t<${fparse charge_time + cooldown_duration}, ${pressure_enclosure_cooldown}, ${pressure_enclosure_desorption}))'
  []
  [solubility_BeO_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${solubility_constant_BeO} * exp(-${solubility_energy_BeO}/T)'
  []
  [solubility_Be_func]
    type = ParsedFunction
    symbol_names = 'T'
    symbol_values = 'temperature_bc_func'
    expression = '${solubility_constant_Be} * exp(-${solubility_energy_Be}/T)'
  []
  [max_time_step_size_func]
    type = ParsedFunction
    expression = 'if(t < ${fparse charge_time-dt_max_large}, ${dt_max_large}, ${dt_max_small})'
  []
[]

[Materials]
  [diffusion_solubility]
    type = ADGenericFunctionMaterial
    prop_names = 'diffusivity_BeO diffusivity_Be solubility_Be solubility_BeO'
    prop_values = 'diffusivity_BeO_func diffusivity_Be_func solubility_Be_func solubility_BeO_func'
    outputs = all
  []
  [converter_to_nonAD]
    type = MaterialADConverter
    ad_props_in = 'diffusivity_Be diffusivity_BeO'
    reg_props_out = 'diffusivity_Be_nonAD diffusivity_BeO_nonAD'
    outputs = all
  []
  [interface_jump]
    type = SolubilityRatioMaterial
    solubility_primary = solubility_BeO
    solubility_secondary = solubility_Be
    boundary = interface
    concentration_primary = deuterium_concentration_BeO
    concentration_secondary = deuterium_concentration_Be
    outputs = all
  []
[]

[Postprocessors]
  [avg_flux_left]
    type = SideDiffusiveFluxAverage
    variable = deuterium_concentration_BeO
    boundary = left
    diffusivity = diffusivity_BeO_nonAD
  []
  [avg_flux_total] # total flux coming out of the sample in atoms/microns^2/s
    type = ScalePostprocessor
    value = avg_flux_left
    scaling_factor = '${fparse 2 * concentration_scaling}'
    # Factor of 2 because symmetry is assumed and only one-half of the specimen is modeled.
    # Thus, the total flux coming out of the specimen (per unit area)
    # is twice the flux calculated at the left side of the domain.
    # The 'concentration_scaling' parameter is used to get a consistent concentration unit
  []
  [temperature]
    type = ElementAverageValue
    block = 0
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [diffusion_Be]
    type = ElementAverageValue
    block = 1
    variable = diffusivity_Be
  []
  [diffusion_BeO]
    type = ElementAverageValue
    block = 0
    variable = diffusivity_BeO
  []
  [solubility_Be]
    type = ElementAverageValue
    block = 1
    variable = solubility_Be
  []
  [solubility_BeO]
    type = ElementAverageValue
    block = 0
    variable = solubility_BeO
  []
  [gold_solubility_ratio]
    type = ParsedPostprocessor
    pp_names = 'solubility_BeO solubility_Be'
    expression = 'solubility_BeO / solubility_Be'
  []
  [BeO_interface]
    type = SideAverageValue
    boundary = interface
    variable = deuterium_concentration_BeO
  []
  [Be_interface]
    type = SideAverageValue
    boundary = interface_other
    variable = deuterium_concentration_Be
  []
  [variable_ratio]
    type = ParsedPostprocessor
    pp_names = 'BeO_interface Be_interface'
    expression = 'BeO_interface / Be_interface'
  []
  [dt]
    type = TimestepSize
  []
  [max_time_step_size_pp]
    type = FunctionValuePostprocessor
    function = max_time_step_size_func
    execute_on = 'INITIAL TIMESTEP_END'
    outputs = none
  []
[]

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

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-12
  end_time = ${endtime}
  automatic_scaling = true
  compute_scaling_once = false
  nl_max_its = 7
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${dt_start_charging}
    optimal_iterations = 5
    growth_factor = 1.1
    cutback_factor_at_failure = .9
    timestep_limiting_postprocessor = max_time_step_size_pp
    time_t = ' 0                    ${charge_time}        ${fparse charge_time + cooldown_duration}'
    time_dt = '${dt_start_charging} ${dt_start_cooldown}  ${dt_start_desorption}'
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  csv = true
  [exodus]
    type = Exodus
    output_material_properties = true
  []
[]

(test/tests/bcs/EquilibriumBC/equilibriumBC.i)

# Simple test to ensure that EquilibriumBC's enclosure_var can take in a field variable

[Mesh]
  type = GeneratedMesh
  dim = 1
[]

[Variables]
  [u]
  []
[]

[AuxVariables]
  [v]
    initial_condition = 1
  []
[]

[Kernels]
  [diff]
    type = MatDiffusion
    variable = u
    diffusivity = 1
  []
  [time]
    type = TimeDerivative
    variable = u
  []
[]

[BCs]
  [right]
    type = DirichletBC
    value = 0
    variable = u
    boundary = 'right'
  []
  [left]
    type = EquilibriumBC
    variable = u
    enclosure_var = v
    boundary = 'left'
    Ko = 1
    temperature = 1
  []
[]

[Postprocessors]
  [flux_surface_left]
    type = SideDiffusiveFluxIntegral
    variable = u
    diffusivity = 1
    boundary = 'left'
    execute_on = 'initial nonlinear linear timestep_end'
  []
[]

[Executioner]
  type = Transient
  dt = .1
  num_steps = 1
  solve_type = PJFNK
  automatic_scaling = true
  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  line_search = 'bt'
  scheme = 'bdf2'
[]

[Outputs]
  file_base = equilibriumBC_out
  perf_graph = true
  [csv]
    type = CSV
    execute_on = 'initial timestep_end'
  []
[]

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

# Verification Problem #1a from TMAP4/TMAP7 V&V document
# Tritium diffusion through SiC layer with depleting source at 2100 C.
# No Soret effect, solubility, or trapping included.

# Physical Constants
# Note that we do NOT use the same number of digits as in TMAP4/TMAP7.
# This is to be consistent with PhysicalConstant.h
kb = '${units 1.380649e-23 J/K}' # Boltzmann constant
R = '${units 8.31446261815324 J/mol/K}' # Gas constant

# Data used in TMAP4/TMAP7 case
temperature = '${units 2373 K}'
initial_pressure = '${units 1e6 Pa}'
volume_enclosure = '${units 5.20e-11 m^3 -> mum^3}'
surface_area = '${units 2.16e-6 m^2 -> mum^2}'
diffusivity_SiC = '${units ${fparse 1.58e-4*exp(-308000.0/(R*temperature))} m^2/s -> mum^2/s}'
solubility_constant = '${units ${fparse 7.244e22 / temperature} at/m^3/Pa -> at/mum^3/Pa}'
slab_thickness = '${units 3.30e-5 m -> mum}'

# Useful equations/conversions
concentration_to_pressure_conversion_factor = '${units ${fparse kb*temperature} Pa*m^3 -> Pa*mum^3}'

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 150
  xmax = '${slab_thickness}'
[]

[Variables]
  # concentration in the SiC layer in atoms/microns^3
  [u]
  []
  # pressure of the enclosure in Pa
  [v]
    family = SCALAR
    order = FIRST
    initial_condition = '${fparse initial_pressure}'
  []
[]

[Kernels]
  [diff]
    type = MatDiffusion
    variable = u
    diffusivity = '${diffusivity_SiC}'
  []
  [time]
    type = TimeDerivative
    variable = u
  []
[]

[ScalarKernels]
  [time]
    type = ODETimeDerivative
    variable = v
  []
  [flux_sink]
    type = EnclosureSinkScalarKernel
    variable = v
    flux = scaled_flux_surface_left
    surface_area = '${surface_area}'
    volume = '${volume_enclosure}'
    concentration_to_pressure_conversion_factor = '${concentration_to_pressure_conversion_factor}'
  []
[]

[BCs]
  # The concentration on the outer boundary of the SiC layer is kept at 0
  [right]
    type = DirichletBC
    value = 0
    variable = u
    boundary = 'right'
  []
  # The surface of the slab in contact with the source is assumed to be in equilibrium with the source enclosure
  [left]
    type = EquilibriumBC
    variable = u
    enclosure_var = v
    boundary = 'left'
    Ko = '${solubility_constant}'
    temperature = ${temperature}
  []
[]

[Postprocessors]
  # flux of tritium through the outer SiC surface - compare to TMAP7
  [flux_surface_right]
    type = SideDiffusiveFluxIntegral
    variable = u
    diffusivity = '${diffusivity_SiC}'
    boundary = 'right'
    execute_on = 'initial nonlinear linear timestep_end'
    outputs = 'console csv exodus'
  []
  # flux of tritium through the surface of SiC layer in contact with enclosure
  [flux_surface_left]
    type = SideDiffusiveFluxIntegral
    variable = u
    diffusivity = '${diffusivity_SiC}'
    boundary = 'left'
    execute_on = 'initial nonlinear linear timestep_end'
    outputs = ''
  []
  # scale the flux to get inward direction
  [scaled_flux_surface_left]
    type = ScalePostprocessor
    scaling_factor = -1
    value = flux_surface_left
    execute_on = 'initial nonlinear linear timestep_end'
    outputs = 'console csv exodus'
  []
  # integral of the tritium flux through outer surface of SiC layer
  [integral_release_flux_right]
    type = PressureReleaseFluxIntegral
    variable = u
    boundary = 'right'
    diffusivity = '${diffusivity_SiC}'
    surface_area = '${surface_area}'
    volume = '${volume_enclosure}'
    concentration_to_pressure_conversion_factor = '${concentration_to_pressure_conversion_factor}'
    outputs = 'console'
  []
  # cumulative sum of integral_release_flux_right over time (i.e. cumulative amount released)
  [cumulative_release_right]
    type = CumulativeValuePostprocessor
    postprocessor = 'integral_release_flux_right'
    outputs = 'console'
  []
  # released fraction based on outer layer flux - compare to TMAP4
  [released_fraction_right]
    type = ScalePostprocessor
    value = 'cumulative_release_right'
    scaling_factor = '${fparse 1./(initial_pressure)}'
    outputs = 'console csv exodus'
  []
  # Make a postprocessor take the value of the scalar value v
  [v_value]
    type = ScalarVariable
    variable = v
  []
  # released fraction based on inner layer flux on v - compare to TMAP7
  [released_fraction_left]
    type = LinearCombinationPostprocessor
    pp_names = 'v_value'
    pp_coefs = '${fparse -1./(initial_pressure)}'
    b = 1
  []
[]

[Executioner]
  type = Transient
  dt = .1
  end_time = 140
  solve_type = PJFNK
  automatic_scaling = true
  dtmin = .1
  l_max_its = 30
  nl_max_its = 5
  petsc_options = '-snes_converged_reason -ksp_monitor_true_residual'
  petsc_options_iname = '-pc_type -mat_mffd_err'
  petsc_options_value = 'lu       1e-5'
  line_search = 'bt'
  scheme = 'bdf2'
  timestep_tolerance = 1e-8
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
  perf_graph = true
  [dof]
    type = DOFMap
    execute_on = 'initial'
  []
  [csv]
    type = CSV
    execute_on = 'initial timestep_end'
  []
[]

(test/tests/val-2e/val-2ed.i)

# Validation Problem #2ed 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

# Solubility data used in TMAP7 case
solubility_exponent = 0.9297 # -
solubility_pre = '${units ${fparse 9.355e22 / 1e18} at/mum^3/Pa^0.9297}'
solubility_energy = '${units ${fparse 5918 * R} J/mol}'

# Modeling data used in current case
file_name = 'val-2ed_out'

!include val-2e_base_three_gases.i

[Kernels]
  # Gas flow kernels
  # Equation for D2 in enclosure upstream
  [MatReaction_upstream_D2_reaction]
    type = ADMatReaction
    variable = D2_pressure_upstream
    v = 'sqrt_PH2_sqrt_PD2_upstream'
    reaction_rate = -1
    extra_vector_tags = 'ref'
  []
  [MatReaction_upstream_D2_re_reaction]
    type = ADMatReaction
    variable = D2_pressure_upstream
    v = 'HD_pressure_upstream'
    reaction_rate = 0.5
    extra_vector_tags = 'ref'
  []
  # Equation for H2 in enclosure upstream
  [MatReaction_upstream_H2_reaction]
    type = ADMatReaction
    variable = H2_pressure_upstream
    v = 'sqrt_PH2_sqrt_PD2_upstream'
    reaction_rate = -1
    extra_vector_tags = 'ref'
  []
  [MatReaction_upstream_H2_re_reaction]
    type = ADMatReaction
    variable = H2_pressure_upstream
    v = 'HD_pressure_upstream'
    reaction_rate = 0.5
    extra_vector_tags = 'ref'
  []
  # Equation for HD in enclosure upstream
  [MatReaction_upstream_HD_reaction]
    type = ADMatReaction
    variable = HD_pressure_upstream
    v = 'sqrt_PH2_sqrt_PD2_upstream'
    reaction_rate = 2
    extra_vector_tags = 'ref'
  []
  [MatReaction_upstream_HD_re_reaction]
    type = ADMatReaction
    variable = HD_pressure_upstream
    v = 'HD_pressure_upstream'
    reaction_rate = -1
    extra_vector_tags = 'ref'
  []
  # Membrane upstream
  [MatReaction_upstream_outflux_membrane_D]
    type = ADMatBodyForce
    variable = D2_pressure_upstream
    material_property = 'membrane_reaction_rate_right_D'
    extra_vector_tags = 'ref'
  []
  [MatReaction_upstream_outflux_membrane_H]
    type = ADMatBodyForce
    variable = H2_pressure_upstream
    material_property = 'membrane_reaction_rate_right_H'
    extra_vector_tags = 'ref'
  []

  # Equation for D2 enclosure downstream
  [MatReaction_downstream_D2_reaction]
    type = ADMatReaction
    variable = D2_pressure_downstream
    v = 'sqrt_PH2_sqrt_PD2_downstream'
    reaction_rate = -1
    extra_vector_tags = 'ref'
  []
  [MatReaction_downstream_D2_re_reaction]
    type = ADMatReaction
    variable = D2_pressure_downstream
    v = 'HD_pressure_downstream'
    reaction_rate = 0.5
    extra_vector_tags = 'ref'
  []
  # Equation for H2 enclosure downstream
  [MatReaction_downstream_H2_reaction]
    type = ADMatReaction
    variable = H2_pressure_downstream
    v = 'sqrt_PH2_sqrt_PD2_downstream'
    reaction_rate = -1
    extra_vector_tags = 'ref'
  []
  [MatReaction_downstream_H2_re_reaction]
    type = ADMatReaction
    variable = H2_pressure_downstream
    v = 'HD_pressure_downstream'
    reaction_rate = 0.5
    extra_vector_tags = 'ref'
  []
  # Equation for HD enclosure downstream
  [MatReaction_downstream_HD_reaction]
    type = ADMatReaction
    variable = HD_pressure_downstream
    v = 'sqrt_PH2_sqrt_PD2_downstream'
    reaction_rate = 2
    extra_vector_tags = 'ref'
  []
  [MatReaction_downstream_HD_re_reaction]
    type = ADMatReaction
    variable = HD_pressure_downstream
    v = 'HD_pressure_downstream'
    reaction_rate = -1
    extra_vector_tags = 'ref'
  []
  # Membrane downstream
  [MatReaction_downstream_influx_membrane_D]
    type = ADMatBodyForce
    variable = D2_pressure_downstream
    material_property = 'membrane_reaction_rate_left_D'
    extra_vector_tags = 'ref'
  []
  [MatReaction_downstream_influx_membrane_H]
    type = ADMatBodyForce
    variable = H2_pressure_downstream
    material_property = 'membrane_reaction_rate_left_H'
    extra_vector_tags = 'ref'
  []
[]

[Materials]
  [membrane_reaction_rate_right_D]
    type = ADParsedMaterial
    property_name = 'membrane_reaction_rate_right_D'
    postprocessor_names = flux_surface_right_D
    expression = 'flux_surface_right_D * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}/2'
  []
  [membrane_reaction_rate_left_D]
    type = ADParsedMaterial
    property_name = 'membrane_reaction_rate_left_D'
    postprocessor_names = flux_surface_left_D
    expression = 'flux_surface_left_D * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}/2'
  []
  [membrane_reaction_rate_right_H]
    type = ADParsedMaterial
    property_name = 'membrane_reaction_rate_right_H'
    postprocessor_names = flux_surface_right_H
    expression = 'flux_surface_right_H * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}/2'
  []
  [membrane_reaction_rate_left_H]
    type = ADParsedMaterial
    property_name = 'membrane_reaction_rate_left_H'
    postprocessor_names = flux_surface_left_H
    expression = 'flux_surface_left_H * ${surface_area} / ${volume_enclosure} * ${concentration_to_pressure_conversion_factor}/2'
  []
[]

[BCs]
  # The surface of the slab in contact with the source is assumed to be in equilibrium with the source enclosure
  [right_diffusion_D]
    type = EquilibriumBC
    variable = D_concentration
    enclosure_var = D2_partial_pressure_upstream
    boundary = 'right'
    Ko = '${solubility_pre}'
    activation_energy = '${solubility_energy}'
    p = '${solubility_exponent}'
    temperature = ${temperature}
  []
  [left_diffusion_D]
    type = EquilibriumBC
    variable = D_concentration
    enclosure_var = D2_partial_pressure_downstream
    boundary = 'left'
    Ko = '${solubility_pre}'
    activation_energy = '${solubility_energy}'
    p = '${solubility_exponent}'
    temperature = ${temperature}
  []
  [right_diffusion_H]
    type = EquilibriumBC
    variable = H_concentration
    enclosure_var = H2_partial_pressure_upstream
    boundary = 'right'
    Ko = '${solubility_pre}'
    activation_energy = '${solubility_energy}'
    p = '${solubility_exponent}'
    temperature = ${temperature}
  []
  [left_diffusion_H]
    type = EquilibriumBC
    variable = H_concentration
    enclosure_var = H2_partial_pressure_downstream
    boundary = 'left'
    Ko = '${solubility_pre}'
    activation_energy = '${solubility_energy}'
    p = '${solubility_exponent}'
    temperature = ${temperature}
  []
[]

Overview
Input Parameters
Input Files