Divertor Monoblock During Pulsed Operation

Contact: Pierre-Clement Simon (pierreclement.simon.at.inl.gov), Masashi Shimada (masashi.shimada.at.inl.gov)

TMAP8 is used to model tritium transport in a divertor monoblock to elucidate the effects of pulsed operation (up to fifty 1600-second plasma discharge and cool-down cycles) on the tritium in-vessel inventory source term and ex-vessel release term (i.e., tritium retention and permeation) for safety analysis. This example reproduces the results presented in Shimada et al. (2024).

General description of the simulation case and corresponding input file

Introduction

In a magnetic confinement fusion system, the divertor components will be subjected to intense particle and heat fluxes from plasma, as well as 14-MeV neutrons stemming from deuterium-tritium (D-T) reactions, thus creating extremely high temperature and hydrogen concentration gradients across tens of millimeters. Furthermore, the thermomechanical properties of the divertor materials will evolve over time in response to displacement damage and gas/solid transmutations. Thus, predictive modeling of tritium transport in the divertor poses enormous challenges, as it requires simulation capabilities capable of (1) multi-material configurations, (2) high thermal cycles, (3) complex 2D and 3D geometries, (4) and microstructural evolution resulting from displacement damage and gas/solid transmutations. For example, the ITER divertor, which consists of tungsten (W) monoblocks bonded to cooling tubes made of copper-chromium-zirconium (CuCrZr) alloy, is designed to withstand high heat fluxes (∼10 MW $var element = document.getElementById("moose-equation-585c4760-094f-46c1-9204-26171630a532");katex.render("m^{\u22122}", element, {displayMode:false,throwOnError:false});$ ) during steady-state plasma operation from intense D-T plasma flux (∼1024 $var element = document.getElementById("moose-equation-b1f57c6f-9b35-43ee-9199-9372410b085e");katex.render("m^{\u22122}s^{\u22121}", element, {displayMode:false,throwOnError:false});$ ). In this example, we simulate tritium retention and permeation as well as thermal transport in an ITER-like 2D W-Cu-CuCrZr monoblock and demonstrate the first three (out of four) aforementioned required simulation capabilities for tritium transport in the divertor using TMAP8, as presented in Shimada et al. (2024). This simulation was also based on the cases published in Hodille et al. (2021).

The sections below describe the simulation details and explain how they translate into the example TMAP8 input file. Not every part of the input file is explained here. If the interested reader has more questions about this example case and how to modify this input file to adapt it to a different case, feel free to reach out to the TMAP8 development team on the TMAP8 GitHub discussion page.

Divertor monoblock geometry and mesh generation

Figure 1 shows the geometry, mesh, and temperature distribution of the 2D monoblock employed in Shimada et al. (2024) and in the present example. The 2D monoblock consisted of three materials:

W monoblock ( $var element = document.getElementById("moose-equation-9b7eb9a5-a045-496e-ac9e-6e869c0e2f22");katex.render("r", element, {displayMode:false,throwOnError:false});$ (mm) > 8.5),
Cu interlayer (7.5 < $var element = document.getElementById("moose-equation-c4b4f460-fda1-4c0f-9309-e0f8fb8a017d");katex.render("r", element, {displayMode:false,throwOnError:false});$ (mm) < 8.5),
CuCrZr tube (6.0 < $var element = document.getElementById("moose-equation-eddbee63-de16-46ee-9f98-7c3200d51c2e");katex.render("r", element, {displayMode:false,throwOnError:false});$ (mm) < 7.5), and
H $var element = document.getElementById("moose-equation-dfec62ad-421f-467b-a330-f5b4bbfab6ed");katex.render("_2", element, {displayMode:false,throwOnError:false});$ O cooling ( $var element = document.getElementById("moose-equation-95c41dbd-5c54-4e24-b6e7-d0afbd807755");katex.render("r", element, {displayMode:false,throwOnError:false});$ (mm) < 6.0),

