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 or , and vacuum boundary conditions are applied at (where is the deformed length of the slab). The S2 neutron transport equation may be written as a neutron diffusion equation under these constraints:

(1)

with the following Marshak boundary condition:

(2)

is the scalar neutron flux (cm ). is the total cross section (cm), which has an inverse dependence on temperature:

(3)

where is the reference temperature of the slab (K) and is the slab temperature (K). is the combined scattering and k-eigenvalue fission source:

(4)

where is the k-eigenvalue, is the average number of neutrons released per fission reaction, is the fractional probability of a fission reaction and is the fractional probability of a scattering reaction. Both reaction probabilities are constant over the length of the slab. is defined as:

(5)

where is Avogadro's number, is the molar mass of the slab (g mol), is the unheated slab density (g cm), and is the reference microscopic total cross section (cm). The neutron flux solution is normalized to a chosen power (eV s):

(6)

where is the energy released per neutron interaction (eV). The temperature distribution in the slab is governed by the following heat conduction equation:

(7)

with the following convection boundary condition:

(8)

is the thermal conductivity of the material (eV s cm K), which has a linear dependence on temperature:

(9)

where is the reference thermal conductivity (eV s cm K). is a convective heat transfer coefficient (eV s cm K). The final set of physics considered in this problem is thermal expansion, where the differential slab length (cm) expands according to the local thermal strain :

(10)(11)

where is the initial slab length, and the linear expansion coefficient (K) is:

(12)

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

(13)

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, and 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)(15)

and

(16)

The neutron flux in the center of the slab () is arbitrary, and is therefore considered an input setting. The slab eigenvalue can be calculated with Eq. (17):

(17)

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

(18)

where is the Gauss hypergeometric function. The manufactured benchmark parameters are calculated with Eq. (19) and Eq. (20):

(19)(20)

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

(21)

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 cm and a slab k-eigenvalue of .

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

ParameterValue
(g cm-3)
(cm)
(g mol-1)
(K)
(eV)
(eV s-1)
(eV s-1 K-2)
(cm-2 s-1)

Cardinal Multi-Group Monte Carlo Results

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 S2 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.

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 cm is used to ensure rays reach a vacuum boundary condition (to remove initialization bias), which is followed by an active length of cm. Temperature feedback is implemented by interpolating a cross section grid with a spacing of 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 . The resulting flux and temperature distributions can be found in Figure 1.

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

Figure 1: 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 ( pcm) and an expanded slab length 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
  []
[]
(test/tests/neutronics/mg/rr_thermal_expansion/openmc.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<<<{"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'
[]
(test/tests/neutronics/mg/rr_thermal_expansion/solid.i)

References

  1. 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.[Export]
  2. 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.[Export]