In-pile, long-lived (or stable) isotopes

note:More Information Available

The material presented here is a subset of the information found in Toptan et al. (2023). See the paper for more information.

A solid sphere of radius $var element = document.getElementById("moose-equation-4737e006-8d25-4622-8ad3-3b69875ce271");katex.render("a", element, {displayMode:false,throwOnError:false});$ is subject to following boundary conditions:

A radial symmetry is considered at $var element = document.getElementById("moose-equation-ebbc750d-b6ec-42e1-899c-2a17e9765993");katex.render("r=0", element, {displayMode:false,throwOnError:false});$ : $var element = document.getElementById("moose-equation-ae21250c-d74f-4548-8416-7146a8ae478f");katex.render("\\frac{\\partial C}{\\partial t}\\big|_{r=0}=0", element, {displayMode:false,throwOnError:false});$ for all $var element = document.getElementById("moose-equation-9ef69d60-1acf-4bdc-b42d-cc6577905c77");katex.render("t\\geq 0", element, {displayMode:false,throwOnError:false});$ .
Concentration is kept at zero at the outer surface: $var element = document.getElementById("moose-equation-9ef6026b-dfc8-487d-a7bf-dcfe6fcd0540");katex.render("C(a,t)=0", element, {displayMode:false,throwOnError:false});$ for all $var element = document.getElementById("moose-equation-81f781a7-302d-4e9c-ba9f-ccf20d6d1255");katex.render("t\\geq 0", element, {displayMode:false,throwOnError:false});$ .

The in-pile conditions for a stable fission product ( $var element = document.getElementById("moose-equation-dde15b1c-3b4a-4266-8e49-0d963115e2aa");katex.render("\\lambda=0", element, {displayMode:false,throwOnError:false});$ ) are considered with the following settings: the initial concentration is zero in the sphere ( $var element = document.getElementById("moose-equation-b906e585-f954-443b-8cfd-7f32deb8c98a");katex.render("C_0=0", element, {displayMode:false,throwOnError:false});$ for $var element = document.getElementById("moose-equation-49d83eb3-9842-470a-962a-53d63b79da58");katex.render("0\\leq r\\leq a", element, {displayMode:false,throwOnError:false});$ ), and a non-zero source ( $var element = document.getElementById("moose-equation-93348aa8-b669-436a-ba96-ce1a90c6bd51");katex.render("p\\neq 0", element, {displayMode:false,throwOnError:false});$ ). In this study, we assume a constant source generation, $var element = document.getElementById("moose-equation-1e83e1cb-2e07-49d1-ad34-3427d31bf946");katex.render("p(t)=p", element, {displayMode:false,throwOnError:false});$ . The mass diffusion solution (or conservation of fission products), one of the Bison's governing equations, is solved in the spherical coordinates. The analytical expressions (Nabielek et al., 1974; Toptan et al., 2023) are given below for concentration, release rate, and total fractional release.

Concentration

The concentration profile, $var element = document.getElementById("moose-equation-13f76189-3e9d-46de-b560-b10094c6d4ad");katex.render("C(r,t)", element, {displayMode:false,throwOnError:false});$ is expressed as below. The steady-state ( $var element = document.getElementById("moose-equation-b57c6b24-ae51-4005-8162-ff10ff31ec7b");katex.render("\\tau\\rightarrow\\infty", element, {displayMode:false,throwOnError:false});$ ) solution for the concentration reduces to the first term on the right-hand side (RHS) of this relation. $(1)var element = document.getElementById("moose-equation-322dd8df-890f-4cf5-ac53-f3b1da4db584");katex.render(" \\frac{C(\\rho,\\tau)}{p/6D^\\prime} = (1-\\rho^2) - \\frac{12}{\\rho} \\sum_{n=1}^\\infty (-1)^{n+1} \\frac{\\exp{\\left(-n^2\\pi^2\\tau\\right)}}{n^3\\pi^3} \\sin{(n\\pi\\rho)}", element, {displayMode:true,throwOnError:false});$ with the following dimensionless numbers:

