Tritium Transport in Li2TiO3 with Voronoi Grain Boundaries

[griffda]

Hi everyone,

I am trying to model tritium transport in a 2D, polycrystalline, grain-level simulation of lithium ceramic (Li2TiO3). I would like to use voronoi tessellation to configure the grain geometry and mesh, however I am having difficulty in getting this to work properly.

The current script generates the polycrystalline mesh (see screenshot), but the T transport data in the output is not physical and no T-trapping data is generated. I see that the example use of the PolyCrystalVoronoi user object tracks grain growth and voids, however in this case these parameters are not needed.

Below is my input file. Any help would be hugely appreciated, thank you in advance.

# TMAP8 Tritium Transport in Li2TiO3 with Voronoi Grain Boundaries
# 
# Features:
# - Voronoi polycrystalline structure (50 grains, 8 order parameters)
# - GBs as fast diffusion pathways (no trapping)
# - Temperature-dependent properties
# - Dirichlet BCs on all boundaries

# ============================================================================
# PARAMETERS
# ============================================================================

# Voronoi parameters
grain_count = 50                     # Number of grains in the domain
op_num = 8                           # Number of order parameters (typically sqrt(grain_count))
avg_grain_size = 2e-6                # Target average grain size (m)
gb_width = 0.11e-6                   # Interface width (7 nm, from example)

# Calculate domain size based on grain count and average size
domain_size = ${fparse sqrt(grain_count) * avg_grain_size * 1.5}

# Mesh discretization
nx_base = 60                         # Base elements (will be refined)
ny_base = 60

# Physical parameters
temperature = 773                    # K
k_B = 8.617e-5                      # eV/K

# Bulk material properties
D_0_bulk = 9.5e-6                   # m^2/s
E_D_bulk = 0.80                     # eV

# GB properties (100× faster)
D_0_gb = 9.3e-7                     # m^2/s
E_D_gb = 0.41                       # eV

# Calculate diffusivities
D_bulk = ${fparse D_0_bulk * exp(-E_D_bulk / (k_B * temperature))}
D_gb = ${fparse D_0_gb * exp(-E_D_gb / (k_B * temperature))}

# Trapping parameters
trapping_prefactor = 1e-16
trapping_energy = 0.52
detrapping_prefactor = 1e13
detrapping_energy = 1.37

# Calculate rate constants (pre-multiply by Arrhenius factors)
alpha_t = ${fparse trapping_prefactor * exp(-trapping_energy / (k_B * temperature))}
alpha_r = ${fparse detrapping_prefactor * exp(-detrapping_energy / (k_B * temperature))}

trapping_site_fraction = 1e-3
li2tio3_density = 3.43e28   
trap_per_free = 1.0

# Normalize
cl = 1e18
N_0 = 1e25
E_formation = 1.5
N_defects = ${fparse N_0 * exp(-E_formation / (k_B * temperature))}
N_normalized = ${fparse N_defects / cl}

# Source terms
source_bulk = 1e13
source_gb = 1e13                    # at/m^3/s

# Simulation time
end_time = 5

# ============================================================================
# MESH
# ============================================================================

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = ${nx_base}
  ny = ${ny_base}
  xmax = ${domain_size}
  ymax = ${domain_size}
  elem_type = QUAD4
  uniform_refine = 1  # Initial refinement
[]

# ============================================================================
# GLOBAL PARAMETERS
# ============================================================================

[GlobalParams]
  op_num = ${op_num}
  var_name_base = gr
  grain_num = ${grain_count}
[]

# ============================================================================
# USER OBJECTS
# ============================================================================

[UserObjects]
  [voronoi]
    type = PolycrystalVoronoi
    grain_num = ${grain_count}
    rand_seed = 12345
    coloring_algorithm = jp
  []
  
  [grain_tracker]
    type = GrainTracker
    threshold = 0.2
    compute_var_to_feature_map = true
    flood_entity_type = elemental
  []
[]

# ============================================================================
# VARIABLES (Tritium transport)
# ============================================================================