where $var element = document.getElementById("moose-equation-9e1e8377-9595-4c57-a856-20a38b53c28a");katex.render("r=\\sqrt{x^2+y^2}", element, {displayMode:false,throwOnError:false});$ in a Cartesian coordinate system, and the center point ( $var element = document.getElementById("moose-equation-dfb972a1-0c92-484e-b286-2089ec2c2f99");katex.render("x", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-ca1d7826-9c9b-48e4-988a-7c3a924a581b");katex.render("y", element, {displayMode:false,throwOnError:false});$ ) = (0, 0) is defined at the center of the CuCrZr tube.

Figure 1: 2D monoblock: (left) geometry and mesh; (right) temperature distribution. This corresponds to Fig. 1 in Ref. Shimada et al. (2024).

MOOSE is equipped with a set of user-friendly, built-in mesh generators for creating meshes based on simple geometries (e.g., a monoblock). The built-in mesh generator ConcentricCircleMeshGenerator is used to create meshes for the 2D monoblock geometry. Only the left-hand (−14.0 < $var element = document.getElementById("moose-equation-31592b50-5ef9-4abf-ac41-ce1d43f6104c");katex.render("x", element, {displayMode:false,throwOnError:false});$ (mm) < 0 & −14.0 < $var element = document.getElementById("moose-equation-9fa68ebd-4970-4897-8da7-0a68a0b8527c");katex.render("y", element, {displayMode:false,throwOnError:false});$ (mm) < 14.0) portion of the monoblock is used, such that the mesh size is halved by assuming symmetry at the $var element = document.getElementById("moose-equation-de024195-a766-4f3b-9ecd-44d82c565bc1");katex.render("x = 0", element, {displayMode:false,throwOnError:false});$ plane. This example uses a total of 9,418 nodes, creating 28,402 degrees of freedom in the non-linear system. The finer mesh size and time step are required to compute the trapping-limited diffusion regime, as described in TMAP8 verification case ver-1d. However, a set of relatively low detrapping energies (< 0.85 eV) is used here to reduce the computing time via the relatively coarse mesh size.

In the input file, the mesh is defined as:

[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
  []
[]

Nomenclature of variables and physical parameters

Table 1 lists the variables and physical parameters used in this example with their units.

Table 1: Nomenclature of variables and physical parameters used in this example.

Variables and governing equations

To simulate tritium and thermal transport, we define two sets of partial differential equations (PDEs). First, the strong form of the mass conservation equation for solute (mobile) T atoms, $var element = document.getElementById("moose-equation-ceb0effc-baf3-40e2-836a-ecd8268164d6");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ , is written as:

$(1)var element = document.getElementById("moose-equation-e50196a2-a3b0-4e6d-8385-e7837f85ef0a");katex.render(" \\frac{\\partial C_s}{\\partial t} + \\nabla \\cdot (-D \\nabla C_s) + \\alpha_t \\frac{C_t^e}{N} C_s - \\alpha_r C_t = 0.", element, {displayMode:true,throwOnError:false});$

We use three sets of mass conservation equations to calculate the behaviors of solute T atoms in three different materials (i.e., W, Cu, CuCrZr). Second, the strong form of the conservation of energy equation is written as:

$(2)var element = document.getElementById("moose-equation-87bb8dad-76d5-4007-9014-24bfd249e507");katex.render(" \\rho c_p \\frac{\\partial T}{\\partial t} - \\nabla \\cdot (k_T \\nabla T) = 0.", element, {displayMode:true,throwOnError:false});$

Table 1 lists all the symbols and definitions used in Eq. (1) and Eq. (2).

The next step is to convert these two strong-form PDEs into their weak forms by multiplying with a test function, $var element = document.getElementById("moose-equation-8315f75b-3073-464c-99e5-39a1524f35e7");katex.render("\\psi", element, {displayMode:false,throwOnError:false});$ , and integrating over a domain, $var element = document.getElementById("moose-equation-5c24dc04-006e-4701-bf28-64992e3e9fe7");katex.render("\\Omega", element, {displayMode:false,throwOnError:false});$ with surface $var element = document.getElementById("moose-equation-82ccc58d-eeb0-46a3-a75b-f454be00d70a");katex.render("\\delta \\Omega", element, {displayMode:false,throwOnError:false});$ and outward-facing normal vector $var element = document.getElementById("moose-equation-504a9e9e-08b4-49f7-aabf-5d335f58740e");katex.render("\\hat{n}", element, {displayMode:false,throwOnError:false});$ . Using the divergence theorem, one can obtain the weak form of the mass conservation equation (Eq. (1)) as follows:

$(3)var element = document.getElementById("moose-equation-aad808ae-73fd-432a-a0c6-f62ce39f21d1");katex.render(" \\int_{\\Omega}\\psi \\frac{\\partial C_s}{\\partial t} -\\int_{\\delta \\Omega} \\psi D \\nabla C_s \\cdot \\hat{n} +\\int_{\\Omega} \\nabla \\psi \\cdot D \\nabla C_s +\\int_{\\Omega} \\psi \\alpha_t \\frac{C_t^e}{N} C_s -\\int_{\\Omega} \\psi \\alpha_r C_t = 0.", element, {displayMode:true,throwOnError:false});$

Similarly, the weak form of the conservation of energy equation can be written by multiplying Eq. (2) with a test function, $var element = document.getElementById("moose-equation-48986f00-89d0-438b-807f-69c251f0821d");katex.render("\\psi", element, {displayMode:false,throwOnError:false});$ , and integrating over a domain, $var element = document.getElementById("moose-equation-c1caf253-9b3e-41e7-98ff-bd5fb18e460d");katex.render("\\Omega", element, {displayMode:false,throwOnError:false});$ with surface $var element = document.getElementById("moose-equation-3aa44c30-b6a1-4a7f-8ca7-5122093929e1");katex.render("\\delta \\Omega", element, {displayMode:false,throwOnError:false});$ and outward-facing normal vector $var element = document.getElementById("moose-equation-136c885f-bc94-4744-aa2c-593caeb97a82");katex.render("\\hat{n}", element, {displayMode:false,throwOnError:false});$ . Using the divergence theorem, one can obtain the weak form of the conservation of energy equation as follows:

$(4)var element = document.getElementById("moose-equation-adc43435-4026-4c2b-98a1-bef528e40f45");katex.render(" \\int_{\\Omega}\\psi c_p \\frac{\\partial T}{\\partial t} -\\int_{\\delta \\Omega} \\psi k_T \\nabla T \\cdot \\hat{n} +\\int_{\\Omega} \\nabla \\psi \\cdot k_T \\nabla T = 0.", element, {displayMode:true,throwOnError:false});$

Then, to solve for the PDEs and physical phenomena, we can select appropriate kernels and boundary conditions (BCs) from MOOSE’s extensive library. The following three subsections describe each kernel, BC, and numerical method utilized in the present work.

In the input file, the variables are defined as:

[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
  []
[]

Note the usage of initial_condition and block parameters in order to set initial conditions and material block applicability for each variable.

Kernels

For the mass conservation equation, we use ADTimeDerivative and ADMatDiffusion to represent the 1 $var element = document.getElementById("moose-equation-56d3a27a-a960-4b36-be44-0c429c0f157a");katex.render("^{\\text{st}}", element, {displayMode:false,throwOnError:false});$ and 3 $var element = document.getElementById("moose-equation-196fd126-4329-4c67-a01a-d51160aa2387");katex.render("^{\\text{rd}}", element, {displayMode:false,throwOnError:false});$ terms of Eq. (3). The TMAP8 kernels TrappingNodalKernel and ReleasingNodalKernel represent the 4 $var element = document.getElementById("moose-equation-7cdb9c42-6dc7-4432-9f17-42857227d847");katex.render("^{\\text{th}}", element, {displayMode:false,throwOnError:false});$ and 5 $var element = document.getElementById("moose-equation-0e2bd974-df17-4feb-b572-7bb50b6fc9cd");katex.render("^{\\text{th}}", element, {displayMode:false,throwOnError:false});$ terms of Eq. (3), respectively, and simulate the trapping/release behavior of hydrogen isotopes in/from trap sites. For the conservation of energy equation, SpecificHeatConductionTimeDerivative and HeatConduction solve the 1 $var element = document.getElementById("moose-equation-eb364e82-3108-4a47-9c2e-621fe4835ec1");katex.render("^{\\text{st}}", element, {displayMode:false,throwOnError:false});$ and 3 $var element = document.getElementById("moose-equation-98ba4e67-9434-46fc-b83d-c67a72f32fae");katex.render("^{\\text{rd}}", element, {displayMode:false,throwOnError:false});$ terms of Eq. (4). ADTimeDerivative and ADMatDiffusion are kernels from the core MOOSE framework, and SpecificHeatConductionTimeDerivative and HeatConduction are kernels from the MOOSE Heat Transfer Module. These pre-made kernels are commonly re-used to represent time derivative, diffusion, and heat conduction terms in a given material within a variety of MOOSE-based simulations.

[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
  []
[]

NodalKernels are used for non-diffusive variables (trapped species).

[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
  []
[]

AuxVariables and AuxKernels

AuxVariables are used to track quantities that are not solved for by the PDEs, but are needed to compute materials properties, or are desirable to obtain as outputs, such as the total concentration of tritium or flux values. AuxKernels provide the expressions that define the AuxVariables.

[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
  []
[]

[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
  []
[]

Boundary conditions and history

For the mass conservation equation, we use FunctionNeumannBC and DirichletBC to solve the 2 $var element = document.getElementById("moose-equation-c0f3ff53-316e-4802-897c-2dcddc7f66bb");katex.render("^{\\text{nd}}", element, {displayMode:false,throwOnError:false});$ term of Eq. (3). FunctionNeumannBC is used to treat the time-dependent plasma exposure at the top (plasma-exposed) surface, and DirichletBC is used to set the BC of the solute T atom concentration to zero at the inner CuCrZr tube (at $var element = document.getElementById("moose-equation-6a451f21-7463-4372-a759-79f5bf8c454f");katex.render("r = 6", element, {displayMode:false,throwOnError:false});$ mm). For the energy conservation equation, we use FunctionNeumannBC and DirichletBC to solve the 2 $var element = document.getElementById("moose-equation-4a2cc95f-6029-4241-b278-eb20c554a306");katex.render("^{\\text{nd}}", element, {displayMode:false,throwOnError:false});$ term of Eq. (4). FunctionNeumannBC is used to treat the time-dependent heat flux at the top (plasma-exposed) surface, and FunctionDirichletBC is used to treat the temperature increase in the cooling tube. FunctionNeumannBC, FunctionDirichletBC, DirichletBC, and FunctionDirichletBC are MOOSE objects commonly used to represent the BCs of the variables to be solved.

We simulate a 20,000-second plasma discharge, with each 1,600-second cycle consisting of a 100-second plasma ramp-up, a 400-second steady-state plasma discharge, a 100-second plasma ramp-down, and a 1,000-second waiting phase. Up to 50 cycles are simulated to achieve the total discharge. Figure 2 shows the integrated (solute, total and trapped) tritium concentration profiles in the monoblock. It shows that implantation fluxes and temperatures vary linearly up to (from) their steady-state values from (up to) their initial values during ramp-up (ramp-down).

Figure 2: Temperature profiles (orange) and integrated tritium concentration profiles (blue) during two 1,600-second-cycle plasma discharges. This corresponds to Fig. 2 in Shimada et al. (2024).

During the steady-state plasma discharge, we set a heat flux of 10 MW $var element = document.getElementById("moose-equation-5b937597-1acc-4127-8468-2cc0a43b996a");katex.render("\\cdot", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-6f6b3796-6067-4c24-849c-5acc21a93d08");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ at the top of the 2D monoblock (at $var element = document.getElementById("moose-equation-9a27d9d6-4c46-4f33-81a5-5f65788bfb20");katex.render("y = 14.0", element, {displayMode:false,throwOnError:false});$ mm) and a cooling temperature of 552 K at the inner CuCrZr tube (at $var element = document.getElementById("moose-equation-ae1817fc-9e27-49ca-8073-bbebe1a862be");katex.render("r = 6.0", element, {displayMode:false,throwOnError:false});$ mm). We assume a 100% T plasma with a 5.0 $var element = document.getElementById("moose-equation-d8149ad8-7b0b-4022-929b-cd31386a6852");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-3a29c97a-b765-42d4-a169-3564e0ce5c6e");katex.render("^{23}", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-c47adb48-6f0f-4a6c-9a79-1ccc384596d3");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-cd3dfbbc-7d36-41ba-b4f4-5ce43f91378a");katex.render("\\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-9e9d8186-a549-40d0-b1de-081e2184622d");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ plasma particle flux (which is half of the full 1.0 $var element = document.getElementById("moose-equation-9032ce33-9f66-485d-829e-16dcffb04d11");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-3ae89f18-c69e-4a2e-856f-4dc1c6f4297e");katex.render("^{24}", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-dc7457e3-2c8c-48af-929b-3d789a464ae6");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-54cc1b26-6bf4-4c8c-bc9c-aa7efcf4cc40");katex.render("\\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-4c70568a-211e-46e1-8ee3-d11c01904547");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ DT plasma particle flux), and only 0.1% of the incident plasma particle flux (5.0 $var element = document.getElementById("moose-equation-87fdf55c-089e-4f7f-adda-97334eb0e568");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-30c160ff-a059-441b-84d3-05268d895e51");katex.render("^{20}", element, {displayMode:false,throwOnError:false});$ m $var element = document.getElementById("moose-equation-7fe9814b-90f3-4cbd-9330-47ce7a5f333a");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-54f8ef3f-eabf-43b6-a82d-79b21264816d");katex.render("\\cdot", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-61f2e38e-ce46-4e4a-8743-afea1e1ee222");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ) entered the first layer of mesh at the exposed surface ( $var element = document.getElementById("moose-equation-eaa9f0e6-55a4-44a3-aaf8-c77cceef5bfe");katex.render("y = 14.0", element, {displayMode:false,throwOnError:false});$ mm) as in Shimada et al. (2024). The solute T atom concentration is set to zero at the inner CuCrZr tube (at $var element = document.getElementById("moose-equation-99e00a3f-1f7d-4970-9b7d-f2eb9277884d");katex.render("r = 6.0", element, {displayMode:false,throwOnError:false});$ mm).

We treated this plasma exposure by setting the flux BC of the solute T atom concentration at the exposed surface ( $var element = document.getElementById("moose-equation-5e1512be-b947-40d3-9e50-5cb8f16f50a9");katex.render("y = 14.0", element, {displayMode:false,throwOnError:false});$ mm) as a simplification of the complex plasma implantation and recombination phenomena, which would require a very fine mesh and increase computational costs.

[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

In this example, Functions are used to define the time-dependent tritium and heat flux at the exposed surface, as well as the coolant temperature. Functions can also be spatially-dependent. The functions are then provided to BCs, as described above.

[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)))))))))'
  []