$var element = document.getElementById("moose-equation-17485eec-7fba-4165-88cc-956b18092dd3");katex.render("\\rho=r/a", element, {displayMode:false,throwOnError:false});$ is the dimensionless radius (unitless),
$var element = document.getElementById("moose-equation-da36a37f-8ce8-4733-8164-f789e6bbb629");katex.render("\\tau=D^\\prime t", element, {displayMode:false,throwOnError:false});$ is the diffusion time (unitless) in terms of the reduced diffusivity coefficient, $var element = document.getElementById("moose-equation-b762763b-288f-46d6-b2d0-bf09ed6e4514");katex.render("D^\\prime=D/a^2", element, {displayMode:false,throwOnError:false});$ (s $var element = document.getElementById("moose-equation-ae53184e-5e52-420b-8a46-41deef05ca10");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ).

The concentration, $var element = document.getElementById("moose-equation-c2e986e2-6add-481c-a3f8-5571d012e9fe");katex.render("C", element, {displayMode:false,throwOnError:false});$ , at a specific point is computed in the code via

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [conc_point_25]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = conc
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '${fparse 0.25*kernel_radius} 0 0'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial timestep_end'
  []
[]

Figure 1 compares analytical and computed concentration, $var element = document.getElementById("moose-equation-0c3b5003-37e8-456e-a2c6-0aacab641de2");katex.render("C/(p/6D^\\prime)", element, {displayMode:false,throwOnError:false});$ , at four different radial locations ( $var element = document.getElementById("moose-equation-08510e17-697a-4729-ac74-042baff21f01");katex.render("\\rho=", element, {displayMode:false,throwOnError:false});$ 0.25, 0.50, 0.75, and 0.95) as a function of dimensionless diffusion time, $var element = document.getElementById("moose-equation-51216d1f-b980-4e35-996b-3f2e4ab694b7");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ . $var element = document.getElementById("moose-equation-aaf9e487-275d-4f40-85b5-7d642477d989");katex.render("\\mu=0.1", element, {displayMode:false,throwOnError:false});$ is arbitrarily chosen; however, the concentration solution is independent of $var element = document.getElementById("moose-equation-c22c159f-b27c-46c3-a87d-e50f33257c85");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ , see Eq. (1).

Figure 1: Comparison of analytical and computed concentration, $var element = document.getElementById("moose-equation-4e0571af-64e4-4413-82ef-c025c65aea15");katex.render("C/(p/6D^\\prime)", element, {displayMode:false,throwOnError:false});$ at four different radial locations ( $var element = document.getElementById("moose-equation-3a665bea-fec5-481b-bb53-3e7367da5b9d");katex.render("\\rho=", element, {displayMode:false,throwOnError:false});$ 0.25, 0.50, 0.75, and 0.95) as a function of dimensionless diffusion time, $var element = document.getElementById("moose-equation-ec1e7ec1-818d-4c9b-92a7-1311c53f5b23");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ (arbitrarily chosen $var element = document.getElementById("moose-equation-44452653-dd08-4c80-9be2-873c99c1b8b5");katex.render("\\mu=0.1", element, {displayMode:false,throwOnError:false});$ ).

Release-rate over birth-rate (R/B)

The release-rate over birth-rate (R/B) is defined as $var element = document.getElementById("moose-equation-996f36fa-422e-42fb-a4ec-41a6318fef1e");katex.render(" (R/B) = \\frac{ R(t) }{ \\frac{4}{3}\\pi a^3 p(t) }", element, {displayMode:true,throwOnError:false});$ in terms of the ratio of release rate $var element = document.getElementById("moose-equation-987558ab-18ce-4028-b128-b1d0adbfd8d7");katex.render("R(t)", element, {displayMode:false,throwOnError:false});$ , the source $var element = document.getElementById("moose-equation-67eecbdd-aff4-4895-84fd-f3675ae9d883");katex.render("p(t)", element, {displayMode:false,throwOnError:false});$ , and the sphere volume $var element = document.getElementById("moose-equation-8ddbb9fb-acdb-4e42-ae21-9a1a36a8e7ca");katex.render("\\frac{4}{3}\\pi a^3", element, {displayMode:false,throwOnError:false});$ . The release rate, $var element = document.getElementById("moose-equation-f45bf046-8668-4ebe-aee2-62c92924948d");katex.render("R(t)", element, {displayMode:false,throwOnError:false});$ is computed in the code via

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [release_inc]
    type = SideIntegralMassFlux<<<{"description": "Computes the integrated mass flux over a given surface.", "href": "../../source/postprocessors/SideIntegralMassFlux.html"}>>>
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = conc
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = exterior
    arrhenius_prpty_name<<<{"description": "Specify the name of the material property"}>>> = diffusion_coef
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial timestep_end'
  []