[Variables]
  [mobile]
    initial_condition = 1e-10
    scaling = 1e-15
  []
  [trapped]
    initial_condition = 0
    scaling = 1e-15
  []
  [PolycrystalVariables]
    var_name_base = gr
    op_num = ${op_num}
  []
[]

# ============================================================================
# INITIAL CONDITIONS
# ============================================================================

[ICs]
  [PolycrystalICs]
    [PolycrystalColoringIC]
      polycrystal_ic_uo = voronoi
    []
  []
[]


# ============================================================================
# AUXVARIABLES
# ============================================================================

[AuxVariables]
  [bnds]
    # Grain boundary indicator (0 in grains, 1 at GBs)
    order = FIRST
    family = LAGRANGE
  []
  
  [unique_grains]
    order = CONSTANT
    family = MONOMIAL
  []
  
  [var_indices]
    order = CONSTANT
    family = MONOMIAL
  []
  
  [temp]
    initial_condition = ${temperature}
  []
  
  [trap_sites]
    initial_condition = ${N_normalized}
  []
[]

# ============================================================================
# AUXKERNELS
# ============================================================================

[AuxKernels]
  [bnds_aux]
    type = BndsCalcAux
    variable = bnds
    execute_on = 'initial timestep_end'
  []
  
  [unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    flood_counter = grain_tracker
    field_display = UNIQUE_REGION
    execute_on = 'initial timestep_end'
  []
  
  [var_indices]
    type = FeatureFloodCountAux
    variable = var_indices
    flood_counter = grain_tracker
    field_display = VARIABLE_COLORING
    execute_on = 'initial timestep_end'
  []
[]

# ============================================================================
# KERNELS (Tritium transport + grain variables)
# ============================================================================

[Kernels]
  # Tritium transport kernels
  [mobile_time]
    type = TimeDerivative
    variable = mobile
  []
  
  [mobile_diffusion]
    type = MatDiffusion
    variable = mobile
    diffusivity = diffusivity
  []
  
  [mobile_source]
    type = BodyForce
    variable = mobile
    value = ${source_bulk}
  []
  
  [mobile_trapped_sink]
    type = ADCoefCoupledTimeDerivative
    variable = mobile
    v = trapped 
    coef = -${trap_per_free}
  []
  
  # Add time derivatives for grain variables to keep them active (frozen)
  [./gr0_time]
    type = TimeDerivative
    variable = gr0
  [../]
  [./gr1_time]
    type = TimeDerivative
    variable = gr1
  [../]
  [./gr2_time]
    type = TimeDerivative
    variable = gr2
  [../]
  [./gr3_time]
    type = TimeDerivative
    variable = gr3
  [../]
  [./gr4_time]
    type = TimeDerivative
    variable = gr4
  [../]
  [./gr5_time]
    type = TimeDerivative
    variable = gr5
  [../]
  [./gr6_time]
    type = TimeDerivative
    variable = gr6
  [../]
  [./gr7_time]
    type = TimeDerivative
    variable = gr7
  [../]
[]

# ============================================================================
# NODAL KERNELS (Trapping/detrapping)
# ============================================================================

[Functions]
  [trap_sites_function]
    type = ParsedFunction
    expression = '${trapping_site_fraction}'  # Base value
  []
[]

[NodalKernels]
  [trapped_time_nodal]
    type = TimeDerivativeNodalKernel
    variable = trapped
  []
  
  [trapping]
    type = TrappingNodalKernel
    variable = trapped
    mobile_concentration = mobile
    alpha_t = ${alpha_t}
    trapping_energy = 0
    N = ${li2tio3_density}
    Ct0 = 'trap_sites_function'
    temperature = temp
    trap_per_free = ${trap_per_free}
  []
  
  [releasing]
    type = ReleasingNodalKernel
    variable = trapped
    alpha_r = ${alpha_r}
    detrapping_energy = 0
    temperature = temp
  []
[]

# ============================================================================
# BOUNDARY CONDITIONS
# ============================================================================

[BCs]
  [left_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'left'
    value = 0
  []
  [right_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'right'
    value = 0
  []
  [top_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'top'
    value = 0
  []
  [bottom_bc]
    type = DirichletBC
    variable = mobile
    boundary = 'bottom'
    value = 0
  []
[]

# ============================================================================
# MATERIALS
# ============================================================================

[Materials]
  [diffusivity]
    type = ParsedMaterial
    property_name = diffusivity
    coupled_variables = 'bnds'
    expression = '${D_bulk} + (${D_gb} - ${D_bulk}) * bnds^2 * (3 - 2*bnds)'
  []
  
  [trap_distribution]
    type = ParsedMaterial
    property_name = trap_distribution_function
    coupled_variables = 'bnds'
    expression = '${trapping_site_fraction} * (1 - bnds^2 * (3 - 2*bnds))'
  []
  
  [mobility]
    type = GenericConstantMaterial
    prop_names = 'L'
    prop_values = '0.0'
  []

  [kappa]
    type = GenericConstantMaterial
    prop_names = 'kappa_op'
    prop_values = '1e-9'
  []
[]

# ============================================================================
# POSTPROCESSORS
# ============================================================================

[Postprocessors]
  [total_mobile]
    type = ElementIntegralVariablePostprocessor
    variable = mobile
  []
  
  [total_trapped]
    type = ElementIntegralVariablePostprocessor
    variable = trapped
  []
  
  [num_grains]
    type = FeatureFloodCount
    variable = bnds
    threshold = 0.5
  []
  
  [dt]
    type = TimestepSize
  []
  [flux_right]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'right'
    diffusivity = diffusivity
  []
  [flux_left]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'left'
    diffusivity = diffusivity
  []
  [flux_top]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'top'
    diffusivity = diffusivity
  []
  [flux_bottom]
    type = SideDiffusiveFluxIntegral
    variable = mobile
    boundary = 'bottom'
    diffusivity = diffusivity
  []

[]

# ============================================================================
# PRECONDITIONING
# ============================================================================

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

# ============================================================================
# EXECUTIONER
# ============================================================================

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = PJFNK
  
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
  
  l_max_its = 50
  l_tol = 1e-3
  nl_max_its = 10
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-3
  
  start_time = 0
  end_time = ${end_time}
  
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.001             # Smaller initial dt
    optimal_iterations = 8
    iteration_window = 2
    growth_factor = 1.1    # Slower growth
    cutback_factor = 0.5   # More aggressive cutback
  []
[]

# ============================================================================
# OUTPUTS
# ============================================================================

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  
  [console]
    type = Console
    fit_mode = 80
  []
[]

[simopier]

@griffda, just to make sure I understand what you want to do, please confirm or update the following statements. Note that I have not yet use your input yet since I am on travel this week.

• you want to generate a polycrystalline microstructure and model tritium generation and transport in it, including trapping.
• You aim to have accelerated grain boundary diffusion.
• You problem is that you are not seeing any trapping happening.
• You also observe grain evolution, which you do not need in your simulation.
• You currently aim to generate tritium in your system, and then have it diffuse to the edges of your domain, where the value is set to 0.

Could also specify what you mean by:

but the T transport data in the output is not physical
What is happening that is not physical?

Some small suggestions:

• for the grains, I would recommend using periodic boundary conditions instead of the Neumann boundary condition imposed by default. It would prevent the formation of elongated grains around the edges that you can see in the image.

[simopier]

Hi @griffda, sorry for the delay here. I have been cleaning up some input files to be able to share them with you.
I'll do my best to finalize this this week and send this over to you.
Do not hesitate to ping me here if you do not hear back by Thursday.

[griffda]

pinging you @simopier

[simopier]

Hi @griffda,

Sorry for the excessive delay.. Below is some information. I am also sharing additional inputs via email.

you want to generate a polycrystalline microstructure and model tritium generation and transport in it, including trapping

Correct. Similar to what is conducted in this paper with Cs in SiC. I notice that this was conducted using the MARMOT code for fission fuel materials which is not open source. An example input file for this may be useful.

Note that in this paper, I did not include trapping explicitly in the model. This paper uses the asymptotic expansion homogenization (AEH) approach to determine the effective Cs transport, rather than model the transport directly. The earlier paper on Ag transport in SiC referenced in the one you listed provides a bit more details about the nature of the simulations. I will share with you via emails the input files relevant to these calculations. It will consist of:

• a first input file used to generate the polycrystal with the desired grain boundary thickness/smoothness. This input file is similar to the one used for the first part of this example case for pore-scale diffusion. It uses outputs of a python script that automatically creates the desired polycrystalls in EBSD formats (desired domain size, resolution, grain dimensions, etc. with the periodic grain structure required for the AEH method).
• a second input file used to perform the AEH calculations. This input file should work with TMAP8's executable (I say should because parts of it might be deprecated and might need to be updated. Error messages should help with this, but let me know if that's the case and you have any questions). Note that this one only accounts for diffusion. If you assume that trapping/detrapping occurs much faster than diffusion, you could simplify the problem to a difffusion-only calculations using D_eff = K_r / (K_t + K_r) with K_t and K_r the trapping and release rates, respectively. If you want to account explicitly for trapping during transients, then the AEH method (which assumes steady-state), might not be appropriate for this study.

The reason I mention grain evolution is because the example input files I have seen using the PolyCrystallineVoronoi user object are also modelling grain growth through the GrainTracker parameter. However, this is not something I want to track in my simulation. My concern is that the grain tracker is a required parameter.

I am sorry if I am missing something, but the grain tracker should help you track the grain boundaries by providing the order parameters for the grains, which you can combine to identify the position of the grain boundaries. That should help you define local properties such as species diffusivity at grain boundaries. However, you do not need to solve for the grain order parameters. In MOOSE language, the grain order parameters are static aux variables. In your input file, they are instead variables and will be solved for.
In that context (not having to solve for it), do you have any additional concerns about the fact that the grain tracker might be a required parameter?

What is happening that is not physical?

The mobile concentration of T in the output file is a constant value in the region e-20 (the source term is constant at 1e13). In other more simple tests I have completed the maximum mobile concentration is much closer to the source term value and the simulation reaches steady state after around 5s.

What do you mean by In other more simple tests? I would suggest trying to run something similar on a domain without grain boundaries (just one single grain), just to check that it indeed does what is expected. Then only 3-4 grains, then a larger simulation. Once we find a case that works, and then a case that doesn't as soon as we add the next step, it'll be easier to debug and move forward.

I hope that helps!
I apologize again for the delay. Thank you for your patience.

[griffda]

Hi @simopier,

Thank you for your response.

If you want to account explicitly for trapping during transients, then the AEH method (which assumes steady-state), might not be appropriate for this study.

I will take a look at your email with the input files. From what I have seen in my other simulations, the system approaches steady state reasonably fast, therefore this might be ok. However it is worth checking this over in the other problems I have been trying to run here: #336

I am sorry if I am missing something, but the grain tracker should help you track the grain boundaries by providing the order parameters for the grains, which you can combine to identify the position of the grain boundaries. That should help you define local properties such as species diffusivity at grain boundaries. However, you do not need to solve for the grain order parameters. In MOOSE language, the grain order parameters are static aux variables. In your input file, they are instead variables and will be solved for.
In that context (not having to solve for it), do you have any additional concerns about the fact that the grain tracker might be a required parameter?

I have probably misunderstood what it does. I had understood that it "tracks" the change in geometry of the grain with time, which is something I am not wanting to capture here. Could you give me an example of how I should use it in this case?

I would suggest trying to run something similar on a domain without grain boundaries (just one single grain), just to check that it indeed does what is expected.

I will do this and get back to you again.

[simopier]

Hi @griffda, I hope you are making good progress here too!

The grain tracker keeps track of grains, but it doesn't mean that grains have to evolve. So the cost can remain relatively low.
You can find more documentation there: https://mooseframework.inl.gov/TMAP8/source/postprocessors/GrainTracker.html

The input file I had sent to you via email should have provided an example.
Let's discuss if it wasn't clear.