GraphiteMatrixThermal

Computes thermal conductivity (W/m/K) and specific heat capacity (J/kg/K) for graphite matrix.

Description

A list of various graphite grades is provided in Table 1 along with the manufacturing descriptions. See Toptan et al. (2021) for more about the model implementations and Toptan et al. (2021) for homogenization techniques applied to the thermal conductivity for non-zero packing fraction.

Table 1: Graphite types.

Grade	Coke source	Binder	Forming process	Vendor
A3-3		Phenolic resin	Quasi-isostatic cold molding	Hobeg
A3-27		Synthesized resin 2	Quasi-isostatic cold molding	Hobeg
H-451	Petroleum	Coal-tar pitch	Extrusion	Great Lakes Carbon Company
IG-110	Petroleum	Coal-tar pitch	Isostatic pressing	ToyoTanso, Japan
PCEA	Petroleum		Extrusion	Graftech, USA
2020	Petroleum		Molded	Mersen

Thermal Conductivity

The thermal conductivity of the matrix material—used in both the fueled and fuel-free regions of TRISO fuel compacts/pebbles—is given for various graphite grades. The thermal conductivity of the fueled zone is based on a homogenization of TRISO particles in the matrix to produce an effective thermal conductivity (Toptan et al., 2021). The mathematical form of the correlation derived for the effective thermal conductivity applies to both the fueled (non-zero packing fraction) and fuel-free regions (zero packing fraction).

Grade H-451 graphite

The thermal conductivity of grade H-451 graphite is formulated according to NEA Nuclear Science Committee (2018) as a function of neutron fluence and temperature, which is considered valid over the temperature range $var element = document.getElementById("moose-equation-eed5c1f4-2538-4d63-8a81-b8e7dfeac5fb");katex.render("[", element, {displayMode:false,throwOnError:false});$ 500 K, 1800 K $var element = document.getElementById("moose-equation-15461123-55b9-444f-9d5e-7a1c3dc255db");katex.render("]", element, {displayMode:false,throwOnError:false});$ . The mathematical formulations for the temperature-dependent thermal conductivity at various neutron fluence values are tabulated in Table 2.

Table 2: Thermal conductivity of grade H-451 graphite as a function of neutron fluence ( $var element = document.getElementById("moose-equation-8610e6eb-263c-49c3-aa44-2f86943f273b");katex.render("\\times 10^{25}", element, {displayMode:false,throwOnError:false});$ n/m $var element = document.getElementById("moose-equation-aa4844f8-9fa9-46ce-8eca-e26cd5823fe6");katex.render("^2", element, {displayMode:false,throwOnError:false});$ ) and temperature (K).

Fluence ( $var element = document.getElementById("moose-equation-b909dfc4-bfb3-4ec8-beaf-ffa43d454375");katex.render("\\times 10^{25}", element, {displayMode:false,throwOnError:false});$ n/m $var element = document.getElementById("moose-equation-715ffdad-99cf-4e3f-a4c8-1f0b23dfa300");katex.render("^2", element, {displayMode:false,throwOnError:false});$ )	Thermal conductivity, $var element = document.getElementById("moose-equation-be97c9fe-64e7-47f8-86ed-23ef12efdbaf");katex.render("k", element, {displayMode:false,throwOnError:false});$ (W/m-K)
0.0 (unirradiated)	$var element = document.getElementById("moose-equation-2bd585bc-bee5-41d3-8851-d3cd1064b7fc");katex.render("k=3.28248\\times 10^{-5} T_K^2 - 1.24890\\times 10^{-1} T_K + 1.692145\\times 10^{2}", element, {displayMode:false,throwOnError:false});$
0.2	$var element = document.getElementById("moose-equation-d70b0fcd-3fb6-49ae-af85-0d4b5af9ceac");katex.render("k= 2.63319\\times 10^{-9} T_K^3 + 1.34998\\times 10^{-5} T_K^2 - 8.03309\\times 10^{-2} T_K + 1.31991\\times 10^{2}", element, {displayMode:false,throwOnError:false});$
0.5	$var element = document.getElementById("moose-equation-04d47137-b592-4db1-a406-5bda4e3fde9b");katex.render("k= 3.60774\\times 10^{-9} T_K^3 + 6.45202\\times 10^{-6} T_K^2 - 6.38157\\times 10^{-2} T_K + 1.17459\\times 10^{2}", element, {displayMode:false,throwOnError:false});$
1.0	$var element = document.getElementById("moose-equation-b5fe456f-16db-4b76-b716-5b21517e8b57");katex.render("k= 4.47613\\times 10^{-9} T_K^3 - 1.86131\\times 10^{-6} T_K^2 - 4.14540\\times 10^{-2} T_K + 9.55394\\times 10^{1}", element, {displayMode:false,throwOnError:false});$
3.0	$var element = document.getElementById("moose-equation-e8b12850-51ea-4403-9126-c611d41d25c6");katex.render("k= 4.19346\\times 10^{-6} T_K^2 - 2.13523\\times 10^{-2} T_K + 5.41993\\times 10^{1}", element, {displayMode:false,throwOnError:false});$
8.0	$var element = document.getElementById("moose-equation-7d907eda-d82a-46cb-afdb-7fd283deced2");katex.render("k= 3.63529\\times 10^{-6} T_K^2 - 1.89383\\times 10^{-2} T_K + 5.04684\\times 10^{1}", element, {displayMode:false,throwOnError:false});$

Grade 2020 graphite

The thermal conductivity of grade 2020 graphite is formulated according to NEA Nuclear Science Committee (2018) as a function of neutron fluence and temperature, which is considered valid over the temperature range $var element = document.getElementById("moose-equation-c81b605b-0629-40e4-bd75-bd66cc58e458");katex.render("[", element, {displayMode:false,throwOnError:false});$ 300 K, 2000 K $var element = document.getElementById("moose-equation-9f3d47c5-90a4-4260-9bf1-d88cf33db5b9");katex.render("]", element, {displayMode:false,throwOnError:false});$ . The mathematical formulations for the temperature-dependent thermal conductivity at various neutron fluence values are tabulated in Table 3.

Table 3: Thermal conductivity of grade 2020 graphite as a function of neutron fluence ( $var element = document.getElementById("moose-equation-916e6190-80eb-4058-9d61-2b058e95157f");katex.render("\\times 10^{22}", element, {displayMode:false,throwOnError:false});$ n/m $var element = document.getElementById("moose-equation-f0c13266-481d-4b69-8d4a-a15c4e755a78");katex.render("^2", element, {displayMode:false,throwOnError:false});$ ) and temperature (K).

Grade A3-3 and A3-27 graphite

The thermal conductivity of the matrix material—used in both the fueled and fuel-free regions of TRISO fuel compacts/pebbles—is given by Miller et al. (2018) and Gontard and Nabielek (1990) for grade A3-3 and A3-27 graphite materials. A correlation is given for fast fluence with neutron energy threshold $var element = document.getElementById("moose-equation-9f084b95-3ca2-44a9-9ff8-0da8c347cb31");katex.render("E>", element, {displayMode:false,throwOnError:false});$ 0.18 MeV.

The thermal conductivity, $var element = document.getElementById("moose-equation-6b85e775-2777-4d80-9afa-58c20f71203a");katex.render("k_{graphite}", element, {displayMode:false,throwOnError:false});$ (W/m-K) of the graphite matrix is $var element = document.getElementById("moose-equation-471ae913-838a-4b29-87a0-c2bc44521b81");katex.render(" k_{graphite} = k_{unirr} \\kappa_{\\phi} \\kappa_{\\rho},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-8f724834-5004-4201-ae35-7933f7bad432");katex.render("k_{unirr}", element, {displayMode:false,throwOnError:false});$ (W/m-K) is the temperature-dependent thermal conductivity of the unirradiated matrix material, $var element = document.getElementById("moose-equation-3c09c44e-8271-47db-b31f-6fc04b1f0689");katex.render("\\kappa_{\\phi}", element, {displayMode:false,throwOnError:false});$ (dimensionless) is a correction for irradiation damage, and $var element = document.getElementById("moose-equation-490bdc76-f687-45b5-9c6a-f9112a80d4b7");katex.render("\\kappa_{\\rho}", element, {displayMode:false,throwOnError:false});$ (dimensionless) is a correction factor for densities other than the density of the reference material (1.7 g/cm $var element = document.getElementById("moose-equation-adbe28a6-8999-429c-8fde-fd275b011114");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ ).

The thermal conductivity $var element = document.getElementById("moose-equation-dd88a4e7-a0da-42cf-bc8e-7c526e7245b4");katex.render("k_{unirr}", element, {displayMode:false,throwOnError:false});$ (W/m-K) of the unirradiated matrix depends on the nature of the matrix material (A3-3 or A3-27) and the temperature of the final heat treatment during fabrication (1800 or 1950 $var element = document.getElementById("moose-equation-1a49595f-8313-4005-9713-c5e6746db8eb");katex.render("\\degree{C}", element, {displayMode:false,throwOnError:false});$ ) (Gontard and Nabielek, 1990):

$var element = document.getElementById("moose-equation-5b4c142c-023a-4a46-8902-82b797d0ce60");katex.render(" k_{unirr} = k_{100} \\left(1 - \\alpha (T-100) e^{\\delta T} \\right),", element, {displayMode:true,throwOnError:false});$ where the $var element = document.getElementById("moose-equation-4547bb4b-ea1c-4c04-9190-86cad9d4d87c");katex.render("k_{100}", element, {displayMode:false,throwOnError:false});$ (W/m-K) is the thermal conductivity of material at 100 $var element = document.getElementById("moose-equation-13d00b24-7249-4995-8cd5-08fd2f959a13");katex.render("\\degree{C}", element, {displayMode:false,throwOnError:false});$ , the non-dimenstional, empirical coefficients, $var element = document.getElementById("moose-equation-63162593-8cc1-4e71-91e9-475731287e70");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-63525237-8bd3-4a25-be30-d72cd87f3fa2");katex.render("\\delta", element, {displayMode:false,throwOnError:false});$ , are tabulated in Table 4, and $var element = document.getElementById("moose-equation-88a6d760-13e5-4143-b453-d67cd0e5be89");katex.render("T", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-2d2d06ad-613e-4c5e-a4a3-4b933b42a193");katex.render("\\degree{C}", element, {displayMode:false,throwOnError:false});$ ) is the temperature of the matrix. The material A3-3 heat-treated at 1800 $var element = document.getElementById("moose-equation-71d861d5-9926-4308-a0c2-4d5013a55b14");katex.render("\\degree{C}", element, {displayMode:false,throwOnError:false});$ is used by default.

Table 4: Coefficients for unirradiated thermal conductivity (Gontard and Nabielek, 1990).