[]

Thus, the (R/B) ratio for a long-lived (or stable; with $var element = document.getElementById("moose-equation-1b4047ad-e190-4255-a3d3-2b967cbab459");katex.render("\\lambda=0", element, {displayMode:false,throwOnError:false});$ ) fission product under the in-pile conditions becomes $var element = document.getElementById("moose-equation-cb8004dd-ca64-4add-a9de-4b5be841affc");katex.render(" (R/B) = 1 - 6 \\sum_{n=1}^\\infty \\frac{\\exp{\\left( -n^2\\pi^2 \\tau \\right)}}{n^2\\pi^2}", element, {displayMode:true,throwOnError:false});$

Figure 2 compares analytical and computed release-rate over birth-rate, (R/B) as a function of the diffusion time, $var element = document.getElementById("moose-equation-f59b3ac7-5d5f-420f-8b7d-246fe8c1f9d9");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ .

Figure 2: Comparison of analytical and computed release-rate over birth-rate.

Fractional release

The analytical expression for the fractional release, $var element = document.getElementById("moose-equation-49c21ff9-364b-44c4-a7b7-57454e6f6389");katex.render("F(t)", element, {displayMode:false,throwOnError:false});$ is expressed as: $var element = document.getElementById("moose-equation-9ccd7c0e-4674-4496-802f-8e0ef78bf053");katex.render(" F(\\tau) = 1 - 6 \\sum_{n=1}^\\infty \\frac{1-\\exp{\\left( -n^2\\pi^2 \\tau \\right)}}{n^4\\pi^4\\tau}", element, {displayMode:true,throwOnError:false});$

The fractional release is computed in the code via

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [x_released]
    type = FractionalRelease<<<{"description": "Compute the fractional release of fission products.", "href": "../../source/postprocessors/FractionalRelease.html"}>>>
    released<<<{"description": "Postprocessor to be used as the released value."}>>> = released
    retained<<<{"description": "Postprocessor to be used as the retained value."}>>> = retained
  []
[]

Figure 3 compares analytical and computed release fractions, $var element = document.getElementById("moose-equation-7e2cdc0d-7867-481c-80f9-1ad805d44f44");katex.render("F(\\tau)", element, {displayMode:false,throwOnError:false});$ as a function of diffusion time, $var element = document.getElementById("moose-equation-b769b152-7ded-45ce-8605-324df870e8e2");katex.render("\\tau", element, {displayMode:false,throwOnError:false});$ .

Figure 3: Comparison of analytical and computed total release fractions.

Summary

The Bison predictions agree well with the analytical results for concentration profiles, release rate, release-rate over birth-rate (R/B), and total fractional release out of a solid sphere under in-pile conditions for a stable fission product.

References

H. Nabielek, H. Hick, M. Wagner-Löffler, and E. H. Voice. Performance limits of coated particle fuel, Part III: Fission product migration in HTR fuel. Technical Report DP-Report-828(Pt.3), O.E.C.D High Temperature Reactor Project Dragon, 6 1974.