[]

Material properties

We use the tritium mass transport properties listed in Table 2 to solve the mass conservation equation (Eq. (3)) for two variables: $var element = document.getElementById("moose-equation-a94dca05-42b0-42aa-8c45-bd65b38188b0");katex.render("C_s", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-0a10fd55-e040-4b7a-a6af-5c5883261602");katex.render("C_t", element, {displayMode:false,throwOnError:false});$ (i.e., solute and trapped T atom concentrations) in three materials and the thermal properties (i.e., density, temperature-dependent specific heat, and thermal conductivity) listed in Table 3 to solve the energy conservation equation (Eq. (4)) for one variable: the temperature $var element = document.getElementById("moose-equation-1aba7aaf-a610-437b-a82c-d490ff8fe1c8");katex.render("T", element, {displayMode:false,throwOnError:false});$ . The diffusivity is defined as $var element = document.getElementById("moose-equation-366c35c9-2d77-4628-b45e-2222e11c556c");katex.render("D=D_0 \\exp\u2061(-E_D/k_B/T)", element, {displayMode:false,throwOnError:false});$ and the solubility is defined as $var element = document.getElementById("moose-equation-05bcb767-529c-40f2-9854-2fad16dcd75b");katex.render("K_s=K_{s,0} \\exp\u2061(-E_s/k_B/T)", element, {displayMode:false,throwOnError:false});$ .

Table 2: Tritium mass transport and trapping properties used in the W, Cu, and CuCrZr layers. $var element = document.getElementById("moose-equation-41f8b0d5-29e8-440e-a4c9-a865ed20eced");katex.render("^{**}", element, {displayMode:false,throwOnError:false});$ NOTE: Two-component solubility in W kept the maximum solubility ratio between W and Cu to 104 at low temperature. The references are provided in Shimada et al. (2024).

Material	$var element = document.getElementById("moose-equation-433dea20-5f55-4630-aaf9-0f989398a364");katex.render("D_0", element, {displayMode:false,throwOnError:false});$ (m $var element = document.getElementById("moose-equation-66a8ade1-5d48-4250-afc7-f7258699e533");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s)	$var element = document.getElementById("moose-equation-f6759740-ea57-49ea-ba44-ef6ee52dbf05");katex.render("E_D", element, {displayMode:false,throwOnError:false});$ (eV)	$var element = document.getElementById("moose-equation-0fb0b5d5-5f5a-4235-9ef8-e0e3827b57cb");katex.render("K_{s,0}", element, {displayMode:false,throwOnError:false});$ (Pa $var element = document.getElementById("moose-equation-b3e3645f-e6f7-4b70-8b0a-d41a08526911");katex.render("^{1/2}", element, {displayMode:false,throwOnError:false});$ )	$var element = document.getElementById("moose-equation-415db12a-2618-485a-b774-9ae899821052");katex.render("E_s", element, {displayMode:false,throwOnError:false});$ (eV)	Detrapping energy: $var element = document.getElementById("moose-equation-c8fe6102-c8fa-4175-8304-c3abf82d80dc");katex.render("E_{dt}", element, {displayMode:false,throwOnError:false});$ (eV)	Trap density: $var element = document.getElementById("moose-equation-01535738-aff8-46e9-88d6-d097724763e9");katex.render("n_{trap}", element, {displayMode:false,throwOnError:false});$ (at.fr.)
W	2.4 $var element = document.getElementById("moose-equation-f7880ebe-aa1a-4792-9e77-19a0f829a277");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-05f7431b-49b8-4cf1-a85d-a63fe62d9ffd");katex.render("^{-7}", element, {displayMode:false,throwOnError:false});$	0.39	1.87 $var element = document.getElementById("moose-equation-38b284c9-5a63-41bc-961c-97fa43ea462e");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-00d071f6-1f80-4e3c-904f-85706bbac60b");katex.render("^{24}", element, {displayMode:false,throwOnError:false});$	1.04	0.85	1.0 $var element = document.getElementById("moose-equation-5020ee9a-b31b-4bf7-86ee-fb1dcba0ab31");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-82c34602-005f-4cfe-952d-d67980e2273f");katex.render("^{-4}", element, {displayMode:false,throwOnError:false});$
W			3.14 $var element = document.getElementById("moose-equation-364ca606-f4ea-43d3-a027-f0ff2f4e7471");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-0d15735b-6651-43be-bd98-6fb56d3fb1a6");katex.render("^{20}", element, {displayMode:false,throwOnError:false});$	0.57 $var element = document.getElementById("moose-equation-0360fd36-4dee-4f9c-a44c-bddd47cf02c1");katex.render("^{**}", element, {displayMode:false,throwOnError:false});$
Cu	6.6 $var element = document.getElementById("moose-equation-68d73c6f-19d8-4a49-b67a-a0b4e2291c6a");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-49788f3f-bf10-409e-9a1c-9bee5eee2926");katex.render("^{-7}", element, {displayMode:false,throwOnError:false});$	0.39	3.14 $var element = document.getElementById("moose-equation-8e2c511f-e63e-47cb-901b-7a66b537990b");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-4b9d4ecf-229e-46f7-9397-fdc938c57c5c");katex.render("^{24}", element, {displayMode:false,throwOnError:false});$	0.57	0.50	5.0 $var element = document.getElementById("moose-equation-1b71421a-6d5b-49d8-b430-abad327eff8d");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-99815ab0-5d73-43f6-8d70-460cd246ae73");katex.render("^{-5}", element, {displayMode:false,throwOnError:false});$
CuCrZr	3.9 $var element = document.getElementById("moose-equation-3c64f4f6-162b-4374-ad4e-6d1045a86735");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-01566592-d2b2-44c2-83a6-6e763641bb1b");katex.render("^{-7}", element, {displayMode:false,throwOnError:false});$	0.42	4.28 $var element = document.getElementById("moose-equation-2feacaea-3d6a-4dd5-9be0-b274802e00ce");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-739ae8fc-3ad1-49a2-aa74-d335fa82c6a3");katex.render("^{23}", element, {displayMode:false,throwOnError:false});$	0.39	0.83	5.0 $var element = document.getElementById("moose-equation-adb7f82e-73ac-43a3-9bdb-d1368973e825");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-c5e18d78-ec19-45f9-bf24-d91c8646f1f8");katex.render("^{-5}", element, {displayMode:false,throwOnError:false});$

