- 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
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:
where is the calculated solubility ratio jump (available as an InterfaceMaterial property, named the solubility_ratio), is the species concentration on the th side of the interface, and is the solubility on the th side of the interface. The subscript corresponds to the primary side and 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
- 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
Input Files
(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/divertor_monoblock/divertor_monoblock.i)
### Nomenclatures
### C_mobile_j mobile H concentration in "j" material, where j = CuCrZr, Cu, W
### C_trapped_j trapped H concentration in "j" material, where j = CuCrZr, Cu, W
### C_total_j total H concentration in "j" material, where j = CuCrZr, Cu, W
###
### S_empty_j empty site concentration in "j" material, where j = CuCrZr, Cu, W
### S_trapped_j trapped site concentration in "j" material, where j = CuCrZr, Cu, W
### S_total_j total site H concentration in "j" material, where j = CuCrZr, Cu, W
###
### F_permeation permeation flux
### F_recombination recombination flux
###
### Sc_ Scaled
### Int_ Integrated
### ScInt_ Scaled and integrated
tungsten_atomic_density = ${units 6.338e28 m^-3}
## Mesh generated by ConcentricCircleMeshGenerator and saved as "2DMonoblock.e"
[Mesh]
[ccmg]
type = ConcentricCircleMeshGenerator
num_sectors = 36
rings = '1 30 20 110'
radii = '${units 6 mm -> m} ${units 7.5 mm -> m} ${units 8.5 mm -> m}'
has_outer_square = on
pitch = ${units 28 mm -> m}
portion = left_half
preserve_volumes = false
smoothing_max_it = 3
[]
[ssbsg1]
type = SideSetsBetweenSubdomainsGenerator
input = ccmg
primary_block = '4' # W
paired_block = '3' # Cu
new_boundary = '4to3'
[]
[ssbsg2]
type = SideSetsBetweenSubdomainsGenerator
input = ssbsg1
primary_block = '3' # Cu
paired_block = '4' # W
new_boundary = '3to4'
[]
[ssbsg3]
type = SideSetsBetweenSubdomainsGenerator
input = ssbsg2
primary_block = '3' # Cu
paired_block = '2' # CuCrZr
new_boundary = '3to2'
[]
[ssbsg4]
type = SideSetsBetweenSubdomainsGenerator
input = ssbsg3
primary_block = '2' # CuCrZr
paired_block = '3' # Cu
new_boundary = '2to3'
[]
[ssbsg5]
type = SideSetsBetweenSubdomainsGenerator
input = ssbsg4
primary_block = '2' # CuCrZr
paired_block = '1' # H2O
new_boundary = '2to1'
[]
[bdg]
type = BlockDeletionGenerator
input = ssbsg5
block = '1' # H2O
[]
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[temperature]
order = FIRST
family = LAGRANGE
initial_condition = ${units 300 K}
[]
######################### Variables for W (block = 4)
[C_mobile_W]
order = FIRST
family = LAGRANGE
initial_condition = ${units 1.0e-20 m^-3}
block = 4
[]
[C_trapped_W]
order = FIRST
family = LAGRANGE
initial_condition = ${units 1.0e-15 m^-3}
block = 4
[]
######################### Variables for Cu (block = 3)
[C_mobile_Cu]
order = FIRST
family = LAGRANGE
initial_condition = ${units 5.0e-17 m^-3}
block = 3
[]
[C_trapped_Cu]
order = FIRST
family = LAGRANGE
initial_condition = ${units 1.0e-15 m^-3}
block = 3
[]
######################### Variables for CuCrZr (block = 2)
[C_mobile_CuCrZr]
order = FIRST
family = LAGRANGE
initial_condition = ${units 1.0e-15 m^-3}
block = 2
[]
[C_trapped_CuCrZr]
order = FIRST
family = LAGRANGE
initial_condition = ${units 1.0e-15 m^-3}
block = 2
[]
[]
[AuxVariables]
[flux_y]
order = FIRST
family = MONOMIAL
[]
############################## AuxVariables for W (block = 4)
[Sc_C_mobile_W]
block = 4
[]
[Sc_C_trapped_W]
block = 4
[]
[C_total_W]
block = 4
[]
[Sc_C_total_W]
block = 4
[]
[S_empty_W]
block = 4
[]
[Sc_S_empty_W]
block = 4
[]
[S_trapped_W]
block = 4
[]
[Sc_S_trapped_W]
block = 4
[]
[S_total_W]
block = 4
[]
[Sc_S_total_W]
block = 4
[]
############################## AuxVariables for Cu (block = 3)
[Sc_C_mobile_Cu]
block = 3
[]
[Sc_C_trapped_Cu]
block = 3
[]
[C_total_Cu]
block = 3
[]
[Sc_C_total_Cu]
block = 3
[]
[S_empty_Cu]
block = 3
[]
[Sc_S_empty_Cu]
block = 3
[]
[S_trapped_Cu]
block = 3
[]
[Sc_S_trapped_Cu]
block = 3
[]
[S_total_Cu]
block = 3
[]
[Sc_S_total_Cu]
block = 3
[]
############################## AuxVariables for CuCrZr (block = 2)
[Sc_C_mobile_CuCrZr]
block = 2
[]
[Sc_C_trapped_CuCrZr]
block = 2
[]
[C_total_CuCrZr]
block = 2
[]
[Sc_C_total_CuCrZr]
block = 2
[]
[S_empty_CuCrZr]
block = 2
[]
[Sc_S_empty_CuCrZr]
block = 2
[]
[S_trapped_CuCrZr]
block = 2
[]
[Sc_S_trapped_CuCrZr]
block = 2
[]
[S_total_CuCrZr]
block = 2
[]
[Sc_S_total_CuCrZr]
block = 2
[]
[]
[Kernels]
############################## Kernels for W (block = 4)
[diff_W]
type = ADMatDiffusion
variable = C_mobile_W
diffusivity = diffusivity_W
block = 4
extra_vector_tags = ref
[]
[time_diff_W]
type = ADTimeDerivative
variable = C_mobile_W
block = 4
extra_vector_tags = ref
[]
[coupled_time_W]
type = ScaledCoupledTimeDerivative
variable = C_mobile_W
v = C_trapped_W
factor = 1e0
block = 4
extra_vector_tags = ref
[]
[heat_conduction_W]
type = HeatConduction
variable = temperature
diffusion_coefficient = thermal_conductivity_W
block = 4
extra_vector_tags = ref
[]
[time_heat_conduction_W]
type = SpecificHeatConductionTimeDerivative
variable = temperature
specific_heat = specific_heat_W
density = density_W
block = 4
extra_vector_tags = ref
[]
############################## Kernels for Cu (block = 3)
[diff_Cu]
type = ADMatDiffusion
variable = C_mobile_Cu
diffusivity = diffusivity_Cu
block = 3
extra_vector_tags = ref
[]
[time_diff_Cu]
type = ADTimeDerivative
variable = C_mobile_Cu
block = 3
extra_vector_tags = ref
[]
[coupled_time_Cu]
type = ScaledCoupledTimeDerivative
variable = C_mobile_Cu
v = C_trapped_Cu
factor = 1e0
block = 3
extra_vector_tags = ref
[]
[heat_conduction_Cu]
type = HeatConduction
variable = temperature
diffusion_coefficient = thermal_conductivity_Cu
block = 3
extra_vector_tags = ref
[]
[time_heat_conduction_Cu]
type = SpecificHeatConductionTimeDerivative
variable = temperature
specific_heat = specific_heat_Cu
density = density_Cu
block = 3
extra_vector_tags = ref
[]
############################## Kernels for CuCrZr (block = 2)
[diff_CuCrZr]
type = ADMatDiffusion
variable = C_mobile_CuCrZr
diffusivity = diffusivity_CuCrZr
block = 2
extra_vector_tags = ref
[]
[time_diff_CuCrZr]
type = ADTimeDerivative
variable = C_mobile_CuCrZr
block = 2
extra_vector_tags = ref
[]
[coupled_time_CuCrZr]
type = ScaledCoupledTimeDerivative
variable = C_mobile_CuCrZr
v = C_trapped_CuCrZr
factor = 1e0
block = 2
extra_vector_tags = ref
[]
[heat_conduction_CuCrZr]
type = HeatConduction
variable = temperature
diffusion_coefficient = thermal_conductivity_CuCrZr
block = 2
extra_vector_tags = ref
[]
[time_heat_conduction_CuCrZr]
type = SpecificHeatConductionTimeDerivative
variable = temperature
specific_heat = specific_heat_CuCrZr
density = density_CuCrZr
block = 2
extra_vector_tags = ref
[]
[]
[AuxKernels]
############################## AuxKernels for W (block = 4)
[Scaled_mobile_W]
variable = Sc_C_mobile_W
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_mobile_W
[]
[Scaled_trapped_W]
variable = Sc_C_trapped_W
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_trapped_W
[]
[total_W]
variable = C_total_W
type = ParsedAux
expression = 'C_mobile_W + C_trapped_W'
coupled_variables = 'C_mobile_W C_trapped_W'
[]
[Scaled_total_W]
variable = Sc_C_total_W
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_total_W
[]
[empty_sites_W]
variable = S_empty_W
type = EmptySitesAux
N = ${units 1.0e0 m^-3} # = ${tungsten_atomic_density} #/m^3 (W lattice density)
# Ct0 = ${units 1.0e-4 m^-3} # E.A. Hodille et al 2021 Nucl. Fusion 61 126003, trap 1
Ct0 = ${units 1.0e-4 m^-3} # E.A. Hodille et al 2021 Nucl. Fusion 61 1260033, trap 2
trap_per_free = 1.0e0 # 1.0e1
trapped_concentration_variables = C_trapped_W
[]
[scaled_empty_W]
variable = Sc_S_empty_W
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_empty_W
[]
[trapped_sites_W]
variable = S_trapped_W
type = NormalizationAux
normal_factor = 1e0
source_variable = C_trapped_W
[]
[scaled_trapped_W]
variable = Sc_S_trapped_W
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_trapped_W
[]
[total_sites_W]
variable = S_total_W
type = ParsedAux
expression = 'S_trapped_W + S_empty_W'
coupled_variables = 'S_trapped_W S_empty_W'
[]
[scaled_total_W]
variable = Sc_S_total_W
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_total_W
[]
############################## AuxKernels for Cu (block = 3)
[Scaled_mobile_Cu]
variable = Sc_C_mobile_Cu
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_mobile_Cu
[]
[Scaled_trapped_Cu]
variable = Sc_C_trapped_Cu
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_trapped_Cu
[]
[total_Cu]
variable = C_total_Cu
type = ParsedAux
expression = 'C_mobile_Cu + C_trapped_Cu'
coupled_variables = 'C_mobile_Cu C_trapped_Cu'
[]
[Scaled_total_Cu]
variable = Sc_C_total_Cu
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_total_Cu
[]
[empty_sites_Cu]
variable = S_empty_Cu
type = EmptySitesAux
N = ${units 1.0e0 m^-3} # = ${tungsten_atomic_density} #/m^3 (W lattice density)
Ct0 = ${units 5.0e-5 m^-3} # R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 3
trap_per_free = 1.0e0 # 1.0e1
trapped_concentration_variables = C_trapped_Cu
[]
[scaled_empty_Cu]
variable = Sc_S_empty_Cu
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_empty_Cu
[]
[trapped_sites_Cu]
variable = S_trapped_Cu
type = NormalizationAux
normal_factor = 1e0
source_variable = C_trapped_Cu
[]
[scaled_trapped_Cu]
variable = Sc_S_trapped_Cu
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_trapped_Cu
[]
[total_sites_Cu]
variable = S_total_Cu
type = ParsedAux
expression = 'S_trapped_Cu + S_empty_Cu'
coupled_variables = 'S_trapped_Cu S_empty_Cu'
[]
[scaled_total_Cu]
variable = Sc_S_total_Cu
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_total_Cu
[]
############################## AuxKernels for CuCrZr (block = 2)
[Scaled_mobile_CuCrZr]
variable = Sc_C_mobile_CuCrZr
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_mobile_CuCrZr
[]
[Scaled_trapped_CuCrZr]
variable = Sc_C_trapped_CuCrZr
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_trapped_CuCrZr
[]
[total_CuCrZr]
variable = C_total_CuCrZr
type = ParsedAux
expression = 'C_mobile_CuCrZr + C_trapped_CuCrZr'
coupled_variables = 'C_mobile_CuCrZr C_trapped_CuCrZr'
[]
[Scaled_total_CuCrZr]
variable = Sc_C_total_CuCrZr
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = C_total_CuCrZr
[]
[empty_sites_CuCrZr]
variable = S_empty_CuCrZr
type = EmptySitesAux
N = ${units 1.0e0 m^-3} # = ${tungsten_atomic_density} #/m^3 (W lattice density)
Ct0 = ${units 5.0e-5 m^-3} # R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 4
# Ct0 = ${units 4.0e-2 m^-3} # R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 5
trap_per_free = 1.0e0 # 1.0e1
trapped_concentration_variables = C_trapped_CuCrZr
[]
[scaled_empty_CuCrZr]
variable = Sc_S_empty_CuCrZr
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_empty_CuCrZr
[]
[trapped_sites_CuCrZr]
variable = S_trapped_CuCrZr
type = NormalizationAux
normal_factor = 1e0
source_variable = C_trapped_CuCrZr
[]
[scaled_trapped_CuCrZr]
variable = Sc_S_trapped_CuCrZr
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_trapped_CuCrZr
[]
[total_sites_CuCrZr]
variable = S_total_CuCrZr
type = ParsedAux
expression = 'S_trapped_CuCrZr + S_empty_CuCrZr'
coupled_variables = 'S_trapped_CuCrZr S_empty_CuCrZr'
[]
[scaled_total_CuCrZr]
variable = Sc_S_total_CuCrZr
type = NormalizationAux
normal_factor = ${tungsten_atomic_density}
source_variable = S_total_CuCrZr
[]
[flux_y_W]
type = DiffusionFluxAux
diffusivity = diffusivity_W
variable = flux_y
diffusion_variable = C_mobile_W
component = y
block = 4
[]
[flux_y_Cu]
type = DiffusionFluxAux
diffusivity = diffusivity_Cu
variable = flux_y
diffusion_variable = C_mobile_Cu
component = y
block = 3
[]
[flux_y_CuCrZr]
type = DiffusionFluxAux
diffusivity = diffusivity_CuCrZr
variable = flux_y
diffusion_variable = C_mobile_CuCrZr
component = y
block = 2
[]
[]
[InterfaceKernels]
[tied_4to3]
type = ADPenaltyInterfaceDiffusion
variable = C_mobile_W
neighbor_var = C_mobile_Cu
penalty = 0.05 # it will not converge with > 0.1, but it creates negative C_mobile _Cu with << 0.1
# jump_prop_name = solubility_ratio_4to3
jump_prop_name = solubility_ratio
boundary = '4to3'
[]
[tied_3to2]
type = ADPenaltyInterfaceDiffusion
variable = C_mobile_Cu
neighbor_var = C_mobile_CuCrZr
penalty = 0.05 # it will not converge with > 0.1, but it creates negative C_mobile _Cu with << 0.1
# jump_prop_name = solubility_ratio_3to2
jump_prop_name = solubility_ratio
boundary = '3to2'
[]
[]
[NodalKernels]
############################## NodalKernels for W (block = 4)
[time_W]
type = TimeDerivativeNodalKernel
variable = C_trapped_W
[]
[trapping_W]
type = TrappingNodalKernel
variable = C_trapped_W
temperature = temperature
alpha_t = 2.75e11 # 1e15
N = 1.0e0 # = (1e0) x (${tungsten_atomic_density} #/m^3)
# Ct0 = 1.0e-4 # E.A. Hodille et al 2021 Nucl. Fusion 61 126003, trap 1
Ct0 = 1.0e-4 # E.A. Hodille et al 2021 Nucl. Fusion 61 1260033, trap 2
trap_per_free = 1.0e0 # 1.0e1
mobile_concentration = 'C_mobile_W'
extra_vector_tags = ref
[]
[release_W]
type = ReleasingNodalKernel
alpha_r = 8.4e12 # 1.0e13
temperature = temperature
# detrapping_energy = 9863.9 # = 0.85 eV E.A. Hodille et al 2021 Nucl. Fusion 61 126003, trap 1
detrapping_energy = 11604.6 # = 1.00 eV E.A. Hodille et al 2021 Nucl. Fusion 61 126003, trap 2
variable = C_trapped_W
[]
############################## NodalKernels for Cu (block = 3)
[time_Cu]
type = TimeDerivativeNodalKernel
variable = C_trapped_Cu
[]
[trapping_Cu]
type = TrappingNodalKernel
variable = C_trapped_Cu
temperature = temperature
alpha_t = 2.75e11 # 1e15
N = 1.0e0 # = ${tungsten_atomic_density} #/m^3 (W lattice density)
Ct0 = 5.0e-5 # R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 3
trap_per_free = 1.0e0 # 1.0e1
mobile_concentration = 'C_mobile_Cu'
extra_vector_tags = ref
[]
[release_Cu]
type = ReleasingNodalKernel
alpha_r = 8.4e12 # 1.0e13
temperature = temperature
detrapping_energy = 5802.3 # = 0.50eV R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 3
variable = C_trapped_Cu
[]
############################## NodalKernels for CuCrZr (block = 2)
[time_CuCrZr]
type = TimeDerivativeNodalKernel
variable = C_trapped_CuCrZr
[]
[trapping_CuCrZr]
type = TrappingNodalKernel
variable = C_trapped_CuCrZr
temperature = temperature
alpha_t = 2.75e11 # 1e15
N = 1.0e0 # = ${tungsten_atomic_density} #/m^3 (W lattice density)
Ct0 = 5.0e-5 # R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 4
# Ct0 = 4.0e-2 # R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 5
trap_per_free = 1.0e0 # 1.0e1
mobile_concentration = 'C_mobile_CuCrZr'
extra_vector_tags = ref
[]
[release_CuCrZr]
type = ReleasingNodalKernel
alpha_r = 8.4e12 # 1.0e13
temperature = temperature
detrapping_energy = 5802.3 # = 0.50eV R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 4
# detrapping_energy = 9631.8 # = 0.83 eV R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038, trap 5
variable = C_trapped_CuCrZr
[]
[]
[BCs]
[C_mob_W_top_flux]
type = FunctionNeumannBC
variable = C_mobile_W
boundary = 'top'
function = mobile_flux_bc_func
[]
[mobile_tube]
type = DirichletBC
variable = C_mobile_CuCrZr
value = 1.0e-18
boundary = '2to1'
[]
[temp_top]
type = FunctionNeumannBC
variable = temperature
boundary = 'top'
function = temp_flux_bc_func
[]
[temp_tube]
type = FunctionDirichletBC
variable = temperature
boundary = '2to1'
function = temp_inner_func
[]
[]
[Functions]
### Maximum mobile flux of 7.90e-13 at the top surface (1.0e-4 [m])
### 1.80e23/m^2-s = (5.0e23/m^2-s) *(1-0.999) = (7.90e-13)*(${tungsten_atomic_density})/(1.0e-4) at steady state
[mobile_flux_bc_func]
type = ParsedFunction
expression = 'if((t % 1600) < 100.0, 0.0 + 7.90e-13*(t % 1600)/100,
if((t % 1600) < 500.0, 7.90e-13,
if((t % 1600) < 600.0, 7.90e-13 - 7.90e-13*((t % 1600)-500)/100, 0.0)))'
[]
### Heat flux of 10MW/m^2 at steady state
[temp_flux_bc_func]
type = ParsedFunction
expression = 'if((t % 1600) < 100.0, 0.0 + 1.0e7*(t % 1600)/100,
if((t % 1600) < 500.0, 1.0e7,
if((t % 1600) < 600.0, 1.0e7 - 1.0e7*((t % 1600)-500)/100, 300)))'
[]
### Maximum coolant temperature of 552K at steady state
[temp_inner_func]
type = ParsedFunction
expression = 'if((t % 1600) < 100.0, 300.0 + (552-300)*(t % 1600)/100,
if((t % 1600) < 500.0, 552,
if((t % 1600) < 600.0, 552.0 - (552-300)*((t % 1600)-500)/100, 300)))'
[]
[timestep_function]
type = ParsedFunction
expression = 'if((t % 1600) < 10.0, 20,
if((t % 1600) < 90.0, 40,
if((t % 1600) < 110.0, 20,
if((t % 1600) < 480.0, 40,
if((t % 1600) < 500.0, 20,
if((t % 1600) < 590.0, 4,
if((t % 1600) < 610.0, 20,
if((t % 1600) < 1500.0, 200,
if((t % 1600) < 1600.0, 40, 2)))))))))'
[]
[]
[Materials]
############################## Materials for W (block = 4)
[diffusivity_W]
type = ADParsedMaterial
property_name = diffusivity_W
coupled_variables = 'temperature'
block = 4
expression = '2.4e-7*exp(-4525.8/temperature)' # H diffusivity in W
outputs = all
[]
[solubility_W]
type = ADParsedMaterial
property_name = solubility_W
coupled_variables = 'temperature'
block = 4
# expression = '2.95e-5 *exp(-12069.0/temperature)' # H solubility in W = (1.87e24)/(${tungsten_atomic_density}) [#/m^3]
expression = '2.95e-5 *exp(-12069.0/temperature) + 4.95e-8 * exp(-6614.6/temperature)' # H solubility in W = (1.87e24)/(${tungsten_atomic_density}) [#/m^3]
outputs = all
[]
[converter_to_regular_W]
type = MaterialADConverter
ad_props_in = 'diffusivity_W'
reg_props_out = 'diffusivity_W_nonAD'
block = 4
[]
[heat_transfer_W]
type = GenericConstantMaterial
prop_names = 'density_W'
prop_values = '19300' # [g/m^3]
block = 4
[]
[specific_heat_W]
type = ParsedMaterial
property_name = specific_heat_W
coupled_variables = 'temperature'
block = 4
expression = '1.16e2 + 7.11e-2 * temperature - 6.58e-5 * temperature^2 + 3.24e-8 * temperature^3 -5.45e-12 * temperature^4' # ~ 132[J/kg-K]
outputs = all
[]
[thermal_conductivity_W]
type = ParsedMaterial
property_name = thermal_conductivity_W
coupled_variables = 'temperature'
block = 4
# expression = '-7.8e-9 * temperature^3 + 5.0e-5 * temperature^2 - 1.1e-1 * temperature + 1.8e2' # ~ 173.0 [ W/m-K] from R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038,
expression = '2.41e2 - 2.90e-1 * temperature + 2.54e-4 * temperature^2 - 1.03e-7 * temperature^3 + 1.52e-11 * temperature^4' # ~ 173.0 [ W/m-K]
outputs = all
[]
############################## Materials for Cu (block = 3)
[diffusivity_Cu]
type = ADParsedMaterial
property_name = diffusivity_Cu
coupled_variables = 'temperature'
block = 3
expression = '6.60e-7*exp(-4525.8/temperature)' # H diffusivity in Cu
outputs = all
[]
[solubility_Cu]
type = ADParsedMaterial
property_name = solubility_Cu
coupled_variables = 'temperature'
block = 3
expression = '4.95e-5*exp(-6614.6/temperature)' # H solubility in Cu = (3.14e24)/(${tungsten_atomic_density}) [#/m^3]
outputs = all
[]
[converter_to_regular_Cu]
type = MaterialADConverter
ad_props_in = 'diffusivity_Cu'
reg_props_out = 'diffusivity_Cu_nonAD'
block = 3
[]
[heat_transfer_Cu]
type = GenericConstantMaterial
prop_names = 'density_Cu'
prop_values = '8960.0' # [g/m^3]
block = 3
[]
[specific_heat_Cu]
type = ParsedMaterial
property_name = specific_heat_Cu
coupled_variables = 'temperature'
block = 3
expression = '3.16e2 + 3.18e-1 * temperature - 3.49e-4 * temperature^2 + 1.66e-7 * temperature^3' # ~ 384 [J/kg-K]
outputs = all
[]
[thermal_conductivity_Cu]
type = ParsedMaterial
property_name = thermal_conductivity_Cu
coupled_variables = 'temperature'
block = 3
# expression = '-3.9e-8 * temperature^3 + 3.8e-5 * temperature^2 - 7.9e-2 * temperature + 4.0e2' # ~ 401.0 [ W/m-K] from R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038,
expression = '4.21e2 - 6.85e-2 * temperature' # ~ 400.0 [ W/m-K]
outputs = all
[]
############################## Materials for CuCrZr (block = 2)
[diffusivity_CuCrZr]
type = ADParsedMaterial
property_name = diffusivity_CuCrZr
coupled_variables = 'temperature'
block = 2
expression = '3.90e-7*exp(-4873.9/temperature)' # H diffusivity in CuCrZr
outputs = all
[]
[solubility_CuCrZr]
type = ADParsedMaterial
property_name = solubility_CuCrZr
coupled_variables = 'temperature'
block = 2
expression = '6.75e-6*exp(-4525.8/temperature)' # H solubility in CuCrZr = (4.28e23)/(${tungsten_atomic_density}) [#/m^3]
outputs = all
[]
[converter_to_regular_CuCrZr]
type = MaterialADConverter
ad_props_in = 'diffusivity_CuCrZr'
reg_props_out = 'diffusivity_CuCrZr_nonAD'
block = 2
[]
[heat_transfer_CuCrZr]
type = GenericConstantMaterial
prop_names = 'density_CuCrZr specific_heat_CuCrZr'
prop_values = '8900.0 390.0' # [g/m^3], [ W/m-K], [J/kg-K]
block = 2
[]
[thermal_conductivity_CuCrZr]
type = ParsedMaterial
property_name = thermal_conductivity_CuCrZr
coupled_variables = 'temperature'
block = 2
# expression = '5.3e-7 * temperature^3 - 6.5e-4 * temperature^2 + 2.6e-1 * temperature + 3.1e2' # ~ 320.0 [ W/m-K] from R. Delaporte-Mathurin et al 2021 Nucl. Fusion 61 036038,
expression = '3.87e2 - 1.28e-1 * temperature' # ~ 349 [ W/m-K]
outputs = all
[]
############################## Materials for others
[interface_jump_4to3]
type = SolubilityRatioMaterial
solubility_primary = solubility_W
solubility_secondary = solubility_Cu
boundary = '4to3'
concentration_primary = C_mobile_W
concentration_secondary = C_mobile_Cu
[]
[interface_jump_3to2]
type = SolubilityRatioMaterial
solubility_primary = solubility_Cu
solubility_secondary = solubility_CuCrZr
boundary = '3to2'
concentration_primary = C_mobile_Cu
concentration_secondary = C_mobile_CuCrZr
[]
[]
[Postprocessors]
############################################################ Postprocessors for W (block = 4)
[F_recombination]
type = SideDiffusiveFluxAverage
boundary = 'top'
diffusivity = 5.01e-24 # (3.01604928)/(6.02e23)/[gram(T)/m^2]
# diffusivity = 5.508e-19 # (1.0e3)*(1.0e3)/(6.02e23)/(3.01604928) [gram(T)/m^2]
variable = Sc_C_total_W
[]
[F_permeation]
type = SideDiffusiveFluxAverage
boundary = '2to1'
diffusivity = 5.01e-24 # (3.01604928)/(6.02e23)/[gram(T)/m^2]
# diffusivity = 5.508e-19 # (1.0e3)*(1.0e3)/(6.02e23)/(3.01604928) [gram(T)/m^2]
variable = Sc_C_total_CuCrZr
[]
[Int_C_mobile_W]
type = ElementIntegralVariablePostprocessor
variable = C_mobile_W
block = 4
[]
[ScInt_C_mobile_W]
type = ScalePostprocessor
value = Int_C_mobile_W
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
[Int_C_trapped_W]
type = ElementIntegralVariablePostprocessor
variable = C_trapped_W
block = 4
[]
[ScInt_C_trapped_W]
type = ScalePostprocessor
value = Int_C_trapped_W
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
[Int_C_total_W]
type = ElementIntegralVariablePostprocessor
variable = C_total_W
block = 4
[]
[ScInt_C_total_W]
type = ScalePostprocessor
value = Int_C_total_W
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
# ############################################################ Postprocessors for Cu (block = 3)
[Int_C_mobile_Cu]
type = ElementIntegralVariablePostprocessor
variable = C_mobile_Cu
block = 3
[]
[ScInt_C_mobile_Cu]
type = ScalePostprocessor
value = Int_C_mobile_Cu
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
[Int_C_trapped_Cu]
type = ElementIntegralVariablePostprocessor
variable = C_trapped_Cu
block = 3
[]
[ScInt_C_trapped_Cu]
type = ScalePostprocessor
value = Int_C_trapped_Cu
scaling_factor = 3.44e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
[Int_C_total_Cu]
type = ElementIntegralVariablePostprocessor
variable = C_total_Cu
block = 3
[]
[ScInt_C_total_Cu]
type = ScalePostprocessor
value = Int_C_total_Cu
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
# ############################################################ Postprocessors for CuCrZr (block = 2)
[Int_C_mobile_CuCrZr]
type = ElementIntegralVariablePostprocessor
variable = C_mobile_CuCrZr
block = 2
[]
[ScInt_C_mobile_CuCrZr]
type = ScalePostprocessor
value = Int_C_mobile_CuCrZr
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
[Int_C_trapped_CuCrZr]
type = ElementIntegralVariablePostprocessor
variable = C_trapped_CuCrZr
block = 2
[]
[ScInt_C_trapped_CuCrZr]
type = ScalePostprocessor
value = Int_C_trapped_CuCrZr
scaling_factor = 3.44e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
[Int_C_total_CuCrZr]
type = ElementIntegralVariablePostprocessor
variable = C_total_CuCrZr
block = 2
[]
[ScInt_C_total_CuCrZr]
type = ScalePostprocessor
value = Int_C_total_CuCrZr
scaling_factor = 3.491e10 # (1.0e3)*(1.0e3)*(${tungsten_atomic_density})/(6.02e23)/(3.01604928) [gram(T)/m^2]
[]
############################################################ Postprocessors for others
[dt]
type = TimestepSize
[]
[temperature_top]
type = PointValue
variable = temperature
point = '0 14.0e-3 0'
[]
[temperature_tube]
type = PointValue
variable = temperature
point = '0 6.0e-3 0'
[]
# limit timestep
[timestep_max_pp] # s
type = FunctionValuePostprocessor
function = timestep_function
[]
[]
[VectorPostprocessors]
[line]
type = LineValueSampler
start_point = '0 14.0e-3 0'
end_point = '0 6.0e-3 0'
num_points = 100
sort_by = 'y'
variable = 'C_total_W C_total_Cu C_total_CuCrZr C_mobile_W C_mobile_Cu C_mobile_CuCrZr C_trapped_W C_trapped_Cu C_trapped_CuCrZr flux_y temperature'
execute_on = timestep_end
[]
[]
[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 # 1e-8 works for 1 cycle
nl_abs_tol = 1e-7 # 1e-11 works for 1 cycle
end_time = 8.0e4 # 50 ITER shots (3.0e4 s plasma, 2.0e4 SSP)
automatic_scaling = true
line_search = 'none'
dtmin = 1e-4
nl_max_its = 18
[TimeStepper]
type = IterationAdaptiveDT
dt = 20
optimal_iterations = 15
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.8
timestep_limiting_postprocessor = timestep_max_pp
[]
[]
[Outputs]
[exodus]
type = Exodus
sync_only = false
# output at key moment in the first two cycles, and then at the end of the simulation
sync_times = '110.0 480.0 590.0 1600.0 1710.0 2080.0 2190.0 3400.0 8.0e4'
[]
csv = true
hide = 'dt
Int_C_mobile_W Int_C_trapped_W Int_C_total_W
Int_C_mobile_Cu Int_C_trapped_Cu Int_C_total_Cu
Int_C_mobile_CuCrZr Int_C_trapped_CuCrZr Int_C_total_CuCrZr'
perf_graph = true
[]
(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
[]
[]