@techreport{NABIELEK1974,
    author = {Nabielek, H. and Hick, H. and Wagner-L\"{o}ffler, M. and Voice, E. H.},
    title = "Performance limits of Coated Particle Fuel, {Part III: F}ission Product Migration in {HTR} Fuel",
    institution = "O.E.C.D High Temperature Reactor Project Dragon",
    number = "DP-Report-828(Pt.3)",
    month = "6",
    place = "A.E.E. Winfrith, Dorchester, Dorset, England",
    year = "1974"
}

A. Toptan, W. Jiang, J. D. Hales, B. W. Spencer, and S. Novascone. Verification of Bison fission product species conservation under TRISO reactor conditions. Journal of Nuclear Materials, 573:154105, 2023. doi:10.1016/j.jnucmat.2022.154105.

@article{TOPTAN2023154105,
    author = "Toptan, A. and Jiang, W. and Hales, J. D. and Spencer, B. W. and Novascone, S.",
    title = "Verification of {B}ison fission product species conservation under {TRISO} reactor conditions",
    journal = "Journal of Nuclear Materials",
    volume = "573",
    pages = "154105",
    doi = "10.1016/j.jnucmat.2022.154105",
    year = "2023"
}

(assessment/verification/TRISO_diffusion/inpile_stable/inpile_stable.i)

#
# TRISO kernel in 1D spherical
#
# In-pile, a stable fission product
# -------------------------------
# Non-zero source, p/=0
# Stable fission product, decay_const=0
# Zero initial concentration, C0=0
#

timestep = 100

p = 1.0
C0 = 0.0
decay_const = 0

kernel_radius = 212.5e-6
Dprime = 1e-4

[GlobalParams]
  order = SECOND
  family = LAGRANGE
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen]
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 ${kernel_radius}'
    mesh_density = 100
    block_names = 'fuel'
    bias = 1
    dual_bias = 0.8
  []
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [conc]
    initial_condition = '${C0}'
  []
[]

[Kernels]
  [mass_dt]
    type = TimeDerivative
    variable = conc
  []
  [mass]
    type = MatDiffusion
    variable = conc
    diffusivity = diffusion_coef
    extra_vector_tags = 'ref'
  []
  [mass_source]
    type = BodyForce
    variable = conc
    value = '${p}'
    extra_vector_tags = 'ref'
  []
  [decay]
    type = Decay
    variable = conc
    radioactive_decay_constant = ${decay_const}
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  [freesurf_conc]
    type = DirichletBC
    variable = conc
    boundary = exterior
    value = 0.0
  []
[]

[Materials]
  [fuel_conc]
    type = GenericConstantMaterial
    block = fuel
    prop_names = diffusion_coef
    prop_values = '${fparse Dprime*kernel_radius^2}'
  []
[]

[Debug]
  show_var_residual_norms = true
  show_var_residual = 'conc'
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'

  line_search = 'none'

  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-20
  nl_max_its = 20

  l_tol = 1e-8
  l_max_its = 30

  start_time = 0.0
  end_time = '${fparse 1/Dprime}'
  num_steps = 100

  dtmax = '${timestep}'
  dtmin = '${timestep}'

  [Predictor]
    type = SimplePredictor
    scale = 1.0
  []
[]

