Griesheimer and Kooreman Analytic Benchmark

Problem Description

This test case models the benchmark problem proposed in Griesheimer and Kooreman (2022), which verifies coupled neutronics, heat conduction, and thermal expansion. The problem consists of a 1D multiplying slab where neutrons are restricted to travel along either $var element = document.getElementById("moose-equation-b6ce04d5-b02d-468f-a067-7288c423bde9");katex.render("\\hat{\\Omega} = (-1, 0, 0)", element, {displayMode:false,throwOnError:false});$ or $var element = document.getElementById("moose-equation-8e422868-9d70-4b5b-89de-b367f790c7c6");katex.render("\\hat{\\Omega} = (1, 0, 0)", element, {displayMode:false,throwOnError:false});$ , and vacuum boundary conditions are applied at $var element = document.getElementById("moose-equation-d2b5ad46-4e22-4b58-8b2c-0a797e3d74c8");katex.render("x = \\pm L/2", element, {displayMode:false,throwOnError:false});$ (where $var element = document.getElementById("moose-equation-10ef438f-8b7f-4f60-8233-69009a539965");katex.render("L", element, {displayMode:false,throwOnError:false});$ is the deformed length of the slab). The S₂ neutron transport equation may be written as a neutron diffusion equation under these constraints:

(1)var element = document.getElementById("moose-equation-d5891c3a-0f02-47c5-896f-0afa82043fd8");katex.render("\\frac{d}{dx}\\left(\\frac{1}{\\Sigma_{t}(x)}\\frac{d\\phi(x)}{dx}\\right) + \\left(\\lambda - 1\\right)\\Sigma_{t}(x)\\phi(x) = 0\\text{,}", element, {displayMode:true,throwOnError:false});

with the following Marshak boundary condition:

(2)var element = document.getElementById("moose-equation-0aca02b5-5760-4b30-9280-0bf21628a6ef");katex.render("\\frac{d\\phi(x)}{dx} + \\Sigma_{t}(x)\\phi(x),\\,\\,\\,\\, x = \\pm \\frac{L}{2}\\text{.}", element, {displayMode:true,throwOnError:false});

