SolubilityRatioMaterial

Calculates the jump in concentration across an interface.

Overview

The SolubilityRatioMaterial object is used to calculate the jump in species concentration across a material interface, given different solubility values for those materials and species concentration on either side of the interface. Solubilities for the "primary" and "secondary" sides can be provided via the Materials System using automatic differentiation (a.k.a. an ADMaterialProperty).

The solubility ratio jump is calculated using the following relationship:

var element = document.getElementById("moose-equation-aa8236bf-fd47-4ff9-bcc6-2beff84740af");katex.render("J = \\frac{c_p}{S_p} - \\frac{c_s}{S_s}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-133ce457-eb1f-415c-bb20-f7bc8afbce63");katex.render("J", element, {displayMode:false,throwOnError:false});$ is the calculated solubility ratio jump (available as an InterfaceMaterial property, named the solubility_ratio), $var element = document.getElementById("moose-equation-de09db2b-50c9-4ad0-bef9-c3e34cd8ae00");katex.render("c_i", element, {displayMode:false,throwOnError:false});$ is the species concentration on the $var element = document.getElementById("moose-equation-d5117411-a27e-4a23-b19a-0c06b79a45b5");katex.render("i", element, {displayMode:false,throwOnError:false});$ th side of the interface, and $var element = document.getElementById("moose-equation-1a726d68-39a7-4dc9-b5fe-2e1ef9bc9213");katex.render("S_i", element, {displayMode:false,throwOnError:false});$ is the solubility on the $var element = document.getElementById("moose-equation-179d431a-8dfe-4e2e-89d4-6d45c60692d4");katex.render("i", element, {displayMode:false,throwOnError:false});$ th side of the interface. The $var element = document.getElementById("moose-equation-0b9e3dd3-872f-4d71-bfa4-b854122401a5");katex.render("p", element, {displayMode:false,throwOnError:false});$ subscript corresponds to the primary side and $var element = document.getElementById("moose-equation-a575f8e4-581b-4287-a2e2-49af45efbf1f");katex.render("s", element, {displayMode:false,throwOnError:false});$ corresponds to the secondary side. The ratio jump material property can then be used in InterfaceKernels, such as ADPenaltyInterfaceDiffusion, to solve for the appropriate concentrations on either side of the interface.

Example Input File Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [interface_jump]
    type = SolubilityRatioMaterial<<<{"description": "Calculates the jump in concentration across an interface.", "href": "SolubilityRatioMaterial.html"}>>>
    solubility_primary<<<{"description": "The material property on the primary side of the interface"}>>> = solubility_BeO
    solubility_secondary<<<{"description": "The material property on the secondary side of the interface"}>>> = solubility_Be
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = interface
    concentration_primary<<<{"description": "Primary side non-linear variable for jump computation"}>>> = deuterium_concentration_BeO
    concentration_secondary<<<{"description": "Secondary side non-linear variable for jump computation"}>>> = deuterium_concentration_Be
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = all
  []
[]

Input Parameters

concentration_primaryPrimary side non-linear variable for jump computation
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Primary side non-linear variable for jump computation
concentration_secondarySecondary side non-linear variable for jump computation
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Secondary side non-linear variable for jump computation
solubility_primaryThe material property on the primary side of the interface
C++ Type:std::string
Controllable:No
Description:The material property on the primary side of the interface
solubility_secondaryThe material property on the secondary side of the interface
C++ Type:std::string
Controllable:No
Description:The material property on the secondary side of the interface

Required 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
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
declare_suffixAn optional suffix parameter that can be appended to any declared 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 declared properties. The suffix will be prepended with a '_' character.

Optional 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

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs 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

(tmap8/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
  []
[]

Overview
Example Input File Syntax
Input Parameters