[Functions]
  [analytic]
    type = ParsedFunction
    expression = '1 - x*x'
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
    execute_on = timestep_end
  []

  # Release rate
  [release_inc]
    type = SideIntegralMassFlux
    variable = conc
    boundary = exterior
    arrhenius_prpty_name = diffusion_coef
    execute_on = 'initial timestep_end'
  []

  [released]
    type = TimeIntegratedPostprocessor
    value = release_inc
    execute_on = 'initial timestep_end'
  []
  [total]
    type = ElementIntegralVariablePostprocessor
    variable = conc
    execute_on = 'initial timestep_end'
  []
  [x_released]
    type = FractionalRelease
    released = released
    retained = retained
  []
  [retained]
    type = ElementIntegralVariablePostprocessor
    variable = conc
  []

  [conc_point_25]
    type = PointValue
    variable = conc
    point = '${fparse 0.25*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_50]
    type = PointValue
    variable = conc
    point = '${fparse 0.50*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_75]
    type = PointValue
    variable = conc
    point = '${fparse 0.75*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_95]
    type = PointValue
    variable = conc
    point = '${fparse 0.95*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_25]
    type = ParsedPostprocessor
    expression = 'conc_point_25 / (p / (6 * Dprime))'
    pp_names = conc_point_25
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_50]
    type = ParsedPostprocessor
    expression = 'conc_point_50 / (p / (6 * Dprime))'
    pp_names = conc_point_50
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_75]
    type = ParsedPostprocessor
    expression = 'conc_point_75 / (p / (6 * Dprime))'
    pp_names = conc_point_75
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_95]
    type = ParsedPostprocessor
    expression = 'conc_point_95 / (p / (6 * Dprime))'
    pp_names = conc_point_95
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_25]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.25 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_50]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.50 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_75]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.75 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_95]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.95 0 0'
    execute_on = 'initial timestep_end'
  []
  [error_25]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_25 analytic_dimensionless_25'
    expression = '(bison_dimensionless_25 - analytic_dimensionless_25)/analytic_dimensionless_25*100'
    execute_on = 'initial timestep_end'
  []
  [error_50]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_50 analytic_dimensionless_50'
    expression = '(bison_dimensionless_50 - analytic_dimensionless_50)/analytic_dimensionless_50*100'
    execute_on = 'initial timestep_end'
  []
  [error_75]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_75 analytic_dimensionless_75'
    expression = '(bison_dimensionless_75 - analytic_dimensionless_75)/analytic_dimensionless_75*100'
    execute_on = 'initial timestep_end'
  []
  [error_95]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_95 analytic_dimensionless_95'
    expression = '(bison_dimensionless_95 - analytic_dimensionless_95)/analytic_dimensionless_95*100'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
  [console]
    type = Console
    show = 'analytic_dimensionless_25 analytic_dimensionless_50 analytic_dimensionless_75 analytic_dimensionless_95
            bison_dimensionless_25 bison_dimensionless_50 bison_dimensionless_75 bison_dimensionless_95
            error_25 error_50 error_75 error_95 released x_released'
  []
  [out]
    type = CSV
    execute_on = final
  []
[]

(assessment/verification/TRISO_diffusion/inpile_stable/inpile_stable.i)

#
# TRISO kernel in 1D spherical
#
# In-pile, a stable fission product
# -------------------------------
# Non-zero source, p/=0
# Stable fission product, decay_const=0
# Zero initial concentration, C0=0
#

timestep = 100

p = 1.0
C0 = 0.0
decay_const = 0

kernel_radius = 212.5e-6
Dprime = 1e-4

[GlobalParams]
  order = SECOND
  family = LAGRANGE
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen]
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 ${kernel_radius}'
    mesh_density = 100
    block_names = 'fuel'
    bias = 1
    dual_bias = 0.8
  []
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [conc]
    initial_condition = '${C0}'
  []
[]

[Kernels]
  [mass_dt]
    type = TimeDerivative
    variable = conc
  []
  [mass]
    type = MatDiffusion
    variable = conc
    diffusivity = diffusion_coef
    extra_vector_tags = 'ref'
  []
  [mass_source]
    type = BodyForce
    variable = conc
    value = '${p}'
    extra_vector_tags = 'ref'
  []
  [decay]
    type = Decay
    variable = conc
    radioactive_decay_constant = ${decay_const}
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  [freesurf_conc]
    type = DirichletBC
    variable = conc
    boundary = exterior
    value = 0.0
  []
[]

[Materials]
  [fuel_conc]
    type = GenericConstantMaterial
    block = fuel
    prop_names = diffusion_coef
    prop_values = '${fparse Dprime*kernel_radius^2}'
  []
[]

[Debug]
  show_var_residual_norms = true
  show_var_residual = 'conc'
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'

  line_search = 'none'

  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-20
  nl_max_its = 20

  l_tol = 1e-8
  l_max_its = 30

  start_time = 0.0
  end_time = '${fparse 1/Dprime}'
  num_steps = 100

  dtmax = '${timestep}'
  dtmin = '${timestep}'

  [Predictor]
    type = SimplePredictor
    scale = 1.0
  []
[]