Table 3: Thermal properties used in the W, Cu, and CuCrZr layers. The references are provided in Shimada et al. (2024).

Material	Density: $var element = document.getElementById("moose-equation-61eaa009-0720-49bb-8f79-ab948f4675a7");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ (g m $var element = document.getElementById("moose-equation-81552d67-4762-4509-ba3c-cf733b19d83f");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$ )	Specific heat: $var element = document.getElementById("moose-equation-bbbca982-15ef-4cab-98cf-41c3b5449c13");katex.render("c_p", element, {displayMode:false,throwOnError:false});$ (J kg $var element = document.getElementById("moose-equation-29836dab-0a18-49d1-a976-4782ce9bc36b");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ K $var element = document.getElementById("moose-equation-f22ed72e-22fc-412f-9bc7-e455d7bc795a");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ )	Thermal conductivity: $var element = document.getElementById("moose-equation-28fecd3a-fc64-436e-bd1f-6b4a3d2d24e1");katex.render("k_T", element, {displayMode:false,throwOnError:false});$ (W m $var element = document.getElementById("moose-equation-92807018-ceea-4204-8711-1c634a16ee91");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ K $var element = document.getElementById("moose-equation-27d5ff88-f7b4-4596-99d1-bf69bd6caa73");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ )
W	19,300	1.16 $var element = document.getElementById("moose-equation-e21d34af-4c8c-43ac-ab82-ed8276f3431d");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-018537ab-980b-4ff4-9f6c-f50b269ed27c");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ +7.11 $var element = document.getElementById("moose-equation-04835872-a30c-493c-9a94-f084a11df055");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-baba1bed-ea37-42aa-a54b-2c54058245b2");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ T –6.58 $var element = document.getElementById("moose-equation-5aff6a82-f550-4737-a37f-b2093b9a7234");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-95b6a8c3-2958-418c-9ff0-5882d8176410");katex.render("^{-5}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-3ca4562a-402b-4d50-b79b-85da27c4bd14");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ +3.24 $var element = document.getElementById("moose-equation-66f03792-eddd-4825-9d97-8e1cf74ce563");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-8c29d7a5-d157-4f06-b140-1221eea53623");katex.render("^{-8}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-ccc9f959-13e4-4708-85bc-f1b371f31d02");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ –5.45 $var element = document.getElementById("moose-equation-bf5dd14d-f07d-4476-b280-d7eb875b4cc3");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-27b6124a-9a6f-47a3-8801-0a81ef9aefd6");katex.render("^{-12}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-5e9e697e-454e-4281-bb48-24af278ae4f9");katex.render("^{4}", element, {displayMode:false,throwOnError:false});$ (293 < T (K) < 2500)	2.41 $var element = document.getElementById("moose-equation-816bdfe5-ae03-4ee0-bb43-673f72c0502c");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-d4fc7fa7-5a7f-49a6-b96e-2ff94e5fda6e");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ –2.90 $var element = document.getElementById("moose-equation-d3ea2f4b-a606-4f01-ad13-062254ee3061");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-c7862e39-bddf-4347-9529-02a3353d9038");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ T + 2.54 $var element = document.getElementById("moose-equation-6c9f4105-ef95-4d4e-889f-73093fac4d2e");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-556867bc-2b16-4a23-81ca-9d80a96daef1");katex.render("^{-4}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-a5c8e735-4ffd-4dcb-8b53-65268ee9b453");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ –1.03 $var element = document.getElementById("moose-equation-b1e1d07f-079d-4206-b19f-9dc41fa40310");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-c1bcc997-1c0a-4468-8c33-7c0782dc0cc6");katex.render("^{-7}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-a5c30659-acb1-4fda-aba7-8a1478e293f2");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ +1.52 $var element = document.getElementById("moose-equation-f8dce809-8264-4036-a08d-8ab0ee832c52");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-dfbd0e94-f172-4b29-81d8-ff21a0852834");katex.render("^{-11}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-c8784603-71b2-4a32-aa0a-47e4689fa3dc");katex.render("^{4}", element, {displayMode:false,throwOnError:false});$ (293 < T (K) < 2500)
Cu	8,960	4.21 $var element = document.getElementById("moose-equation-2c62b3f2-407f-426c-a043-5bba267ce03c");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-4c861e04-0ccc-4ef6-829d-551a00cd1763");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ –6.85 $var element = document.getElementById("moose-equation-388ea8fb-b18b-4375-81c6-c1b67094ab1a");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-7eac5726-d6db-498f-a57f-46253c002a8c");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ T (293 < T (K) < 873)	3.16 $var element = document.getElementById("moose-equation-1afc75fe-7d6c-4c86-8b29-a722fe09d050");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-5dc02e5b-0c49-4512-ac77-e152df5c38a2");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ +3.18 $var element = document.getElementById("moose-equation-565ce86c-4775-4e82-b44a-c3eafcaabb94");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-29e2fb65-f4aa-4d4e-8e58-d565a9c4a4af");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ T –3.49 $var element = document.getElementById("moose-equation-16e2e014-e390-47b9-a6bd-5e310dd47d6c");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-bb5c9968-5864-459e-ae6a-19696bb1f797");katex.render("^{-4}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-0a963f15-36e4-4fff-8279-4afb875a6fa7");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ +1.66 $var element = document.getElementById("moose-equation-3c12e799-c487-40dc-9442-206b3bf76d85");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-95a435d0-10da-47a1-b2e3-1922b191fb81");katex.render("^{-7}", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-cf5eb4ae-0e45-4623-84d5-7c7a2360edac");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ (293 < T (K) < 873)
CuCrZr	8,900	390	3.87 $var element = document.getElementById("moose-equation-01ef0ea2-e6c9-4160-9e79-752210bd489d");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-165bb4d2-807a-42cd-aec1-7d54e008210f");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ –1.28 $var element = document.getElementById("moose-equation-a1113d1f-4789-4b3a-9dce-b1f272045acb");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-224377bc-814f-4fe7-baaf-1f3f7a0faf2d");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ T (293 < T (K) < 927)

