Links to other pages for this example : - Overview - Additional theory - Preparing the simulation
Input file preparation
This calculation is carried out in 3 stages
Implantation: the material sample is bombarded with deuterium ions at constant temperature
Resting: the source term is switched off, slow off-gassing occurs due to the low temperature
Desorption: Temperature is raised at a constant rate and the rate of desorption is measured from each surface
Multiapp
For this example a multiapp structure has been used to combine the three simulation steps. The final, desorption step is the dominant process and each simulation step launches the step which precedes it before beginning its own calculation. The transfer block is used to copy the results from variables in sub-app calculations back to the calling app.
[MultiApps]
[sub_app]
type = FullSolveMultiApp
#app_type = BlueKiteApp
input_files = 'resting_multi.i'
# positions = '0 0.029 0
# 0.015 0.0.029 0
# 0.015 0.014 0
# 0.0 0.014 0'
execute_on = INITIAL
clone_master_mesh = true
[]
[]
[Transfers]
# [to_sub]
# type = MultiAppCopyTransfer
# multi_app = sub_app
# source_variable = 'surface_temp'
# variable = 'surface_temp'
# direction = to_multiapp
# []
[mobile_transfer]
type = MultiAppCopyTransfer #MultiAppScalarToAuxScalarTransfer
multi_app = sub_app
source_variable = 'Mobile'
variable = 'Mobile'
direction = from_multiapp
[]
[trap_1_transfer]
type = MultiAppCopyTransfer #MultiAppScalarToAuxScalarTransfer
multi_app = sub_app
source_variable = 'Trapped_1'
variable = 'Trapped_1'
direction = from_multiapp
[]
[trap_2_transfer]
type = MultiAppCopyTransfer #MultiAppScalarToAuxScalarTransfer
multi_app = sub_app
source_variable = 'Trapped_2'
variable = 'Trapped_2'
direction = from_multiapp
[]
[trap_3_transfer]
type = MultiAppCopyTransfer #MultiAppScalarToAuxScalarTransfer
multi_app = sub_app
source_variable = 'Trapped_3'
variable = 'Trapped_3'
direction = from_multiapp
[]
[]
Common parameters
The parameters used for this example have been taken from Delaporte-Mathurin et al. (2019), the experimental set-up follows Ogorodnikova et al. (2003)
Table 1: Example innput parameters
| parameter name | symbol | value | Unit |
|---|---|---|---|
| Density | |||
| Lattice sites | |||
| Trap Density (low energy intrinsic) | |||
| Trap Density (high energy intrinsic) | |||
| Maximum extrinsic trap Density type a | |||
| Maximum extrinsic trap Density type b | |||
| Diffusion rate pre-exponential factor | |||
| Diffusion Energy | 0.39 | ||
| Lattice constant | |||
| De-trapping pre-exponential factor (all traps) | |||
| De-trapping Energy, trap 1 | |||
| De-trapping Energy, trap 2 | |||
| De-trapping Energy, trap 3 | |||
| Temperature (implantation and resting) | |||
| Temperature ramping rate | |||
| Simulation time | |||
| Extrinsic trap creation rate | |||
| Extrinsic trap creation rate | |||
| Plastic deformation region |
Also common between all calculation stages are:
A zero valued dirichlet condition is applied for the mobile phase at all external surfaces
The core set of kernels are identical for all Achlys calculations. Though here the external source term is only used for the implantation phase.
Material properties as described in the table above
The problem geometry and mesh
Stage 1: Implantation
In this simulation:
A gaussian implantation profile is used for the source term
A dirichlet boundary condition with a zero value for the surface concentration of the mobile phase is used for all external surfaces
Constant values are used for the densities of two instrisic trap types, labelled trap 1 and trap2. These are uniformly distributed throughout the entire material domain.
A third trap is simulated whose concentration increases over time as the material is irradiated. This trap only appears in a narrow region near the surface exposed to particle bombardment.
The temperature is kept at a constant 300 K
The total simulation time is 400 s
Mesh
[Mesh]
[./unlabelled]
file = /Users/sdixon/projects/blue_kite/problems/thermal_desorption/optimisation/ogorodnikova_1k.msh
type = FileMeshGenerator
construct_side_list_from_node_list = true
[../]
[./block_1]
type = AllSideSetsByNormalsGenerator
input = unlabelled
#fixed_normal = false
[../]
[]
Source term: the implantation profile
The implantation profile for the experiment is a function of the composition, energy, and direction of the incident particles. This is generally determined by an external tool such as SRIM. Here a gaussian profile with mean 4.5 nm and standard deviation of 4.5 nm has been used following Delaporte-Mathurin et al. (2019). The scale factor applied means units of are obtained for the source term when the value of the gaussian function is multiplied by the particle flux in .
It is also possible to define other profiles (exonweib shown below) or alternatively to read in values from a CSV file.
[Functions]
[./dts]
type = PiecewiseLinear
x = '0 1e-4 5e-4 1e-3'
y = '1e-4 1e-4 5e-4 5e-4'
[../]
[./Gaussian_implant]
type = ParsedFunction
value = 'scale * exp( -0.5 * ((x - mean) / sd)^2)'
vars = 'scale mean sd'
vals = '0.93e8 4.5e-9 4.5e-9'
[../]
[./Exponweib_implant]
type = ParsedFunction
value = 'peak * (a * c / scale) * (( 1 - exp(-1 * (((x*1e9)-loc)/scale)^c))^(a-1) ) * exp(-1 * (((x*1e9)-loc)/scale)^c)*(((x*1e9)-loc)/scale)^(c-1)'
vars = 'peak a c loc scale'
vals = '5.82e8 0.90 1.96 -0.013 4.909'
[../]
[]
Kernels
The kernels shown below represent the standard set for any Achlys calculation. Note only that the source term has been selected to use the gaussian implantation profile defined above. The "value" field here is the value of the particle flux which has been scaled by the density of solute atoms. For the case of deuterium in tungsten the scale factor is and the units are .
[Kernels]
[./H3_source]
type = ADBodyForce
variable = Mobile
value = 4e-10
function = Gaussian_implant
[]
[./H3_diffusion_eq1]
type = ADMatDiffusion
variable = Mobile
Diffusivity = D
[../]
[./mobile_time_deriv]
type = ADTimeDerivative
variable = Mobile
[../]
[./trapping_equilibrium_equation1]
type = ADTrappingEquilibriumEquation
variable = Trapped_1
v = Mobile
n_traps = n1
vi = V1
[../]
[./trapping_equilibrium_equation2]
type = ADTrappingEquilibriumEquation
variable = Trapped_2
v = Mobile
n_traps = n2
vi = V2
[../]
[./trapping_equilibrium_equation3]
type = ADTrappingEquilibriumEquation
variable = Trapped_3
v = Mobile
n_traps = n3
vi = V3
[../]
[./trapped_time_deriv_couple]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_1
[../]
[./trapped_time_deriv_couple2]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_2
[../]
[./trapped_time_deriv_couple3]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_3
[../]
[./trapped_time_deriv]
type = ADTimeDerivative
variable = Trapped_1
[../]
[./trapped_time_deriv2]
type = ADTimeDerivative
variable = Trapped_2
[../]
[./trapped_time_deriv3]
type = ADTimeDerivative
variable = Trapped_3
[../]
[]
Material Properties
A specific material class has been developed to model the transient evolution of extrinsic traps that is used in this calculation. The example inputs for the ExtrinsicTransientTrappingMaterial class are shown below.
Using this class means that the density of availble type-3 traps increases over time within this simulation stage. The rate of increase depends on the flux, whose shape is defined by the "function" parameter and whose magnitude is given by the "flux" field.
[Materials]
[./implant]
type = ExtrinsicTransientTrappingMaterial2
# Energies
E_diff = 0.39
E1 = 0.87
E2 = 1.0
E3 = 1.5
k_boltz = 8.617333E-5
# pre-exponential rate constants
v0 = 1e13
D0 = 4.1e-7
lambda = 1.1E-10
# unused
rho = 1 # unecessary scaling factor, do not use
# site densities
n_sol = 6
n1 = 1e-3
n2 = 4e-4
n3a_max = 1e-1
n3b_max = 1e-2
# trap creation rates
eta_a = 6e-4
eta_b = 2e-4
# flux distribution parameters
flux = 4e-10
function = Gaussian_implant
xp = 1e-6
#Temperature
const_T = 300
block = 'Tungsten'
# thermal properties
conductivity = 150 # W/K
Cp = 137 # J/(kg K)
density = 19300 # kg/m3
[../]
[]
Boundary Conditions
[BCs]
[./PFC]
type = ADDirichletBC
variable = Mobile
boundary = 1
value = 0
[../]
[]
Executioner
There are 3 things to note within this input block:
The total simulation time is set by the end_time parameter. This is where we specify the duration of the implantation phase
Adaptive timestepping is used to achieve the most efficient time-steps.
PetSc allows any numerical scheme to be substituted here to carry out the solve. Here a bi-conjugate gradient linear method with LU preconditioner and second-order time integration has been used.
[Executioner]
type = Transient
solve_type = NEWTON
scheme = bdf2
automatic_scaling = True
# compute_scaling_once=False
# scaling_group_variables = 'Trapped_3 Trapped_2 Trapped_1 Mobile'
resid_vs_jac_scaling_param = 0.8
l_max_its = 100
nl_max_funcs = 7
# nl_rel_tol = 1e-7
nl_abs_tol = 1e-29
# Set PETSc parameters to optimize solver efficiency
timestep_tolerance = 0.01
petsc_options_iname = '-ksp_type -pc_type -pc_factor_shift_type'
petsc_options_value = 'bcgs lu NONZERO'
end_time = 400
[TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 4
cutback_factor = 0.8
growth_factor = 1.2
dt = 5e-4
[]
[]
Visualising implantation
The implantation step of this calculation can be run as a standalone simulation. The input file can be run by passing it's path using the "-i" flag as shown below; outputs will be written to the same folder in which the input is located.
achlys_exe=~/projects/achlys/achlys-opt
input_file=~/projects/achlys/problems/thermal_desorption/ogorodnikova/tds_multiapp/implant_sub.i
./achlys-opt -i $input_file
Paraview can be used to visualise the resulting exodus file. in addition the total quantities of deuterium in each phase will be written to a csv file at each timestep.
Stage 2: Resting
In this stage:
The source term is switched off so no new particles are implanted
The creation of new type-3 traps also stops when the implantation flux stops meaning the concentration and distribution of traps remains static from the end of implantation
Desorption continues to occur from the external surfaces, though this is limited due to the low temperature; the boundary condition remains a zero-valued dirichlet condition
The temperature is a constant 300 K
The total resting time is 50 s
Kernels
The Kernels here are the standard set for Achlys calculations. Note only that the source term has now been removed for the resting phase.
[Kernels]
[./H3_diffusion_eq1]
type = ADMatDiffusion
variable = Mobile
Diffusivity = D
[../]
[./mobile_time_deriv]
type = ADTimeDerivative
variable = Mobile
[../]
[./trapping_equilibrium_equation1]
type = ADTrappingEquilibriumEquation
variable = Trapped_1
v = Mobile
n_traps = n1
vi = V1
[../]
[./trapping_equilibrium_equation2]
type = ADTrappingEquilibriumEquation
variable = Trapped_2
v = Mobile
n_traps = n2
vi = V2
[../]
[./trapping_equilibrium_equation3]
type = ADTrappingEquilibriumEquation
variable = Trapped_3
v = Mobile
n_traps = n3
vi = V3
[../]
[./trapped_time_deriv_couple]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_1
[../]
[./trapped_time_deriv_couple2]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_2
[../]
[./trapped_time_deriv_couple3]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_3
[../]
[./trapped_time_deriv]
type = ADTimeDerivative
variable = Trapped_1
[../]
[./trapped_time_deriv2]
type = ADTimeDerivative
variable = Trapped_2
[../]
[./trapped_time_deriv3]
type = ADTimeDerivative
variable = Trapped_3
[../]
[]
Material Properties
[Materials]
[./implant]
type = ExtrinsicStaticTrappingMaterial
# Energies
E_diff = 0.39
E1 = 0.87
E2 = 1.0
E3 = 1.5
k_boltz = 8.617333E-5
# pre-exponential rate constants
v0 = 1.0E13
D0 = 4.1E-7
lambda = 1.1E-10
# unused
rho = 1 # unecessary scaling factor, do not use
# site densities
n_sol = 6
n1 = 1E-3 #1.0E25
n2 = 4e-4
n3a_max = 1e-1
n3b_max = 1e-2
# trap creation rates
eta_a = 6e-4
eta_b = 2e-4
trap_evolution_time = 400
# flux distribution parameters
flux = 4e-10
function = Gaussian_implant
xp = 1e-6
#Temperature
const_T = 300
block = 'Tungsten' #'Implantation_region'
# thermal properties
conductivity = 150 # W/K
Cp = 137 # J/(kg K)
density = 19300 # kg/m3
[../]
[]
Boundary Conditions
During resting the mobile concentration at the surface is considered to be zero. Dirichlet conditions are applied at both surfaces.
[BCs]
[./desorption]
type = ADDirichletBC
variable = Mobile
boundary = 1
value = 0
[../]
[./back_desorption]
type = ADDirichletBC
variable = Mobile
boundary = 2
value = 0
[../]
[]
Stage 3: Desorption
In this stage:
The temperature is raised at a constant rate of 8 K/s.
The boundary condition remains a zero-valued dirichlet condition.
The concentration and distribution of traps remains static. This state is unchanged from the end of the implantation calculation.
The rate of desorption from each surface is measured using the ADSideFluxIntegral postprocessor and results are written to a CSV file.
Kernels
[Kernels]
[./H3_diffusion_eq1]
type = ADMatDiffusion
variable = Mobile
Diffusivity = D
[../]
[./mobile_time_deriv]
type = ADTimeDerivative
variable = Mobile
[../]
[./trapping_equilibrium_equation1]
type = ADTrappingEquilibriumEquation
variable = Trapped_1
v = Mobile
n_traps = n1
vi = V1
[../]
[./trapping_equilibrium_equation2]
type = ADTrappingEquilibriumEquation
variable = Trapped_2
v = Mobile
n_traps = n2
vi = V2
[../]
[./trapping_equilibrium_equation3]
type = ADTrappingEquilibriumEquation
variable = Trapped_3
v = Mobile
n_traps = n3
vi = V3
[../]
[./trapped_time_deriv_couple]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_1
[../]
[./trapped_time_deriv_couple2]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_2
[../]
[./trapped_time_deriv_couple3]
type = ADCoupledTimeDerivative
variable = Mobile
v = Trapped_3
[../]
[./trapped_time_deriv]
type = ADTimeDerivative
variable = Trapped_1
[../]
[./trapped_time_deriv2]
type = ADTimeDerivative
variable = Trapped_2
[../]
[./trapped_time_deriv3]
type = ADTimeDerivative
variable = Trapped_3
[../]
[]
Material Properties
[Materials]
[./implant]
type = ExtrinsicStaticTrappingMaterialRampingT
# Energies
E_diff = 0.39
E1 = 0.87
E2 = 1.0
E3 = 1.50
k_boltz = 8.617333E-5
# pre-exponential rate constants
v0 = 1.0E13
D0 = 4.1E-7
lambda = 1.1E-10
# unused
rho = 1 # unecessary scaling factor, do not use
# site densities
n_sol = 6
n1 = 1E-3 #1.0E25
n2 = 4e-4
n3a_max = 1e-1
n3b_max = 1e-2
# trap creation rates
eta_a = 6e-4
eta_b = 2e-4
trap_evolution_time = 400
# flux distribution parameters
flux = 4e-10
function = Gaussian_implant
xp = 1e-6
#Temperature
initial_T = 300
beta = 8
block = 'Tungsten' # 'Stopping_zone' # Plastic_region
# thermal properties
conductivity = 150 # W/K
Cp = 137 # J/(kg K)
density = 19300 # kg/m3
[../]
[]
Boundary Conditions
For the desorption phase the mobile concentration at the surface is considered to be zero. Dirichlet conditions are applied at both surfaces.
[BCs]
[./desorption]
type = ADDirichletBC
variable = Mobile
boundary = 1
value = 0
[../]
[./back_desorption]
type = ADDirichletBC
variable = Mobile
boundary = 2
value = 0
[../]
[]
Outputs
[Postprocessors]
[./total_mobile]
type = VariableIntegral
variable = Mobile
[../]
[./total_trap_1]
type = VariableIntegral
variable = Trapped_1
[../]
[./total_trap_2]
type = VariableIntegral
variable = Trapped_2
[../]
[./total_trap_3]
type = VariableIntegral
variable = Trapped_3
[../]
[/pfc_flux]
type = ADSideFluxIntegral
variable = Mobile
boundary = 1
diffusivity = D
[../]
[/back_flux]
type = ADSideFluxIntegral
variable = Mobile
boundary = 2
diffusivity = D
[../]
[]
[Outputs]
#execute_on = 'timestep_end'
exodus = false # Output Exodus format
csv = true
[]
References
- Rémi Delaporte-Mathurin, Etienne A. Hodille, Jonathan Mougenot, Yann Charles, and Christian Grisolia.
Finite element analysis of hydrogen retention in iter plasma facing components using festim.
Nuclear Materials and Energy, 21:100709, 2019.
URL: https://www.sciencedirect.com/science/article/pii/S2352179119300547, doi:https://doi.org/10.1016/j.nme.2019.100709.[BibTeX]
- O.V Ogorodnikova, J Roth, and M Mayer.
Deuterium retention in tungsten in dependence of the surface conditions.
Journal of Nuclear Materials, 313-316:469–477, 2003.
Plasma-Surface Interactions in Controlled Fusion Devices 15.
URL: https://www.sciencedirect.com/science/article/pii/S0022311502013752, doi:https://doi.org/10.1016/S0022-3115(02)01375-2.[BibTeX]