Material	$var element = document.getElementById("moose-equation-f2147699-fc2a-4f06-a91f-e5809d5c1f4b");katex.render("T (\\degree{C})", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-8f834147-da84-404d-8da1-70590d3498d2");katex.render("k_{100}", element, {displayMode:false,throwOnError:false});$ (W/m-K)	$var element = document.getElementById("moose-equation-11a63d81-e00d-4403-afa2-431cdd250b0b");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-aa21cbf0-5c3a-4856-ac9e-45ebb4f074c8");katex.render("\\delta", element, {displayMode:false,throwOnError:false});$
A3-3	1800	50.8	1.1810 $var element = document.getElementById("moose-equation-b9d01d7b-2006-47f5-b038-f25d4f7507a4");katex.render("\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	-7.8453 $var element = document.getElementById("moose-equation-db1eba1d-05d9-4d5b-b5a2-7d35db1af8ae");katex.render("\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$
A3-3	1950	64.6	1.4079 $var element = document.getElementById("moose-equation-ffd6ffe9-e4d3-4bfe-a974-f057f18b00c8");katex.render("\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	-9.0739 $var element = document.getElementById("moose-equation-2ff95d1f-61ce-4241-9836-d4a95076a0a6");katex.render("\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$
A3-27	1800	47.4	9.7556 $var element = document.getElementById("moose-equation-ae7e209f-3161-4105-b743-2261a8339321");katex.render("\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$	-6.0360 $var element = document.getElementById("moose-equation-0f3beea4-0847-4331-ba9f-f29af10cdbb4");katex.render("\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$
A3-27	1950	62.2	1.4621 $var element = document.getElementById("moose-equation-c0f8b74f-eeec-42a0-9349-6c9aec32ed47");katex.render("\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$	-9.6050 $var element = document.getElementById("moose-equation-77fc0c55-0d9c-4724-9766-90335a3ea31c");katex.render("\\times 10^{-4}", element, {displayMode:false,throwOnError:false});$

Neutron fluence decreases the thermal conductivity of the matrix. A correction factor $var element = document.getElementById("moose-equation-28fae780-c108-4354-9c09-d1323b3e1f28");katex.render("\\kappa_{\\phi}", element, {displayMode:false,throwOnError:false});$ (dimensionless) based on experimental data is applied to take the effect of irradiation damage into account (Gontard and Nabielek, 1990), which is given by

$(1)var element = document.getElementById("moose-equation-ad0c18a1-efbd-4a9f-ac1f-1a89e5931ff5");katex.render("\\kappa_{\\phi} = 1.0 - \\left(0.94 - 0.604t_C\\right) \\left(1 - \\exp\\left[-(2.96 - 1.955t_C) \\frac{\\phi}{1.52}\\right] \\right) - \\left(0.043t_C- 0.008t_C^8 \\right) \\frac{\\phi}{1.52}, ", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-a36b921f-8349-42be-b280-03aa98fb4108");katex.render("t_C=T/1000", element, {displayMode:false,throwOnError:false});$ in terms of the temperature of the matrix, $var element = document.getElementById("moose-equation-0685daba-541c-4436-8c72-209fc2114c9a");katex.render("T", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-89b1865e-a502-47e1-a37c-66e51f7b7f61");katex.render("\\degree{C}", element, {displayMode:false,throwOnError:false});$ ), and $var element = document.getElementById("moose-equation-c6b7f56a-605b-425e-b52a-477732e9a9e0");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ (10 $var element = document.getElementById("moose-equation-1ec0bbf2-6d09-4d14-b61d-cac412cf2125");katex.render("^{25}", element, {displayMode:false,throwOnError:false});$ n/m $var element = document.getElementById("moose-equation-dce0cedd-3d87-49bb-8dab-d2f70813985a");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-e993f39e-eefd-4549-96c9-c4350e661ba9");katex.render("E>", element, {displayMode:false,throwOnError:false});$ 0.18 MeV) is the fast neutron fluence. Note that in the original reference (Gontard and Nabielek, 1990), $var element = document.getElementById("moose-equation-8f2c9657-3ff3-4ea9-bf00-1260c0b3b22e");katex.render("T", element, {displayMode:false,throwOnError:false});$ is mistakenly used instead of $var element = document.getElementById("moose-equation-450a6f60-7271-4e51-aa0f-1e6f001534f0");katex.render("\\Gamma", element, {displayMode:false,throwOnError:false});$ in the first instance of $var element = document.getElementById("moose-equation-85a660b9-40d5-45a7-b80b-4b7bc8421fda");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ in Eq. (1).

Finally, a correction factor $var element = document.getElementById("moose-equation-556abb2e-cde8-4df0-bc5c-22c30026a704");katex.render("\\kappa_{\\rho}", element, {displayMode:false,throwOnError:false});$ (dimensionless) is applied for densities other than the density of the reference material (1700 kg/m $var element = document.getElementById("moose-equation-13d5b32c-6347-499c-9f43-c09bec41b636");katex.render("^3", element, {displayMode:false,throwOnError:false});$ ) used to obtain the experimental data on unirradiated material. It is given by Miller et al. (2018):

$var element = document.getElementById("moose-equation-1ac8fcf0-f0c2-423c-a7ff-4d0be079bc42");katex.render(" \\kappa_{\\rho} = \\frac{\\rho}{1700},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-6ef0a86b-4501-4637-acfc-3ca93a71c88b");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ (kg/m $var element = document.getElementById("moose-equation-926c5fcb-625b-4b2c-bc22-1b7a823efd66");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ ) is the density of the matrix (initial_matrix_density). The parameter initial_matrix_density may differ greatly from the compact density tracked by the "density" material property in fueled regions of the compact.

Grade IG-110 graphite

note:IG-110 thermal conductivity

The original reference, JAEA report (Shibata et al., 2009) provides a model for the irradiated IG-110 thermal conductivity; however, the model given in the report has lack of information to regenerate the model behavior as shown in Figure 1.1.2 of the report. Therefore, the model was not integrated into the code. Additionally, the report doesn't provide a model for the unirradiated IG-110 thermal conductivity, a plot provided instead (see Figure 1.3.11 in the JAEA report). For this reason, the model behavior is captured via a polynomial function that is fit to the digitized data from Figure 1.3.11.

The unirradiated thermal conductivity of IG-110 graphite is digitized from Figure 1.3.11 in the JAEA report (Shibata et al., 2009) and fitted to a third-order polynomial function as follows: $var element = document.getElementById("moose-equation-bd0dfe0b-4bf7-4dfb-9778-bc517f6c0e1c");katex.render(" k = -2.13817804\\times 10^{-8} T_C^3 + 9.70154119\\times 10^{-5} T_C^2 - 0.154796557 T_C + 127.996024", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-3601a6ae-7d20-49b6-8359-b5b2ca1979cb");katex.render("k", element, {displayMode:false,throwOnError:false});$ is the thermal conductivity in W/m-K and $var element = document.getElementById("moose-equation-c4f96604-c542-4578-98e6-2d4b317c878d");katex.render("T_C", element, {displayMode:false,throwOnError:false});$ is the temperature in $var element = document.getElementById("moose-equation-acbc515c-136a-422c-8efc-b4e379376cbb");katex.render("\\degree", element, {displayMode:false,throwOnError:false});$ C.

Grade PCEA graphite

The irradiated thermal conductivity of PCEA graphite, $var element = document.getElementById("moose-equation-2bf233a3-3ee4-46f9-916c-9c3faed51333");katex.render("k_{irr}", element, {displayMode:false,throwOnError:false});$ (W/m-K) is from Hawkes et al. (2015) and Snead and Burchell (1995) as: $var element = document.getElementById("moose-equation-e899a524-a176-49d8-940c-b52b874ca1b2");katex.render(" \\frac{k_{irr}}{k_0} = \\left( 0.25 - 0.00017 T_C \\right) A \\log{(dpa)} + 0.000683 T_C", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-cf1daf86-e6ae-4c96-b3fa-8067e7a5730c");katex.render("k_0", element, {displayMode:false,throwOnError:false});$ is the unirradiated PCEA graphite, $var element = document.getElementById("moose-equation-76db3592-dfeb-40dd-a769-270876acb98d");katex.render("T_C", element, {displayMode:false,throwOnError:false});$ is the temperature in $var element = document.getElementById("moose-equation-8f28bd6a-2597-4186-825f-77cfe09e8192");katex.render("\\degree", element, {displayMode:false,throwOnError:false});$ C, and $var element = document.getElementById("moose-equation-38b9320d-850c-401a-a02d-409bf307237c");katex.render("dpa", element, {displayMode:false,throwOnError:false});$ is the number of displacements per carbon atom from fast neutrons. The multiplier used to convert fast fluence ( $var element = document.getElementById("moose-equation-4ae8af4d-98e2-48e5-9312-779882463c98");katex.render(">", element, {displayMode:false,throwOnError:false});$ 0.18 MeV) to dpa is 8.23 $var element = document.getElementById("moose-equation-84246b92-08dc-474e-9345-befebb28d78a");katex.render("\\times 10^{-26}", element, {displayMode:false,throwOnError:false});$ dpa/(n/m $var element = document.getElementById("moose-equation-6ae56c7d-f4a3-4f9d-ab33-1345e9e9be8c");katex.render("^2", element, {displayMode:false,throwOnError:false});$ ) (Hawkes et al., 2015). Here, the thermal conductivity of unirradiated PCEA graphite, $var element = document.getElementById("moose-equation-91b434a5-43dc-4f1e-b682-003601b9f036");katex.render("k_0", element, {displayMode:false,throwOnError:false});$ (W/m-K) is taken from Srinivasan and others (2021) and Swank et al. (2018), which is expressed as: $var element = document.getElementById("moose-equation-2cbcbfc5-16a6-49fb-8e6b-eec4a444d173");katex.render(" k_0 = 183.81 - 0.204 T_C + 1.0\\times 10^{-4} T_C^2", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-ef586489-14b1-484a-afab-8d5859ffc50b");katex.render("T_C", element, {displayMode:false,throwOnError:false});$ is the temperature in $var element = document.getElementById("moose-equation-0ea42bdd-deb6-40ad-b52f-7c4b0a0f0fa5");katex.render("\\degree", element, {displayMode:false,throwOnError:false});$ C.

Homogenization

The presence of TRISO particles in the matrix modifies its thermal conductivity. The homogenized (or effective) thermal conductivity of the graphite matrix containing TRISO particles ( $var element = document.getElementById("moose-equation-7a3f1884-e66b-4c56-b3c2-0040f629c196");katex.render("k_{homogenized}", element, {displayMode:false,throwOnError:false});$ ) is computed according to seven modeling options available in the code (see Table 5). The models are formulated in terms of the graphite thermal conductivity ( $var element = document.getElementById("moose-equation-d7e1bd57-5adb-40e0-a0b7-f86d70428540");katex.render("k_{graphite}", element, {displayMode:false,throwOnError:false});$ ), the thermal conductivity of a TRISO particle ( $var element = document.getElementById("moose-equation-6bc7de76-5f60-4724-95eb-2d9b951ad12c");katex.render("k_{particle}", element, {displayMode:false,throwOnError:false});$ ), and the packing fraction of TRISO particles within the graphite matrix ( $var element = document.getElementById("moose-equation-13208288-d14c-47e7-ab8a-4911316a1617");katex.render("PF", element, {displayMode:false,throwOnError:false});$ ). The packing fraction is defined as the ratio of the volume of the TRISO particles to the volume of the embedding matrix: $var element = document.getElementById("moose-equation-b4d50dc6-dae1-468e-8ad5-98116d3dd929");katex.render(" PF = \\frac{N_p V_p}{V_{fr}},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-5b27a055-e383-4d74-ae81-dd3ddc7e2c99");katex.render("N_p", element, {displayMode:false,throwOnError:false});$ (dimensionless) is the number of TRISO particles embedded in the matrix fueled region, $var element = document.getElementById("moose-equation-02f398a5-675b-4c7a-991b-e4f07be80445");katex.render("V_p", element, {displayMode:false,throwOnError:false});$ (cm $var element = document.getElementById("moose-equation-890d0b08-db39-4795-b81d-eb182dbdff96");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ ) is the average volume of a TRISO particle, and $var element = document.getElementById("moose-equation-d64fc1aa-9353-4199-90c6-877a5922acd0");katex.render("V_{fr}", element, {displayMode:false,throwOnError:false});$ (cm $var element = document.getElementById("moose-equation-6026cf40-f290-4751-aa8c-a602c04a175f");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ ) is the volume of the matrix fueled region (i.e., the volume containing TRISO particles). If the fuel element (compact or pebble) contains a fuel-free region, the volume of the fuel-free region is not included in the calculation of the packing fraction. The packing fraction is a fuel specification (i.e., a user input).

Table 5: Homogenization modeling options. $var element = document.getElementById("moose-equation-7dbba9b0-4fff-41b7-b32a-4e4a1702af26");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-16a3b415-044b-40ac-a0fe-ef720409ba68");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ depend on the thermal conductivity of the TRISO particles ( $var element = document.getElementById("moose-equation-5340bd82-ba41-4e55-be62-5f04132d7f2a");katex.render("k_{particle}", element, {displayMode:false,throwOnError:false});$ ), and the graphite matrix ( $var element = document.getElementById("moose-equation-7a482edd-de00-4d03-a43e-d26c9f042d90");katex.render("k_{graphite}", element, {displayMode:false,throwOnError:false});$ ), are given by $var element = document.getElementById("moose-equation-b43f978d-e27b-4278-bfe3-e700b79fc6d0");katex.render("\\alpha = k_{particle} / k_{graphite}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-05c02eb2-f571-44ff-ae5a-86e1339e14ec");katex.render("\\beta = (\\alpha -1)/(\\alpha + 2)", element, {displayMode:false,throwOnError:false});$ .

Model	Reference	Formulation
D-EMT	McLaughlin (1977)	$var element = document.getElementById("moose-equation-eab14254-035f-4efe-a615-86f24ba5d61e");katex.render("\\left( \\dfrac{k_{homogenized}-k_{particle}}{k_{graphite}-k_{particle}} \\right)^3 \\dfrac{k_{graphite}}{k_{homogenized}} = (1-PF)^3", element, {displayMode:false,throwOnError:false});$
Chiew–Glandt	Chiew and Glandt (1983)	$var element = document.getElementById("moose-equation-54429ced-6284-45fa-8e10-235604a2fce7");katex.render("\\dfrac{k_{homogenized}}{k_{graphite}} = \\dfrac{1 + 2 \\beta PF + (2 \\beta^3 - 0.1 \\beta) PF^2 + 0.05 PF^3 e^{4.5\\beta}}{1 - \\beta PF}", element, {displayMode:false,throwOnError:false});$
Maxwell	Maxwell (1954)	$var element = document.getElementById("moose-equation-26adb854-70c0-4443-a8c8-7ab6295ae4db");katex.render("\\dfrac{k_{homogenized}}{k_{graphite}} = 1 + \\dfrac{3 PF}{\\left( \\dfrac{\\alpha+2}{\\alpha-1} \\right) - PF }", element, {displayMode:false,throwOnError:false});$
Reduced Maxwell		$var element = document.getElementById("moose-equation-a1104f89-c29c-489c-8759-52ce0ceee0ca");katex.render("\\dfrac{k_{homogenized}}{k_{graphite}} = \\dfrac{1-PF}{1+0.5PF}", element, {displayMode:false,throwOnError:false});$
Bruggeman (or EMT)	Bruggeman (1935)	$var element = document.getElementById("moose-equation-8ed679f1-ec55-4f59-b694-cc67244740de");katex.render("\\dfrac{k_{homogenized}}{k_{graphite}} = \\alpha A + \\sqrt{\\alpha^2 A^2 + \\dfrac{\\alpha}{2}}\\quad \\text{with}\\quad A = \\dfrac{1}{4} \\left( 3PF -1 + \\dfrac{1}{\\alpha} \\left[2-3PF\\right] \\right)", element, {displayMode:false,throwOnError:false});$
Series (or harmonic mean)		$var element = document.getElementById("moose-equation-e046b034-d0c4-4096-9457-7f9714cc6b2b");katex.render("\\dfrac{k_{homogenized}}{k_{graphite}} = \\dfrac{1}{1+PF \\left( \\dfrac{1}{\\alpha} -1 \\right)}", element, {displayMode:false,throwOnError:false});$
Parallel (or volume average)		$var element = document.getElementById("moose-equation-b133ecae-e2e7-48c2-b288-aaf638f867a5");katex.render("\\dfrac{k_{homogenized}}{k_{graphite}} = 1+PF \\left( \\alpha -1 \\right)", element, {displayMode:false,throwOnError:false});$

The default modeling option, the differential effective medium theory (D-EMT) method, is recommended based on our analyses in Toptan et al. (2021). The D-EMT yields a cubic equation and we numerically solve this cubic equation according to the subroutine given in Toptan et al. (2020). The thermal conductivity of a TRISO particle was estimated based on its geometry and on the thermal conductivities of all its constituent layers (Folsom et al., 2015; Stainsby and others, 2009; Toptan et al., 2021). A value of 4.13 W/m-K has been adopted.

Specific Heat Capacity

A list of specific heat capacity correlations for graphite and its various grades are provided in Table 6. The specific heat capacity of grade 2020 graphite is considered to be the same as of grade H-451 graphite.

Table 6: Specific heat capacity modeling options. $var element = document.getElementById("moose-equation-24052813-490b-4287-bdf6-547722a2f548");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ (kg/m $var element = document.getElementById("moose-equation-590366de-b3ea-4e70-a0a3-fa841fd74c84");katex.render("^3", element, {displayMode:false,throwOnError:false});$ ) is the density of the graphite, $var element = document.getElementById("moose-equation-78971fd4-dfad-4545-a875-6534ef123345");katex.render("T_C", element, {displayMode:false,throwOnError:false});$ is the temperature in Celcius, and $var element = document.getElementById("moose-equation-4d617c49-39fa-449c-abf9-c5fe11dff7de");katex.render("T_K", element, {displayMode:false,throwOnError:false});$ is the temperature in Kelvin.

Model	Reference	Formulation, $var element = document.getElementById("moose-equation-6b0459b3-cb6b-48e8-b411-c401ee5e789b");katex.render("\\bm{c_P}", element, {displayMode:false,throwOnError:false});$ (J/kg-K)	Range of applicability
Carbon brick	IAEA (2000)	$var element = document.getElementById("moose-equation-2cafe8ae-a2ba-45f1-8d3d-cc2c4ee78820");katex.render("c_P= \\dfrac{1.75\\times 10^6}{\\rho} \\left(0.645 + 3.14 \\left[ \\dfrac{T_C}{1000} \\right] - 2.809 \\left[ \\dfrac{T_C}{1000} \\right]^2 + 0.959 \\left[ \\dfrac{T_C}{1000} \\right]^3 \\right)", element, {displayMode:false,throwOnError:false});$	up to 1200 $var element = document.getElementById("moose-equation-5645dbf2-6b4b-4482-a70b-cc7c77fdefe7");katex.render("\\degree{C}", element, {displayMode:false,throwOnError:false});$
H-451	Butland and Maddison (1973)	$var element = document.getElementById("moose-equation-b3fe4434-2ee9-470e-a97d-37e17cce1eed");katex.render("c_P= 4184 \\left( 0.54212 - 2.42667\\times 10^{-6}T_K -\\dfrac{90.2725}{T_K} -\\dfrac{43449.3}{T_K^2} + \\dfrac{1.59309\\times 10^7}{T_K^3} - \\dfrac{1.43688\\times 10^9}{T_K^4} \\right)", element, {displayMode:false,throwOnError:false});$	250K to 3000K
IG-110	Srinivasan and others (2021) and Maruyama et al. (1995)	$var element = document.getElementById("moose-equation-057d16d4-7306-4749-a38c-45aa7d34ba44");katex.render("c_P= \\left( 0.6572 + 0.0026 T_C - 2.0\\times 10^{-6} T_C^2 \\right) / 0.001", element, {displayMode:false,throwOnError:false});$	up to 800 $var element = document.getElementById("moose-equation-c4a6dd45-0479-4ceb-99d1-f4b77005a7e2");katex.render("\\degree", element, {displayMode:false,throwOnError:false});$ C

If the temperature is out of valid temperature boundaries, the specific heat capacity value at the temperature boundary is used. Note that the input parameter specific_heat_scale_factor can be used to scale the specific heat capacity.

Homogenization

The presence of TRISO particles affects the specific heat capacity of a fueled matrix. A simple volumetric-average homogenization method is applied first to the particle and then to the fueled matrix. The TRISO particle homogenized specific heat may be supplied. The default value of 748.72 J/kg-K is provided from Cho et al. (2009) where the particle characteristics in Table 7 are volumetric-average homogenized.

Table 7: Default TRISO Particle Characteristics for Homogenization

Material	Radius (cm)	Specific Heat (J/kg-K)
UO $var element = document.getElementById("moose-equation-d9158df9-99d7-45aa-878c-f227beb2d185");katex.render("_2", element, {displayMode:false,throwOnError:false});$ Kernel	0.0251	312
Buffer	0.03425	710
Inner PyC	0.03824	710
SiC	0.04177	1300
Outer PyC	0.04577	710

The packing fraction is used to homogenize the TRISO particle specific heat with the calculated specific heat of the matrix as

$var element = document.getElementById("moose-equation-50bef1ba-7994-4a1b-9688-3ca184de2145");katex.render(" c_{P-homogeneous} = PF c_{P-particle} + (1 - PF) c_P.", element, {displayMode:true,throwOnError:false});$

If the packing fraction is zero, then an un-fueled graphite matrix specific heat is calculated.

Example Input Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [matrix_thermal_props]
    type = GraphiteMatrixThermal<<<{"description": "Computes thermal conductivity (W/m/K) and specific heat capacity (J/kg/K) for graphite matrix.", "href": "GraphiteMatrixThermal.html"}>>>
    flux_conversion_factor<<<{"description": "Convert fast neutron flux E>0.10 to E>0.18 MeV"}>>> = 0.8
    packing_fraction<<<{"description": "Volumetric packing fraction (dimensionless) of particles"}>>> = 0.35
    etc_model<<<{"description": "Option for the correlation used to compute effective thermal conductivity (ETC)"}>>> = DEMT
    temperature<<<{"description": "Coupled temperature"}>>> = temp
    initial_matrix_density<<<{"description": "Initial pure matrix density value for A3-3 or A3-27 graphite (kg/m^3)."}>>> = 1200.0
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = all
  []
[]

Input Parameters

flux_conversion_factorConvert fast neutron flux E>0.10 to E>0.18 MeV
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Convert fast neutron flux E>0.10 to E>0.18 MeV
graphite_gradeA3_3_1800The grade of graphite A3_3_1800 A3_3_1950 A3_27_1800 A3_27_1950 H_451 IG_110 NG2020 PCEA NuclearGrade
Default:A3_3_1800
C++ Type:MooseEnum
Options:A3_3_1800, A3_3_1950, A3_27_1800, A3_27_1950, H_451, IG_110, NG2020, PCEA, NuclearGrade
Controllable:No
Description:The grade of graphite A3_3_1800 A3_3_1950 A3_27_1800 A3_27_1950 H_451 IG_110 NG2020 PCEA NuclearGrade
packing_fractionVolumetric packing fraction (dimensionless) of particles
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Volumetric packing fraction (dimensionless) of particles
temperatureCoupled temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature

Required Parameters

TRISO_conductivity4.13Particle thermal conductivity (W/m-K)
Default:4.13
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Particle thermal conductivity (W/m-K)
TRISO_specific_heat748.72Particle specific heat (J/kg-K)
Default:748.72
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Particle specific heat (J/kg-K)
base_nameOptional parameter that allows the user to define multiple material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple material systems on the same block, i.e. for multiple phases
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
densitydensityName of density material property.
Default:density
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Name of density material property.
etc_modelDEMTOption for the correlation used to compute effective thermal conductivity (ETC)
Default:DEMT
C++ Type:MooseEnum
Options:BRUGGEMAN, CHIEWGLANDT, DEMT, MAXWELL, PARALLEL, REDUCEDMAXWELL, SERIES
Controllable:No
Description:Option for the correlation used to compute effective thermal conductivity (ETC)
fast_neutron_fluencefast_neutron_fluenceCoupled fast (E>0.10 MeV) neutron fluence (n/m^2)
Default:fast_neutron_fluence
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Coupled fast (E>0.10 MeV) neutron fluence (n/m^2)
initial_matrix_density1700Initial pure matrix density value for A3-3 or A3-27 graphite (kg/m^3).
Default:1700
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Initial pure matrix density value for A3-3 or A3-27 graphite (kg/m^3).

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

specific_heat_scale_factor1Optional scaling factor applied to the overall specific heat.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Optional scaling factor applied to the overall specific heat.
thermal_conductivity_scale_factor1Optional scaling factor applied to the overall thermal conductivity.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Optional scaling factor applied to the overall thermal conductivity.

Advanced: Scaling Factors Parameters

Input Files

(test/tests/triso_pebble/1D_pebble_from_csv.i)
(test/tests/solid_mechanics/graphite_grade_elasticity_tensor/test.i)
(test/tests/triso/graphite_matrix_thermal/bruggeman.i)
(test/tests/triso/graphite_matrix_thermal/chiew_glandt.i)
(test/tests/triso/graphite_matrix_thermal/demt.i)
(test/tests/triso/graphite_matrix_thermal/reduced_maxwell.i)
(test/tests/triso/graphite_matrix_thermal/parallel.i)
(test/tests/triso_pebble/3D_pebble.i)
(examples/TRISO/pebble/3D_pebble_with_failed_particles.i)
(test/tests/triso/graphite_matrix_thermal/grade_pcea.i)
(test/tests/triso_pebble/1D_pebble_picard.i)
(test/tests/triso_pebble/3D_pebble_location_from_file.i)
(test/tests/triso/graphite_matrix_thermal/maxwell.i)
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_base.i)
(test/tests/triso/graphite_matrix_thermal/grade_h_451.i)
(examples/TRISO/pebble/3D_pebble.i)
(test/tests/triso/graphite_matrix_thermal/series.i)
(test/tests/triso/graphite_matrix_thermal/exact.i)
(test/tests/triso_pebble/1D_pebble.i)
(test/tests/triso/graphite_matrix_thermal/grade_ig_110.i)
(test/tests/triso_pebble/3D_pebble_from_csv.i)

References

D.A.G. Bruggeman. Berechnung verschiedener physikalischer konstanten von heterogenen substanzen. i. dielektrizitätskonstanten und leitfähigkeiten der mischköorper aus isotropen substanzen. Annalen der physik, 416(7):636–664, 1935. doi:10.1002/andp.19354160705.

@article{Bruggeman1935,
    author = "Bruggeman, {D.A.G.}",
    title = {Berechnung verschiedener physikalischer konstanten von heterogenen substanzen. I. Dielektrizit\"{a}tskonstanten und leitf\"{a}higkeiten der mischk\"{o}orper aus isotropen substanzen},
    journal = "Annalen der physik",
    volume = "416",
    number = "7",
    pages = "636-664",
    doi = "10.1002/andp.19354160705",
    year = "1935"
}

A.T.D. Butland and R.J. Maddison. The specific heat of graphite: An evaluation of measurements. Journal of Nuclear Materials, 49(1):45–56, 1973. doi:10.1016/0022-3115(73)90060-3.

@article{Butland1973,
    author = "Butland, A.T.D. and Maddison, R.J.",
    title = "The specific heat of graphite: {A}n evaluation of measurements",
    journal = "Journal of Nuclear Materials",
    volume = "49",
    number = "1",
    pages = "45-56",
    doi = "10.1016/0022-3115(73)90060-3",
    year = "1973"
}

Y. C. Chiew and E. D. Glandt. The effect of structure on the conductivity of a dispersion. Journal of Colloid and Interface Science, 91(1):90–104, 1983. doi:10.1016/0021-9797(83)90238-2.

@article{chiew1983,
    author = "Chiew, Y. C. and Glandt, E. D.",
    title = "The effect of structure on the conductivity of a dispersion",
    journal = "Journal of Colloid and Interface Science",
    volume = "91",
    number = "1",
    year = "1983",
    doi = "10.1016/0021-9797(83)90238-2",
    pages = "90-104"
}

Nam Zin Cho, Hui Yu, and Jong Woon Kim. Two-temperature homogenized model for steady-state and transient thermal analyses of a pebble with distributed fuel particles. Annals of Nuclear Energy, 36:448–457, May 2009. URL: https://www.sciencedirect.com/science/article/pii/S030645490800323X, doi:10.1016/j.anucene.2008.11.031.

@article{cho2009,
    author = "Cho, Nam Zin and Yu, Hui and Kim, Jong Woon",
    title = "Two-temperature homogenized model for steady-state and transient thermal analyses of a pebble with distributed fuel particles",
    journal = "Annals of Nuclear Energy",
    volume = "36",
    issue = "4",
    pages = "448-457",
    year = "2009",
    month = "May",
    doi = "10.1016/j.anucene.2008.11.031",
    url = "https://www.sciencedirect.com/science/article/pii/S030645490800323X"
}

Charles Folsom, Changhu Xing, Colby Jensen, Heng Ban, and Douglas W. Marshall. Experimental measurement and numerical modeling of the effective thermal conductivity of TRISO fuel compacts. Journal of Nuclear Materials, 458:198–205, 2015. doi:10.1016/j.jnucmat.2014.12.042.

@article{FOLSOM2015198,
    author = "Folsom, Charles and Xing, Changhu and Jensen, Colby and Ban, Heng and Marshall, Douglas W.",
    title = "Experimental measurement and numerical modeling of the effective thermal conductivity of {TRISO} fuel compacts",
    journal = "Journal of Nuclear Materials",
    volume = "458",
    pages = "198-205",
    year = "2015",
    issn = "0022-3115",
    doi = "10.1016/j.jnucmat.2014.12.042"
}

R. Gontard and H. Nabielek. Performance Evaluation of Modern HTR TRISO Fuels. Report HTA-1B-05/90, Forschungszentrums Juelich, July 1990.

@techreport{Gontard1990,
    author = "Gontard, R. and Nabielek, H.",
    title = "{Performance Evaluation of Modern HTR TRISO Fuels}",
    institution = "Forschungszentrums Juelich",
    year = "1990",
    month = "July",
    number = "HTA-1B-05/90",
    type = "Report"
}

G. L. Hawkes, J. W. Sterbentz, and B. T. Pham. Sensitivity evaluation of the AGR 3-4 experiment thermal model irradiated in the advanced test reactor. In ASME International Mechanical Engineering Congress and Exposition, volume 6B: Energy. 11 2015. doi:10.1115/IMECE2015-53544.

@inproceedings{HAWKES20156B,
    author = "Hawkes, G. L. and Sterbentz, J. W. and Pham, B. T.",
    title = "Sensitivity Evaluation of the {AGR 3-4} Experiment Thermal Model Irradiated in the Advanced Test Reactor",
    volume = "6B: Energy",
    booktitle = "ASME International Mechanical Engineering Congress and Exposition",
    doi = "10.1115/IMECE2015-53544",
    month = "11",
    year = "2015"
}

T. Maruyama, T. Kaito, S. Onose, and I. Shibahara. Change in physical properties of high density isotropic graphites irradiated in the ‘JOYO’ fast reactor. Journal of Nuclear Materials, 225:267–272, 1995. doi:10.1016/0022-3115(95)00051-8.

@article{MARUYAMA1995267,
    author = "Maruyama, T. and Kaito, T. and Onose, S. and Shibahara, I.",
    title = "Change in Physical Properties of High Density Isotropic Graphites Irradiated in the {‘JOYO’} Fast Reactor",
    journal = "Journal of Nuclear Materials",
    volume = "225",
    pages = "267--272",
    doi = "10.1016/0022-3115(95)00051-8",
    year = "1995"
}

J. C. Maxwell. A Treatise on Electricity and Magnetism. Oxford University Press, Dover, New York, 3rd edition, 1954.

@book{Maxwell1954,
    author = "Maxwell, J. C.",
    title = "A Treatise on Electricity and Magnetism",
    chapter = "Ch. IX. Conduction Through Heterogeneous Media",
    edition = "3rd",
    publisher = "Oxford University Press",
    address = "Dover, New York",
    year = "1954"
}

R. McLaughlin. A study of the differential scheme for composite materials. Int. J. Eng. Sci., 15(4):237–244, 1977. doi:10.1016/ 0020- 7225(77)90058- 1.

@article{McLaughlin1977,
    author = "McLaughlin, R.",
    title = "A study of the differential scheme for composite materials",
    journal = "Int. J. Eng. Sci.",
    volume = "15",
    number = "4",
    pages = "237-244",
    doi = "10.1016/ 0020- 7225(77)90058- 1",
    year = "1977"
}

G.K. Miller, D.A. Petti, J.T. Maki, D.L. Knudson, and W.F. Skerjanc. PARFUME Theory and Model Basis Report. Report INL/EXT-08-14497 (Rev.1), Idaho National Laboratory, September 2018.

@techreport{Miller2018,
    author = "Miller, G.K. and Petti, D.A. and Maki, J.T. and Knudson, D.L. and Skerjanc, W.F.",
    title = "{PARFUME Theory and Model Basis Report}",
    institution = "Idaho National Laboratory",
    year = "2018",
    month = "September",
    number = "INL/EXT-08-14497 (Rev.1)",
    type = "Report"
}

T. Shibata, M. Eto, E. Kunimoto, S. Shiozawa, K. Sawa, Oku T., and T. Maruyama. Draft of standard for graphite core components in high temperature gas-cooled reactor. Technical Report JAEAResearch2009-042, JAEA, 2009. URL: https://jopss.jaea.go.jp/pdfdata/JAEA-Research-2009-042.pdf.

@techreport{JAEAResearch2009042,
    author = "Shibata, T. and Eto, M. and Kunimoto, E. and Shiozawa, S. and Sawa, K. and T., Oku and Maruyama, T.",
    title = "Draft Of Standard For Graphite Core Components In High Temperature Gas-Cooled Reactor",
    institution = "JAEA",
    number = "JAEAResearch2009-042",
    url = "https://jopss.jaea.go.jp/pdfdata/JAEA-Research-2009-042.pdf",
    year = "2009"
}

L. L. Snead and T. D. Burchell. Reduction in thermal conductivity due to neutron irradiation. In 22nd Biennial Conference on Carbon, 774–775. 1995.

@inproceedings{SNEAD1995774,
    author = "Snead, L. L. and Burchell, T. D.",
    title = "Reduction in Thermal Conductivity Due to Neutron Irradiation",
    booktitle = "22nd Biennial Conference on Carbon",
    pages = "774--775",
    year = "1995"
}

M. Srinivasan and others. Appendices to the assessment of graphite properties and degradation, including source dependence. Technical Report TLR/RES/DE/REB-2021-08, United States Nuclear Regulatory Commision, 8 2021. URL: https://www.nrc.gov/docs/ML2121/ML21215A346.pdf.

@techreport{ML21215A346,
    author = "Srinivasan, M. and others",
    title = "Appendices to the Assessment of Graphite Properties and Degradation, Including Source Dependence",
    institution = "United States Nuclear Regulatory Commision",
    number = "TLR/RES/DE/REB-2021-08",
    month = "8",
    url = "https://www.nrc.gov/docs/ML2121/ML21215A346.pdf",
    year = "2021"
}

R. Stainsby and others. Investigation of local heat transfer phenomena in a pebble bed HTGR core. Technical Report NR001/RP/002 R01, AMEC, May 2009.

@techreport{Stainsby2009,
    author = "Stainsby, R. and others",
    title = "Investigation of local heat transfer phenomena in a pebble bed {HTGR} core",
    institution = "AMEC",
    number = "NR001/RP/002 R01",
    month = "May",
    year = "2009"
}

W. D. Swank, D. T. Rohrbaugh, D. L. Cottle, and W. E. Windes. High temperature annealing of irradiated graphite. Technical Report INL/CON-18-45742 (Rev. 0), Idaho National Laboratory, 2018. URL: https://www.osti.gov/biblio/1484477.

@techreport{INL_CON-18-45742,
    author = "Swank, W. D. and Rohrbaugh, D. T. and Cottle, D. L. and Windes, W. E.",
    title = "High Temperature Annealing of Irradiated Graphite",
    number = "INL/CON-18-45742 (Rev. 0)",
    institution = "Idaho National Laboratory",
    url = "https://www.osti.gov/biblio/1484477",
    year = "2018"
}

A. Toptan, W. Jiang, S. R. Novascone, and J. D. Hales. Matrix graphite material models in pebbles and compacts for Bison. Tech. Rep. INL/EXT-21-64643, Idaho National Laboratory, Idaho Falls, ID United States, 9 2021. URL: https://www.osti.gov/biblio/1825268, doi:10.2172/1825268.

@techreport{INL/EXT-21-64643,
    author = "Toptan, A. and Jiang, W. and Novascone, S. R. and Hales, J. D.",
    title = "Matrix graphite material models in pebbles and compacts for {Bison}",
    month = "9",
    type = "Tech. Rep.",
    number = "INL/EXT-21-64643",
    institution = "Idaho National Laboratory",
    url = "https://www.osti.gov/biblio/1825268",
    address = "Idaho Falls, ID United States",
    doi = "10.2172/1825268",
    year = "2021"
}

Aysenur Toptan, Jason D. Hales, Richard L. Williamson, Stephen R. Novascone, Giovanni Pastore, and David J. Kropaczek. Modeling of gap conductance for LWR fuel rods applied in the BISON code. Journal of Nuclear Science and Technology, 57(8):963–974, 2020. doi:10.1080/00223131.2020.1740808.

@article{TOPTAN2020963,
    author = "Toptan, Aysenur and Hales, Jason D. and Williamson, Richard L. and Novascone, Stephen R. and Pastore, Giovanni and Kropaczek, David J.",
    title = "Modeling of gap conductance for {LWR} fuel rods applied in the {BISON} code",
    journal = "Journal of Nuclear Science and Technology",
    year = "2020",
    volume = "57",
    number = "8",
    pages = "963--974",
    issn = "0022-3131",
    doi = "10.1080/00223131.2020.1740808"
}

Aysenur Toptan, Wen Jiang, Jason D. Hales, Benjamin W. Spencer, Albert Casagranda, and Stephen R. Novascone. FEA-aided investigation of the effective thermal conductivity in a medium with embedded spheres. Nuclear Engineering and Design, 381:111355, 2021. doi:10.1016/j.nucengdes.2021.111355.

@article{TOPTAN2021111355,
    author = "Toptan, Aysenur and Jiang, Wen and Hales, Jason D. and Spencer, Benjamin W. and Casagranda, Albert and Novascone, Stephen R.",
    title = "{FEA}-aided investigation of the effective thermal conductivity in a medium with embedded spheres",
    journal = "Nuclear Engineering and Design",
    pages = "111355",
    year = "2021",
    doi = "10.1016/j.nucengdes.2021.111355",
    volume = "381"
}

IAEA. Heat Transport and Afterheat Removal for Gas Cooled Reactors Under Accident Conditions. Technical Report TECDOC-1163, International Atomic Energy Agency, 2000.

@techreport{IAEA2000,
    author = "{IAEA}",
    title = "{Heat Transport and Afterheat Removal for Gas Cooled Reactors Under Accident Conditions}",
    institution = "International Atomic Energy Agency",
    year = "2000",
    number = "TECDOC-1163"
}

NEA Nuclear Science Committee. Nea benchmark of the modular high-temperature gas-cooled reactor-350 mw core design: volumes i and ii. Technical Report NEA/NSC/R(2017)4, Organisation for Economic Co-operation and Development, 2018. URL: https://www.oecd.org/officialdocuments/publicdisplaydocumentpdf/?cote=NEA/NSC/R(2017)4&docLanguage=En.

@techreport{OECDNEA2018,
    author = "{NEA Nuclear Science Committee}",
    title = "NEA Benchmark of the Modular High-Temperature Gas-Cooled Reactor-350 MW Core Design: Volumes I and II",
    institution = "Organisation for Economic Co-operation and Development",
    number = "NEA/NSC/R(2017)4",
    url = "https://www.oecd.org/officialdocuments/publicdisplaydocumentpdf/?cote=NEA/NSC/R(2017)4\&docLanguage=En",
    year = "2018"
}

(test/tests/triso/graphite_matrix_thermal/demt.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The DEMT model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = DEMT
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
[]

[Outputs]
  csv = true
  file_base = 'demt_out'
[]

(test/tests/triso_pebble/1D_pebble_from_csv.i)

[GlobalParams]
  order = SECOND
  family = LAGRANGE
  flux_conversion_factor = 0.85
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen]
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 0.0138 0.018 0.02'
    mesh_density = '5    10   10'
    block_names = 'core fuel shell'
  []
[]

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[DiracKernels]
  [vpp_point_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = heat_source
    value_name = heat_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    fuel_element_blocks = fuel
    execute_on = 'INITIAL'
    num_particles = 100
  []
[]

[VectorPostprocessors]
  [heat_source]
    type = TRISOMonteCarloPointSourceFromFile
    file_base_time = monte_carlo_out_time
    file_base_source = monte_carlo_out_release_heat_inc
    x_header_in_file = 'time:time'
    y_header_in_file = 'release_heat_inc:release_heat_inc'
    source_name = heat_inc
    num_files = 2
    num_particles = 100
    format = columns
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
  [point_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = heat_source
    value_name = heat_inc
    point_source_location = point_source_location
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 773.15
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = 'core fuel shell'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'core fuel shell'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 'core fuel shell'
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6'
    y = '1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 5.75e19
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 6.28e17
  []
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = fuel
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fuel_region_density]
    type = ParsedMaterial
    block = fuel
    property_name = density
    expression = 1750.0
  []
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = shell
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = shell
    property_name = density
    expression = 1750.0
  []
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = core
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = core
    property_name = density
    expression = 1400.0
  []
[]

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

[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-9
  nl_abs_tol = 1e-9
  nl_max_its = 15

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 10
  num_steps = 2
[]

[Outputs]
  print_linear_residuals = true
  time_step_interval =  1
  exodus = true
  perf_graph = true
[]

(test/tests/solid_mechanics/graphite_grade_elasticity_tensor/test.i)

# This test case is prepared to test the elastic properties of
# grade H-451 grade material model.

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [mesh]
    type = CircularCrossSectionMeshGenerator
    num_sectors = 30
    offset = '0.0 0.0'
    elements_per_ring = '0 4 0'
    block_names = 'null shell null2'
    coordinates = '0.039 0.064 0.065'
  []
[]

[AuxVariables]
  [temperature]
    initial_condition = 673.15
  []
  [fast_neutron_fluence]
    order = CONSTANT
    family = MONOMIAL
  []
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 0.075 * t'
  []
  [flux_history]
    type = PiecewiseLinear
    x = '0   1e4'
    y = '0 16e21'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [perm]
    add_variables = true
    strain = SMALL
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx strain_yy strain_zz'
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temperature
    function = temp_fcn
  []
  [fast_neutron_fluence]
    type = MaterialRealAux
    variable = fast_neutron_fluence
    property = fast_neutron_fluence
    execute_on = timestep_end
  []
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    execute_on = 'initial linear'
  []
[]

[BCs]
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0
  []
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [Pressure]
    [inside_pressure]
      boundary = 1001
      factor = 1e4
    []
  []
[]

[Materials]
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = flux_history
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [normal_vectors_triso]
    type = NormalVectorsTRISO
  []
  [elasticity_tensor]
    type = GraphiteGradeElasticityTensor
    temperature = temperature
    graphite_grade = H_451
    fast_neutron_fluence = fast_neutron_fluence
    flux_conversion_factor = 0.85
  []
  [thermal_props]
    type = GraphiteMatrixThermal
    graphite_grade = H_451
    temperature = temperature
    packing_fraction = 0.0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    flux_conversion_factor = 0.85
  []
  [density]
    type = StrainAdjustedDensity
    strain_free_density = 1900.0
  []
[]

[Postprocessors]
  [temperature]
    type = ElementExtremeValue
    value_type = 'max'
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [fluence]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fast_neutron_fluence
    execute_on = 'initial timestep_end'
  []
  [sigma_x]
    type = ElementExtremeValue
    value_type = 'max'
    variable = stress_xx
    execute_on = 'initial timestep_end'
  []
  [sigma_y]
    type = ElementExtremeValue
    value_type = 'max'
    variable = stress_yy
    execute_on = 'initial timestep_end'
  []
  [sigma_z]
    type = ElementExtremeValue
    value_type = 'max'
    variable = stress_zz
    execute_on = 'initial timestep_end'
  []
  [strain_x]
    type = ElementExtremeValue
    value_type = 'max'
    variable = strain_xx
    execute_on = 'initial timestep_end'
  []
  [strain_y]
    type = ElementExtremeValue
    value_type = 'max'
    variable = strain_yy
    execute_on = 'initial timestep_end'
  []
  [strain_z]
    type = ElementExtremeValue
    value_type = 'max'
    variable = strain_zz
    execute_on = 'initial timestep_end'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'
  line_search = 'none'

  l_tol = 1e-5
  nl_max_its = 10
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-10
  start_time = 0.0
  num_steps = 10
  dt = 1e3
[]

[Outputs]
  csv = true
[]

(test/tests/triso/graphite_matrix_thermal/bruggeman.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The BRUGGEMAN (or EMT) model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf    thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.35  1.0'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             alpha := thcond_particle/k_graphite;
             coeff_A := 0.25 * (3.0 * pf - 1.0 + (2.0 - 3.0 * pf) / alpha);
             f_packing_rhs := alpha * coeff_A + sqrt((alpha * coeff_A)^2 + 0.5 * alpha);
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = BRUGGEMAN
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = 'bruggeman_out'
[]

(test/tests/triso/graphite_matrix_thermal/chiew_glandt.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix
# for each graphite type: A3-3-1800, A3-3-1950, A3-27-1800, and A3-27-1950.
#
# The CHIEWGLANDT model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp spheat_scale_factor'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72 1.1'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             spheat_scale_factor * (cp * (1 - pf) + part_cp * pf)'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf    thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.35  0.9'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             p_beta := (thcond_particle - k_graphite) / (thcond_particle + 2 * k_graphite);
             f_packing_rhs := (1 + 2 * p_beta * pf +
                           (2 * p_beta^3 - 0.1 * p_beta) * pf^2 +
                           0.05 * exp(4.5 * p_beta) * pf^3) / (1.0 - p_beta * pf);
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = CHIEWGLANDT
    specific_heat_scale_factor = 1.1
    thermal_conductivity_scale_factor = 0.9
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = exact_A3_3_1800_out
[]

(test/tests/triso/graphite_matrix_thermal/demt.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The DEMT model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = DEMT
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
[]

[Outputs]
  csv = true
  file_base = 'demt_out'
[]

(test/tests/triso/graphite_matrix_thermal/reduced_maxwell.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The reduced form of MAXWELL (REDUCEDMAXWELL) model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.35  1.0'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             f_packing_rhs := (1.0 - pf) / (1.0 + 0.5 * pf);
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = REDUCEDMAXWELL
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = 'reduced_maxwell_out'
[]

(test/tests/triso/graphite_matrix_thermal/parallel.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The PARALLEL model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf    thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.35  1.0'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             alpha := thcond_particle/k_graphite;
             f_packing_rhs := 1.0 + pf * (alpha -1.0);
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = PARALLEL
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = 'parallel_out'
[]

(test/tests/triso_pebble/3D_pebble.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  flux_conversion_factor = 0.85
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = pebble.e
  []
[]

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 10
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  [sub]
    type = SamplerTransientMultiApp
    input_files = triso_1d.i
    sampler = sample
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
[]

[DiracKernels]
  [vpp_point_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    sampler = sample
    fuel_element_blocks = fuel
    execute_on = 'INITIAL'
  []
[]

[Transfers]
  [released_heat_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = released_heat_inc
    from_postprocessor = released_heat_inc
  []
  [temp_bc]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub
    vector_postprocessor = sample_points
    postprocessor = temp_bc
    vector_name = temperature
  []
[]

[Controls]
  [cmdline]
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
[]

[VectorPostprocessors]
  [sample_points]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location
  []
  [sampler_data]
    type = SamplerData
    execute_on = 'TIMESTEP_BEGIN'
    sampler = sample
  []
  [released_heat_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [point_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 773.15
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = 'core fuel shell'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'core fuel shell'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 'core fuel shell'
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6'
    y = '1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 5.75e19
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 6.28e17
  []
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = fuel
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fuel_region_density]
    type = ParsedMaterial
    block = fuel
    property_name = density
    expression = 1750.0
  []
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = shell
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = shell
    property_name = density
    expression = 1750.0
  []
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = core
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = core
    property_name = density
    expression = 1400.0
  []
[]

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

[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-9
  nl_abs_tol = 1e-9
  nl_max_its = 20

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 20000
  num_steps = 2
[]

[Outputs]
  print_linear_residuals = true
  time_step_interval =  1
  exodus = true
  csv = true
  perf_graph = true
[]

(examples/TRISO/pebble/3D_pebble_with_failed_particles.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  flux_conversion_factor = 0.85
  energy_per_fission = 3.204e-11 #[J/fission]
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = pebble.e
  []
[]

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

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 100
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
  [sample_failed]
    type = MonteCarlo
    num_rows = 10
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  [sub]
    type = SamplerTransientMultiApp
    input_files = triso_1d.i
    sampler = sample
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
  [sub_failed]
    type = SamplerTransientMultiApp
    input_files = triso_1d_failed.i
    sampler = sample_failed
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
[]

[DiracKernels]
  [heat_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = release_heat_inc
    value_name = release_heat_inc:release_heat_inc
    point_source_location = point_source_location
  []
  [Cs_source]
    type = TRISOMonteCarloPointSource
    variable = conc_Cs
    point_source_values = release_Cs_inc
    value_name = release_Cs_inc:release_Cs_inc
    point_source_location = point_source_location
  []

  [heat_source_failed]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = release_heat_inc_failed
    value_name = release_heat_inc_failed:release_heat_inc
    point_source_location = point_source_location_failed
  []
  [Cs_source_failed]
    type = TRISOMonteCarloPointSource
    variable = conc_Cs
    point_source_values = release_Cs_inc_failed
    value_name = release_Cs_inc_failed:release_Cs_inc
    point_source_location = point_source_location_failed
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    sampler = sample
    fuel_element_blocks = 2
    particle_space = 0.0002
    execute_on = 'INITIAL'
  []
  [point_source_location_failed]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    sampler = sample_failed
    fuel_element_blocks = 2
    particle_space = 0.0002
    execute_on = 'INITIAL'
  []
[]

[Transfers]
  [release_heat_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = release_heat_inc
    from_postprocessor = release_heat_inc
  []
  [release_Cs_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = release_Cs_inc
    from_postprocessor = release_Cs_inc
  []
  [release_heat_inc_failed]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub_failed
    sampler = sample_failed
    to_vector_postprocessor = release_heat_inc_failed
    from_postprocessor = release_heat_inc
  []
  [release_Cs_inc_failed]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub_failed
    sampler = sample_failed
    to_vector_postprocessor = release_Cs_inc_failed
    from_postprocessor = release_Cs_inc
  []
  [temp_bc]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub
    vector_postprocessor = sample_points
    postprocessor = temp_bc
    vector_name = temperature
  []
  [temp_bc_failed]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub_failed
    vector_postprocessor = sample_points_failed
    postprocessor = temp_bc
    vector_name = temperature
  []
[]

[Controls]
  [cmdline]
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
  [cmdline_failed]
    type = MultiAppSamplerControl
    multi_app = sub_failed
    sampler = sample_failed
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
[]

[VectorPostprocessors]
  [sample_points]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location
  []
  [sample_points_failed]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location_failed
  []
  [release_heat_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [release_Cs_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [heat_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = release_heat_inc
    value_name = release_heat_inc:release_heat_inc
    point_source_location = point_source_location
    execute_on = 'TIMESTEP_END'
  []
  [Cs_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = release_Cs_inc
    value_name = release_Cs_inc:release_Cs_inc
    point_source_location = point_source_location
    execute_on = 'TIMESTEP_END'
  []

  [release_heat_inc_failed]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [release_Cs_inc_failed]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [heat_source_output_failed]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = release_heat_inc_failed
    value_name = release_heat_inc_failed:release_heat_inc
    point_source_location = point_source_location_failed
    execute_on = 'TIMESTEP_END'
  []
  [Cs_source_output_failed]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = release_Cs_inc_failed
    value_name = release_Cs_inc_failed:release_Cs_inc
    point_source_location = point_source_location_failed
    execute_on = 'TIMESTEP_END'
  []
[]

[Variables]
  [temperature]
    initial_condition = 773.15
  []
  [conc_Cs]
    scaling = 1e14
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
  [Cs_diff_coef]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e3 4.51008e7' #change time (s) for desired EFPD
    y = '0 1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 2.927e18
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [mass_Cs_dt]
    type = TimeDerivative
    variable = conc_Cs
  []
  [mass_Cs]
    type = ArrheniusDiffusion
    variable = conc_Cs
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = '1 2 3'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = '1 2 3'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = '1 2 3'
    execute_on = timestep_end
  []
  [Cs_diff_coef]
    type = MaterialRealAux
    variable = Cs_diff_coef
    property = arrhenius_diffusion_coef_Cs
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
  [freesurf_conc_Cs]
    type = DirichletBC
    variable = conc_Cs
    boundary = exterior
    value = 0.0
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 1.109e18
  []

  ##### fuel region propereties
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = 2
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fueled_region_density]
    type = ParsedMaterial
    block = 2
    property_name = density
    expression = 1750.0
  []
  [fuel_conc_Cs]
    type = ArrheniusDiffusionCoef
    block = 2
    d1 = 1e-8 # m^2/s
    q1 = 0.0 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []

  ##### Shell properties
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = 3
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = 3
    property_name = density
    expression = 1750.0
  []

  [shell_conc_Cs]
    type = ArrheniusDiffusionCoef
    block = 3
    d1 = 1e-8 # m^2/s
    q1 = 0.0 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []

  ##### pebble core
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = 1
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = 1
    property_name = density
    expression = 1400.0
  []
  [core_conc_Cs]
    type = ArrheniusDiffusionCoef
    block = 1
    d1 = 1e-8 # m^2/s
    q1 = 0.0 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []
[]

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

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  line_search = 'none'

  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-9
  nl_max_its = 15

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 50000
  num_steps = 10
[]

[Postprocessors]

  ### irradiation conditions
  [pebble_power]
    type = ElementIntegralPower
    variable = temperature
    use_material_fission_rate = true
    fission_rate_material = fission_rate
    block = 2
    execute_on = 'initial timestep_end'
  []
  [heat_from_pebble]
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = exterior
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [fast_neutron_maximum_fluence_MP]
    type = ElementExtremeMaterialProperty
    mat_prop = fast_neutron_fluence
    block = '1 2 3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []
  [fast_neutron_maximum_flux_MP]
    type = ElementExtremeMaterialProperty
    mat_prop = fast_neutron_flux
    block = '1 2 3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []
  [fission_rate_maximum]
    type = ElementExtremeMaterialProperty
    mat_prop = fission_rate
    block = '1 2 3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  ### temperatures
  [fuel_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '2'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [fuel_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = '2'
    execute_on = 'initial timestep_end'
  []
  [fuel_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '2'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [shell_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '3'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [shell_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = '3'
    execute_on = 'initial timestep_end'
  []
  [shell_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [core_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '1'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [core_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = '1'
    execute_on = 'initial timestep_end'
  []
  [core_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '1'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  print_linear_residuals = true
  time_step_interval =  1
  exodus = true
  csv = true
  perf_graph = true
[]

(test/tests/triso/graphite_matrix_thermal/grade_pcea.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix
# for grade PCEA graphite.
# Assuming zero packing of TRISO particles in the graphite matrix.
#
# The exact solution is performed through Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 100 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'tk              pf  part_cp rho        '
    symbol_values = 'temperature_avg 0.0 748.72 density_avg'
    expression = 'cp := 1.75e6 * (0.645 + 3.14 * ((tk-273.15)/1000) - 2.809 * ((tk-273.15)/1000)^2 + 0.959 * ((tk-273.15)/1000)^3) / rho;
             cp * (1 - pf) + part_cp * pf'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    graphite_grade = PCEA
    packing_fraction = 0.0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    temperature = temp
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [avg_cP]
    type = ElementAverageValue
    variable = specific_heat
  []
  [exact_cP]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [diff_cP]
    type = DifferencePostprocessor
    value1 = avg_cP
    value2 = exact_cP
    outputs = none
  []
  [max_diff_cP]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = diff_cP
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
[]

[Outputs]
  csv = true
  file_base = 'exact_PCEA_out'
[]

(test/tests/triso_pebble/1D_pebble_picard.i)

[GlobalParams]
  order = SECOND
  family = LAGRANGE
  flux_conversion_factor = 0.85
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen]
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 0.0138 0.018 0.02'
    mesh_density = '5    10   10'
    block_names = 'core fuel shell'
  []
[]

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 10
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  [sub]
    type = SamplerTransientMultiApp
    input_files = triso_1d.i
    sampler = sample
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
[]

[DiracKernels]
  [vpp_point_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    sampler = sample
    fuel_element_blocks = fuel
    execute_on = 'INITIAL'
  []
[]

[Transfers]
  [released_heat_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = released_heat_inc
    from_postprocessor = released_heat_inc
  []
  [sub_temp_bc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = temp_bc
    from_postprocessor = aver_temp_exterior
  []
  [temp_bc]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub
    vector_postprocessor = sample_points
    postprocessor = temp_bc
    vector_name = temperature
  []
[]

[Controls]
  [cmdline]
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
[]

[VectorPostprocessors]
  [sample_points]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location
  []
  [sampler_data]
    type = SamplerData
    execute_on = 'TIMESTEP_BEGIN'
    sampler = sample
  []
  [released_heat_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [temp_bc]
    type = StochasticResults
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
  [point_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
  [temp_bc_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = temp_bc
    value_name = sub_temp_bc:aver_temp_exterior
    point_source_location = point_source_location
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 773.15
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = FissionRateHeatSource
    variable = temperature
    block = fuel
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = 'core fuel shell'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'core fuel shell'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 'core fuel shell'
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6'
    y = '1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 5.75e19
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 6.28e17
  []
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = fuel
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fuel_region_density]
    type = ParsedMaterial
    block = fuel
    property_name = density
    expression = 1750.0
  []
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = shell
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = shell
    property_name = density
    expression = 1750.0
  []
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = core
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = core
    property_name = density
    expression = 1400.0
  []
[]

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

[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-12
  nl_abs_tol = 1e-12
  nl_max_its = 15

  nl_forced_its = 2

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 20000
  num_steps = 2

  fixed_point_max_its = 20
  accept_on_max_fixed_point_iteration = true
  fixed_point_rel_tol = 1e-8
  fixed_point_abs_tol = 1e-10
[]

[Outputs]
  exodus = true
[]

(test/tests/triso_pebble/3D_pebble_location_from_file.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  flux_conversion_factor = 0.85
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = pebble.e
  []
[]

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 10
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  [sub]
    type = SamplerTransientMultiApp
    input_files = triso_1d.i
    sampler = sample
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
[]

[DiracKernels]
  [vpp_point_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocationFromFile
    sampler = sample
    particle_location_file_name = location.txt
    execute_on = 'INITIAL'
  []
[]

[Transfers]
  [released_heat_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = released_heat_inc
    from_postprocessor = released_heat_inc
  []
  [temp_bc]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub
    vector_postprocessor = sample_points
    postprocessor = temp_bc
    vector_name = temperature
  []
[]

[Controls]
  [cmdline]
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
[]

[VectorPostprocessors]
  [sample_points]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location
  []
  [sampler_data]
    type = SamplerData
    execute_on = 'TIMESTEP_BEGIN'
    sampler = sample
  []
  [released_heat_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [point_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 773.15
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = 'core fuel shell'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'core fuel shell'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 'core fuel shell'
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6'
    y = '1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 5.75e19
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 6.28e17
  []
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = fuel
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fuel_region_density]
    type = ParsedMaterial
    block = fuel
    property_name = density
    expression = 1750.0
  []
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = shell
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = shell
    property_name = density
    expression = 1750.0
  []
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = core
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = core
    property_name = density
    expression = 1400.0
  []
[]

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

[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-9
  nl_abs_tol = 1e-9
  nl_max_its = 20

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 20000
  num_steps = 2
[]

[Outputs]
  print_linear_residuals = true
  time_step_interval =  1
  exodus = true
  csv = true
  perf_graph = true
[]

(test/tests/triso/graphite_matrix_thermal/maxwell.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The MAXWELL model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf    thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.35  1.0'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             alpha := thcond_particle/k_graphite;
             f_packing_rhs := ( 1.0 + 3.0 * pf / ((alpha+2.0)/(alpha-1.0) - pf) );
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = MAXWELL
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = 'maxwell_out'
[]

(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_base.i)

# rings:                       DTF  |  DRV  |  IR  |  OR  |  SR
# sorption interfaces:           DRV_IR   IR_OR  OR_SR

center= 0.0
dtf_outer= 0.000175
driver_outer= 0.006275
ir_outer= 0.011805
or_outer= 0.016290
sr_outer= 0.026480


initial_fuel_density = 1700

[GlobalParams]
  family = LAGRANGE
  order = SECOND
  elem_type = QUAD8
  flux_conversion_factor = 0.85
  end_point = '${sr_outer} 0 0.0'
[]

[Mesh]
  coord_type = RZ
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 726
    ny = 1
    xmax = ${sr_outer}
    # This is quasi-2D case, so dimension in y direction does not matter.
    ymax = 1e-4
  []

  # If you have different radia for rings, you need to modify the bounding box range below.
  [subdomain0]
    input = gen
    type = SubdomainBoundingBoxGenerator
    bottom_left = '${center} 0 0'
    block_id = 0
    block_name = 'DTF'
    top_right = '${dtf_outer} 12.5e-3 0'
  []
  [subdomain1]
    input = subdomain0
    type = SubdomainBoundingBoxGenerator
    bottom_left = '${dtf_outer} 0 0'
    block_id = 1
    block_name = 'DRV'
    top_right = '${driver_outer} 12.5e-3 0'
  []
  [subdomain2]
    input = subdomain1
    type = SubdomainBoundingBoxGenerator
    bottom_left = '${driver_outer} 0 0'
    block_id = 2
    block_name = 'IR'
    top_right = '${ir_outer} 12.5e-3 0'
  []
  [subdomain3]
    input = subdomain2
    type = SubdomainBoundingBoxGenerator
    bottom_left = '${ir_outer} 0 0'
    block_id = 3
    block_name = 'OR'
    top_right = '${or_outer} 12.5e-3 0'
  []
  [subdomain4]
    input = subdomain3
    type = SubdomainBoundingBoxGenerator
    bottom_left = '${or_outer} 0 0'
    block_id = 4
    block_name = 'SR'
    top_right = '${sr_outer} 12.5e-3 0'
  []
  [breakmesh]
    input = subdomain4
    type = BreakMeshByBlockGenerator
    block_pairs = '1 2; 2 3; 3 4'
    add_interface_on_two_sides = true
    split_interface = true
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 1200
  []
[]

[AuxVariables]
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative #HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source_dtf]
    type = HeatSource
    function = DTF_volumetric_heat #volumetric_heat_DTF
    block = 'DTF'
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source_drv]
    type = HeatSource
    function = DRV_volumetric_heat #volumetric_heat_DRV
    block = 'DRV'
    variable = temperature
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'DTF DRV IR OR SR'
    execute_on = timestep_end
  []
[]

[Functions]
  [DTF_volumetric_heat]
    type = PiecewiseLinear
    data_file = capsule1_dtf_out.csv
    x_index_in_file = 0
    y_index_in_file = 9
    xy_in_file_only = false
    format = columns
    scale_factor = 16630067523
    # we need to convert the point source to volumetric source
    # 20/(pi*${dtf_outer}^2*12.5e-3)
    # 20: number of DTF particles
    # ${dtf_outer} DTF fuel outer radius
    # 12.5e-3 Compact height
  []
  [DTF_MassIsotope]
    type = PiecewiseLinear
    data_file = capsule1_dtf_out.csv
    x_index_in_file = 0
    # y_index_in_file = 7 #This number should not be hardcoded
    xy_in_file_only = false
    format = columns
    scale_factor = 16630067523
  []
  [DRV_volumetric_heat]
    type = PiecewiseLinear
    data_file = capsule1_driver_out.csv
    x_index_in_file = 0
    y_index_in_file = 9
    xy_in_file_only = false
    format = columns
    scale_factor = 1198178193
    # 1872/(pi*((0.006155)^2-(${dtf_outer})^2)*12.5e-3)
    # 1872: number of driver particles
    # 0.006155 compact outer radius
    # ${dtf_outer} DTF fuel outer radius
    # 12.5e-3 Compact height
  []
  [DRV_MassIsotope]
    type = PiecewiseLinear
    data_file = capsule1_driver_out.csv
    x_index_in_file = 0
    # y_index_in_file = 7 #This number should not be hardcoded
    xy_in_file_only = false
    format = columns
    scale_factor = 1198178193
  []

  # DRV to IR sorption partial pressures
  [pressure_DRV_IR_error]
    type = ParsedFunction
    symbol_names = 'pressure_DRV_outer pressure_IR_inner'
    symbol_values = 'pressure_DRV_outer pressure_IR_inner'
    expression = 'pressure_DRV_outer - pressure_IR_inner'
  []

  # IR to OR sorption partial pressures
  [pressure_IR_OR_error]
    type = ParsedFunction
    symbol_names = 'pressure_IR_outer pressure_OR_inner'
    symbol_values = 'pressure_IR_outer pressure_OR_inner'
    expression = 'pressure_IR_outer - pressure_OR_inner'
  []

  # OR to SR sorption partial pressures
  [pressure_OR_SR_error]
    type = ParsedFunction
    symbol_names = 'pressure_OR_outer pressure_SR_inner'
    symbol_values = 'pressure_OR_outer pressure_SR_inner'
    expression = 'pressure_OR_outer - pressure_SR_inner'
  []

  # DRV to IR flux
  [flux_DRV_IR_error]
    type = ParsedFunction
    symbol_names = 'flux_DRV_outer flux_IR_inner'
    symbol_values = 'flux_DRV_outer flux_IR_inner'
    expression = 'flux_DRV_outer + flux_IR_inner'
  []

  # IR to OR flux
  [flux_IR_OR_error]
    type = ParsedFunction
    symbol_names = 'flux_IR_outer flux_OR_inner'
    symbol_values = 'flux_IR_outer flux_OR_inner'
    expression = 'flux_IR_outer + flux_OR_inner'
  []

  # OR to SR flux
  [flux_OR_SR_error]
    type = ParsedFunction
    symbol_names = 'flux_OR_outer flux_SR_inner'
    symbol_values = 'flux_OR_outer flux_SR_inner'
    expression = 'flux_OR_outer + flux_SR_inner'
  []
[]

[BCs]
  [fuel_compact_outer]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'DRV_IR'
    function = temp_compact_outer
    extra_vector_tags = 'ref'
  []
  [matrix_inner]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'IR_DRV'
    function = temp_matrix_inner
    extra_vector_tags = 'ref'
  []
  [matrix_outer]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'IR_OR'
    function = temp_matrix_outer
    extra_vector_tags = 'ref'
  []
  [graphite_inner]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'OR_IR'
    function = temp_graphite_inner
    extra_vector_tags = 'ref'
  []
  [graphite_outer]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'OR_SR'
    function = temp_graphite_outer
    extra_vector_tags = 'ref'
  []
  [sink_inner]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'SR_OR'
    function = temp_sink_inner
    extra_vector_tags = 'ref'
  []
  [sink_outer]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'right'
    function = temp_sink_outer
    extra_vector_tags = 'ref'
  []
[]

[Functions]
  [temp_compact_outer]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 3
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [temp_matrix_inner]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 4
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [temp_matrix_outer]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 5
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [temp_graphite_inner]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 6
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [temp_graphite_outer]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 7
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [temp_sink_inner]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 8
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [temp_sink_outer]
    type = PiecewiseLinear
    x_index_in_file = 0
    y_index_in_file = 9
    xy_in_file_only = false
    format = columns
    data_file= AGR-34_capsule_daily_data/Cap1Temps.csv
  []
  [power_history]
    type = PiecewiseLinear
    x = '0 1e3 4.51008e7' # those are used for some of the graphite thermal properties.
    y = '0 1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 4.8156e+19 #those are used for some of the graphite thermal properties.
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  # FastNeutronFlux - factor
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 0.5519e+18
  []
  [matrix_thermal]
    type = GraphiteMatrixThermal
    block = 'DTF DRV IR'
    graphite_grade = H_451  # subject to change
    packing_fraction = 0
    temperature = temperature
  []

  [matrix_density]
    type = StrainAdjustedDensity
    block = 'DTF DRV IR'
    density = ${initial_fuel_density}
  []
  [graphite_thermal]
    type = GraphiteMatrixThermal
    block = 'OR SR'
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1603.0
  []
  [graphite_density]
    type = StrainAdjustedDensity
    block = 'OR SR'
    density = 1603.0
  []
[]

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

[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 = 5e-4
  nl_abs_tol = 1e-7 # 1e-10
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 50

  start_time = 0.0
  end_time = 31890240
  dt = 100

  dtmax = 2e6
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
    optimal_iterations = 6
    growth_factor = 1.5
    linear_iteration_ratio = 100
  []

  [Predictor]
    type = SimplePredictor
    scale = 1
  []
[]

[Postprocessors]
  [compact_DTF_volume]
    type = ElementIntegralMaterialProperty
    mat_prop = 1
    block = DTF
  []

  ### temperatures
  [fuel_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'DRV'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [fuel_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'DRV'
    execute_on = 'initial timestep_end'
  []
  [fuel_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'DRV'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [matrix_ring_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'IR'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [matrix_ring_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'IR'
    execute_on = 'initial timestep_end'
  []
  [matrix_ring_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'IR'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [graphite_ring_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'OR'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [graphite_ring_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'OR'
    execute_on = 'initial timestep_end'
  []
  [graphite_ring_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'OR'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [sink_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'SR'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [sink_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = 'SR'
    execute_on = 'initial timestep_end'
  []
  [sink_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = 'SR'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  # check partial pressure equality and integral flux conservation from DRV to IR
  [pressure_DRV_outer]
    type = SideAverageMaterialProperty
    boundary = 'DRV_IR'
    property = partial_pressure
  []
  [pressure_IR_inner]
    type = SideAverageMaterialProperty
    boundary = 'IR_DRV'
    property = partial_pressure
  []
  [pressure_DRV_IR_error]
    type = FunctionValuePostprocessor
    function = pressure_DRV_IR_error
  []

  # check partial pressure equality and integral flux conservation from IR to OR
  [pressure_IR_outer]
    type = SideAverageMaterialProperty
    boundary = 'IR_OR'
    property = partial_pressure
  []
  [pressure_OR_inner]
    type = SideAverageMaterialProperty
    boundary = 'OR_IR'
    property = partial_pressure
  []
  [pressure_IR_OR_error]
    type = FunctionValuePostprocessor
    function = pressure_IR_OR_error
  []

  # check partial pressure equality and integral flux conservation from OR to SR
  [pressure_OR_outer]
    type = SideAverageMaterialProperty
    boundary = 'OR_SR'
    property = partial_pressure
  []
  [pressure_SR_inner]
    type = SideAverageMaterialProperty
    boundary = 'SR_OR'
    property = partial_pressure
  []
  [pressure_OR_SR_error]
    type = FunctionValuePostprocessor
    function = pressure_OR_SR_error
  []
[]

[VectorPostprocessors]
  [temp]
    type = LineValueSampler
    variable = temperature
    start_point = '0.0 0 0.0'
    num_points = 1000
    use_displaced_mesh = false
    execute_on = final
    sort_by = x
    outputs = line_plot_temp
  []
[]

[Outputs]
  csv= true
  [line_plot_temp]
    type = CSV
    execute_on = 'FINAL'
    time_step_interval =  1
    create_final_symlink = true
    file_base = capsule_base_radial
  []
  [chkfile]
    type = CSV
    execute_on = 'FINAL'
  []
[]

(test/tests/triso/graphite_matrix_thermal/grade_h_451.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix
# for grade H-451 graphite.
# Assuming zero packing of TRISO particles in the graphite matrix.
#
# The exact solution is performed through Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'tk              pf  part_cp'
    symbol_values = 'temperature_avg 0.0 748.72'
    expression = 'cp := (0.54212 - 2.42667e-6 * tk - 90.2725 / tk - 43449.3 / tk^2 + 1.59309e7 / tk^3 - 1.43688e9 / tk^4) * 4814;
             cp * (1 - pf) + part_cp * pf'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    graphite_grade = H_451
    packing_fraction = 0.0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    temperature = temp
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [avg_cP]
    type = ElementAverageValue
    variable = specific_heat
  []
  [exact_cP]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [diff_cP]
    type = DifferencePostprocessor
    value1 = avg_cP
    value2 = exact_cP
    outputs = none
  []
  [max_diff_cP]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = diff_cP
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
[]

[Outputs]
  csv = true
  file_base = 'exact_H_451_out'
[]

(examples/TRISO/pebble/3D_pebble.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  flux_conversion_factor = 0.85
  energy_per_fission = 3.204e-11 #[J/fission]
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = pebble.e
  []
[]

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

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 100
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  [sub]
    type = SamplerTransientMultiApp
    input_files = triso_1d.i
    sampler = sample
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
[]

[DiracKernels]
  [heat_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = release_heat_inc
    value_name = release_heat_inc:release_heat_inc
    point_source_location = point_source_location
  []
  [Cs_source]
    type = TRISOMonteCarloPointSource
    variable = conc_Cs
    point_source_values = release_Cs_inc
    value_name = release_Cs_inc:release_Cs_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    sampler = sample
    fuel_element_blocks = 2
    particle_space = 0.0002
    execute_on = 'INITIAL'
  []
[]

[Transfers]
  [release_heat_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = release_heat_inc
    from_postprocessor = release_heat_inc
  []
  [release_Cs_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = release_Cs_inc
    from_postprocessor = release_Cs_inc
  []
  [sub_temp_bc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = temp_bc
    from_postprocessor = aver_temp_exterior
  []
  [temp_bc]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub
    vector_postprocessor = sample_points
    postprocessor = temp_bc
    vector_name = temperature
  []
[]

[Controls]
  [cmdline]
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
[]

[VectorPostprocessors]
  [sample_points]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location
  []
  [release_heat_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [release_Cs_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [heat_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = release_heat_inc
    value_name = release_heat_inc:release_heat_inc
    point_source_location = point_source_location
  []
  [Cs_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = release_Cs_inc
    value_name = release_Cs_inc:release_Cs_inc
    point_source_location = point_source_location
  []
  [temp_bc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [temp_bc_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = temp_bc
    value_name = sub_temp_bc:aver_temp_exterior
    point_source_location = point_source_location
  []
[]

[Variables]
  [temperature]
    initial_condition = 773.15
  []
  [conc_Cs]
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
  [Cs_diff_coef]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e3 4.51008e7' #change time (s) for desired EFPD
    y = '0 1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 2.927e18
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [mass_Cs_dt]
    type = TimeDerivative
    variable = conc_Cs
  []
  [mass_Cs]
    type = ArrheniusDiffusion
    variable = conc_Cs
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = '1 2 3'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = '1 2 3'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = '1 2 3'
    execute_on = timestep_end
  []
  [Cs_diff_coef]
    type = MaterialRealAux
    variable = Cs_diff_coef
    property = arrhenius_diffusion_coef_Cs
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
  [freesurf_conc_Cs]
    type = DirichletBC
    variable = conc_Cs
    boundary = exterior
    value = 0.0
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 1.109e18
  []

  ##### fuel region propereties
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = 2
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fueled_region_density]
    type = ParsedMaterial
    block = 2
    property_name = density
    expression = 1750.0
  []
  [fuel_conc_Cs]
    type = ArrheniusDiffusionCoef
    block = 2
    d1 = 1e-8 # m^2/s
    q1 = 0.0 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []

  ##### Shell properties
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = 3
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = 3
    property_name = density
    expression = 1750.0
  []

  [shell_conc_Cs]
    type = ArrheniusDiffusionCoef
    block = 3
    d1 = 1e-8 # m^2/s
    q1 = 0.0 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []

  ##### pebble core
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = 1
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = 1
    property_name = density
    expression = 1400.0
  []
  [core_conc_Cs]
    type = ArrheniusDiffusionCoef
    block = 1
    d1 = 1e-8 # m^2/s
    q1 = 0.0 # J/mol
    temperature = temperature
    arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
  []
[]

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

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  line_search = 'none'

  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-9
  nl_max_its = 15

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 500000
  num_steps = 10
[]

[Postprocessors]

  ### irradiation conditions
  [pebble_power]
    type = ElementIntegralPower
    variable = temperature
    use_material_fission_rate = true
    fission_rate_material = fission_rate
    block = 2
    execute_on = 'initial timestep_end'
  []
  [heat_from_pebble]
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = exterior
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [fast_neutron_maximum_fluence_MP]
    type = ElementExtremeMaterialProperty
    mat_prop = fast_neutron_fluence
    block = '1 2 3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []
  [fast_neutron_maximum_flux_MP]
    type = ElementExtremeMaterialProperty
    mat_prop = fast_neutron_flux
    block = '1 2 3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []
  [fission_rate_maximum]
    type = ElementExtremeMaterialProperty
    mat_prop = fission_rate
    block = '1 2 3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  ### temperatures
  [fuel_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '2'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [fuel_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = '2'
    execute_on = 'initial timestep_end'
  []
  [fuel_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '2'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [shell_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '3'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [shell_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = '3'
    execute_on = 'initial timestep_end'
  []
  [shell_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '3'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []

  [core_minimum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '1'
    value_type = 'min'
    execute_on = 'initial timestep_end'
  []
  [core_average_temperature]
    type = AverageNodalVariableValue
    variable = temperature
    block = '1'
    execute_on = 'initial timestep_end'
  []
  [core_maximum_temperature]
    type = NodalExtremeValue
    variable = temperature
    block = '1'
    value_type = 'max'
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  print_linear_residuals = true
  time_step_interval =  1
  exodus = true
  csv = true
  perf_graph = true
[]

(test/tests/triso/graphite_matrix_thermal/series.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix.
#
# The SERIES model
# is employed to compute the homogenized thermal conductivity of
# the graphite matrix containing TRISO particles.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.35 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf    thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.35  1.0'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             alpha := thcond_particle/k_graphite;
             f_packing_rhs := 1.0 / ( 1.0 + pf * (1.0/alpha - 1.0) );
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.35
    etc_model = SERIES
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = 'series_out'
[]

(test/tests/triso/graphite_matrix_thermal/exact.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix
# for each graphite type: A3-3-1800, A3-3-1950, A3-27-1800, and A3-27-1950.
# Assuming zero packing of TRISO particles in the graphite matrix.
#
# The exact solution is performed through Functions/thcond_exact for
# the homogenized thermal conductivity and Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg density_avg pf   part_cp'
    symbol_values = 'temperature_avg fluence_avg density_avg 0.0 748.72'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             cp := 1.1 * 1.75e6 * (0.645 + 3.14 * thouC - 2.809 * thouC^2 + 0.959 * thouC^3) / density_avg;
             cp * (1 - pf) + part_cp * pf'
  []
  [thcond_unirradiated]
    type = ParsedFunction
    symbol_names = 'temperature_avg coef0 coef1 coef2'
    symbol_values = 'temperature_avg 50.8 1.1810e-3 -7.8453e-4'
    expression = 'coef0 * (1.0 - coef1 * (temperature_avg - 373.15) * exp(coef2 * (temperature_avg - 273.15)))'
  []
  [thcond_exact]
    type = ParsedFunction
    symbol_names = 'temperature_avg fluence_avg init_density thcond_unirradiated thcond_particle pf    thcond_scale_factor'
    symbol_values = 'temperature_avg fluence_avg 1200 thcond_unirradiated 4.13            0.0   0.9'
    expression = 'thouC := (temperature_avg - 273.15) / 1000;
             gamma := 0.94 - 0.604 * thouC;
             fluence_cut := if (fluence_avg * 0.8 * 1e-25 > 10.5, 10.5, fluence_avg * 0.8 * 1e-25);
             fluence := fluence_cut / 1.52;
             beta := 2.96 - 1.955 * thouC;
             epsilon := 0.043 * thouC - 0.008 * thouC^8;
             f_fluence := 1 - gamma * (1 - exp(-beta * fluence)) - epsilon * fluence;
             f_matrix := init_density / 1700;
             k_graphite := thcond_unirradiated * f_fluence * f_matrix;
             f_packing_rhs := 1.0;
             thcond_scale_factor * k_graphite * f_packing_rhs'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    packing_fraction = 0.0
    specific_heat_scale_factor = 1.1
    thermal_conductivity_scale_factor = 0.9
    temperature = temp
    initial_matrix_density = 1200.0
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [spheat_avg]
    type = ElementAverageValue
    variable = specific_heat
  []
  [spheat_exact]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [spheat_diff]
    type = DifferencePostprocessor
    value1 = spheat_avg
    value2 = spheat_exact
    outputs = none
  []
  [spheat_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = spheat_diff
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
  [thcond_exact]
    type = FunctionValuePostprocessor
    function = thcond_exact
    outputs = none
  []
  [thcond_diff]
    type = DifferencePostprocessor
    value1 = thcond_avg
    value2 = thcond_exact
    outputs = none
  []
  [thcond_max_diff]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = thcond_diff
  []
[]

[Outputs]
  csv = true
  file_base = exact_A3_3_1800_out
[]

(test/tests/triso_pebble/1D_pebble.i)

[GlobalParams]
  order = SECOND
  family = LAGRANGE
  flux_conversion_factor = 0.85
[]

[Mesh]
  coord_type = RSPHERICAL
  [gen]
    type = TRISO1DMeshGenerator
    elem_type = EDGE3
    coordinates = '0 0.0138 0.018 0.02'
    mesh_density = '5    10   10'
    block_names = 'core fuel shell'
  []
[]

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[Samplers]
  [sample]
    type = MonteCarlo
    num_rows = 10
    distributions = 'normal_kernel_r normal_buffer_t normal_ipyc_t normal_sic_t normal_opyc_t'
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  [sub]
    type = SamplerTransientMultiApp
    input_files = triso_1d.i
    sampler = sample
    execute_on = 'TIMESTEP_BEGIN'
    mode = normal
  []
[]

[DiracKernels]
  [vpp_point_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    sampler = sample
    fuel_element_blocks = fuel
    execute_on = 'INITIAL'
  []
[]

[Transfers]
  [released_heat_inc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = released_heat_inc
    from_postprocessor = released_heat_inc
  []
  [sub_temp_bc]
    type = SamplerPostprocessorTransfer
    from_multi_app = sub
    sampler = sample
    to_vector_postprocessor = temp_bc
    from_postprocessor = aver_temp_exterior
  []
  [temp_bc]
    type = MultiAppVectorPostprocessorTransfer
    to_multi_app = sub
    vector_postprocessor = sample_points
    postprocessor = temp_bc
    vector_name = temperature
  []
[]

[Controls]
  [cmdline]
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'kernel_radius buffer_thickness IPyC_thickness SiC_thickness OPyC_thickness'
  []
[]

[VectorPostprocessors]
  [sample_points]
    type = TRISOMonteCarloPointValueSampler
    variable = temperature
    execute_on = 'INITIAL TIMESTEP_BEGIN'
    point_source_location = point_source_location
  []
  [sampler_data]
    type = SamplerData
    execute_on = 'TIMESTEP_BEGIN'
    sampler = sample
  []
  [released_heat_inc]
    type = StochasticResults
    execute_on = 'TIMESTEP_BEGIN'
  []
  [temp_bc]
    type = StochasticResults
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
  [point_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = released_heat_inc
    value_name = released_heat_inc:released_heat_inc
    point_source_location = point_source_location
  []
  [temp_bc_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = temp_bc
    value_name = sub_temp_bc:aver_temp_exterior
    point_source_location = point_source_location
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 773.15
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = FissionRateHeatSource
    variable = temperature
    block = fuel
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = 'core fuel shell'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'core fuel shell'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 'core fuel shell'
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6'
    y = '1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 5.75e19
  []
[]

[UserObjects]
  [particle_geometry]
    type = TRISOGeometry
    outer_core = core_outer_boundary
    outer_fueled_region = fuel_outer_boundary
    outer_shell = shell_outer_boundary
    include_particle = false
    include_pebble = true
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 6.28e17
  []
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = fuel
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fuel_region_density]
    type = ParsedMaterial
    block = fuel
    property_name = density
    expression = 1750.0
  []
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = shell
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = shell
    property_name = density
    expression = 1750.0
  []
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = core
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = core
    property_name = density
    expression = 1400.0
  []
[]

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

[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-12
  nl_abs_tol = 1e-12
  nl_max_its = 15

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 20000
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

(test/tests/triso/graphite_matrix_thermal/grade_ig_110.i)

# This regression test is designed to test the thermal properties
# (specific heat capacity and thermal conductivity) of the graphite matrix
# for grade IG-110 graphite.
# Assuming zero packing of TRISO particles in the graphite matrix.
#
# The exact solution is performed through Functions/spheat_exact for the
# specific heat capacity.
#
# Differences between the BISON predictions and exact solutions are
# computed in Postprocessors. The exact and BISON predictions agree well (<1e-10).

[GlobalParams]
  displacements = 'disp_x'
[]

[Mesh]
  [unit_line]
    type = GeneratedMeshGenerator
    dim = 1
    xmax = 1e-2
  []
[]

[Variables]
  [disp_x]
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 600
  []
[]

[AuxKernels]
  [temp_ramp]
    type = FunctionAux
    variable = temp
    function = temp_fcn
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [unit_line]
    strain = FINITE
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []
  [strain_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '-t * 1e-4'
  []
[]

[Functions]
  [temp_fcn]
    type = ParsedFunction
    expression = '600 + 25 * t'
  []
  [fluence_fcn]
    type = ParsedFunction
    expression = '1e26 * t / 5'
  []
  [spheat_exact]
    type = ParsedFunction
    symbol_names = 'tk              pf  part_cp'
    symbol_values = 'temperature_avg 0.0 748.72'
    expression = 'cp := (0.6572 + 0.0026*(tk-273.15) -2.0e-6*(tk-273.15)^2) / 0.001;
             cp * (1 - pf) + part_cp * pf'
  []
[]

[Materials]
  [fluence]
    type = GenericFunctionMaterial
    prop_names = fast_neutron_fluence
    prop_values = fluence_fcn
    outputs = all
  []
  [matrix_stress]
    type = ComputeFiniteStrainElasticStress
  []
  [matrix_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.25
  []
  [matrix_thermal_props]
    type = GraphiteMatrixThermal
    flux_conversion_factor = 0.8
    graphite_grade = IG_110
    packing_fraction = 0.0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    temperature = temp
    outputs = all
  []
  [matrix_density]
    type = StrainAdjustedDensity
    strain_free_density = 1200.0
    outputs = all
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Postprocessors]
  [temperature_avg]
    type = ElementAverageValue
    variable = temp
  []
  [density_avg]
    type = ElementAverageValue
    variable = density
  []
  [fluence_avg]
    type = ElementAverageValue
    variable = fast_neutron_fluence
  []

  [avg_cP]
    type = ElementAverageValue
    variable = specific_heat
  []
  [exact_cP]
    type = FunctionValuePostprocessor
    function = spheat_exact
    outputs = none
  []
  [diff_cP]
    type = DifferencePostprocessor
    value1 = avg_cP
    value2 = exact_cP
    outputs = none
  []
  [max_diff_cP]
    type = TimeExtremeValue
    value_type = abs_max
    postprocessor = diff_cP
  []

  [thcond_avg]
    type = ElementAverageValue
    variable = thermal_conductivity
  []
[]

[Outputs]
  csv = true
  file_base = 'exact_IG_110_out'
[]

(test/tests/triso_pebble/3D_pebble_from_csv.i)

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  flux_conversion_factor = 0.85
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = pebble.e
  []
[]

[Distributions]
  [normal_kernel_r]
    type = TruncatedNormal
    mean = 213.35e-6
    standard_deviation = 4.4e-6
    lower_bound = 1.9575e-04
    upper_bound = 2.3095e-04
  []
  [normal_buffer_t]
    type = TruncatedNormal
    mean = 98.9e-6
    standard_deviation = 8.4e-6
    lower_bound = 6.53e-05
    upper_bound = 1.325e-04
  []
  [normal_ipyc_t]
    type = TruncatedNormal
    mean = 40.4e-6
    standard_deviation = 2.5e-6
    lower_bound = 3.0400e-05
    upper_bound = 5.0400e-05
  []
  [normal_sic_t]
    type = TruncatedNormal
    mean = 35.2e-6
    standard_deviation = 1.2e-6
    lower_bound = 3.0400e-05
    upper_bound = 4.0000e-05
  []
  [normal_opyc_t]
    type = TruncatedNormal
    mean = 43.4e-6
    standard_deviation = 2.9e-6
    lower_bound = 3.1800e-05
    upper_bound = 5.5000e-05
  []
[]

[DiracKernels]
  [vpp_point_source]
    type = TRISOMonteCarloPointSource
    variable = temperature
    point_source_values = heat_source
    value_name = heat_inc
    point_source_location = point_source_location
  []
[]

[UserObjects]
  [point_source_location]
    type = TRISOMonteCarloPointSourceLocation
    min_radius = 0.0138
    max_radius = 0.018
    geometry = SPHERE
    fuel_element_blocks = fuel
    execute_on = 'INITIAL'
    num_particles = 100
  []
[]

[VectorPostprocessors]
  [heat_source]
    type = TRISOMonteCarloPointSourceFromFile
    file_base_time = monte_carlo_out_time
    file_base_source = monte_carlo_out_release_heat_inc
    x_header_in_file = 'time:time'
    y_header_in_file = 'release_heat_inc:release_heat_inc'
    source_name = heat_inc
    num_files = 2
    num_particles = 100
    format = columns
    execute_on = 'INITIAL TIMESTEP_BEGIN'
  []
  [point_source_output]
    type = TRISOMonteCarloPointSourceOutput
    vector_postprocessor = heat_source
    value_name = heat_inc
    point_source_location = point_source_location
  []
[]

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

[Variables]
  [temperature]
    initial_condition = 773.15
  []
[]

[AuxVariables]
  [density]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [density]
    type = MaterialRealAux
    variable = density
    property = density
    block = 'core fuel shell'
    execute_on = 'initial linear'
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 'core fuel shell'
    execute_on = timestep_end
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 'core fuel shell'
    execute_on = timestep_end
  []
[]

[BCs]
  [coolant_temp]
    type = DirichletBC
    variable = temperature
    boundary = exterior
    value = 773.15
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 76e6'
    y = '1 1'
  []
  [fission_rate]
    type = LinearCombinationFunction
    functions = power_history
    w = 5.75e19
  []
[]

[Materials]
  [fission_rate]
    type = GenericFunctionMaterial
    prop_names = fission_rate
    prop_values = fission_rate
  []
  [fast_neutron_flux]
    type = FastNeutronFlux
    calculate_fluence = true
    flux_function = power_history
    factor = 6.28e17
  []
  [fueled_region_thermal]
    type = GraphiteMatrixThermal
    block = fuel
    graphite_grade = A3_27_1800
    packing_fraction = 0.22
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [fuel_region_density]
    type = ParsedMaterial
    block = fuel
    property_name = density
    expression = 1750.0
  []
  [shell_region_thermal]
    type = GraphiteMatrixThermal
    block = shell
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1750.0
  []
  [shell_region_density]
    type = ParsedMaterial
    block = shell
    property_name = density
    expression = 1750.0
  []
  [core_region_thermal]
    type = GraphiteMatrixThermal
    block = core
    graphite_grade = A3_27_1800
    packing_fraction = 0.0
    temperature = temperature
    initial_matrix_density = 1400.0
  []
  [core_region_density]
    type = ParsedMaterial
    block = core
    property_name = density
    expression = 1400.0
  []
[]

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

[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-9
  nl_abs_tol = 1e-9
  nl_max_its = 15

  l_tol = 1e-5
  l_max_its = 30

  start_time = 0.0
  dt = 10
  num_steps = 2
[]

[Outputs]
  print_linear_residuals = true
  time_step_interval =  1
  exodus = true
  perf_graph = true
[]

Description
Example Input Syntax
Input Parameters
Input Files
References