[Functions]
  [analytic]
    type = ParsedFunction
    expression = '1 - x*x'
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
    execute_on = timestep_end
  []

  # Release rate
  [release_inc]
    type = SideIntegralMassFlux
    variable = conc
    boundary = exterior
    arrhenius_prpty_name = diffusion_coef
    execute_on = 'initial timestep_end'
  []

  [released]
    type = TimeIntegratedPostprocessor
    value = release_inc
    execute_on = 'initial timestep_end'
  []
  [total]
    type = ElementIntegralVariablePostprocessor
    variable = conc
    execute_on = 'initial timestep_end'
  []
  [x_released]
    type = FractionalRelease
    released = released
    retained = retained
  []
  [retained]
    type = ElementIntegralVariablePostprocessor
    variable = conc
  []

  [conc_point_25]
    type = PointValue
    variable = conc
    point = '${fparse 0.25*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_50]
    type = PointValue
    variable = conc
    point = '${fparse 0.50*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_75]
    type = PointValue
    variable = conc
    point = '${fparse 0.75*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_95]
    type = PointValue
    variable = conc
    point = '${fparse 0.95*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_25]
    type = ParsedPostprocessor
    expression = 'conc_point_25 / (p / (6 * Dprime))'
    pp_names = conc_point_25
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_50]
    type = ParsedPostprocessor
    expression = 'conc_point_50 / (p / (6 * Dprime))'
    pp_names = conc_point_50
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_75]
    type = ParsedPostprocessor
    expression = 'conc_point_75 / (p / (6 * Dprime))'
    pp_names = conc_point_75
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_95]
    type = ParsedPostprocessor
    expression = 'conc_point_95 / (p / (6 * Dprime))'
    pp_names = conc_point_95
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_25]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.25 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_50]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.50 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_75]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.75 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_95]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.95 0 0'
    execute_on = 'initial timestep_end'
  []
  [error_25]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_25 analytic_dimensionless_25'
    expression = '(bison_dimensionless_25 - analytic_dimensionless_25)/analytic_dimensionless_25*100'
    execute_on = 'initial timestep_end'
  []
  [error_50]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_50 analytic_dimensionless_50'
    expression = '(bison_dimensionless_50 - analytic_dimensionless_50)/analytic_dimensionless_50*100'
    execute_on = 'initial timestep_end'
  []
  [error_75]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_75 analytic_dimensionless_75'
    expression = '(bison_dimensionless_75 - analytic_dimensionless_75)/analytic_dimensionless_75*100'
    execute_on = 'initial timestep_end'
  []
  [error_95]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_95 analytic_dimensionless_95'
    expression = '(bison_dimensionless_95 - analytic_dimensionless_95)/analytic_dimensionless_95*100'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
  [console]
    type = Console
    show = 'analytic_dimensionless_25 analytic_dimensionless_50 analytic_dimensionless_75 analytic_dimensionless_95
            bison_dimensionless_25 bison_dimensionless_50 bison_dimensionless_75 bison_dimensionless_95
            error_25 error_50 error_75 error_95 released x_released'
  []
  [out]
    type = CSV
    execute_on = final
  []
[]

(assessment/verification/TRISO_diffusion/inpile_stable/inpile_stable.i)

#
# TRISO kernel in 1D spherical
#
# In-pile, a stable fission product
# -------------------------------
# Non-zero source, p/=0
# Stable fission product, decay_const=0
# Zero initial concentration, C0=0
#

timestep = 100

p = 1.0
C0 = 0.0
decay_const = 0

kernel_radius = 212.5e-6
Dprime = 1e-4

[GlobalParams]
  order = SECOND
  family = LAGRANGE
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen]
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 ${kernel_radius}'
    mesh_density = 100
    block_names = 'fuel'
    bias = 1
    dual_bias = 0.8
  []
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
[]

[Variables]
  [conc]
    initial_condition = '${C0}'
  []
[]

[Kernels]
  [mass_dt]
    type = TimeDerivative
    variable = conc
  []
  [mass]
    type = MatDiffusion
    variable = conc
    diffusivity = diffusion_coef
    extra_vector_tags = 'ref'
  []
  [mass_source]
    type = BodyForce
    variable = conc
    value = '${p}'
    extra_vector_tags = 'ref'
  []
  [decay]
    type = Decay
    variable = conc
    radioactive_decay_constant = ${decay_const}
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  [freesurf_conc]
    type = DirichletBC
    variable = conc
    boundary = exterior
    value = 0.0
  []
[]

[Materials]
  [fuel_conc]
    type = GenericConstantMaterial
    block = fuel
    prop_names = diffusion_coef
    prop_values = '${fparse Dprime*kernel_radius^2}'
  []