[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
  []
[]

Interface Kernels

TMAP8 assumes equilibrium between the chemical potentials of the diffusing species similar to how TMAP4 and TMAP7 treated the diffusing species across two different materials with different chemical potentials. SolubilityRatioMaterial` is used to treat solute concentration differences stemming from a difference in T solubilities across the interface. The solubility ratio jump is calculated via the following: $(5)var element = document.getElementById("moose-equation-f9b119e8-782a-4581-9acb-7516ef343c76");katex.render(" J = \\frac{C_{s,1}}{K_{s,1}} - \\frac{C_{s,2}}{K_{s,2}}.", element, {displayMode:true,throwOnError:false});$

ADPenaltyInterfaceDiffusion is used to conserve the particle flux at this interface between two different solubilities. The extremely low tritium solubility in W at low temperature leads to an extremely low solute concentration in W, creating a large solute concentration difference between W and Cu. Two-component solubility in W is used to keep the maximum solubility ratio between W and Cu to 104 at low temperature to avoid the convergence issue associated with calculating two significantly different solute concentration in W and Cu.

[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'
  []
[]

Numerical method

We use a standard preconditioner: the “single matrix preconditioner.” The Newton method is used to model the transient tritium and thermal transport in a 2D monoblock. It is important to note that MOOSE is equipped with the built-in Message Passing Interface (MPI) protocol, as tritium and thermal transport analysis of fifty 1,600-second cycle plasma discharges in the 2D monoblock is performed in under 2 hours using a single device/computer (3.5 GHz Apple M2 Pro, 10-Core CPU/16-Core GPU) with this MPI feature.

[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
  []
[]

Results

The simulation results from this example are shown in Figure 3 and Figure 4. For more results, information, and discussion about the results for this example case and their significance, the reader is referred to Shimada et al. (2024).

Figure 3: Tritium concentration profile in W (left), Cu (center), and CuCrZr (right) after ten 1,600-second cycles ( $var element = document.getElementById("moose-equation-39459795-cec7-4fde-8dea-5aaa9a2e5eea");katex.render("t = 14912", element, {displayMode:false,throwOnError:false});$ s). This corresponds to Fig. 4A in Shimada et al. (2024).

Figure 4: Tritium concentration profile in W (left), Cu (center), and CuCrZr (right) after fifty 1,600-second cycles ( $var element = document.getElementById("moose-equation-eb0f166a-8be2-427d-a016-1d03c4379c5b");katex.render("t = 78912", element, {displayMode:false,throwOnError:false});$ s). This corresponds to Fig. 4B in Shimada et al. (2024).

Complete input file

Below is the complete input file, which can be run reliably with approximately 4 processor cores. Note that this input file has not been optimized for computational costs.

## Tokamak divertor monoblock
## Application: TMAP8
## POC: Masashi Shimada (masashi.shimada at inl.gov)
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

### 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
[]

References

E. A. Hodille, R. Delaporte-Mathurin, J. Denis, M. Pecovnik, E. Bernard, Y. Ferro, R. Sakamoto, Y. Charles, J. Mougenot, A. De Backer, C. S. Becquart, S. Markelj, and C. Grisolia. Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under iter-like conditions. Nuclear Fusion, 61:126003, 10 2021. URL: https://iopscience.iop.org/article/10.1088/1741-4326/ac2abc, doi:10.1088/1741-4326/AC2ABC.

@article{Hodille2021126003,
    author = "Hodille, E. A. and Delaporte-Mathurin, R. and Denis, J. and Pecovnik, M. and Bernard, E. and Ferro, Y. and Sakamoto, R. and Charles, Y. and Mougenot, J. and Backer, A. De and Becquart, C. S. and Markelj, S. and Grisolia, C.",
    doi = "10.1088/1741-4326/AC2ABC",
    issn = "0029-5515",
    issue = "12",
    journal = "Nuclear Fusion",
    month = "10",
    pages = "126003",
    publisher = "IOP Publishing",
    title = "Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under ITER-like conditions",
    volume = "61",
    url = "https://iopscience.iop.org/article/10.1088/1741-4326/ac2abc",
    year = "2021"
}

Masashi Shimada, Pierre-Clément A. Simon, Casey T. Icenhour, and Gyanender Singh. Toward a high-fidelity tritium transport modeling for retention and permeation experiments. Fusion Engineering and Design, 203:114438, 6 2024. URL: https://linkinghub.elsevier.com/retrieve/pii/S0920379624002916, doi:10.1016/J.FUSENGDES.2024.114438.

@article{Shimada2024114438,
    author = "Shimada, Masashi and Simon, Pierre-Clément A. and Icenhour, Casey T. and Singh, Gyanender",
    doi = "10.1016/J.FUSENGDES.2024.114438",
    issn = "0920-3796",
    journal = "Fusion Engineering and Design",
    month = "6",
    pages = "114438",
    publisher = "North-Holland",
    title = "Toward a high-fidelity tritium transport modeling for retention and permeation experiments",
    volume = "203",
    url = "https://linkinghub.elsevier.com/retrieve/pii/S0920379624002916",
    year = "2024"
}

(fusion/mcf/divertor_monoblock/divertor_monoblock.i)

## Tokamak divertor monoblock
## Application: TMAP8
## POC: Masashi Shimada (masashi.shimada at inl.gov)
## If using or referring to this model, please cite as explained in
## https://mooseframework.inl.gov/virtual_test_bed/citing.html

### 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
[]

General description of the simulation case and corresponding input file
Results
Complete input file
References

Symbol	Variable or Physical Property	Unit
$var element = document.getElementById("moose-equation-c57596d9-f3f3-49f3-9dd2-f376aa9e46ab");katex.render("C_s", element, {displayMode:false,throwOnError:false});$	Concentration of solute (mobile) species	m $var element = document.getElementById("moose-equation-a1607c02-e95c-4180-a332-712096dd8546");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-60ade4c7-264a-4d60-afda-c91cb863d721");katex.render("C_t", element, {displayMode:false,throwOnError:false});$	Concentration of trapped species	m $var element = document.getElementById("moose-equation-7b608275-6c6c-4c0c-806a-7855a8b4baec");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-43cdbf35-034a-49ae-a6bb-03ea90a5e9eb");katex.render("T", element, {displayMode:false,throwOnError:false});$	Temperature	K
$var element = document.getElementById("moose-equation-de8e29c7-7b16-46b3-81d6-ed067b12ff48");katex.render("D", element, {displayMode:false,throwOnError:false});$	Diffusivity of solute (mobile) species	m $var element = document.getElementById("moose-equation-f35d3f15-4886-4215-b027-1bab3e55d8a7");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ s $var element = document.getElementById("moose-equation-bbf971b2-3a08-40c8-b6b8-a9516023ed4a");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-e2ed5350-8ec3-4a50-9329-47fca327be9c");katex.render("K_s", element, {displayMode:false,throwOnError:false});$	Solubility of solute (mobile) species	m $var element = document.getElementById("moose-equation-57c31a80-1531-4eda-9630-dab7e08fd428");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$ Pa $var element = document.getElementById("moose-equation-6d305d3d-844d-485f-9d00-759087edc748");katex.render("^{-1/2}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-2c497213-8bc5-4109-9d1b-555cc9ff7f40");katex.render("C_{total}", element, {displayMode:false,throwOnError:false});$	Total concentration of species $var element = document.getElementById("moose-equation-bf31cb82-a825-4f06-a78c-5cb92b05a769");katex.render("C_{total}=C_s + C_t", element, {displayMode:false,throwOnError:false});$	m $var element = document.getElementById("moose-equation-dccb711f-8a02-4da0-bfa4-d3f836c15174");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-73b4fe33-3783-487e-bbbd-22c57821a16c");katex.render("C_t^e", element, {displayMode:false,throwOnError:false});$	Concentration of empty trapping sites	m $var element = document.getElementById("moose-equation-337d6b20-c0d9-43af-95f7-796853ebe086");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-e70df9aa-e269-4d2f-bb2a-acbce2bf0140");katex.render("C_t^0", element, {displayMode:false,throwOnError:false});$	Concentration of trapping sites $var element = document.getElementById("moose-equation-5a801d00-818a-4bbd-981c-2a9cf4a1fc77");katex.render("C_t^e=C_t^0-C_t", element, {displayMode:false,throwOnError:false});$	m $var element = document.getElementById("moose-equation-8a0c2c6f-dce6-476a-a296-305eea8d9d22");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-f20d63bd-8621-4e6c-b90a-8c3de8c143bd");katex.render("\\alpha_t", element, {displayMode:false,throwOnError:false});$	Trapping rate coefficient, $var element = document.getElementById("moose-equation-e83d68b7-bd92-43be-a627-b8a763fc027f");katex.render("2.75 \\times 10^{11}", element, {displayMode:false,throwOnError:false});$	s $var element = document.getElementById("moose-equation-bdbdf361-e810-497f-bae4-9cb3049906ca");katex.render("^{\u22121}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-2bae83ff-bd9f-4543-982f-0744a12a40fb");katex.render("\\alpha_r", element, {displayMode:false,throwOnError:false});$	Release rate coefficient $var element = document.getElementById("moose-equation-0f5fde1d-a33e-4839-bef5-70ede0cbb4dc");katex.render("\\alpha_r=\\alpha_{r0} \\exp\u2061((-E_{dt})/(k_b T))", element, {displayMode:false,throwOnError:false});$	s $var element = document.getElementById("moose-equation-3cee6a40-9262-44d2-82a7-fd4c1e4d528f");katex.render("^{\u22121}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-308d2bdd-65e3-4f8c-ba02-7c6f5ad25aa5");katex.render("\\alpha_{r0}", element, {displayMode:false,throwOnError:false});$	Pre-exponential factor, $var element = document.getElementById("moose-equation-769eb284-a630-4c10-9379-d0a7aac9abee");katex.render("8.40 \\times 10^{12}", element, {displayMode:false,throwOnError:false});$	s $var element = document.getElementById("moose-equation-9a4b2f96-17d0-4ad0-9e93-ebe126c80f5f");katex.render("^{\u22121}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-3efba1f5-a3c7-4395-8c57-5591023939f3");katex.render("E_{dt}", element, {displayMode:false,throwOnError:false});$	Detrapping energy	eV
$var element = document.getElementById("moose-equation-973378af-dade-4db4-b372-5fd5ab79be77");katex.render("N", element, {displayMode:false,throwOnError:false});$	Atomic number density	m $var element = document.getElementById("moose-equation-c92d51f6-f0b9-43d5-b98f-691df048445d");katex.render("^{\u22123}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-5c5a8a03-7812-4bcb-a0c8-663071d655d5");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$	Density	g m $var element = document.getElementById("moose-equation-13f3cc9a-32fb-40dd-8a94-8b9f42a6a1b1");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-ec00a3a5-2de4-4975-a564-c5fc648e65c5");katex.render("c_p", element, {displayMode:false,throwOnError:false});$	Specific heat	J kg $var element = document.getElementById("moose-equation-3fbe7bb8-b12f-4520-8b3d-4f6262ac941c");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ K $var element = document.getElementById("moose-equation-deaa27eb-7417-402c-8a44-c2b7dc5696cb");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-1d89aa12-23a0-489e-ba28-d18ba59ef2f9");katex.render("k_T", element, {displayMode:false,throwOnError:false});$	Thermal conductivity	W m $var element = document.getElementById("moose-equation-efaa7183-2d07-4537-87ef-2798355bba26");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ K $var element = document.getElementById("moose-equation-7510ac13-b283-4fce-a87c-37270b5b774f");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$