$var element = document.getElementById("moose-equation-a8c727e5-fcbc-4c1d-a289-c77c350e564a");katex.render("\\phi(x)", element, {displayMode:false,throwOnError:false});$ is the scalar neutron flux (cm $var element = document.getElementById("moose-equation-a129432e-bdd3-4940-b9fd-006355b6aa8c");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-f88869b5-30b7-48b6-bf3b-ae3776453228");katex.render("s^{-1}", element, {displayMode:false,throwOnError:false});$ ). $var element = document.getElementById("moose-equation-a1a352d4-276f-4332-98b0-a0db0cf6e19f");katex.render("\\Sigma_{t}(x)", element, {displayMode:false,throwOnError:false});$ is the total cross section (cm $var element = document.getElementById("moose-equation-dafc2673-0bef-4be4-8d78-8a9fc5470e99");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ), which has an inverse dependence on temperature:

(3)var element = document.getElementById("moose-equation-f18fcdc8-fe04-41dd-92f8-360dc0f4f756");katex.render("\\Sigma_{t}(x) = \\Sigma_{t,0}\\frac{T_{0}}{T(x)}\\text{,}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-d002a360-75d5-4749-8b72-3863222eb378");katex.render("T_{0}", element, {displayMode:false,throwOnError:false});$ is the reference temperature of the slab (K) and $var element = document.getElementById("moose-equation-2fec800a-85bf-4a5d-a717-6c8d92171328");katex.render("T(x)", element, {displayMode:false,throwOnError:false});$ is the slab temperature (K). $var element = document.getElementById("moose-equation-e308d0f2-a7a0-4a2d-a4a2-9452cba4fbc3");katex.render("\\lambda", element, {displayMode:false,throwOnError:false});$ is the combined scattering and k-eigenvalue fission source:

(4)var element = document.getElementById("moose-equation-040dea85-3c7a-43a3-8ac8-d36e75380e7e");katex.render("\\lambda = \\frac{1}{k_{eff}}\\nu\\frac{\\Sigma_{f}}{\\Sigma_{t}} + \\frac{\\Sigma_{s}}{\\Sigma_{t}}\\text{,}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-a14fa50c-7a66-4ba6-8a6d-044773a23679");katex.render("k_{eff}", element, {displayMode:false,throwOnError:false});$ is the k-eigenvalue, $var element = document.getElementById("moose-equation-ea801026-d150-4a67-88d2-cc1a77bd4122");katex.render("\\nu", element, {displayMode:false,throwOnError:false});$ is the average number of neutrons released per fission reaction, $var element = document.getElementById("moose-equation-c8225d57-d911-4bef-b417-26dd5d07efca");katex.render("\\Sigma_{f}/\\Sigma_{t}", element, {displayMode:false,throwOnError:false});$ is the fractional probability of a fission reaction and $var element = document.getElementById("moose-equation-bd8f0a68-353a-4fb6-aa04-e1b9db90ce77");katex.render("\\Sigma_{s}/\\Sigma_{t}", element, {displayMode:false,throwOnError:false});$ is the fractional probability of a scattering reaction. Both reaction probabilities are constant over the length of the slab. $var element = document.getElementById("moose-equation-3108df3f-db35-4730-9866-c79bccd1b0d7");katex.render("\\Sigma_{t,0}", element, {displayMode:false,throwOnError:false});$ is defined as:

(5)var element = document.getElementById("moose-equation-f59e21d4-88a6-42eb-b2e6-57afec9cdf27");katex.render("\\Sigma_{t,0} = \\frac{N_{A}}{A}\\rho_{0}\\sigma_{t,0}\\text{,}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-9d3dd1da-333f-46dd-b59b-dc5238f70aa0");katex.render("N_{A}", element, {displayMode:false,throwOnError:false});$ is Avogadro's number, $var element = document.getElementById("moose-equation-155e8033-045c-4716-a0ef-99c6bca391df");katex.render("A", element, {displayMode:false,throwOnError:false});$ is the molar mass of the slab (g mol $var element = document.getElementById("moose-equation-a904fd4c-be1f-4f5a-9a3d-e654d2793a0e");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ), $var element = document.getElementById("moose-equation-10e69aa5-3705-4a56-b6bc-944063f59dea");katex.render("\\rho_{0}", element, {displayMode:false,throwOnError:false});$ is the unheated slab density (g cm $var element = document.getElementById("moose-equation-2a4e8388-0c11-4de0-b3a4-bfd10599203b");katex.render("^{-3}", element, {displayMode:false,throwOnError:false});$ ), and $var element = document.getElementById("moose-equation-f8989ed2-3454-450d-a0df-60e1d652dec7");katex.render("\\sigma_{t,0}", element, {displayMode:false,throwOnError:false});$ is the reference microscopic total cross section (cm $var element = document.getElementById("moose-equation-3f05a9aa-4c54-48c9-9c97-6ed8b1bb2d60");katex.render("^2", element, {displayMode:false,throwOnError:false});$ ). The neutron flux solution is normalized to a chosen power $var element = document.getElementById("moose-equation-fa386f45-87cf-4f11-b3ce-565c1c01ebaa");katex.render("P", element, {displayMode:false,throwOnError:false});$ (eV s $var element = document.getElementById("moose-equation-d3a8497b-ddf3-4080-af08-7867021fb066");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ):

(6)var element = document.getElementById("moose-equation-d88b3b3a-c338-4880-9e5e-729ac18b1d09");katex.render("P = \\int_{-L/2}^{L/2}q\\Sigma_{t}(x)\\phi(x)\\,dx", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-348609ee-f9bd-49ed-9fb4-33906fed3a8f");katex.render("q", element, {displayMode:false,throwOnError:false});$ is the energy released per neutron interaction (eV). The temperature distribution in the slab is governed by the following heat conduction equation:

(7)var element = document.getElementById("moose-equation-a3257d9b-3056-465e-b5d8-33dbf58061e5");katex.render("\\frac{d}{dx}\\left(\\kappa(x)\\frac{d T(x)}{dx}\\right) + q\\Sigma_{t}(x)\\phi(x) = 0\\text{,}", element, {displayMode:true,throwOnError:false});

with the following convection boundary condition:

(8)var element = document.getElementById("moose-equation-f5bbd307-9b3f-44b3-839a-34d7e6105fdf");katex.render("\\kappa(x)\\frac{dT(x)}{dx} + h\\left(T(x) - T_{0}\\right),\\,\\,\\,\\, x = \\pm \\frac{L}{2}\\text{.}", element, {displayMode:true,throwOnError:false});

$var element = document.getElementById("moose-equation-32432214-9c27-4cc4-8c02-267be1a74bdc");katex.render("\\kappa(x)", element, {displayMode:false,throwOnError:false});$ is the thermal conductivity of the material (eV s $var element = document.getElementById("moose-equation-fe7fc6d2-a0fe-40d9-90b0-e80d14f59711");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ cm^{$var element = document.getElementById("moose-equation-84daaa39-cc8a-41b4-b3f0-62c9769c4275");katex.render("{-1}", element, {displayMode:false,throwOnError:false});$ K} $var element = document.getElementById("moose-equation-9901bb63-d465-402e-b220-a546c8b2a985");katex.render("{-1}", element, {displayMode:false,throwOnError:false});$ ), which has a linear dependence on temperature:

(9)var element = document.getElementById("moose-equation-5adae1f4-3aeb-4dfe-bf4b-7868e2a7a38b");katex.render("\\kappa(x) = \\kappa_{0}T(x)", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-0d5278ce-f153-4345-be72-1d23ec181229");katex.render("\\kappa_{0}", element, {displayMode:false,throwOnError:false});$ is the reference thermal conductivity (eV s $var element = document.getElementById("moose-equation-ea69205b-eded-43b6-857a-83d013035f51");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ cm $var element = document.getElementById("moose-equation-332376e0-c64a-49d7-b4ab-117a045a4f05");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ K $var element = document.getElementById("moose-equation-d56c6e9d-5cd3-41ec-8a60-82728b21d0ef");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ ). $var element = document.getElementById("moose-equation-d9c9538d-e002-4e04-ae40-8d9c36347fde");katex.render("h", element, {displayMode:false,throwOnError:false});$ is a convective heat transfer coefficient (eV s $var element = document.getElementById("moose-equation-a10d7a57-8cc6-4812-9e75-aacfd18af62f");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ cm $var element = document.getElementById("moose-equation-0ca9ba43-5f1a-4e19-8e84-09a50fe91489");katex.render("^{-2}", element, {displayMode:false,throwOnError:false});$ K $var element = document.getElementById("moose-equation-061e3bda-c784-4dc2-9f9b-5a0db0132f1e");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ). The final set of physics considered in this problem is thermal expansion, where the differential slab length $var element = document.getElementById("moose-equation-13f16f96-a688-4559-b216-2a8afa278c9d");katex.render("\\ell(x)", element, {displayMode:false,throwOnError:false});$ (cm) expands according to the local thermal strain $var element = document.getElementById("moose-equation-c681923e-1f93-4275-9dc0-cc7dfa3d49e0");katex.render("\\epsilon_{x}(x)", element, {displayMode:false,throwOnError:false});$ :

(10)var element = document.getElementById("moose-equation-9d8ef5fd-28ec-4fc3-bccb-6c24bbeb992f");katex.render("\\ell(x) = \\ell_{0}(1 + \\epsilon_{x}(x))\\text{,}", element, {displayMode:true,throwOnError:false});

(11)var element = document.getElementById("moose-equation-9ddda018-ce00-4f29-a8a6-6e3aa203c7fd");katex.render("\\epsilon_{x}(x) = \\int_{T_{0}}^{T(x)}\\alpha(\\tau)\\, d\\tau\\text{,}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-e5b80a4a-55cb-45df-a6c3-5163d28b768f");katex.render("\\ell_{0}", element, {displayMode:false,throwOnError:false});$ is the initial slab length, and the linear expansion coefficient $var element = document.getElementById("moose-equation-8c6b714f-d253-4026-afb8-01e3c3626ed3");katex.render("\\alpha(T)", element, {displayMode:false,throwOnError:false});$ (K $var element = document.getElementById("moose-equation-31237b16-2a01-4b08-914f-4f4fb10d217a");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ ) is:

(12)var element = document.getElementById("moose-equation-6cfb80b0-ce23-4658-a8f0-177a953a8812");katex.render("\\alpha(T) = \\frac{1}{2\\sqrt{T_{0}T(x)}}\\text{.}", element, {displayMode:true,throwOnError:false});

Mass is conserved by ensuring that the following expression is satisfied:

(13)var element = document.getElementById("moose-equation-e446629c-5b95-487e-8a76-214f40e63c3e");katex.render("\\rho(x)\\ell(x) = \\rho_{0}\\ell_{0}\\text{.}", element, {displayMode:true,throwOnError:false});

Analytical Solutions

The analytical solutions to the coupled problem were derived by Griesheimer and Kooreman (2022) under the ansatz that the temperature and neutron flux fields are proportional to each other. To ensure this ansatz is not violated, $var element = document.getElementById("moose-equation-dfb6ea31-dd7d-4eae-9d1e-544fce86c442");katex.render("\\sigma_{t,0}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-c4401b03-f3f9-4fb5-b9d4-558dbcd16ba5");katex.render("h", element, {displayMode:false,throwOnError:false});$ are manufactured based on other benchmark parameters. For details regarding the derivation of the benchmark solution, readers are referred to Griesheimer and Kooreman (2022).

The analytic solutions for the neutron flux (Eq. (14)), temperature (Eq. (15)), and power density (Eq. (16)) can be found below:

(14)var element = document.getElementById("moose-equation-ce81aa86-b7b9-4fbd-8411-7f7f0fe96d0e");katex.render("\\phi(x) = \\phi(0)\\sqrt{1 - \\frac{(\\lambda - 1)P^2 x^2}{L^2 q^2 \\phi^2(0)}}\\text{,}", element, {displayMode:true,throwOnError:false});

(15)var element = document.getElementById("moose-equation-83a7cd9f-f5d6-4435-b1ab-f8f216c83643");katex.render("T(x) = \\frac{q \\Sigma_{t,0} T_{0} L}{P}\\phi(x)\\text{,}", element, {displayMode:true,throwOnError:false});

and

(16)var element = document.getElementById("moose-equation-1a4b8f87-edd3-4344-ae57-55f91aa21e14");katex.render("\\Sigma_{t}(x)\\phi(x)q = \\frac{P}{L}\\text{.}", element, {displayMode:true,throwOnError:false});

The neutron flux in the center of the slab ( $var element = document.getElementById("moose-equation-cf4b5563-1902-4f78-865b-bc41028658b8");katex.render("\\phi(0)", element, {displayMode:false,throwOnError:false});$ ) is arbitrary, and is therefore considered an input setting. The slab eigenvalue can be calculated with Eq. (17):

(17)var element = document.getElementById("moose-equation-8ca3689f-72b7-4b77-bc34-4220b809b7eb");katex.render("\\lambda = \\frac{1}{2}\\left(1 + \\sqrt{1 + \\frac{16 q^2 \\phi^2(0)}{P^2}}\\right)\\text{,}", element, {displayMode:true,throwOnError:false});

and the expanded slab length can be calculated with Eq. (18):

(18)var element = document.getElementById("moose-equation-69b4121c-b581-456e-a9d0-0d47c162f05f");katex.render("L = \\sqrt[3]{\\left(\\frac{L_{0}}{_2F_1\\left(1/4,1/2,3/2,1/\\lambda\\right)}\\right)^{4}\\left(\\frac{\\frac{q^2 \\phi^2(0)}{P(\\lambda - 1)} - \\frac{P}{4}}{\\kappa_{0}T_{0}^2(1 - 1/\\lambda)}\\right)}\\text{,}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-18265d39-1742-4311-9c87-550197634aae");katex.render("_2F_1(a,b,c,d)", element, {displayMode:false,throwOnError:false});$ is the Gauss hypergeometric function. The manufactured benchmark parameters are calculated with Eq. (19) and Eq. (20):

(19)var element = document.getElementById("moose-equation-5859da00-5355-4e76-968d-be153dff667e");katex.render("\\sigma_{t, 0} = \\frac{1}{M\\rho_{0}T_{0}}\\sqrt{\\frac{P}{\\kappa_{0}L(\\lambda - 1)}}\\text{;}", element, {displayMode:true,throwOnError:false});

(20)var element = document.getElementById("moose-equation-7b256703-1092-47d7-9bf7-ef84cf7e5e4b");katex.render("h = \\left(\\sqrt{\\frac{L(\\lambda - 1)}{\\kappa_{0}P}} - \\frac{2T_{0}}{P}\\right)^{-1}\\text{.}", element, {displayMode:true,throwOnError:false});

To ensure positivity of the solution, the power must be chosen such that:

(21)var element = document.getElementById("moose-equation-7e034e19-9243-4b22-85be-57005c5f1b71");katex.render("P \\geq \\frac{4 \\kappa_{0}T_{0}^2}{L(\\lambda - 1)}\\text{.}", element, {displayMode:true,throwOnError:false});

The canonical parameters of the benchmark proposed by Griesheimer and Kooreman (2022) can be found in Table 1. These parameters result in an expanded slab length of $var element = document.getElementById("moose-equation-6ed739c9-c3fb-44f2-af63-060608f8c0d5");katex.render("L = 106.47", element, {displayMode:false,throwOnError:false});$ cm and a slab k-eigenvalue of $var element = document.getElementById("moose-equation-2500a616-15f0-45f4-8a8b-8af1de428ed1");katex.render("k_{eff} = 0.29557", element, {displayMode:false,throwOnError:false});$ .

Table 1: Thermal expansion benchmark parameters proposed in Griesheimer and Kooreman (2022).

Parameter	Value
$var element = document.getElementById("moose-equation-f978c82d-f7b3-406a-afef-a93ddf0b43d4");katex.render("\\rho_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-675aca02-65ac-4b1e-b91c-d28c3b17ae64");katex.render("1.2", element, {displayMode:false,throwOnError:false});$ (g cm^-3)
$var element = document.getElementById("moose-equation-18f5dc1b-9346-4630-befd-26d742c2b4bb");katex.render("L_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-3b3864a5-861f-4408-a0a8-0fcec65e6cfa");katex.render("100", element, {displayMode:false,throwOnError:false});$ (cm)
$var element = document.getElementById("moose-equation-38a29251-c88c-4a42-807d-df7062f8f025");katex.render("A", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-7505ceed-54a2-4b1e-abec-49c1d4320ccc");katex.render("180", element, {displayMode:false,throwOnError:false});$ (g mol^-1)
$var element = document.getElementById("moose-equation-0e8aebbc-d599-43ad-bfea-65c797ca5fbb");katex.render("T_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-2d8263d4-86bd-4e2a-88a6-381d8084ab34");katex.render("293", element, {displayMode:false,throwOnError:false});$ (K)
$var element = document.getElementById("moose-equation-59d560a2-f051-4c23-a5df-06d769a95547");katex.render("q", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-aa894e62-143b-4846-8482-32c4b39c5481");katex.render("1\\times 10^{8}", element, {displayMode:false,throwOnError:false});$ (eV)
$var element = document.getElementById("moose-equation-da8eddd7-bef4-4eb8-89bf-1765897086cf");katex.render("P", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-0930fdf6-e227-48d4-a066-2f56b8977f18");katex.render("1\\times 10^{22}", element, {displayMode:false,throwOnError:false});$ (eV s^-1)
$var element = document.getElementById("moose-equation-10f140bf-2da2-4fe1-a68a-1a6cde6cdef5");katex.render("\\kappa_{0}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-f289f524-c760-44cb-bef1-a967907d2e0b");katex.render("1.25\\times 10^{19}", element, {displayMode:false,throwOnError:false});$ (eV s^-1 K^-2)
$var element = document.getElementById("moose-equation-92a7c5a7-0795-4647-9de5-e075aca083a0");katex.render("\\phi(0)", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-c9129ba8-f7a4-427f-bf54-8736660ea6b0");katex.render("2.5\\times 10^{14}", element, {displayMode:false,throwOnError:false});$ (cm^-2 s^-1)
$var element = document.getElementById("moose-equation-42e155e1-88e6-4482-9038-f3f84476c17d");katex.render("\\nu\\Sigma_{f}/\\Sigma_{t}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-3446a144-caa7-48cd-a87b-c4aa97f8de49");katex.render("1.5", element, {displayMode:false,throwOnError:false});$
$var element = document.getElementById("moose-equation-3f78de36-06fa-42fc-9ef7-0ce4d8c8dd54");katex.render("\\Sigma_{s}/\\Sigma_{t}", element, {displayMode:false,throwOnError:false});$	$var element = document.getElementById("moose-equation-9b59d0c2-c0e8-4aa6-b217-ca0e3ef0b5c9");katex.render("0.45", element, {displayMode:false,throwOnError:false});$

Cardinal Multi-Group Monte Carlo Model

Eltawila et al. (2026) used the Griesheimer and Kooreman benchmark problem to verify Cardinal's thermomechanical-neutronics coupling to OpenMC. Those benchmark results required several small modifications to OpenMC's multi-group Monte Carlo solver to enable S₂ neutron transport, precluding the inclusion of input files in Cardinal's regression suite. Readers interested in these results are referred to Eltawila et al. (2026) for more information. Cardinal was found to converge to the analytic solutions with a linear rate of convergence as cells, mesh elements, and temperature grid points were refined. Excellent agreement between Cardinal and the analytic solution was obtained in the case of maximal refinement.

construction:To run this model you need to build your application with the S

var element = document.getElementById("moose-equation-a7def132-58cd-4f6f-94fb-9ed5947bb185");katex.render("_2", element, {displayMode:false,throwOnError:false});

OpenMC patch

This model is using S $var element = document.getElementById("moose-equation-8d377247-e110-4a12-92e0-ca97da904aec");katex.render("_2", element, {displayMode:false,throwOnError:false});$ neutron transport to compare with the 1D analytic benchmark. To replicate the results, you have to build your application with this patch of OpenMC. You can do this simply with:

cd cardinal/contrib/openmc
git remote add s2patch https://github.com/meltawila/openmc_S2transport.git
git fetch s2patch
git cherry-pick d34189845b394895ed210008218c40bc754bce25

However, we no do not frequently rebase this branch of OpenMC off the latest version of OpenMC, so you may encounter other compatibility issues (usually small).

Cardinal is used to couple OpenMC with MOOSE via Picard iteration. OpenMC tallies neutron flux, heating, and the multiplication factor. A conceptual depiction of the meshes used as well as the data transfers that occur on each Picard iteration is shown in Figure 1. Three different meshes/geometries are involved. The thermomechanics mesh is a simple 1D mesh of $var element = document.getElementById("moose-equation-8cdd6cc8-5111-4b15-8c3a-4efbc9af0f2e");katex.render("N", element, {displayMode:false,throwOnError:false});$ EDGE2 elements; the thermomechanical physics do not require higher dimensional (e.g. 2D or 3D) meshes for this problem. The OpenMC geometry is represented using DAGMC, where cells are bounded by a triangulated surface (TRI3 elements). To facilitate data transfers between OpenMC and MOOSE thermomechanics, there is also an intermediate volumetric data transfer mesh composed of TET4 elements; no physics solves happen on this intermediate mesh. Figure 1 highlights one element (for MOOSE and the Cardinal intermediate data transfer mesh) and one OpenMC cell each in blue, with adjacent portions of the geometry shown in gray for context.

Figure 1: Data transfers between OpenMC and MOOSE and conceptual depiction of the different meshes used.

We show the input files for a coarse mesh here with more details on results convergence shown in the next section.

The script used to generate the OpenMC model is as follows:

import openmc
import openmc.mgxs as mgxs
import numpy as np
import h5py

T0 = 293 # K
Tmax = 390 # K
Tmin = 290 # K
L0 = 100 # cm
L = 106.47 # cm
rho = 1.2 # g/cc
N_A = 6.022e23 # Avagadro's number
A = 180 # mass number for slab material
num_dens = rho*N_A/A
P = 1.0e22 # eV/s
q = 1e8 # eV
k0= 1.25e19 # eV/(s-cm-K^2)
phi0 = 2.5e14 # 1/s-cm^2
s = 0.45 # Sigma_s/Sigma_t
f = 1.5 # nu Sigma_f/Sigma_t
nu = f/(1-s)
lam = 0.5*(1+np.sqrt(1+(16*q*q*phi0*phi0)/(P*P)))
Sig_t0 = np.sqrt(P/((lam-1)*k0*L))/(T0) # 1/cm
sig_t0 = Sig_t0/num_dens # cm^2
N = 20
infdim = 0.5

model = openmc.Model()

dagmc_univ = openmc.DAGMCUniverse(filename="slab20v.h5m")
model.geometry = openmc.Geometry(root=dagmc_univ)
cells = dagmc_univ.get_all_cells().values()

slab = openmc.Material(1, "slab")
slab.set_density('macro',1.)
slab.add_macroscopic('slab_xs')

model.materials = openmc.Materials([slab])
model.materials.cross_sections = 'one_group.h5'

groups = mgxs.EnergyGroups(group_edges=[0.0, 20.0e6])
NT = 100
temps = np.linspace(Tmin,Tmax,NT)
xsdata = openmc.XSdata('slab_xs', energy_groups=groups, temperatures=temps, num_delayed_groups=0)
xsdata.order = 0

for T in temps:
    Sig_t = (Sig_t0 * T0) / T
    Sig_s = s*Sig_t
    nu_Sig_f = f*Sig_t
    Sig_f = nu_Sig_f / nu
    Sig_a = Sig_f
    xsdata.set_total(np.array([Sig_t]),temperature=T)
    xsdata.set_scatter_matrix(np.array([[[Sig_s]]]),temperature=T)
    xsdata.set_absorption(np.array([Sig_a]),temperature=T)
    xsdata.set_fission(np.array([Sig_f]),temperature=T)
    xsdata.set_nu_fission(np.array([nu_Sig_f]),temperature=T)
    xsdata.set_kappa_fission(np.array([q*Sig_t]),temperature=T)

one_g_XS_file = openmc.MGXSLibrary(groups)
one_g_XS_file.add_xsdata(xsdata)
one_g_XS_file.export_to_hdf5('one_group.h5')

model.settings = openmc.Settings()
batches = 150
inactive = 50
particles = 50000
model.settings.energy_mode = 'multi-group'
model.settings.batches = batches
model.settings.inactive = inactive
model.settings.particles = particles
model.settings.output = {'tallies': True,'summary':False}

bounds = [-60, -infdim, -infdim, 60, infdim, infdim]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
mu = openmc.stats.Discrete([0.0],[1.0])
phi = openmc.stats.Discrete([0,np.pi],[0.5,0.5])
bidirectional_x = openmc.stats.PolarAzimuthal(mu=mu,phi=phi)
model.settings.source = openmc.IndependentSource(space=uniform_dist,angle=bidirectional_x)
model.settings.temperature = {'default': (Tmin+Tmax)/2,
                        'method': 'nearest',
                        'tolerance': 50,
                        'range': (Tmin, Tmax)}

model.export_to_model_xml()

The neutronics input file is as follows:

P = ${units 1.0e22 eV/s}
eV_to_J = ${units 1.602176634e-19 J/eV}
initial_temperature = ${units 293 K}

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [file]
    type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../source/meshgenerators/FileMeshGenerator.html"}>>>
    file<<<{"description": "The filename to read."}>>> = slab20v.e
  []
[]

[GlobalParams<<<{"href": "../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  use_displaced_mesh = true
[]

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.0
  []
  [disp_z]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.0
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [bins]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [bins]
    type = SkinnedBins<<<{"description": "Bins created by a skinner", "href": "../source/auxkernels/SkinnedBins.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bins
    skinner<<<{"description": "MOAB mesh skinner"}>>> = moab
  []
[]

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [temperature]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temperature
    function<<<{"description": "The initial condition function."}>>> = ${initial_temperature}
  []
[]

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = OpenMCCellAverageProblem
  verbose = true
  source_rate_normalization = 'kappa_fission'
  inactive_batches = 150
  batches = 300
  particles = 10000000
  power = '${fparse P*eV_to_J}'
  cell_level = 0
  temperature_variables = 'temperature'
  temperature_blocks = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
  volume_calculation = vol
  skinner = moab
  [Tallies<<<{"href": "../syntax/Problem/Tallies/index.html"}>>>]
    [kappa_fission]
      type = CellTally<<<{"description": "A class which implements distributed cell tallies.", "href": "../source/tallies/CellTally.html"}>>>
      estimator<<<{"description": "Type of tally estimator to use in OpenMC"}>>> = collision
      block<<<{"description": "Subdomains for which to add tallies in OpenMC. If not provided, tallies will be applied over the entire domain corresponding to the [Mesh] block."}>>> = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
      name<<<{"description": "Auxiliary variable name(s) to use for OpenMC tallies. If not specified, defaults to the names of the scores"}>>> = kappa_fission
      output<<<{"description": "UNRELAXED field(s) to output from OpenMC for each tally score. unrelaxed_tally_std_dev will write the standard deviation of each tally into auxiliary variables named *_std_dev. unrelaxed_tally_rel_error will write the relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) of each tally into auxiliary variables named *_rel_error. unrelaxed_tally will write the raw unrelaxed tally into auxiliary variables named *_raw (replace * with 'name')."}>>> = unrelaxed_tally_std_dev
    []
    [flux]
      type = CellTally<<<{"description": "A class which implements distributed cell tallies.", "href": "../source/tallies/CellTally.html"}>>>
      estimator<<<{"description": "Type of tally estimator to use in OpenMC"}>>> = collision
      block<<<{"description": "Subdomains for which to add tallies in OpenMC. If not provided, tallies will be applied over the entire domain corresponding to the [Mesh] block."}>>> = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
      score<<<{"description": "Score(s) to use in the OpenMC tallies. If not specified, defaults to 'kappa_fission'"}>>> = flux
      output<<<{"description": "UNRELAXED field(s) to output from OpenMC for each tally score. unrelaxed_tally_std_dev will write the standard deviation of each tally into auxiliary variables named *_std_dev. unrelaxed_tally_rel_error will write the relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) of each tally into auxiliary variables named *_rel_error. unrelaxed_tally will write the raw unrelaxed tally into auxiliary variables named *_raw (replace * with 'name')."}>>> = unrelaxed_tally_std_dev
    []
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  start_time = 0.0
  end_time = 1e5
  dt = 50
[]

[MultiApps<<<{"href": "../syntax/MultiApps/index.html"}>>>]
  [solid]
    type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../source/multiapps/TransientMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'solid.i'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Transfers<<<{"href": "../syntax/Transfers/index.html"}>>>]
  [kappa_fission_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer<<<{"description": "Transfers field data at the MultiApp position using the finite element shape functions from the origin application.", "href": "../source/transfers/MultiAppGeneralFieldShapeEvaluationTransfer.html"}>>>
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = solid
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = kappa_fission
    source_variable<<<{"description": "The variable to transfer from."}>>> = kappa_fission
    from_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the from-app to evaluate an adjusting factor."}>>> = heat_source
    to_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the to-app to evaluate an adjusting factor."}>>> = source_integral
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    displaced_source_mesh<<<{"description": "Whether or not to use the displaced mesh for the source mesh."}>>> = true
    displaced_target_mesh<<<{"description": "Whether or not to use the displaced mesh for the target mesh."}>>> = true
  []
  [temperature_from_solid]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = solid
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = temperature
    source_variable<<<{"description": "The variable to transfer from."}>>> = temperature_praux
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    displaced_source_mesh<<<{"description": "Whether or not to use the displaced mesh for the source mesh."}>>> = true
    displaced_target_mesh<<<{"description": "Whether or not to use the displaced mesh for the target mesh."}>>> = true
  []
  [dispx]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = disp_praux
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = disp_x
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = solid
    displaced_source_mesh<<<{"description": "Whether or not to use the displaced mesh for the source mesh."}>>> = true
    displaced_target_mesh<<<{"description": "Whether or not to use the displaced mesh for the target mesh."}>>> = true
  []
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [heat_source]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = kappa_fission
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'transfer initial timestep_end'
  []
  [k]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
  []
  [k_std_dev]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
    output<<<{"description": "The value to output. Options are $k_{eff}$ (mean), the standard deviation of $k_{eff}$ (std_dev), or the relative error of $k_{eff}$ (rel_err)."}>>> = 'std_dev'
  []
[]

[VectorPostprocessors<<<{"href": "../syntax/VectorPostprocessors/index.html"}>>>]
  [temperature]
    type = ElementValueSampler<<<{"description": "Samples values of variables on elements.", "href": "../source/vectorpostprocessors/ElementValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'temperature'
    sort_by<<<{"description": "What to sort the samples by. Options include 'x', 'y', 'z', 'id', and the name of any of the sampled quantities (which each create a vector of the same name)."}>>> = x
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'FINAL'
  []
[]

[UserObjects<<<{"href": "../syntax/UserObjects/index.html"}>>>]
  [vol]
    type = OpenMCVolumeCalculation<<<{"description": "Stochastic volume calculation of the OpenMC cells which map to MOOSE", "href": "../source/userobjects/OpenMCVolumeCalculation.html"}>>>
    n_samples<<<{"description": "Number of samples to use for the stochastic volume calculation"}>>> = 200000
  []
  [moab]
    type = MoabSkinner<<<{"description": "Re-generate the OpenMC geometry on-the-fly according to changes in the mesh geometry and/or contours in temperature and density", "href": "../source/userobjects/MoabSkinner.html"}>>>
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin'
    verbose<<<{"description": "Whether to print diagnostic information"}>>> = true
    temperature<<<{"description": "Temperature variable by which to bin elements"}>>> = temperature
    n_temperature_bins<<<{"description": "Number of temperature bins"}>>> = 20.0
    temperature_min<<<{"description": "Lower bound of temperature bins"}>>> = 0
    temperature_max<<<{"description": "Upper bound of temperature bins"}>>> = 1000
    build_graveyard<<<{"description": "Whether to build a graveyard around the geometry"}>>> = true
    graveyard_scale_inner<<<{"description": "Multiplier on mesh bounding box to form inner graveyard surface"}>>> = 1.5
    graveyard_scale_outer<<<{"description": "Multiplier on mesh bounding box to form outer graveyard surface"}>>> = 1.6
  []
[]

The thermomechanics input file is as follows:

T0 = ${units 293 K}
L0 = ${units 100 cm}
L = ${units 106.47 cm}
P = ${units 1.0e22 eV/s}
q = ${units 1e8 eV}
k0 = ${units 1.25e19 eV/s/cm/K^2} # k(T) = k0 T(x)
phi0 = ${units 2.5e14 1/s/cm^2} # (flux at the origin)
eV_to_J = ${units 1.602176634e-19 J/eV}
lam = '${fparse 0.5*(1+sqrt(1+(16*q*q*phi0*phi0)/(P*P)))}' # eigenvalue solution
h = '${fparse 1/(sqrt(L*(lam-1)/(k0*P)) - (2*T0)/(P))*eV_to_J }' # W/cm^2-K

E = ${units 1e6 N/m^2}
nu = 0.0
rho = ${units 20 g/cm^3}
c_p = ${units 0.3 J/g/K}

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [generated]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 1
    nx<<<{"description": "Number of elements in the X direction"}>>> = 500
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = -50
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 50
  []
  [pin]
    type = ExtraNodesetGenerator<<<{"description": "Creates a new node set and a new boundary made with the nodes the user provides.", "href": "../source/meshgenerators/ExtraNodesetGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = generated
    new_boundary<<<{"description": "The names of the boundaries to create"}>>> = pin
    coord<<<{"description": "The nodes with coordinates you want to be in the nodeset. Separate multple coords with ';' (Either this parameter or \"nodes\" must be supplied)."}>>> = '0 0 0'
  []
[]

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = ReferenceResidualProblem
  extra_tag_vectors = 'ref'
  reference_vector = 'ref'
[]

[GlobalParams<<<{"href": "../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x'
  use_displaced_mesh = true
[]

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [temperature]
  []
  [disp_x]
  []
[]

[Physics<<<{"href": "../syntax/Physics/index.html"}>>>]
  [SolidMechanics<<<{"href": "../syntax/Physics/SolidMechanics/index.html"}>>>]
    [QuasiStatic<<<{"href": "../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
      [block1]
        strain<<<{"description": "Strain formulation"}>>> = FINITE
        add_variables<<<{"description": "Add the displacement variables"}>>> = true
        eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = eigenstrain
        generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = 'strain_xx'
        temperature<<<{"description": "The temperature"}>>> = temperature
      []
    []
  []
[]

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [temperature]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temperature
    function<<<{"description": "The initial condition function."}>>> = ${T0}
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [heat_conduction]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
  [kappa_fission]
    type = CoupledForce<<<{"description": "Implements a source term proportional to the value of a coupled variable. Weak form: $(\\psi_i, -\\sigma v)$.", "href": "../source/kernels/CoupledForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    v<<<{"description": "The coupled variable which provides the force"}>>> = kappa_fission
  []
  [dt]
    type = HeatConductionTimeDerivative<<<{"description": "Time derivative term $\\rho c_p \\frac{\\partial T}{\\partial t}$ of the thermal energy conservation equation.", "href": "../source/kernels/HeatConductionTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [kappa_fission]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = '${fparse P*eV_to_J/L0 }'
  []
  [temperature_praux]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [disp_praux]
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [temperature_projection]
    type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../source/auxkernels/ProjectionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = temperature_praux
    v<<<{"description": "Variable to take the value of."}>>> = temperature
  []
  [disp_projection]
    type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../source/auxkernels/ProjectionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = disp_praux
    v<<<{"description": "Variable to take the value of."}>>> = disp_x
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [thermal_parameters]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    temperature<<<{"description": "Coupled Temperature"}>>> = temperature
    thermal_conductivity_temperature_function<<<{"description": "Thermal conductivity as a function of temperature."}>>> = conductivity
    specific_heat<<<{"description": "The specific heat value"}>>> = ${c_p}
  []
  [density]
    type = ParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../source/materials/ParsedMaterial.html"}>>>
    property_name<<<{"description": "Name of the parsed material property"}>>> = density
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = ${rho}
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = ${E}
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = ${nu}
  []
  [stress]
    type = ComputeFiniteStrainElasticStress<<<{"description": "Compute stress using elasticity for finite strains", "href": "../source/materials/ComputeFiniteStrainElasticStress.html"}>>>
  []
  [thermal_expansion_strain]
    type = ComputeMeanThermalExpansionFunctionEigenstrain<<<{"description": "Computes eigenstrain due to thermal expansion using a function that describes the mean thermal expansion as a function of temperature", "href": "../source/materials/ComputeMeanThermalExpansionFunctionEigenstrain.html"}>>>
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = ${T0}
    thermal_expansion_function_reference_temperature<<<{"description": "Reference temperature for thermal_exansion_function (IMPORTANT: this is different in general from the stress_free_temperature)"}>>> = ${T0}
    thermal_expansion_function<<<{"description": "Function describing the mean thermal expansion as a function of temperature"}>>> = cte_func_mean
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = eigenstrain
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [sides_convective_BC]
    type = ConvectiveFluxFunction<<<{"description": "Determines boundary value by fluid heat transfer coefficient and far-field temperature", "href": "../source/bcs/ConvectiveFluxFunction.html"}>>>
    T_infinity<<<{"description": "Function describing far-field temperature"}>>> = ${T0}
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left right'
    coefficient<<<{"description": "Function describing heat transfer coefficient"}>>> = ${h}
  []
  [pin_x]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'pin'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [conductivity]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${fparse k0*eV_to_J} * t'
  []
  [cte_func_mean]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '1 / (2*sqrt(${T0}*t))'
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = 'PJFNK'
  verbose = true
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-7
  petsc_options = '-ksp_snes_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  snesmf_reuse_base = false
  line_search = 'none'
  nl_max_its = 100
  [TimeStepper<<<{"href": "../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 1e-1
    optimal_iterations = 10
    growth_factor = 1.25
    cutback_factor = 0.8
  []
[]

[Preconditioning<<<{"href": "../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [source_integral]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = kappa_fission
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'transfer initial timestep_end'
  []
  [strain_xx]
    type = ElementAverageValue<<<{"description": "Computes the volumetric average of a variable", "href": "../source/postprocessors/ElementAverageValue.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = strain_xx
  []
  [total_length]
    type = FunctionElementIntegral<<<{"description": "Integrates a function over elements", "href": "../source/postprocessors/FunctionElementIntegral.html"}>>>
    function<<<{"description": "Name of function to integrate"}>>> = 1
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_BEGIN LINEAR TIMESTEP_END'
  []
[]

It is important to note that incremental (hypoelastic) formulations of finite strain behavior inherently introduce inaccuracies arising from the time integration error. This limitation is inherent to this formulation and therefore exists in the MOOSE implementation and the implementations of other major commercial codes. The magnitude of these errors grows with increasing strain, so this typically becomes a practical issue only for very large strains, for which hyperelastic models are more appropriate.

The main sources of error in this coupled model are: (i) the statistical error from use of finite particles in OpenMC; (ii) the user-input cross-section data, which is defined with 1 K $var element = document.getElementById("moose-equation-4c87f7a9-a6c4-4ad1-a6ff-b45f15ff43ac");katex.render("\\Delta", element, {displayMode:false,throwOnError:false});$ T spacings used in all cases while taking the nearest temperature point cross section; and (iii) discretization error from a finite spatial mesh in both OpenMC (temperature feedback resolution) and MOOSE thermomechanics. Additional sources of error are present, but anticipated to be negligible due to the use of fine settings for these simulation parameters: (i) finite Picard iterations, and (ii) nonzero nonlinear tolerances in the MOOSE thermomechanics model.

Results

We see the heating result in Figure 2 compared against the analytic solution, as well as flux and temperature solutions in Figure 3 and Figure 4. Figures are only included for $var element = document.getElementById("moose-equation-7bd87510-d6bf-4b9f-b70c-b8fe45ed4665");katex.render("P=1.0 \\times 10^{22}", element, {displayMode:false,throwOnError:false});$ eV/s and on a fine mesh where the number of OpenMC and mesh elements is $var element = document.getElementById("moose-equation-58728443-59a0-46ca-878d-aaf7948c03d0");katex.render("N = 100", element, {displayMode:false,throwOnError:false});$ , since all results have a similar shape. Both the heating and flux tally everywhere agree with the analytic solution within 3 $var element = document.getElementById("moose-equation-ece58ab4-bc32-4478-b6f7-50f68b9f2dff");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ .

Figure 2: Heating result from OpenMC for N = 100 mesh elements compared against the analytic solution.

Figure 3: Flux result from OpenMC for N = 100 mesh elements compared against the analytic solution.

Figure 4: Temperature result from MOOSE for N = 100 mesh elements compared against the analytic solution.

Table 2 shows results convergence with mesh refinement.

Table 2: Results and corresponding errors with different mesh sizes (Eltawila et al., 2026).

Resolution	Heated length (cm)	Error (μm)	$var element = document.getElementById("moose-equation-7f9a0c3b-71e3-48c9-b98d-b3c81a580db1");katex.render("k_{\\text{eff}}", element, {displayMode:false,throwOnError:false});$	Error (pcm)
Analytic solution	106.47	–	0.29557	–
5	106.3298	-1402	0.29608 ± 0.00001	51 ± 1
10	106.4382	-318	0.29619 ± 0.00001	62 ± 1
20	106.4671	-29	0.29557 ± 0.00001	0 ± 1
50	106.4730	30	0.29536 ± 0.00001	-21 ± 1
100	106.4745	45	0.29549 ± 0.00001	-8 ± 1
200	106.4748	48	0.29556 ± 0.00001	-1 ± 1

In Figure 5 and Figure 6, convergence of flux and temperature relative L $var element = document.getElementById("moose-equation-7e3132a7-35eb-43ae-adc9-e94e3d479d19");katex.render("^2", element, {displayMode:false,throwOnError:false});$ norms, defined as

$var element = document.getElementById("moose-equation-5c986418-1d49-4290-825b-f4e09052fe62");katex.render(" \\varepsilon_\\phi=\\frac{\\|\\phi_a-\\phi_{\\mathrm{sln}}\\|_2}{\\|\\phi_a\\|_2}", element, {displayMode:true,throwOnError:false});$

and

$var element = document.getElementById("moose-equation-58f70c0a-6262-4192-a269-c767f14a32b7");katex.render(" \\varepsilon_T=\\frac{\\|T_a-T_{\\mathrm{sln}}\\|_2}{\\|T_a\\|_2}", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-11f5cf4c-5da8-490f-ac3c-d215814585b2");katex.render("\\phi_a", element, {displayMode:false,throwOnError:false});$ is the analytic flux solution, $var element = document.getElementById("moose-equation-7acf3ee2-af81-43e3-a2e5-1772217d462b");katex.render("\\phi_{sln.}", element, {displayMode:false,throwOnError:false});$ is the Cardinal flux solution, $var element = document.getElementById("moose-equation-23550907-4a9d-4010-b7b7-89d7845862a0");katex.render("T_a", element, {displayMode:false,throwOnError:false});$ is the analytic temperature solution, and $var element = document.getElementById("moose-equation-96adcac8-614c-4afe-ad5c-8986030cc7eb");katex.render("T_{sln.}", element, {displayMode:false,throwOnError:false});$ is the Cardinal temperature solution.

Convergence of flux relative L$^2$ norm.

Figure 5: Convergence of flux relative L $var element = document.getElementById("moose-equation-4bf5a6a3-b9bf-4271-87b9-91cc092accd6");katex.render("^2", element, {displayMode:false,throwOnError:false});$ norm.

Convergence of temperature relative L$^2$ norm.

Figure 6: Convergence of temperature relative L $var element = document.getElementById("moose-equation-d3012f01-d5ab-4dfb-bc0a-928947ad8a4a");katex.render("^2", element, {displayMode:false,throwOnError:false});$ norm.

A first-order convergence of the temperature is achieved with mesh refinement, while the flux starts to converge following a first-order trend then stagnates where other error components dominate as discussed earlier. The stochastic noise error influence is highlighted by keeping the number of particles constant for each line, where we see further convergence when there is less stochastic noise.

We hypothesize that the observed first-order convergence arises from the use of cell tallies in OpenMC, where heating results are tallied as constant monomials. This may be limiting the convergence rate of the coupled solve.

Cardinal Random Ray Solutions

The Griesheimer and Kooreman problem is solved in Cardinal using The Random Ray Method (TRRM) in a similar manner to the multi-group Monte Carlo approach presented in Eltawila et al. (2026). The neutronics domain is discretized into 200 regions, which are further subdivided into 31 tetrahedral elements per subregion (for a total of 6200 flat source regions). The conversion to tetrahedral elements is necessary to use on-the-fly geometry skinning in Cardinal for material deformation. Unfortunately, this conversion introduces non-physical source region gradients in the y and z directions which bias the TRRM solution. It is important to note that this bias only exists for solutions using TRRM; multi-group Monte Carlo solutions in Cardinal do not experience any bias as the fission source is continuous. 100 inactive batches are used to converge the scattering and fission sources, and 100 active batches are used to accumulate tally results. An inactive ray length of $var element = document.getElementById("moose-equation-3445fa9c-4df8-4a65-b8d4-e562afb70ee6");katex.render("150", element, {displayMode:false,throwOnError:false});$ cm is used to ensure rays reach a vacuum boundary condition (to remove initialization bias), which is followed by an active length of $var element = document.getElementById("moose-equation-efb6e4a9-5afd-4c53-b14e-0535a5648aa2");katex.render("7500", element, {displayMode:false,throwOnError:false});$ cm. Temperature feedback is implemented by interpolating a cross section grid with a spacing of $var element = document.getElementById("moose-equation-1d928805-b48e-468a-8b87-89533ae455d7");katex.render("1", element, {displayMode:false,throwOnError:false});$ K.

The neutronics model is coupled with a heat conduction and thermal expansion solver from the MOOSE Heat Transfer module in an identical manner to the approach presented in Eltawila et al. (2026). Temperatures and displacements are solved on a 1D mesh with linear Lagrange basis functions. A quasi-static finite strain formulation is used, and a time derivative is included in the temperature model to act as physics-based relaxation. The solution is iterated between neutronics and thermomechanics with Picard iterations; it is considered to be converged when the relative difference between neutronics field variables between Picard iterations is less than $var element = document.getElementById("moose-equation-fbc16e45-f6d7-49d4-8f39-e6b1a42da0ed");katex.render("10^{-8}", element, {displayMode:false,throwOnError:false});$ . The resulting flux and temperature distributions can be found in Figure 7.

Scalar flux and temperature distributions predicted with [!ac](TRRM) in Cardinal.

Figure 7: Scalar flux and temperature distributions predicted with TRRM in Cardinal.

The impact of the y-z induced source region gradients can be seen in the flux solution, resulting in reasonable (but not perfect) agreement with Eq. (14). The thermomechanics problem uses a true 1D mesh, resulting in better agreement between the numerical temperature predictions and Eq. (15). The TRRM solution predicts $var element = document.getElementById("moose-equation-49f1f431-9916-4e62-80a7-23b0ec447c0c");katex.render("k_{eff} = 0.29520", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-be61fa5f-404f-473c-bac9-9f6ab2e614e9");katex.render("\\Delta k_{eff} = \u221237", element, {displayMode:false,throwOnError:false});$ pcm) and an expanded slab length $var element = document.getElementById("moose-equation-b164a29f-39d6-4e82-8ccf-5e6eff09632a");katex.render("L = 106.47", element, {displayMode:false,throwOnError:false});$ cm. It is expected that the agreement in both the flux and eigenvalue solutions will improve when hexahedral element support is added to the MoabSkinner in issue #1352,as this will eliminate the y-z source region gradients in the TRRM solution.

The input files for the TRRM solutions to the Griesheimer and Kooreman problem can be found in /test/tests/neutronics/mg/rr_thermal_expansion, where they are used to regression test on-the-fly skinning with TRRM. Please note that these input files have been modified to run a 20 cell mesh, with fewer Picard iterations and rays per batch to reduce the computational requirements of the test.

# Cardinal input file for slab benchmark OpenMC model to couple with thermomechanics

P = 1.0e22 # eV/s
eV_to_J = 1.602e-19 # J/eV

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [file]
    type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../source/meshgenerators/FileMeshGenerator.html"}>>>
    file<<<{"description": "The filename to read."}>>> = slab20v.e
  []
[]

[GlobalParams<<<{"href": "../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  use_displaced_mesh = true
[]

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [disp_x]
  []
  [disp_y]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.0
  []
  [disp_z]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.0
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [bins]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [bins]
    type = SkinnedBins<<<{"description": "Bins created by a skinner", "href": "../source/auxkernels/SkinnedBins.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = bins
    skinner<<<{"description": "MOAB mesh skinner"}>>> = moab
  []
[]

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [temperature]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temp
    function<<<{"description": "The initial condition function."}>>> = 293
  []
[]

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = OpenMCCellAverageProblem
  verbose = true
  source_rate_normalization = 'kappa_fission'

  inactive_batches = 100
  batches = 200
  particles = 10

  power = ${fparse P*eV_to_J}
  cell_level = 0
  temperature_blocks = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
  volume_calculation = vol
  skinner = moab
  [Tallies<<<{"href": "../syntax/Problem/Tallies/index.html"}>>>]
    [kappa_fission]
      type = CellTally<<<{"description": "A class which implements distributed cell tallies.", "href": "../source/tallies/CellTally.html"}>>>
      name<<<{"description": "Auxiliary variable name(s) to use for OpenMC tallies. If not specified, defaults to the names of the scores"}>>> = kappa_fission
      output<<<{"description": "UNRELAXED field(s) to output from OpenMC for each tally score. unrelaxed_tally_std_dev will write the standard deviation of each tally into auxiliary variables named *_std_dev. unrelaxed_tally_rel_error will write the relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) of each tally into auxiliary variables named *_rel_error. unrelaxed_tally will write the raw unrelaxed tally into auxiliary variables named *_raw (replace * with 'name')."}>>> = unrelaxed_tally_std_dev
      normalize_by_global_tally<<<{"description": "Whether to normalize local tallies by a global tally (true) or else by the sum of the local tally (false)"}>>> = false
    []
    [flux]
      type = CellTally<<<{"description": "A class which implements distributed cell tallies.", "href": "../source/tallies/CellTally.html"}>>>
      score<<<{"description": "Score(s) to use in the OpenMC tallies. If not specified, defaults to 'kappa_fission'"}>>> = flux
      output<<<{"description": "UNRELAXED field(s) to output from OpenMC for each tally score. unrelaxed_tally_std_dev will write the standard deviation of each tally into auxiliary variables named *_std_dev. unrelaxed_tally_rel_error will write the relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) of each tally into auxiliary variables named *_rel_error. unrelaxed_tally will write the raw unrelaxed tally into auxiliary variables named *_raw (replace * with 'name')."}>>> = unrelaxed_tally_std_dev
      normalize_by_global_tally<<<{"description": "Whether to normalize local tallies by a global tally (true) or else by the sum of the local tally (false)"}>>> = false
    []
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  start_time = 0.0
  num_steps = 5
  steady_state_detection = true
[]

[MultiApps<<<{"href": "../syntax/MultiApps/index.html"}>>>]
  [solid]
    type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../source/multiapps/TransientMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'solid.i'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Transfers<<<{"href": "../syntax/Transfers/index.html"}>>>]
  [kappa_fission_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer<<<{"description": "Transfers field data at the MultiApp position using the finite element shape functions from the origin application.", "href": "../source/transfers/MultiAppGeneralFieldShapeEvaluationTransfer.html"}>>>
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = solid
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = kappa_fission
    source_variable<<<{"description": "The variable to transfer from."}>>> = kappa_fission
    from_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the from-app to evaluate an adjusting factor."}>>> = heat_source
    to_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the to-app to evaluate an adjusting factor."}>>> = source_integral
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    displaced_source_mesh<<<{"description": "Whether or not to use the displaced mesh for the source mesh."}>>> = true
    displaced_target_mesh<<<{"description": "Whether or not to use the displaced mesh for the target mesh."}>>> = true
  []
  [temp_from_solid]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = solid
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = temp
    source_variable<<<{"description": "The variable to transfer from."}>>> = temperature_praux
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    displaced_source_mesh<<<{"description": "Whether or not to use the displaced mesh for the source mesh."}>>> = true
    displaced_target_mesh<<<{"description": "Whether or not to use the displaced mesh for the target mesh."}>>> = true
  []
  [dispx]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = disp_praux
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = disp_x
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = solid
    displaced_source_mesh<<<{"description": "Whether or not to use the displaced mesh for the source mesh."}>>> = true
    displaced_target_mesh<<<{"description": "Whether or not to use the displaced mesh for the target mesh."}>>> = true
  []
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
  execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_END'
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [heat_source]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = kappa_fission
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'transfer initial timestep_end'
  []
  [k]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
    value_type<<<{"description": "Type of eigenvalue global tally to report"}>>> = 'tracklength'
  []
  [k_std_dev]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
    value_type<<<{"description": "Type of eigenvalue global tally to report"}>>> = 'tracklength'
    output<<<{"description": "The value to output. Options are $k_{eff}$ (mean), the standard deviation of $k_{eff}$ (std_dev), or the relative error of $k_{eff}$ (rel_err)."}>>> = 'std_dev'
  []
[]

[UserObjects<<<{"href": "../syntax/UserObjects/index.html"}>>>]
  [vol]
    type = OpenMCVolumeCalculation<<<{"description": "Stochastic volume calculation of the OpenMC cells which map to MOOSE", "href": "../source/userobjects/OpenMCVolumeCalculation.html"}>>>
    n_samples<<<{"description": "Number of samples to use for the stochastic volume calculation"}>>> = 200000
  []
  [moab]
    type = MoabSkinner<<<{"description": "Re-generate the OpenMC geometry on-the-fly according to changes in the mesh geometry and/or contours in temperature and density", "href": "../source/userobjects/MoabSkinner.html"}>>>
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin'
    verbose<<<{"description": "Whether to print diagnostic information"}>>> = true
    temperature<<<{"description": "Temperature variable by which to bin elements"}>>> = temp
    n_temperature_bins<<<{"description": "Number of temperature bins"}>>> = 1.0
    temperature_min<<<{"description": "Lower bound of temperature bins"}>>> = 0.0
    temperature_max<<<{"description": "Upper bound of temperature bins"}>>> = 1000
    build_graveyard<<<{"description": "Whether to build a graveyard around the geometry"}>>> = true
    graveyard_scale_inner<<<{"description": "Multiplier on mesh bounding box to form inner graveyard surface"}>>> = 1.5
    graveyard_scale_outer<<<{"description": "Multiplier on mesh bounding box to form outer graveyard surface"}>>> = 1.6
  []
[]

# MOOSE Thermomechanics model input file for the slab benchmark problem

T0 = 293 # K
L0 = 100 # cm
L = 106.47 # cm
P = 1.0e22 # eV/s
q = 1e8 # eV
k0 = 1.25e19 # eV/(s-cm-K^2) k(T) = k0 T(x)
phi0 = 2.5e14 # 1/s-cm^2 (flux at the origin)
eV_to_J = 1.602e-19 # J/eV
lam = ${fparse 0.5*(1+sqrt(1+(16*q*q*phi0*phi0)/(P*P)))} # eigenvalue solution
h = ${fparse 1/(sqrt(L*(lam-1)/(k0*P)) - (2*T0)/(P))*eV_to_J } # W/cm^2-K

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [generated]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 1
    nx<<<{"description": "Number of elements in the X direction"}>>> = 20
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = -50
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 50
  []
  [pin]
    type = ExtraNodesetGenerator<<<{"description": "Creates a new node set and a new boundary made with the nodes the user provides.", "href": "../source/meshgenerators/ExtraNodesetGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = generated
    new_boundary<<<{"description": "The names of the boundaries to create"}>>> = pin
    coord<<<{"description": "The nodes with coordinates you want to be in the nodeset. Separate multple coords with ';' (Either this parameter or \"nodes\" must be supplied)."}>>> = '0 0 0'
  []
[]

[GlobalParams<<<{"href": "../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x'
  use_displaced_mesh = true
[]

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [temperature]
  []
  [disp_x]
  []
[]

[Physics<<<{"href": "../syntax/Physics/index.html"}>>>]
  [SolidMechanics<<<{"href": "../syntax/Physics/SolidMechanics/index.html"}>>>]
    [QuasiStatic<<<{"href": "../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
      [block1]
        strain<<<{"description": "Strain formulation"}>>> = FINITE
        add_variables<<<{"description": "Add the displacement variables"}>>> = true
        eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = eigenstrain
        generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = 'strain_xx'
        temperature<<<{"description": "The temperature"}>>> = temperature
      []
    []
  []
[]

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [temperature]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = temperature
    function<<<{"description": "The initial condition function."}>>> = 293
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [heat_conduction]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
  [kappa_fission]
    type = CoupledForce<<<{"description": "Implements a source term proportional to the value of a coupled variable. Weak form: $(\\psi_i, -\\sigma v)$.", "href": "../source/kernels/CoupledForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    v<<<{"description": "The coupled variable which provides the force"}>>> = kappa_fission
  []
  [dt]
    type = HeatConductionTimeDerivative<<<{"description": "Time derivative term $\\rho c_p \\frac{\\partial T}{\\partial t}$ of the thermal energy conservation equation.", "href": "../source/kernels/HeatConductionTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [kappa_fission]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = '${fparse P*eV_to_J/L0 }'
  []
  [temperature_praux]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [disp_praux]
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
   [temperature_projection]
     type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../source/auxkernels/ProjectionAux.html"}>>>
     variable<<<{"description": "The name of the variable that this object applies to"}>>> = temperature_praux
     v<<<{"description": "Variable to take the value of."}>>> = temperature
   []
   [disp_projection]
     type = ProjectionAux<<<{"description": "Returns the specified variable as an auxiliary variable with a projection of the source variable. If they are the same type, this amounts to a simple copy.", "href": "../source/auxkernels/ProjectionAux.html"}>>>
     variable<<<{"description": "The name of the variable that this object applies to"}>>> = disp_praux
     v<<<{"description": "Variable to take the value of."}>>> = disp_x
   []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [thermal_parameters]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    temp = temperature
    thermal_conductivity_temperature_function<<<{"description": "Thermal conductivity as a function of temperature."}>>> = conductivity
    specific_heat<<<{"description": "The specific heat value"}>>> = 0.3
  []
  [density]
    type = Density<<<{"description": "Creates density material property. This class is deprecated, and its functionalityis replaced by StrainAdjustedDensity for cases when the density should be adjustedto account for material deformation. If it is not desired to adjust the density fordeformation, a variety of general-purpose Materials, such as GenericConstantMaterialor ParsedMaterial can be used to define the density. A StrainAdjustedDensity can also be used with '0 0 0' for the displacements."}>>>
    density<<<{"description": "Density"}>>> = 20.0
    use_displaced_mesh<<<{"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."}>>> = false
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 1e6
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.0
  []
  [stress]
    type = ComputeFiniteStrainElasticStress<<<{"description": "Compute stress using elasticity for finite strains", "href": "../source/materials/ComputeFiniteStrainElasticStress.html"}>>>
  []
  [thermal_expansion_strain]
    type = ComputeMeanThermalExpansionFunctionEigenstrain<<<{"description": "Computes eigenstrain due to thermal expansion using a function that describes the mean thermal expansion as a function of temperature", "href": "../source/materials/ComputeMeanThermalExpansionFunctionEigenstrain.html"}>>>
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = ${T0}
    thermal_expansion_function_reference_temperature<<<{"description": "Reference temperature for thermal_exansion_function (IMPORTANT: this is different in general from the stress_free_temperature)"}>>> = ${T0}
    thermal_expansion_function<<<{"description": "Function describing the mean thermal expansion as a function of temperature"}>>> = cte_func_mean
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = eigenstrain
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [sides_convective_BC]
    type = ConvectiveFluxFunction<<<{"description": "Determines boundary value by fluid heat transfer coefficient and far-field temperature", "href": "../source/bcs/ConvectiveFluxFunction.html"}>>>
    T_infinity<<<{"description": "Function describing far-field temperature"}>>> = ${T0}
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left right'
    coefficient<<<{"description": "Function describing heat transfer coefficient"}>>> = ${h}
  []
  [pin_x]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'pin'
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [conductivity]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${fparse k0*eV_to_J} * t'
  []
  [cte_func_mean]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the symbol_values vector."}>>> = 'T0'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = '293'
    expression<<<{"description": "The user defined function."}>>> = '1 / (2*sqrt(T0*t))'
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = 'PJFNK'
  verbose = true
  nl_abs_tol = 5e-6
  nl_rel_tol = 5e-6
  petsc_options = '-ksp_snes_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  snesmf_reuse_base = false
  line_search = 'none'
  nl_max_its = 100
  [TimeStepper<<<{"href": "../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 1e1
    optimal_iterations = 10
    growth_factor = 1.25
    cutback_factor = 0.8
  []
[]

[Preconditioning<<<{"href": "../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [source_integral]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = kappa_fission
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'transfer initial timestep_end'
  []
  [strain_xx]
    type = ElementAverageValue<<<{"description": "Computes the volumetric average of a variable", "href": "../source/postprocessors/ElementAverageValue.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = strain_xx
  []
  [total_length]
    type = FunctionElementIntegral<<<{"description": "Integrates a function over elements", "href": "../source/postprocessors/FunctionElementIntegral.html"}>>>
    function<<<{"description": "Name of function to integrate"}>>> = 1
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_BEGIN LINEAR TIMESTEP_END'
  []
[]

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = ReferenceResidualProblem
  extra_tag_vectors = 'ref'
  reference_vector = 'ref'
[]

References

Mahmoud Eltawila, Pierre-Clément A. Simon, Guillaume L. Giudicelli, Helen Brooks, Nicholas Wozniak, and April J. Novak. Thermomechanics coupling to monte carlo particle transport on unstructured mesh geometries using cardinal. Annals of Nuclear Energy, 229:112068, 2026. doi:https://doi.org/10.1016/j.anucene.2025.112068.

BibTeX

@article{eltawila_thermo,
    author = "Eltawila, Mahmoud and Simon, Pierre-Clément A. and Giudicelli, Guillaume L. and Brooks, Helen and Wozniak, Nicholas and Novak, April J.",
    title = "Thermomechanics coupling to Monte Carlo particle transport on unstructured mesh geometries using Cardinal",
    journal = "Annals of Nuclear Energy",
    volume = "229",
    pages = "112068",
    year = "2026",
    doi = "https://doi.org/10.1016/j.anucene.2025.112068"
}

RIS

TY  - JOUR
AU  - Eltawila, Mahmoud
AU  - Simon, Pierre-Clément A.
AU  - Giudicelli, Guillaume L.
AU  - Brooks, Helen
AU  - Wozniak, Nicholas
AU  - Novak, April J.
TI  - Thermomechanics coupling to Monte Carlo particle transport on unstructured mesh geometries using Cardinal
JO  - Annals of Nuclear Energy
PY  - 2026
VL  - 229
DO  - https://doi.org/10.1016/j.anucene.2025.112068
SP  - 112068
ER  -

Plain Text

Mahmoud Eltawila, Pierre-Clément A. Simon, Guillaume L. Giudicelli, Helen Brooks, Nicholas Wozniak, and April J. Novak. Thermomechanics coupling to monte carlo particle transport on unstructured mesh geometries using cardinal. Annals of Nuclear Energy, 229:112068, 2026. doi:https://doi.org/10.1016/j.anucene.2025.112068.

David P. Griesheimer and Gabriel Kooreman. Analytical Benchmark Solution for 1-D Neutron Transport Coupled with Thermal Conduction and Material Expansion. In Proceedings of PHYSOR. 2022.

BibTeX

@InProceedings{s2_thermal_expansion,
    author = "Griesheimer, David P. and Kooreman, Gabriel",
    title = "{Analytical Benchmark Solution for 1-D Neutron Transport Coupled with Thermal Conduction and Material Expansion}",
    booktitle = "{Proceedings of PHYSOR}",
    year = "2022"
}

RIS

TY  - CPAPER
AU  - Griesheimer, David P.
AU  - Kooreman, Gabriel
TI  - Analytical Benchmark Solution for 1-D Neutron Transport Coupled with Thermal Conduction and Material Expansion
T2  - Proceedings of PHYSOR
PY  - 2022
ER  -

Plain Text

David P. Griesheimer and Gabriel Kooreman. Analytical Benchmark Solution for 1-D Neutron Transport Coupled with Thermal Conduction and Material Expansion. In Proceedings of PHYSOR. 2022.

(doc/content/vv/slab_benchmark_model/model.py)

import openmc
import openmc.mgxs as mgxs
import numpy as np
import h5py

T0 = 293 # K
Tmax = 390 # K
Tmin = 290 # K
L0 = 100 # cm
L = 106.47 # cm
rho = 1.2 # g/cc
N_A = 6.022e23 # Avagadro's number
A = 180 # mass number for slab material
num_dens = rho*N_A/A
P = 1.0e22 # eV/s
q = 1e8 # eV
k0= 1.25e19 # eV/(s-cm-K^2)
phi0 = 2.5e14 # 1/s-cm^2
s = 0.45 # Sigma_s/Sigma_t
f = 1.5 # nu Sigma_f/Sigma_t
nu = f/(1-s)
lam = 0.5*(1+np.sqrt(1+(16*q*q*phi0*phi0)/(P*P)))
Sig_t0 = np.sqrt(P/((lam-1)*k0*L))/(T0) # 1/cm
sig_t0 = Sig_t0/num_dens # cm^2
N = 20
infdim = 0.5

model = openmc.Model()

dagmc_univ = openmc.DAGMCUniverse(filename="slab20v.h5m")
model.geometry = openmc.Geometry(root=dagmc_univ)
cells = dagmc_univ.get_all_cells().values()

slab = openmc.Material(1, "slab")
slab.set_density('macro',1.)
slab.add_macroscopic('slab_xs')

model.materials = openmc.Materials([slab])
model.materials.cross_sections = 'one_group.h5'

groups = mgxs.EnergyGroups(group_edges=[0.0, 20.0e6])
NT = 100
temps = np.linspace(Tmin,Tmax,NT)
xsdata = openmc.XSdata('slab_xs', energy_groups=groups, temperatures=temps, num_delayed_groups=0)
xsdata.order = 0

for T in temps:
    Sig_t = (Sig_t0 * T0) / T
    Sig_s = s*Sig_t
    nu_Sig_f = f*Sig_t
    Sig_f = nu_Sig_f / nu
    Sig_a = Sig_f
    xsdata.set_total(np.array([Sig_t]),temperature=T)
    xsdata.set_scatter_matrix(np.array([[[Sig_s]]]),temperature=T)
    xsdata.set_absorption(np.array([Sig_a]),temperature=T)
    xsdata.set_fission(np.array([Sig_f]),temperature=T)
    xsdata.set_nu_fission(np.array([nu_Sig_f]),temperature=T)
    xsdata.set_kappa_fission(np.array([q*Sig_t]),temperature=T)

one_g_XS_file = openmc.MGXSLibrary(groups)
one_g_XS_file.add_xsdata(xsdata)
one_g_XS_file.export_to_hdf5('one_group.h5')

model.settings = openmc.Settings()
batches = 150
inactive = 50
particles = 50000
model.settings.energy_mode = 'multi-group'
model.settings.batches = batches
model.settings.inactive = inactive
model.settings.particles = particles
model.settings.output = {'tallies': True,'summary':False}

bounds = [-60, -infdim, -infdim, 60, infdim, infdim]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
mu = openmc.stats.Discrete([0.0],[1.0])
phi = openmc.stats.Discrete([0,np.pi],[0.5,0.5])
bidirectional_x = openmc.stats.PolarAzimuthal(mu=mu,phi=phi)
model.settings.source = openmc.IndependentSource(space=uniform_dist,angle=bidirectional_x)
model.settings.temperature = {'default': (Tmin+Tmax)/2,
                        'method': 'nearest',
                        'tolerance': 50,
                        'range': (Tmin, Tmax)}

model.export_to_model_xml()

(doc/content/vv/slab_benchmark_model/openmc.i)

P = ${units 1.0e22 eV/s}
eV_to_J = ${units 1.602176634e-19 J/eV}
initial_temperature = ${units 293 K}

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

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  use_displaced_mesh = true
[]

[Variables]
  [disp_x]
  []
  [disp_y]
    initial_condition = 0.0
  []
  [disp_z]
    initial_condition = 0.0
  []
[]

[AuxVariables]
  [bins]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [bins]
    type = SkinnedBins
    variable = bins
    skinner = moab
  []
[]

[ICs]
  [temperature]
    type = FunctionIC
    variable = temperature
    function = ${initial_temperature}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  source_rate_normalization = 'kappa_fission'
  inactive_batches = 150
  batches = 300
  particles = 10000000
  power = '${fparse P*eV_to_J}'
  cell_level = 0
  temperature_variables = 'temperature'
  temperature_blocks = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
  volume_calculation = vol
  skinner = moab
  [Tallies]
    [kappa_fission]
      type = CellTally
      estimator = collision
      block = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
      name = kappa_fission
      output = unrelaxed_tally_std_dev
    []
    [flux]
      type = CellTally
      estimator = collision
      block = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
      score = flux
      output = unrelaxed_tally_std_dev
    []
  []
[]

[Executioner]
  type = Transient
  start_time = 0.0
  end_time = 1e5
  dt = 50
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    input_files = 'solid.i'
    execute_on = timestep_end
  []
[]

[Transfers]
  [kappa_fission_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = solid
    variable = kappa_fission
    source_variable = kappa_fission
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = source_integral
    execute_on = timestep_end
    displaced_source_mesh = true
    displaced_target_mesh = true
  []
  [temperature_from_solid]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = solid
    variable = temperature
    source_variable = temperature_praux
    execute_on = timestep_end
    displaced_source_mesh = true
    displaced_target_mesh = true
  []
  [dispx]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = disp_praux
    variable = disp_x
    from_multi_app = solid
    displaced_source_mesh = true
    displaced_target_mesh = true
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = kappa_fission
    execute_on = 'transfer initial timestep_end'
  []
  [k]
    type = KEigenvalue
  []
  [k_std_dev]
    type = KEigenvalue
    output = 'std_dev'
  []
[]

[VectorPostprocessors]
  [temperature]
    type = ElementValueSampler
    variable = 'temperature'
    sort_by = x
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [vol]
    type = OpenMCVolumeCalculation
    n_samples = 200000
  []
  [moab]
    type = MoabSkinner
    execute_on = 'timestep_begin'
    verbose = true
    temperature = temperature
    n_temperature_bins = 20.0
    temperature_min = 0
    temperature_max = 1000
    build_graveyard = true
    graveyard_scale_inner = 1.5
    graveyard_scale_outer = 1.6
  []
[]

(doc/content/vv/slab_benchmark_model/solid.i)

T0 = ${units 293 K}
L0 = ${units 100 cm}
L = ${units 106.47 cm}
P = ${units 1.0e22 eV/s}
q = ${units 1e8 eV}
k0 = ${units 1.25e19 eV/s/cm/K^2} # k(T) = k0 T(x)
phi0 = ${units 2.5e14 1/s/cm^2} # (flux at the origin)
eV_to_J = ${units 1.602176634e-19 J/eV}
lam = '${fparse 0.5*(1+sqrt(1+(16*q*q*phi0*phi0)/(P*P)))}' # eigenvalue solution
h = '${fparse 1/(sqrt(L*(lam-1)/(k0*P)) - (2*T0)/(P))*eV_to_J }' # W/cm^2-K

E = ${units 1e6 N/m^2}
nu = 0.0
rho = ${units 20 g/cm^3}
c_p = ${units 0.3 J/g/K}

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 500
    xmin = -50
    xmax = 50
  []
  [pin]
    type = ExtraNodesetGenerator
    input = generated
    new_boundary = pin
    coord = '0 0 0'
  []
[]

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

[GlobalParams]
  displacements = 'disp_x'
  use_displaced_mesh = true
[]

[Variables]
  [temperature]
  []
  [disp_x]
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [block1]
        strain = FINITE
        add_variables = true
        eigenstrain_names = eigenstrain
        generate_output = 'strain_xx'
        temperature = temperature
      []
    []
  []
[]

[ICs]
  [temperature]
    type = FunctionIC
    variable = temperature
    function = ${T0}
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = temperature
  []
  [kappa_fission]
    type = CoupledForce
    variable = temperature
    v = kappa_fission
  []
  [dt]
    type = HeatConductionTimeDerivative
    variable = temperature
  []
[]

[AuxVariables]
  [kappa_fission]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = '${fparse P*eV_to_J/L0 }'
  []
  [temperature_praux]
    family = MONOMIAL
    order = CONSTANT
  []
  [disp_praux]
  []
[]

[AuxKernels]
  [temperature_projection]
    type = ProjectionAux
    variable = temperature_praux
    v = temperature
  []
  [disp_projection]
    type = ProjectionAux
    variable = disp_praux
    v = disp_x
  []
[]

[Materials]
  [thermal_parameters]
    type = HeatConductionMaterial
    temperature = temperature
    thermal_conductivity_temperature_function = conductivity
    specific_heat = ${c_p}
  []
  [density]
    type = ParsedMaterial
    property_name = density
    expression = ${rho}
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = ${E}
    poissons_ratio = ${nu}
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [thermal_expansion_strain]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    stress_free_temperature = ${T0}
    thermal_expansion_function_reference_temperature = ${T0}
    thermal_expansion_function = cte_func_mean
    temperature = temperature
    eigenstrain_name = eigenstrain
  []
[]

[BCs]
  [sides_convective_BC]
    type = ConvectiveFluxFunction
    T_infinity = ${T0}
    variable = temperature
    boundary = 'left right'
    coefficient = ${h}
  []
  [pin_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'pin'
    value = 0.0
  []
[]

[Functions]
  [conductivity]
    type = ParsedFunction
    expression = '${fparse k0*eV_to_J} * t'
  []
  [cte_func_mean]
    type = ParsedFunction
    expression = '1 / (2*sqrt(${T0}*t))'
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  verbose = true
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-7
  petsc_options = '-ksp_snes_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  snesmf_reuse_base = false
  line_search = 'none'
  nl_max_its = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-1
    optimal_iterations = 10
    growth_factor = 1.25
    cutback_factor = 0.8
  []
[]

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

[Outputs]
  exodus = true
  csv = true
[]

[Postprocessors]
  [source_integral]
    type = ElementIntegralVariablePostprocessor
    variable = kappa_fission
    execute_on = 'transfer initial timestep_end'
  []
  [strain_xx]
    type = ElementAverageValue
    variable = strain_xx
  []
  [total_length]
    type = FunctionElementIntegral
    function = 1
    execute_on = 'TIMESTEP_BEGIN LINEAR TIMESTEP_END'
  []
[]

(test/tests/neutronics/mg/rr_thermal_expansion/openmc.i)

# Cardinal input file for slab benchmark OpenMC model to couple with thermomechanics

P = 1.0e22 # eV/s
eV_to_J = 1.602e-19 # J/eV

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

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  use_displaced_mesh = true
[]

[Variables]
  [disp_x]
  []
  [disp_y]
    initial_condition = 0.0
  []
  [disp_z]
    initial_condition = 0.0
  []
[]

[AuxVariables]
  [bins]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [bins]
    type = SkinnedBins
    variable = bins
    skinner = moab
  []
[]

[ICs]
  [temperature]
    type = FunctionIC
    variable = temp
    function = 293
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  source_rate_normalization = 'kappa_fission'

  inactive_batches = 100
  batches = 200
  particles = 10

  power = ${fparse P*eV_to_J}
  cell_level = 0
  temperature_blocks = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20'
  volume_calculation = vol
  skinner = moab
  [Tallies]
    [kappa_fission]
      type = CellTally
      name = kappa_fission
      output = unrelaxed_tally_std_dev
      normalize_by_global_tally = false
    []
    [flux]
      type = CellTally
      score = flux
      output = unrelaxed_tally_std_dev
      normalize_by_global_tally = false
    []
  []
[]

[Executioner]
  type = Transient
  start_time = 0.0
  num_steps = 5
  steady_state_detection = true
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    input_files = 'solid.i'
    execute_on = timestep_end
  []
[]

[Transfers]
  [kappa_fission_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = solid
    variable = kappa_fission
    source_variable = kappa_fission
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = source_integral
    execute_on = timestep_end
    displaced_source_mesh = true
    displaced_target_mesh = true
  []
  [temp_from_solid]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = solid
    variable = temp
    source_variable = temperature_praux
    execute_on = timestep_end
    displaced_source_mesh = true
    displaced_target_mesh = true
  []
  [dispx]
    type = MultiAppGeneralFieldNearestLocationTransfer
    source_variable = disp_praux
    variable = disp_x
    from_multi_app = solid
    displaced_source_mesh = true
    displaced_target_mesh = true
  []
[]

[Outputs]
  csv = true
  execute_on = 'TIMESTEP_END'
[]

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = kappa_fission
    execute_on = 'transfer initial timestep_end'
  []
  [k]
    type = KEigenvalue
    value_type = 'tracklength'
  []
  [k_std_dev]
    type = KEigenvalue
    value_type = 'tracklength'
    output = 'std_dev'
  []
[]

[UserObjects]
  [vol]
    type = OpenMCVolumeCalculation
    n_samples = 200000
  []
  [moab]
    type = MoabSkinner
    execute_on = 'timestep_begin'
    verbose = true
    temperature = temp
    n_temperature_bins = 1.0
    temperature_min = 0.0
    temperature_max = 1000
    build_graveyard = true
    graveyard_scale_inner = 1.5
    graveyard_scale_outer = 1.6
  []
[]

(test/tests/neutronics/mg/rr_thermal_expansion/solid.i)

# MOOSE Thermomechanics model input file for the slab benchmark problem

T0 = 293 # K
L0 = 100 # cm
L = 106.47 # cm
P = 1.0e22 # eV/s
q = 1e8 # eV
k0 = 1.25e19 # eV/(s-cm-K^2) k(T) = k0 T(x)
phi0 = 2.5e14 # 1/s-cm^2 (flux at the origin)
eV_to_J = 1.602e-19 # J/eV
lam = ${fparse 0.5*(1+sqrt(1+(16*q*q*phi0*phi0)/(P*P)))} # eigenvalue solution
h = ${fparse 1/(sqrt(L*(lam-1)/(k0*P)) - (2*T0)/(P))*eV_to_J } # W/cm^2-K

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 20
    xmin = -50
    xmax = 50
  []
  [pin]
    type = ExtraNodesetGenerator
    input = generated
    new_boundary = pin
    coord = '0 0 0'
  []
[]

[GlobalParams]
  displacements = 'disp_x'
  use_displaced_mesh = true
[]

[Variables]
  [temperature]
  []
  [disp_x]
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [block1]
        strain = FINITE
        add_variables = true
        eigenstrain_names = eigenstrain
        generate_output = 'strain_xx'
        temperature = temperature
      []
    []
  []
[]

[ICs]
  [temperature]
    type = FunctionIC
    variable = temperature
    function = 293
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = temperature
  []
  [kappa_fission]
    type = CoupledForce
    variable = temperature
    v = kappa_fission
  []
  [dt]
    type = HeatConductionTimeDerivative
    variable = temperature
  []
[]

[AuxVariables]
  [kappa_fission]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = '${fparse P*eV_to_J/L0 }'
  []
  [temperature_praux]
    family = MONOMIAL
    order = CONSTANT
  []
  [disp_praux]
  []
[]

[AuxKernels]
   [temperature_projection]
     type = ProjectionAux
     variable = temperature_praux
     v = temperature
   []
   [disp_projection]
     type = ProjectionAux
     variable = disp_praux
     v = disp_x
   []
[]

[Materials]
  [thermal_parameters]
    type = HeatConductionMaterial
    temp = temperature
    thermal_conductivity_temperature_function = conductivity
    specific_heat = 0.3
  []
  [density]
    type = Density
    density = 20.0
    use_displaced_mesh = false
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.0
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [thermal_expansion_strain]
    type = ComputeMeanThermalExpansionFunctionEigenstrain
    stress_free_temperature = ${T0}
    thermal_expansion_function_reference_temperature = ${T0}
    thermal_expansion_function = cte_func_mean
    temperature = temperature
    eigenstrain_name = eigenstrain
  []
[]

[BCs]
  [sides_convective_BC]
    type = ConvectiveFluxFunction
    T_infinity = ${T0}
    variable = temperature
    boundary = 'left right'
    coefficient = ${h}
  []
  [pin_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'pin'
    value = 0.0
  []
[]

[Functions]
  [conductivity]
    type = ParsedFunction
    expression = '${fparse k0*eV_to_J} * t'
  []
  [cte_func_mean]
    type = ParsedFunction
    symbol_names = 'T0'
    symbol_values = '293'
    expression = '1 / (2*sqrt(T0*t))'
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'
  verbose = true
  nl_abs_tol = 5e-6
  nl_rel_tol = 5e-6
  petsc_options = '-ksp_snes_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  snesmf_reuse_base = false
  line_search = 'none'
  nl_max_its = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e1
    optimal_iterations = 10
    growth_factor = 1.25
    cutback_factor = 0.8
  []
[]

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

[Postprocessors]
  [source_integral]
    type = ElementIntegralVariablePostprocessor
    variable = kappa_fission
    execute_on = 'transfer initial timestep_end'
  []
  [strain_xx]
    type = ElementAverageValue
    variable = strain_xx
  []
  [total_length]
    type = FunctionElementIntegral
    function = 1
    execute_on = 'TIMESTEP_BEGIN LINEAR TIMESTEP_END'
  []
[]

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

Problem Description
Analytical Solutions
Cardinal Multi - Group Monte Carlo Model
Cardinal Random Ray Solutions
References