[]

[Debug]
  show_var_residual_norms = true
  show_var_residual = 'conc'
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'

  line_search = 'none'

  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-20
  nl_max_its = 20

  l_tol = 1e-8
  l_max_its = 30

  start_time = 0.0
  end_time = '${fparse 1/Dprime}'
  num_steps = 100

  dtmax = '${timestep}'
  dtmin = '${timestep}'

  [Predictor]
    type = SimplePredictor
    scale = 1.0
  []
[]

[Functions]
  [analytic]
    type = ParsedFunction
    expression = '1 - x*x'
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
    execute_on = timestep_end
  []

  # Release rate
  [release_inc]
    type = SideIntegralMassFlux
    variable = conc
    boundary = exterior
    arrhenius_prpty_name = diffusion_coef
    execute_on = 'initial timestep_end'
  []

  [released]
    type = TimeIntegratedPostprocessor
    value = release_inc
    execute_on = 'initial timestep_end'
  []
  [total]
    type = ElementIntegralVariablePostprocessor
    variable = conc
    execute_on = 'initial timestep_end'
  []
  [x_released]
    type = FractionalRelease
    released = released
    retained = retained
  []
  [retained]
    type = ElementIntegralVariablePostprocessor
    variable = conc
  []

  [conc_point_25]
    type = PointValue
    variable = conc
    point = '${fparse 0.25*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_50]
    type = PointValue
    variable = conc
    point = '${fparse 0.50*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_75]
    type = PointValue
    variable = conc
    point = '${fparse 0.75*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [conc_point_95]
    type = PointValue
    variable = conc
    point = '${fparse 0.95*kernel_radius} 0 0'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_25]
    type = ParsedPostprocessor
    expression = 'conc_point_25 / (p / (6 * Dprime))'
    pp_names = conc_point_25
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_50]
    type = ParsedPostprocessor
    expression = 'conc_point_50 / (p / (6 * Dprime))'
    pp_names = conc_point_50
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_75]
    type = ParsedPostprocessor
    expression = 'conc_point_75 / (p / (6 * Dprime))'
    pp_names = conc_point_75
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [bison_dimensionless_95]
    type = ParsedPostprocessor
    expression = 'conc_point_95 / (p / (6 * Dprime))'
    pp_names = conc_point_95
    constant_names = 'p Dprime'
    constant_expressions = '${p} ${Dprime}'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_25]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.25 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_50]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.50 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_75]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.75 0 0'
    execute_on = 'initial timestep_end'
  []
  [analytic_dimensionless_95]
    type = FunctionValuePostprocessor
    function = analytic
    point = '0.95 0 0'
    execute_on = 'initial timestep_end'
  []
  [error_25]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_25 analytic_dimensionless_25'
    expression = '(bison_dimensionless_25 - analytic_dimensionless_25)/analytic_dimensionless_25*100'
    execute_on = 'initial timestep_end'
  []
  [error_50]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_50 analytic_dimensionless_50'
    expression = '(bison_dimensionless_50 - analytic_dimensionless_50)/analytic_dimensionless_50*100'
    execute_on = 'initial timestep_end'
  []
  [error_75]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_75 analytic_dimensionless_75'
    expression = '(bison_dimensionless_75 - analytic_dimensionless_75)/analytic_dimensionless_75*100'
    execute_on = 'initial timestep_end'
  []
  [error_95]
    type = ParsedPostprocessor
    pp_names = 'bison_dimensionless_95 analytic_dimensionless_95'
    expression = '(bison_dimensionless_95 - analytic_dimensionless_95)/analytic_dimensionless_95*100'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  time_step_interval =  1
  exodus = false
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
  [console]
    type = Console
    show = 'analytic_dimensionless_25 analytic_dimensionless_50 analytic_dimensionless_75 analytic_dimensionless_95
            bison_dimensionless_25 bison_dimensionless_50 bison_dimensionless_75 bison_dimensionless_95
            error_25 error_50 error_75 error_95 released x_released'
  []
  [out]
    type = CSV
    execute_on = final
  []
[]

Concentration
Release - rate over birth - rate ( R / B )
Fractional release
Summary
References