- Ct0The fraction of host sites that can contribute to trapping as a function (-)
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:The fraction of host sites that can contribute to trapping as a function (-)
- NThe atomic number density of the host material (1/m^3)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The atomic number density of the host material (1/m^3)
- alpha_tThe trapping rate coefficient. This has units of 1/s (e.g. no number densities are involved)
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The trapping rate coefficient. This has units of 1/s (e.g. no number densities are involved)
- mobile_concentrationThe variable representing the mobile concentration of solute particles (1/m^3)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The variable representing the mobile concentration of solute particles (1/m^3)
- temperatureThe temperature (K)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The temperature (K)
- 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
TrappingNodalKernel
Implements a residual describing the trapping of a species in a material.
Overview
This object implements the residual
where is the trapping rate coefficient, which has dimensions of 1/s, is the atomic number density of the host material, is the concentration of empty trapping sites, and is the concentration of the mobile species. is defined as with being the pre-exponential factor in units of 1/s, being the trapping energy in units of Kelvin, and being the temperature in units of Kelvin.
As outlined in Scaling, it is important to scale different specie numerical concentrations to the same order of magnitude in order to have robust (non)linear solves. Unfortunately, incorporation of scaling tends to obfuscate residual-computing objects a bit, including TrappingNodalKernel. Hopefully we can make things clear here, however. Let's imagine that the trapped concentration is 1000 times greater than the mobile concentration. We can bring the trapped concentration to the same numerical level of the mobile concentration by changing the trapped concentration units from #/volume to k#/volume where 'k' denotes kilo. Now let's figure out how this transformation fits into TrappingNodalKernel. TrappingNodalKernel is producing trapped species, so consequently its residual must have units (in our example) of k#/(volume*s). Trapped concentrations have units of k#/volume and the mobile concentration and host density () have units of #/volume. TrappingNodalKernel computes the empty trapping sites concentration, in the following way:
auto empty_trapping_sites = _Ct0.value(_t, (*_current_node)) * _N;
for (const auto & trap_conc : _trapped_concentrations)
empty_trapping_sites -= (*trap_conc)[_qp] * _trap_per_free;
The trapping concentration, in units of k#/volume, is converted to units of #/volume by multiplying by trap_per_free which, in this example of k#/volume trapping concentration and #/volume mobile concentration, has a value of 1000 #/(k#). We then compute the residual with the code
return -_alpha_t * std::exp(-_detrapping_energy / _temperature[_qp]) * empty_trapping_sites * _mobile_conc[_qp] / (_N * _trap_per_free);
Let's carry through the units: 1/s * #/volume * #/volume / (#/volume * 1000 #/(k#)) -> k#/(s*volume) which is exactly the units that we needed.
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
- 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)
- other_trapped_concentration_variablesOther variables representing trapped particle concentrations.
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Other variables representing trapped particle concentrations.
- trap_per_free1An estimate for the ratio of the concentration magnitude of trapped species to free species. Setting a value for this can be helpful in producing a well-scaled matrix (-)
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:An estimate for the ratio of the concentration magnitude of trapped species to free species. Setting a value for this can be helpful in producing a well-scaled matrix (-)
- trapping_energy0The trapping energy (K)
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The trapping energy (K)
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill
Contribution To Tagged Field Data Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- diag_save_inThe name of auxiliary variables to save this 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.)
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- save_inThe name of auxiliary variables to save this 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.)
- 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
Input Files
- (test/tests/divertor_monoblock/divertor_monoblock.i)
- (test/tests/val-2f/val-2f_trapping_3.i)
- (test/tests/val-2f/val-2f_trapping_4.i)
- (test/tests/ver-1jb/ver-1jb.i)
- (test/tests/val-2f/val-2f_trapping_5.i)
- (test/tests/ver-1d/ver-1d-trapping.i)
- (test/tests/ver-1d/ver-1d-diffusion.i)
- (test/tests/val-2f/val-2f_trapping_1.i)
- (test/tests/val-2d/val-2d.i)
- (test/tests/ver-1dc/ver-1dc_base.i)
- (test/tests/val-2f/val-2f_trapping_intrinsic.i)
- (test/tests/val-2f/val-2f_trapping_2.i)
- (test/tests/ver-1d/ver-1d-trapping_profile.i)
(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-2f/val-2f_trapping_3.i)
# This input file adds the trapping sites 3 for validation case val-2f.
# It is included in val-2f.i
[Variables]
[trapped_3]
order = FIRST
family = LAGRANGE
[]
[]
[Bounds]
[trapped_3_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_3
bound_type = lower
bound_value = '${fparse -1e-20}'
[]
[]
[Kernels]
# trapping 3 kernel
[coupled_time_trap_3]
type = ADCoefCoupledTimeDerivative
variable = deuterium_concentration_W
v = trapped_3
coef = ${trap_per_free_3}
[]
[]
[NodalKernels]
[time_3]
type = TimeDerivativeNodalKernel
variable = trapped_3
[]
[trapping_3]
type = TrappingNodalKernel
variable = trapped_3
mobile_concentration = deuterium_concentration_W
alpha_t = '${trapping_prefactor_3}'
trapping_energy = '${trapping_energy_3}'
N = '${tungsten_density}'
Ct0 = 'trap_distribution_function_3'
temperature = 'temperature'
trap_per_free = ${trap_per_free_3}
[]
[release_3]
type = ReleasingNodalKernel
variable = trapped_3
alpha_r = '${detrapping_prefactor_3}'
detrapping_energy = '${detrapping_energy_3}'
temperature = 'temperature'
[]
[]
[Functions]
[trap_distribution_function_3]
type = ParsedFunction
expression = '${trapping_site_fraction_3} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[trap_distribution_function_3_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_3_inf} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[]
[Materials]
[trap_distribution_function_3]
type = GenericFunctionMaterial
prop_names = trap_distribution_function_3
prop_values = trap_distribution_function_3
[]
[]
[Postprocessors]
[integral_trapped_concentration_3]
type = ElementIntegralVariablePostprocessor
variable = trapped_3
outputs = none
[]
[scaled_trapped_deuterium_3]
type = ScalePostprocessor
scaling_factor = '${fparse ${trap_per_free_3} * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_3
[]
[]
(test/tests/val-2f/val-2f_trapping_4.i)
# This input file adds the trapping sites 4 for validation case val-2f.
# It is included in val-2f.i
[Variables]
[trapped_4]
order = FIRST
family = LAGRANGE
[]
[]
[Kernels]
# trapping 4 kernel
[coupled_time_trap_4]
type = ADCoefCoupledTimeDerivative
variable = deuterium_concentration_W
v = trapped_4
coef = ${trap_per_free_4}
[]
[]
[NodalKernels]
[time_4]
type = TimeDerivativeNodalKernel
variable = trapped_4
[]
[trapping_4]
type = TrappingNodalKernel
variable = trapped_4
mobile_concentration = deuterium_concentration_W
alpha_t = '${trapping_prefactor_4}'
trapping_energy = '${trapping_energy_4}'
N = '${tungsten_density}'
Ct0 = 'trap_distribution_function_4'
temperature = 'temperature'
trap_per_free = ${trap_per_free_4}
[]
[release_4]
type = ReleasingNodalKernel
variable = trapped_4
alpha_r = '${detrapping_prefactor_4}'
detrapping_energy = '${detrapping_energy_4}'
temperature = 'temperature'
[]
[]
[Functions]
[trap_distribution_function_4]
type = ParsedFunction
expression = '${trapping_site_fraction_4} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[trap_distribution_function_4_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_4_inf} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[]
[Materials]
[trap_distribution_function_4]
type = GenericFunctionMaterial
prop_names = trap_distribution_function_4
prop_values = trap_distribution_function_4
[]
[]
[Postprocessors]
[integral_trapped_concentration_4]
type = ElementIntegralVariablePostprocessor
variable = trapped_4
outputs = none
[]
[scaled_trapped_deuterium_4]
type = ScalePostprocessor
scaling_factor = '${fparse ${trap_per_free_4} * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_4
[]
[]
(test/tests/ver-1jb/ver-1jb.i)
# Verification Problem #1jb from TMAP7 V&V document
# Radioactive Decay of Tritium in a Distributed Trap
# Physical Constants
ideal_gas_constant = ${units 8.31446261815324 J/K/mol} # from PhysicalConstants.h
boltzmann_constant = ${units 1.380649e-23 J/K -> eV/K } # from PhysicalConstants.h
# Case and model parameters (adapted from TMAP7)
slab_length = ${units 1.5 m }
tritium_mobile_concentration_initial = ${units 1 atoms/m3}
trapping_sites_atomic_fraction_max = 0.001 # (-)
trapping_sites_fraction_occupied_initial = 0.5 # (-)
normal_center_position = ${fparse slab_length/2} # m
normal_standard_deviation = ${fparse slab_length/4} # m
density_material = ${units 6.34e28 atoms/m^3} # for tungsten
density_scalar = ${units 1e25 atoms/m^3} # used to scale variables and use a dimensionless system
temperature = ${units 300 K} # assumed (TMAP7's input file lists 273 K)
tritium_diffusivity_prefactor = ${units 1.58e-4 m^2/s} # from TMAP7 V&V input file
tritium_diffusivity_energy = ${units 308000.0 J/K} # from TMAP7 V&V input file
tritium_release_prefactor = ${units 1.0e13 1/s} # from TMAP7 V&V input file
tritium_release_energy = ${units 4.2 eV}
tritium_trapping_prefactor = ${units 2.096e15 1/s} # from TMAP7 V&V input file
tritium_trapping_energy = ${tritium_diffusivity_energy} # J/K - from TMAP7 V&V input file
trap_per_free = 1.e-25 # (-)
half_life_s = ${units 12.3232 year -> s}
decay_rate_constant = ${fparse 0.693/half_life_s} # 1/s
# Simulation parameters
num_mesh_element = 50
end_time = ${units 100 year -> s} # s
dt_start = ${fparse end_time/250} # s
dt_max = ${fparse end_time/100} # s
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${num_mesh_element}
xmax = ${slab_length}
[]
[Variables]
# tritium mobile concentration in atoms/m^3 / density_scalar = (-)
[tritium_mobile_concentration_scaled]
initial_condition = ${fparse tritium_mobile_concentration_initial / density_scalar}
[]
# tritium trapped concentration in (atoms/m^3) / density_scalar = (-)
[tritium_trapped_concentration_scaled]
[]
# helium concentration in (atoms/m^3) / density_scalar = (-)
[helium_concentration_scaled]
[]
[]
[Functions]
[trapping_sites_fraction_function] # (atomic fraction)
type = ParsedFunction
expression = '${trapping_sites_atomic_fraction_max} * exp(-1/2*((x-${normal_center_position})/${normal_standard_deviation})^2)'
[]
[density_material_function] # (atoms/m^3)
type = ConstantFunction
value = ${density_material}
[]
[trapping_sites_concentration_function] # (atoms/m^3)
type = CompositeFunction
functions = 'density_material_function trapping_sites_fraction_function'
[]
[initial_trapping_sites_occupied_function] # (-)
type = ConstantFunction
value = ${trapping_sites_fraction_occupied_initial}
[]
[tritium_trapped_concentration_initial_function] # (atoms/m^3)
type = CompositeFunction
functions = 'initial_trapping_sites_occupied_function trapping_sites_concentration_function'
[]
[density_scalar_inverse_function] # (atoms/m^3)^(-1)
type = ConstantFunction
value = ${fparse 1/density_scalar}
[]
[tritium_trapped_concentration_initial_function_scaled] # (-)
type = CompositeFunction
functions = 'density_scalar_inverse_function tritium_trapped_concentration_initial_function'
[]
[]
[ICs]
[tritium_trapped_concentration_IC]
type = FunctionIC
variable = tritium_trapped_concentration_scaled
function = 'tritium_trapped_concentration_initial_function_scaled'
[]
[]
[AuxVariables]
[empty_sites] # (-)
[]
[scaled_empty_sites] # atoms/m^3
[]
[trapped_sites] # (atoms/m^3)
[]
[total_sites] # (atoms/m^3)
[]
[temperature] # (K)
initial_condition = ${temperature}
[]
[tritium_mobile_concentration] # (atoms/m^3)
[]
[tritium_trapped_concentration] # (atoms/m^3)
[]
[helium_concentration] # (atoms/m^3)
[]
[]
[AuxKernels]
[empty_sites]
variable = empty_sites
type = EmptySitesAux
N = ${fparse density_material / density_scalar} # (-)
Ct0 = 'trapping_sites_fraction_function' # atomic fraction
trap_per_free = ${trap_per_free}
trapped_concentration_variables = tritium_trapped_concentration_scaled
[]
[scaled_empty_sites]
variable = scaled_empty_sites
type = NormalizationAux
normal_factor = ${density_scalar}
source_variable = empty_sites
[]
[trapped_sites]
variable = trapped_sites
type = NormalizationAux
normal_factor = ${density_scalar}
source_variable = tritium_trapped_concentration_scaled
[]
[total_sites]
variable = total_sites
type = ParsedAux
coupled_variables = 'trapped_sites scaled_empty_sites'
expression = 'trapped_sites + scaled_empty_sites'
[]
[tritium_mobile_concentration]
type = NormalizationAux
variable = tritium_mobile_concentration
normal_factor = ${density_scalar}
source_variable = tritium_mobile_concentration_scaled
[]
[tritium_trapped_concentration]
type = NormalizationAux
variable = tritium_trapped_concentration
normal_factor = ${density_scalar}
source_variable = tritium_trapped_concentration_scaled
[]
[helium_concentration]
type = NormalizationAux
variable = helium_concentration
normal_factor = ${density_scalar}
source_variable = helium_concentration_scaled
[]
[]
[Kernels]
# kernels for the tritium concentration equation
[time_tritium]
type = TimeDerivative
variable = tritium_mobile_concentration_scaled
[]
[diffusion]
type = MatDiffusion
variable = tritium_mobile_concentration_scaled
diffusivity = diffusivity
[]
[decay_tritium]
type = MatReaction
variable = tritium_mobile_concentration_scaled
v = tritium_mobile_concentration_scaled
reaction_rate = '${fparse -decay_rate_constant}'
[]
[coupled_time_tritium]
type = ScaledCoupledTimeDerivative
variable = tritium_mobile_concentration_scaled
v = tritium_trapped_concentration_scaled
factor = ${trap_per_free}
[]
# re-adding it to the equation of mobile tritium because it is accounted for in coupled_time_tritium, and needs to be removed
[decay_tritium_trapped]
type = MatReaction
variable = tritium_mobile_concentration_scaled
v = tritium_trapped_concentration_scaled
reaction_rate = '${fparse - decay_rate_constant * trap_per_free}'
[]
# kernels for the helium concentration equation
[time_helium]
type = TimeDerivative
variable = helium_concentration_scaled
[]
[decay_helium_mobile]
type = MatReaction
variable = helium_concentration_scaled
v = tritium_mobile_concentration_scaled
reaction_rate = '${fparse decay_rate_constant}'
[]
[decay_helium_trapped]
type = MatReaction
variable = helium_concentration_scaled
v = tritium_trapped_concentration_scaled
reaction_rate = '${fparse decay_rate_constant}'
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = tritium_trapped_concentration_scaled # (atoms/m^3) / density_scalar = (-)
[]
[trapping]
type = TrappingNodalKernel
variable = tritium_trapped_concentration_scaled # (-)
alpha_t = ${tritium_trapping_prefactor} # (1/s)
trapping_energy = ${fparse tritium_trapping_energy/ideal_gas_constant} # (K)
N = ${fparse density_material / density_scalar} # (-)
Ct0 = 'trapping_sites_fraction_function' # (atomic fraction)
mobile_concentration = 'tritium_mobile_concentration_scaled' # (-)
temperature = temperature # (K)
trap_per_free = ${trap_per_free}
[]
[release]
type = ReleasingNodalKernel
variable = tritium_trapped_concentration_scaled # (atoms/m^3) / density_scalar = (-)
alpha_r = ${fparse tritium_release_prefactor} # (1/s)
detrapping_energy = ${fparse tritium_release_energy / boltzmann_constant} # (K)
temperature = temperature # (K)
[]
[decay]
type = ReleasingNodalKernel
variable = tritium_trapped_concentration_scaled # (atoms/m^3) / density_scalar = (-)
alpha_r = ${decay_rate_constant} # (1/s) - decay rate
detrapping_energy = 0 # (K) - The decay rate is independent of temperature
temperature = temperature # (K)
[]
[]
[Materials]
[diffusivity] # (m2/s) tritium diffusivity
type = DerivativeParsedMaterial
property_name = 'diffusivity'
coupled_variables = 'temperature'
expression = '${tritium_diffusivity_prefactor}*exp(-${tritium_diffusivity_energy}/${ideal_gas_constant}/temperature)'
[]
[alpha_t_tot] # (1/s) - trapping rate
type = ParsedMaterial
property_name = 'alpha_t_tot'
coupled_variables = 'temperature'
expression = '${tritium_trapping_prefactor} * exp(- ${tritium_trapping_energy}/${ideal_gas_constant}/temperature)'
outputs = 'all'
[]
[alpha_r_tot] # (1/s) - detrapping rate
type = ParsedMaterial
property_name = 'alpha_r_tot'
coupled_variables = 'temperature'
expression = '${tritium_release_prefactor} * exp(- ${tritium_release_energy} / ${boltzmann_constant}/temperature)'
outputs = 'all'
[]
[]
[Postprocessors]
# Amount of mobile tritium in the sample in (atoms/m^3) / density_scalar * m = (m)
[tritium_mobile_inventory_scaled]
type = ElementIntegralVariablePostprocessor
variable = tritium_mobile_concentration_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[tritium_mobile_inventory] # (atoms/m^2)
type = ScalePostprocessor
value = tritium_mobile_inventory_scaled
scaling_factor = ${density_scalar}
execute_on = 'INITIAL TIMESTEP_END'
[]
# Amount of trapped tritium in the sample (atoms/m^3) / density_scalar * m = (m)
[tritium_trapped_inventory_scaled]
type = ElementIntegralVariablePostprocessor
variable = tritium_trapped_concentration_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[tritium_trapped_inventory] # (atoms/m^2)
type = ScalePostprocessor
value = tritium_trapped_inventory_scaled
scaling_factor = ${density_scalar}
execute_on = 'INITIAL TIMESTEP_END'
[]
# Amount of helium in the sample in atoms/m^3 / density_scalar * m = (m)
[helium_inventory_scaled]
type = ElementIntegralVariablePostprocessor
variable = helium_concentration_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[helium_inventory] # (atoms/m^2)
type = ScalePostprocessor
value = helium_inventory_scaled
scaling_factor = ${density_scalar}
execute_on = 'INITIAL TIMESTEP_END'
[]
# check mass conservation - this should remain constant - (atoms/m^2)
[total_inventory] # (atoms/m^2)
type = SumPostprocessor
values = 'tritium_mobile_inventory tritium_trapped_inventory helium_inventory'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[VectorPostprocessors]
[line]
type = LineValueSampler
start_point = '0 0 0'
end_point = '${slab_length} 0 0'
num_points = ${num_mesh_element}
sort_by = 'x'
variable = 'tritium_mobile_concentration tritium_trapped_concentration helium_concentration'
execute_on = 'timestep_end'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${end_time}
solve_type = NEWTON
scheme = 'bdf2'
dtmin = 1
dtmax = ${dt_max}
petsc_options = '-snes_converged_reason'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt_start}
optimal_iterations = 9
growth_factor = 1.2
cutback_factor = 0.9
[]
[]
[Outputs]
perf_graph = true
file_base = ver-1jb_out
[time_dependent_out]
type = CSV
execute_vector_postprocessors_on = NONE
file_base = ver-1jb_time_dependent_out
[]
[profile_out]
type = CSV
sync_only = true
sync_times = '${units 45 year -> s}'
execute_postprocessors_on = NONE
file_base=ver-1jb_profile_out
[]
[]
(test/tests/val-2f/val-2f_trapping_5.i)
# This input file adds the trapping sites 5 for validation case val-2f.
# It is included in val-2f.i
[Variables]
[trapped_5]
order = FIRST
family = LAGRANGE
[]
[]
[Kernels]
# trapping 5 kernel
[coupled_time_trap_5]
type = ADCoefCoupledTimeDerivative
variable = deuterium_concentration_W
v = trapped_5
coef = ${trap_per_free_5}
[]
[]
[NodalKernels]
[time_5]
type = TimeDerivativeNodalKernel
variable = trapped_5
[]
[trapping_5]
type = TrappingNodalKernel
variable = trapped_5
mobile_concentration = deuterium_concentration_W
alpha_t = '${trapping_prefactor_5}'
trapping_energy = '${trapping_energy_5}'
N = '${tungsten_density}'
Ct0 = 'trap_distribution_function_5'
temperature = 'temperature'
trap_per_free = ${trap_per_free_5}
[]
[release_5]
type = ReleasingNodalKernel
variable = trapped_5
alpha_r = '${detrapping_prefactor_5}'
detrapping_energy = '${detrapping_energy_5}'
temperature = 'temperature'
[]
[]
[Functions]
[trap_distribution_function_5]
type = ParsedFunction
expression = '${trapping_site_fraction_5} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[trap_distribution_function_5_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_5_inf} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[]
[Materials]
[trap_distribution_function_5]
type = GenericFunctionMaterial
prop_names = trap_distribution_function_5
prop_values = trap_distribution_function_5
[]
[]
[Postprocessors]
[integral_trapped_concentration_5]
type = ElementIntegralVariablePostprocessor
variable = trapped_5
outputs = none
[]
[scaled_trapped_deuterium_5]
type = ScalePostprocessor
scaling_factor = '${fparse ${trap_per_free_5} * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_5
[]
[]
(test/tests/ver-1d/ver-1d-trapping.i)
cl = 3.1622e18
trap_per_free = 1e3
N = 3.1622e22
time_scaling = 1
epsilon=10000
temperature = 1000
[Mesh]
type = GeneratedMesh
dim = 1
nx = 1000
xmax = 1
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile]
[]
[trapped]
[]
[]
[AuxVariables]
[empty_sites]
[]
[scaled_empty_sites]
[]
[trapped_sites]
[]
[total_sites]
[]
[]
[AuxKernels]
[empty_sites]
variable = empty_sites
type = EmptySitesAux
N = '${fparse N / cl}'
Ct0 = .1
trap_per_free = ${trap_per_free}
trapped_concentration_variables = trapped
[]
[scaled_empty]
variable = scaled_empty_sites
type = NormalizationAux
normal_factor = ${cl}
source_variable = empty_sites
[]
[trapped_sites]
variable = trapped_sites
type = NormalizationAux
normal_factor = ${trap_per_free}
source_variable = trapped
[]
[total_sites]
variable = total_sites
type = ParsedAux
expression = 'trapped_sites + empty_sites'
coupled_variables = 'trapped_sites empty_sites'
[]
[]
[Kernels]
[diff]
type = MatDiffusion
variable = mobile
diffusivity = '${fparse 1. / time_scaling}'
extra_vector_tags = ref
[]
[time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[]
[coupled_time]
type = ScaledCoupledTimeDerivative
variable = mobile
v = trapped
factor = ${trap_per_free}
extra_vector_tags = ref
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
alpha_t = '${fparse 1e15 / time_scaling}'
N = '${fparse N / cl}'
Ct0 = 0.1
mobile_concentration = 'mobile'
temperature = ${temperature}
trap_per_free = ${trap_per_free}
extra_vector_tags = ref
[]
[release]
type = ReleasingNodalKernel
alpha_r = ${fparse 1e13 / time_scaling}
temperature = ${temperature}
detrapping_energy = ${epsilon}
variable = trapped
[]
[]
[BCs]
[left]
type = FunctionDirichletBC
variable = mobile
function = 'BC_func'
boundary = left
[]
[right]
type = DirichletBC
variable = mobile
value = 0
boundary = right
[]
[]
[Functions]
[BC_func]
type = ParsedFunction
expression = '${fparse cl / cl}*tanh( 3 * t )'
[]
[]
[Postprocessors]
[outflux]
type = SideDiffusiveFluxAverage
boundary = 'right'
diffusivity = 1
variable = mobile
[]
[scaled_outflux]
type = ScalePostprocessor
value = outflux
scaling_factor = ${cl}
[]
[min_trapped]
type = NodalExtremeValue
value_type = MIN
variable = trapped
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = 1000
dtmax = 5
solve_type = NEWTON
scheme = BDF2
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
line_search = 'none'
[TimeStepper]
type = IterationAdaptiveDT
dt = 1e-6
optimal_iterations = 9
growth_factor = 1.1
cutback_factor = 0.909
[]
[]
[Outputs]
exodus = true
csv = true
[dof]
type = DOFMap
execute_on = initial
[]
perf_graph = true
[]
(test/tests/ver-1d/ver-1d-diffusion.i)
cl=3.1622e18
temperature = 1000
[Mesh]
type = GeneratedMesh
dim = 1
nx = 200
xmax = 1
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile][]
[trapped][]
[]
[Kernels]
[./diff]
type = Diffusion
variable = mobile
extra_vector_tags = ref
[../]
[./time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[../]
[coupled_time]
type = CoupledTimeDerivative
variable = mobile
v = trapped
extra_vector_tags = ref
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
alpha_t = 1e15
N = ${fparse 3.1622e22 / cl}
Ct0 = 0.1
mobile_concentration = 'mobile'
temperature = ${temperature}
extra_vector_tags = ref
[]
[release]
type = ReleasingNodalKernel
alpha_r = 1e13
temperature = ${temperature}
detrapping_energy = 100
variable = trapped
[]
[]
[BCs]
[left]
type = DirichletBC
variable = mobile
value = ${fparse 3.1622e18 / cl}
boundary = left
[]
[right]
type = DirichletBC
variable = mobile
value = 0
boundary = right
[]
[]
[Postprocessors]
[outflux]
type = SideDiffusiveFluxAverage
boundary = 'right'
diffusivity = 1
variable = mobile
[]
[scaled_outflux]
type = ScalePostprocessor
value = outflux
scaling_factor = ${cl}
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = 3
dt = .01
dtmin = .01
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
verbose = true
compute_scaling_once = false
[]
[Outputs]
csv = true
exodus = true
[dof]
type = DOFMap
execute_on = initial
[]
perf_graph = true
[]
(test/tests/val-2f/val-2f_trapping_1.i)
# This input file adds the trapping sites 1 for validation case val-2f.
# It is included in val-2f.i
[Variables]
[trapped_1]
order = FIRST
family = LAGRANGE
[]
[]
[Bounds]
[trapped_1_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_1
bound_type = lower
bound_value = '${fparse -1e-20}'
[]
[]
[Kernels]
# trapping 1 kernel
[coupled_time_trap_1]
type = ADCoefCoupledTimeDerivative
variable = deuterium_concentration_W
v = trapped_1
coef = ${trap_per_free_1}
[]
[]
[NodalKernels]
# First traps
[time_1]
type = TimeDerivativeNodalKernel
variable = trapped_1
[]
[trapping_1]
type = TrappingNodalKernel
variable = trapped_1
mobile_concentration = deuterium_concentration_W
alpha_t = '${trapping_prefactor_1}'
trapping_energy = '${trapping_energy_1}'
N = '${tungsten_density}'
Ct0 = 'trap_distribution_function_1'
temperature = 'temperature'
trap_per_free = ${trap_per_free_1}
[]
[release_1]
type = ReleasingNodalKernel
variable = trapped_1
alpha_r = '${detrapping_prefactor_1}'
detrapping_energy = '${detrapping_energy_1}'
temperature = 'temperature'
[]
[]
[Functions]
[trap_distribution_function_1]
type = ParsedFunction
expression = '${trapping_site_fraction_1} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[trap_distribution_function_1_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_1_inf} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[]
[Materials]
[trap_distribution_function_1]
type = GenericFunctionMaterial
prop_names = trap_distribution_function_1
prop_values = trap_distribution_function_1
[]
[]
[Postprocessors]
[integral_trapped_concentration_1]
type = ElementIntegralVariablePostprocessor
variable = trapped_1
outputs = none
[]
[scaled_trapped_deuterium_1]
type = ScalePostprocessor
scaling_factor = '${fparse ${trap_per_free_1} * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_1
[]
[]
(test/tests/val-2d/val-2d.i)
# General parameters
kB = '${units 1.380649e-23 J/K}' # Boltzmann constant (from PhysicalConstants.h - https://physics.nist.gov/cgi-bin/cuu/Value?r)
# Model parameters
TDS_initial_time = '${units 5e3 s}'
TDS_critial_time_1 = '${units 5400 s}'
TDS_critial_time_2 = '${units 5404 s}'
simulation_time = '${units 6.8e3 s}'
outputs_initial_time = '${units 4000 s}'
step_interval_max = 50 # (-)
step_interval_mid = 15 # (-)
step_interval_min = 6 # (-)
bound_value_max = '${units 2e4 at/mum^3}'
bound_value_min = '${units -1e-10 at/mum^3}'
# Diffusion parameters
flux_high = '${units 1e19 at/m^2/s -> at/mum^2/s}'
flux_low = '${units 0 at/mum^2/s}'
diffusivity_coefficient = '${units 4.1e-7 m^2/s -> mum^2/s}'
E_D = '${units 0.39 eV -> J}'
initial_concentration = '${units 1e-10 at/m^3 -> at/mum^3}'
width_source = '${units 3e-9 m -> mum}'
depth_source = '${units 4.6e-9 m -> mum}'
# Traps parameters
N = '${units 6.25e28 at/m^3 -> at/mum^3}'
initial_concentration_trap_2 = 4.4e-10 # (-)
initial_concentration_trap_3 = 1.4e-10 # (-)
trapping_energy = '${fparse ${units 0.39 eV -> J} / ${kB}}'
detrapping_energy_1 = '${fparse ${units 1.2 eV -> J} / ${kB}}'
detrapping_energy_2 = '${fparse ${units 1.6 eV -> J} / ${kB}}'
detrapping_energy_3 = '${fparse ${units 3.1 eV -> J} / ${kB}}'
trapping_site_fraction_1 = 0.002156 # (-)
trapping_site_fraction_2 = 0.00175 # (-)
trapping_site_fraction_3 = 0.0020 # (-)
trapping_rate_prefactor = '${units 9.1316e12 1/s}'
release_rate_profactor = '${units 8.4e12 1/s}'
trap_per_free_1 = 1e6 # (-)
trap_per_free_2 = 1e4 # (-)
trap_per_free_3 = 1e4 # (-)
width_trap1 = '${units 10e-9 m -> mum}'
# thermal parameters
temperature_low = '${units 300 K}'
temperature_high = '${units 1273 K}'
temperature_rate = '${units ${fparse 50 / 60} K/s}'
[Mesh]
active = 'cartesian_mesh'
[cartesian_mesh]
nx_scale = 2
type = CartesianMeshGenerator
dim = 1
dx = '${fparse 10 * ${units 1.5e-9 m -> mum}}
${units 1e-9 m -> mum} ${units 1e-8 m -> mum} ${units 1e-7 m -> mum}
${units 4e-6 m -> mum} ${units 4.407e-6 m -> mum} ${fparse 11 * ${units 7.407e-6 m -> mum}}'
ix = '${fparse 10 * ${nx_scale}}
${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}}
${fparse 50 * ${nx_scale}} ${fparse 2} ${fparse 1}'
subdomain_id = '0 1 1 1 1 1 1'
[]
[cartesian_mesh_coarse]
nx_scale = 1
type = CartesianMeshGenerator
dim = 1
dx = '${fparse 10 * ${units 1.5e-9 m -> mum}}
${units 1e-9 m -> mum} ${units 1e-8 m -> mum} ${units 1e-7 m -> mum}
${units 4e-6 m -> mum} ${units 4.407e-6 m -> mum} ${fparse 11 * ${units 7.407e-6 m -> mum}}'
ix = '${fparse 10 * ${nx_scale}}
${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}} ${fparse 1 * ${nx_scale}}
${fparse 6 * ${nx_scale}} ${fparse 1} ${fparse 1}'
subdomain_id = '0 1 1 1 1 1 1'
[]
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Bounds]
[concentration_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = concentration
bound_type = lower
bound_value = ${bound_value_min}
[]
[trapped_1_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_1
bound_type = lower
bound_value = ${bound_value_min}
[]
[trapped_2_upper_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_2
bound_type = upper
bound_value = ${bound_value_max}
[]
[trapped_2_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_2
bound_type = lower
bound_value = ${bound_value_min}
[]
[trapped_3_upper_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_3
bound_type = upper
bound_value = ${bound_value_max}
[]
[trapped_3_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_3
bound_type = lower
bound_value = ${bound_value_min}
[]
[]
[Variables]
[concentration]
order = FIRST
family = LAGRANGE
initial_condition = ${initial_concentration}
[]
[trapped_1]
order = FIRST
family = LAGRANGE
block = 0
outputs = none
[]
[trapped_2]
order = FIRST
family = LAGRANGE
initial_condition = '${fparse initial_concentration_trap_2 * trapping_site_fraction_2 * N}'
outputs = none
[]
[trapped_3]
order = FIRST
family = LAGRANGE
initial_condition = '${fparse initial_concentration_trap_3 * trapping_site_fraction_3 * N}'
outputs = none
[]
[]
[AuxVariables]
[temperature]
[]
[bounds_dummy]
order = FIRST
family = LAGRANGE
[]
[]
[AuxKernels]
[temperature_Aux]
type = FunctionAux
variable = temperature
function = Temperature_function
execute_on = 'initial timestep_end linear'
[]
[]
[Kernels]
[diffusion_implantation]
type = ADMatDiffusion
variable = concentration
diffusivity = Diffusivity
extra_vector_tags = ref
[]
[time_diffusion_implantation]
type = ADTimeDerivative
variable = concentration
extra_vector_tags = ref
[]
[source]
type = ADBodyForce
variable = concentration
function = concentration_source_norm_function
[]
# trapping kernel
[coupled_time_trap_1]
type = ADCoefCoupledTimeDerivative
variable = concentration
v = trapped_1
coef = ${trap_per_free_1}
block = 0
extra_vector_tags = ref
[]
[coupled_time_trap_2_implantation]
type = ADCoefCoupledTimeDerivative
variable = concentration
v = trapped_2
coef = ${trap_per_free_2}
extra_vector_tags = ref
[]
[coupled_time_trap_3_implantation]
type = ADCoefCoupledTimeDerivative
variable = concentration
v = trapped_3
coef = ${trap_per_free_3}
extra_vector_tags = ref
[]
[]
[NodalKernels]
# First traps
[time_1]
type = TimeDerivativeNodalKernel
variable = trapped_1
[]
[trapping_1]
type = TrappingNodalKernel
variable = trapped_1
mobile_concentration = concentration
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = 'trap_1_distribution_function'
temperature = 'temperature'
trap_per_free = ${trap_per_free_1}
extra_vector_tags = ref
[]
[release_1]
type = ReleasingNodalKernel
variable = trapped_1
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_1}'
temperature = 'temperature'
[]
# Second traps
[time_2]
type = TimeDerivativeNodalKernel
variable = trapped_2
[]
[trapping_2_implantation]
type = TrappingNodalKernel
variable = trapped_2
mobile_concentration = concentration
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = '${trapping_site_fraction_2}'
temperature = 'temperature'
trap_per_free = ${trap_per_free_2}
extra_vector_tags = ref
[]
[release_2_implantation]
type = ReleasingNodalKernel
variable = trapped_2
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_2}'
temperature = 'temperature'
[]
# Third traps
[time_3]
type = TimeDerivativeNodalKernel
variable = trapped_3
[]
[trapping_3_implantation]
type = TrappingNodalKernel
variable = trapped_3
mobile_concentration = concentration
alpha_t = '${trapping_rate_prefactor}'
trapping_energy = '${trapping_energy}'
N = '${N}'
Ct0 = '${trapping_site_fraction_3}'
temperature = 'temperature'
trap_per_free = ${trap_per_free_3}
extra_vector_tags = ref
[]
[release_3_implantation]
type = ReleasingNodalKernel
variable = trapped_3
alpha_r = '${release_rate_profactor}'
detrapping_energy = '${detrapping_energy_3}'
temperature = 'temperature'
[]
[]
[BCs]
[left]
type = ADDirichletBC
variable = concentration
boundary = left
value = 0
[]
[right]
type = ADDirichletBC
variable = concentration
boundary = right
value = 0
[]
[]
[Materials]
[Diffusivity_implantation]
type = ADDerivativeParsedMaterial
property_name = 'Diffusivity'
functor_names = 'Temperature_function'
functor_symbols = 'Temperature_function'
expression = '${diffusivity_coefficient} * exp(- ${E_D} / ${kB} / Temperature_function)'
block = 0
output_properties = 'Diffusivity'
[]
[Diffusivity_Tungsten]
type = ADDerivativeParsedMaterial
property_name = 'Diffusivity'
functor_names = 'Temperature_function'
functor_symbols = 'Temperature_function'
expression = '${diffusivity_coefficient} * exp(- ${E_D} / ${kB} / Temperature_function) * 10'
block = 1
output_properties = 'Diffusivity'
[]
[converter_to_regular]
type = MaterialADConverter
ad_props_in = 'Diffusivity'
reg_props_out = 'Diffusivity_nonAD'
outputs = none
[]
[]
[Functions]
[Temperature_function]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time}, ${temperature_low},
if(t<${TDS_initial_time} + (${temperature_high} - ${temperature_low}) /
${temperature_rate}, ${temperature_low} + ${temperature_rate} * (t - ${TDS_initial_time}),
${temperature_high}))'
[]
[surface_flux_function]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time}, ${flux_high}, ${flux_low})'
[]
[source_distribution_function]
type = ParsedFunction
expression = '1 / ( ${width_source} * sqrt(2 * pi) ) * exp(-0.5 * ((x - ${depth_source}) / ${width_source} ) ^ 2)'
[]
[trap_1_distribution_function]
type = ParsedFunction
expression = ' ${trapping_site_fraction_1} / ( ${width_trap1} * sqrt(2 * pi) ) * exp(-0.5 * ((x - ${depth_source}) / ${width_trap1}) ^ 2)'
[]
[concentration_source_norm_function]
type = ParsedFunction
symbol_names = 'source_distribution_function surface_flux_function'
symbol_values = 'source_distribution_function surface_flux_function'
expression = 'source_distribution_function * surface_flux_function'
[]
[max_dt_size_function]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time} , ${step_interval_mid}, ${step_interval_min})'
[]
[max_dt_size_function_coarse]
type = ParsedFunction
expression = 'if(t<${TDS_initial_time} , ${step_interval_mid},
if(t<${TDS_critial_time_1}, ${step_interval_max},
if(t<${TDS_critial_time_2}, ${step_interval_min}, ${step_interval_max})))'
[]
[]
[Postprocessors]
[flux_surface_left]
type = SideDiffusiveFluxIntegral
variable = concentration
diffusivity = 'Diffusivity_nonAD'
boundary = 'left'
outputs = none
[]
[scaled_flux_surface_left]
type = ScalePostprocessor
scaling_factor = '${units 1 m^2 -> mum^2}'
value = flux_surface_left
execute_on = 'initial nonlinear linear timestep_end'
outputs = 'console csv exodus'
[]
[flux_surface_right]
type = SideDiffusiveFluxIntegral
variable = concentration
diffusivity = 'Diffusivity_nonAD'
boundary = 'right'
outputs = none
[]
[scaled_flux_surface_right]
type = ScalePostprocessor
scaling_factor = '${units 1 m^2 -> mum^2}'
value = flux_surface_right
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[max_time_step_size]
type = FunctionValuePostprocessor
function = max_dt_size_function
execute_on = 'initial nonlinear linear timestep_end'
outputs = none
[]
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = NEWTON
petsc_options_iname = '-pc_type -snes_type'
petsc_options_value = 'lu vinewtonrsls'
end_time = ${simulation_time}
line_search = 'none'
automatic_scaling = true
nl_rel_tol = 1e-10
nl_max_its = 34
[TimeStepper]
type = IterationAdaptiveDT
dt = 1.0
iteration_window = 5
optimal_iterations = 26
growth_factor = 1.1
cutback_factor = 0.9
cutback_factor_at_failure = 0.9
timestep_limiting_postprocessor = max_time_step_size
[]
[]
[Outputs]
file_base = 'val-2d_out'
[csv]
type = CSV
start_time = ${outputs_initial_time}
[]
[exodus]
type = Exodus
start_time = ${outputs_initial_time}
output_material_properties = true
time_step_interval = 20
[]
[]
(test/tests/ver-1dc/ver-1dc_base.i)
# This is the base input file for the ver-1dc case.
# This input file is meant to be included within both the ver-1dc.i
# and ver-1dc_mms.i input files to be complete.
# Its purpose is to centralize the capability common to the two cases
# within one file to minimize redundancy and ease maintenance.
# It is not meant to be run on its own.
epsilon_1 = ${units 100 K}
epsilon_2 = ${units 500 K}
epsilon_3 = ${units 800 K}
temperature = ${units 1000 K}
trapping_site_fraction_1 = 0.10 # (-)
trapping_site_fraction_2 = 0.15 # (-)
trapping_site_fraction_3 = 0.20 # (-)
diffusivity = 1 # m^2/s
[Mesh]
type = GeneratedMesh
dim = 1
nx = ${nx_num}
xmax = 1
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile]
[]
[trapped_1]
[]
[trapped_2]
[]
[trapped_3]
[]
[]
[Kernels]
[diff]
type = Diffusion
variable = mobile
extra_vector_tags = ref
[]
[time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[]
[coupled_time_1]
type = CoupledTimeDerivative
variable = mobile
v = trapped_1
extra_vector_tags = ref
[]
[coupled_time_2]
type = CoupledTimeDerivative
variable = mobile
v = trapped_2
extra_vector_tags = ref
[]
[coupled_time_3]
type = CoupledTimeDerivative
variable = mobile
v = trapped_3
extra_vector_tags = ref
[]
[]
[NodalKernels]
# For first traps
[time_1]
type = TimeDerivativeNodalKernel
variable = trapped_1
[]
[trapping_1]
type = TrappingNodalKernel
variable = trapped_1
alpha_t = '${trapping_rate_coefficient}'
N = '${fparse N / cl}'
Ct0 = '${trapping_site_fraction_1}'
mobile_concentration = 'mobile'
temperature = '${temperature}'
extra_vector_tags = ref
[]
[release_1]
type = ReleasingNodalKernel
alpha_r = '${release_rate_coefficient}'
temperature = '${temperature}'
detrapping_energy = '${epsilon_1}'
variable = trapped_1
[]
# For second traps
[time_2]
type = TimeDerivativeNodalKernel
variable = trapped_2
[]
[trapping_2]
type = TrappingNodalKernel
variable = trapped_2
alpha_t = '${trapping_rate_coefficient}'
N = '${fparse N / cl}'
Ct0 = '${trapping_site_fraction_2}'
mobile_concentration = 'mobile'
temperature = '${temperature}'
extra_vector_tags = ref
[]
[release_2]
type = ReleasingNodalKernel
alpha_r = '${release_rate_coefficient}'
temperature = '${temperature}'
detrapping_energy = '${epsilon_2}'
variable = trapped_2
[]
# For third traps
[time_3]
type = TimeDerivativeNodalKernel
variable = trapped_3
[]
[trapping_3]
type = TrappingNodalKernel
variable = trapped_3
alpha_t = '${trapping_rate_coefficient}'
N = '${fparse N / cl}'
Ct0 = '${trapping_site_fraction_3}'
mobile_concentration = 'mobile'
temperature = '${temperature}'
extra_vector_tags = ref
[]
[release_3]
type = ReleasingNodalKernel
alpha_r = '${release_rate_coefficient}'
temperature = '${temperature}'
detrapping_energy = '${epsilon_3}'
variable = trapped_3
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = ${simulation_time}
dtmax = ${time_interval_max}
solve_type = NEWTON
scheme = ${scheme}
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
line_search = 'none'
[]
(test/tests/val-2f/val-2f_trapping_intrinsic.i)
# This input file adds the intrinsic trapping sites for validation case val-2f.
# It is included in val-2f.i
[Variables]
[trapped_intrinsic]
order = FIRST
family = LAGRANGE
[]
[]
[Bounds]
[trapped_intrinsic_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_intrinsic
bound_type = lower
bound_value = '${fparse -1e-20}'
[]
[]
[Kernels]
# trapping intrinsic kernel
[coupled_time_trap_intrinsic]
type = ADCoefCoupledTimeDerivative
variable = deuterium_concentration_W
v = trapped_intrinsic
coef = ${trap_per_free_intrinsic}
[]
[]
[NodalKernels]
[time_intrinsic]
type = TimeDerivativeNodalKernel
variable = trapped_intrinsic
[]
[trapping_intrinsic]
type = TrappingNodalKernel
variable = trapped_intrinsic
mobile_concentration = deuterium_concentration_W
alpha_t = '${trapping_prefactor_intrinsic}'
trapping_energy = '${trapping_energy_intrinsic}'
N = '${tungsten_density}'
Ct0 = '${trapping_site_fraction_intrinsic}'
temperature = 'temperature'
trap_per_free = ${trap_per_free_intrinsic}
[]
[release_intrinsic]
type = ReleasingNodalKernel
variable = trapped_intrinsic
alpha_r = '${detrapping_prefactor_intrinsic}'
detrapping_energy = '${detrapping_energy_intrinsic}'
temperature = 'temperature'
[]
[]
[Postprocessors]
[integral_trapped_concentration_intrinsic]
type = ElementIntegralVariablePostprocessor
variable = trapped_intrinsic
outputs = none
[]
[scaled_trapped_deuterium_intrinsic]
type = ScalePostprocessor
scaling_factor = '${fparse ${trap_per_free_intrinsic} * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_intrinsic
[]
[]
(test/tests/val-2f/val-2f_trapping_2.i)
# This input file adds the trapping sites 2 for validation case val-2f.
# It is included in val-2f.i
[Variables]
[trapped_2]
order = FIRST
family = LAGRANGE
[]
[]
[Bounds]
[trapped_2_lower_bound]
type = ConstantBounds
variable = bounds_dummy
bounded_variable = trapped_2
bound_type = lower
bound_value = '${fparse -1e-20}'
[]
[]
[Kernels]
# trapping 2 kernel
[coupled_time_trap_2]
type = ADCoefCoupledTimeDerivative
variable = deuterium_concentration_W
v = trapped_2
coef = ${trap_per_free_2}
[]
[]
[NodalKernels]
# Second traps
[time_2]
type = TimeDerivativeNodalKernel
variable = trapped_2
[]
[trapping_2]
type = TrappingNodalKernel
variable = trapped_2
mobile_concentration = deuterium_concentration_W
alpha_t = '${trapping_prefactor_2}'
trapping_energy = '${trapping_energy_2}'
N = '${tungsten_density}'
Ct0 = 'trap_distribution_function_2'
temperature = 'temperature'
trap_per_free = ${trap_per_free_2}
[]
[release_2]
type = ReleasingNodalKernel
variable = trapped_2
alpha_r = '${detrapping_prefactor_2}'
detrapping_energy = '${detrapping_energy_2}'
temperature = 'temperature'
[]
[]
[Functions]
[trap_distribution_function_2]
type = ParsedFunction
expression = '${trapping_site_fraction_2} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[trap_distribution_function_2_inf]
type = ParsedFunction
expression = '${trapping_site_fraction_2_inf} / (1 + exp((x - ${depth_center}) / ${depth_width}))'
[]
[]
[Materials]
[trap_distribution_function_2]
type = GenericFunctionMaterial
prop_names = trap_distribution_function_2
prop_values = trap_distribution_function_2
[]
[]
[Postprocessors]
[integral_trapped_concentration_2]
type = ElementIntegralVariablePostprocessor
variable = trapped_2
outputs = none
[]
[scaled_trapped_deuterium_2]
type = ScalePostprocessor
scaling_factor = '${fparse ${trap_per_free_2} * ${units 1 m^2 -> mum^2}}'
value = integral_trapped_concentration_2
[]
[]
(test/tests/ver-1d/ver-1d-trapping_profile.i)
# This input file is slightly adapted from ver-1d-trapping.i to include a spatially-dependent trap site density
cl=3.1622e18
trap_per_free=1e3
N=3.1622e22
time_scaling=1
epsilon=10000
temperature = 1000
Ct0_surface = .1 # trapping site density close to the sample surface
Ct0_bulk = .01 # trapping site density in the sample bulk
trap_profile_depth = 0.2 # position of the transition from Ct0_surface to Ct0_bulk
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 1
[]
[Functions]
[trapping_sites_density_function]
type = ParsedFunction
expression = 'if(x<${trap_profile_depth}, ${Ct0_surface}, ${Ct0_bulk})'
execute_on = 'INITIAL LINEAR NONLINEAR TIMESTEP_END'
[]
[]
[Problem]
type = ReferenceResidualProblem
extra_tag_vectors = 'ref'
reference_vector = 'ref'
[]
[Variables]
[mobile][]
[trapped][]
[]
[AuxVariables]
[empty_sites][]
[scaled_empty_sites][]
[trapped_sites][]
[total_sites][]
[]
[AuxKernels]
[empty_sites]
variable = empty_sites
type = EmptySitesAux
N = ${fparse N / cl}
Ct0 = trapping_sites_density_function
trap_per_free = ${trap_per_free}
trapped_concentration_variables = trapped
[]
[scaled_empty]
variable = scaled_empty_sites
type = NormalizationAux
normal_factor = ${cl}
source_variable = empty_sites
[]
[trapped_sites]
variable = trapped_sites
type = NormalizationAux
normal_factor = ${trap_per_free}
source_variable = trapped
[]
[total_sites]
variable = total_sites
type = ParsedAux
expression = 'trapped_sites + empty_sites'
coupled_variables = 'trapped_sites empty_sites'
[]
[]
[Kernels]
[./diff]
type = MatDiffusion
variable = mobile
diffusivity = ${fparse 1. / time_scaling}
extra_vector_tags = ref
[../]
[./time]
type = TimeDerivative
variable = mobile
extra_vector_tags = ref
[../]
[coupled_time]
type = ScaledCoupledTimeDerivative
variable = mobile
v = trapped
factor = ${trap_per_free}
extra_vector_tags = ref
[]
[]
[NodalKernels]
[time]
type = TimeDerivativeNodalKernel
variable = trapped
[]
[trapping]
type = TrappingNodalKernel
variable = trapped
alpha_t = ${fparse 1e15 / time_scaling}
N = ${fparse 3.1622e22 / cl}
Ct0 = trapping_sites_density_function
mobile_concentration = 'mobile'
temperature = ${temperature}
trap_per_free = ${trap_per_free}
extra_vector_tags = ref
[]
[release]
type = ReleasingNodalKernel
alpha_r = ${fparse 1e13 / time_scaling}
temperature = ${temperature}
detrapping_energy = ${epsilon}
variable = trapped
[]
[]
[BCs]
[left]
type = DirichletBC
variable = mobile
value = ${fparse 3.1622e18 / cl}
boundary = left
[]
[right]
type = DirichletBC
variable = mobile
value = 0
boundary = right
[]
[]
[Postprocessors]
[outflux]
type = SideDiffusiveFluxAverage
boundary = 'right'
diffusivity = 1
variable = mobile
[]
[scaled_outflux]
type = ScalePostprocessor
value = outflux
scaling_factor = ${cl}
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
end_time = 1000
dt = 1
solve_type = NEWTON
line_search = 'none'
automatic_scaling = true
verbose = true
compute_scaling_once = false
l_tol = 1e-11
[]
[Outputs]
exodus = true
csv = true
[dof]
type = DOFMap
execute_on = initial
[]
perf_graph = true
[]