Model Inputs
There are two VTR models stored on the virtual test bed; a standalone neutronics model and multi-physics model. There are additional sub-app input files that are called by the multi-physics model that can also be ran standalone with slight modification. In the preceding sections, we will walk through the input files for the Griffin standalone model and multi-physics model.
Griffin Standalone
The complete input file for the Griffin neutronics model is shown below.
The mesh for the Griffin model and the cross sections and equivalence factors are hosted on LFS. Please refer to LFS instructions to download them.
# ==============================================================================
# Model description
# Application : Griffin
# ------------------------------------------------------------------------------
# Idaho Falls, INL, June 23rd, 2022
# Author(s): Nicolas Martin
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================
# - VTR GRIFFIN standalone neutronics input
# ==============================================================================
# - The Model has been built based on [1].
# ------------------------------------------------------------------------------
# [1] Heidet, F. and Roglans-Ribas, J. Core design activities of the
# versatile test reactor -- conceptual phase. EPJ Web Conf. 247(2021)
# 01010. doi:10.1051/epjconf/202124701010.
# URL https://doi.org/10.1051/epjconf/202124701010.
# ==============================================================================
# MODEL PARAMETERS
# ==============================================================================
# State ------------------------------------------------------------------------
initial_fuel_temperature = 600. # (K)
initial_salt_temperature = 595. # (K)
initial_cr_fraction = 0.0
# Power ------------------------------------------------------------------------
tpow = 300e6 #(300 MW)
# ==============================================================================
# GLOBAL PARAMETERS
# ==============================================================================
[GlobalParams]
is_meter = true
grid_names = 'tfuel tcool cr'
grid_variables = 'tfuel tcool cr'
library_file = 'xs/vtr_xs_verification.xml'
library_name = 'vtr_xs'
isotopes = 'pseudo'
densities = '1.0'
[]
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Mesh]
[fmg]
type = FileMeshGenerator
file = mesh/vtr_core.e
exodus_extra_element_integers = 'equivalence_id material_id'
[]
[eqvid]
type = ExtraElementIDCopyGenerator
input = fmg
source_extra_element_id = equivalence_id
target_extra_element_ids = 'equivalence_id'
[]
[]
[Equivalence]
type = SPH
transport_system = CFEM-Diffusion
equivalence_library = 'sph/vtr_sph_verification.xml'
library_name = 'vtr_xs'
compute_factors = false
equivalence_grid_names = 'tfuel tcool cr'
equivalence_grid_variables = 'tfuel tcool cr'
[]
# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================
[Variables]
[]
[Kernels]
[]
# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
[tfuel]
initial_condition = ${initial_fuel_temperature}
[]
[tcool]
initial_condition = ${initial_salt_temperature}
[]
[cr]
family = MONOMIAL
order = CONSTANT
initial_condition = ${initial_cr_fraction}
[]
[]
[AuxKernels]
[]
# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[ICs]
[]
[Functions]
[cr_withdrawal]
type = ConstantFunction
value = 1.92 # 1.92 = fully out, 0.88 = fully in
[]
[]
# ==============================================================================
# MATERIALS AND USER OBJECTS
# ==============================================================================
[Materials]
[fuel]
type = CoupledFeedbackMatIDNeutronicsMaterial
block = 1
plus = 1
[]
[refl_shield_sr] # reflector, shield, SR assemblies
type = CoupledFeedbackMatIDNeutronicsMaterial
block = 2
plus = 0
[]
[control_rods] # control rods modeled via RoddedNeutronicsMaterial
type = CoupledFeedbackRoddedNeutronicsMaterial
block = 3
segment_material_ids = '115 116 115' # 115 = non-abs, 116=abs
rod_segment_length = 1.0424
front_position_function = cr_withdrawal
rod_withdrawn_direction = y # Y-axis is the one vertical
isotopes = 'pseudo; pseudo; pseudo'
densities = '1.0 1.0 1.0'
plus = 0
[]
[]
[PowerDensity]
power = ${tpow}
power_density_variable = power_density
integrated_power_postprocessor = integrated_power
power_scaling_postprocessor = power_scaling
family = MONOMIAL
order = CONSTANT
[]
# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[BCs]
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[TransportSystems]
particle = neutron
equation_type = eigenvalue
G = 6
VacuumBoundary = '1 2 3'
[CFEM-Diffusion]
scheme = CFEM-Diffusion
n_delay_groups = 6
family = LAGRANGE
order = FIRST
assemble_scattering_jacobian = true
assemble_fission_jacobian = true
[]
[]
[Executioner]
type = Eigenvalue
solve_type = 'PJFNKMO'
constant_matrices = true
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 100'
[Quadrature]
order = FOURTH
[]
[]
# ==============================================================================
# MULTIAPPS AND TRANSFER
# ==============================================================================
[MultiApps]
[]
[Transfers]
[]
# ==============================================================================
# RESTART
# ==============================================================================
[RestartVariables]
[]
# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
[power]
type = ElementIntegralVariablePostprocessor
variable = power_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Debug]
[]
[Outputs]
file_base = out_griffin_only
exodus = true
csv = true
perf_graph = true
[console]
type = Console
[]
[]Model Parameters
The first block is for model parameters. This could include the reactor state point (fuel temperature, coolant temperature/density, etc.), reactor power level, or any user-desired information. For consistency, the initial state point is set with the state on the cross section file. The cross sections for this model are not parameterized, instead, there only exists one state point. We'll see that the standalone model is used to verify the neutronics solution and multi-group cross section library to Serpent. The power for the VTR is 300 MW and set with tpow.
Global Parameters
Global parameters are parameters that may be referenced throughout the input file. This block is user specific; the purpose is to simplify repeated variable entries. However, be careful to not override a default parameter option in a later block without meaning to.
For this model, we store cross section information here. Notice that the library name and file location are defined. For simplicity, all of the cross section files are stored in the xs directory. Also, the cross section parameterization (names and variables) on the xml file are set. Even though the library is not parameterized, we still note the state variables at which the cross sections are generated at. The unit of the cross sections are specified with the logical, is_meter, either meters or centimeters. We also define the materials as a "pseudo" mixture with a density of 1.0, which refers to the homogenized cross sections generated with Serpent and processed with ISOXML.
If these definitions do not make sense now, we will walk-through their uses and meanings in the subsequent blocks.
[GlobalParams]
is_meter = true
grid_names = 'tfuel tcool cr'
grid_variables = 'tfuel tcool cr'
library_file = 'xs/vtr_xs_verification.xml'
library_name = 'vtr_xs'
isotopes = 'pseudo'
densities = '1.0'
[]Geometry and Mesh
In this section, we will cover the mesh inputs. The full mesh block can be found below.
[Mesh]
[fmg]
type = FileMeshGenerator
file = mesh/vtr_core.e
exodus_extra_element_integers = 'equivalence_id material_id'
[]
[eqvid]
type = ExtraElementIDCopyGenerator
input = fmg
source_extra_element_id = equivalence_id
target_extra_element_ids = 'equivalence_id'
[]
[]The computational mesh is generated with CUBIT, an external code developed at Sandia National Laboratory, from the VTR geometry. An internal INL tool is used to post-process the mesh to ensure consistency with the Serpent generated multi-group cross section library (i.e. material region and equivalent region identification). The resulting output is an Exodus file that is read by MOOSE using the FileMeshGenerator type. We need to directly specify the ids on the mesh using exodus_extra_element_integers which will then read in the material ids and equivalent ids (for the SPH factors). This is performed with the ExtraElementIDCopyGenerator type by defining the "equivalence_id" as the source (from exodus_extra_element_integers) and "equivalence_id" as the target. A simple way to identify the material and equivalent ids is to open the mesh file vtr_core.e in an Exodus supported visualization tool (i.e. ParaView). All of the mesh files are stored in the mesh directory.
Equivalence
Equivalent cross sections are handled in the [Equivalent] block. Super homogenization (SPH) factors are generated from a Serpent solution to correct homogenization errors and preserve the reaction rates of the heterogenous solution. For convenience, all of the equivalence files are stored in the sph directory. The equivalence library is first defined with equivalence_library which is an xml file generated by the ISOXML pre-processing code. The equivalence_grid_names and equivalence_grid_variables are then specified with the state point variables on the equivalence library. These variables define the SPH factor parameterization. In this model, the cross section library consists of only one state point, made up of fuel temperature, coolant temperature, and control rod fraction state parameters. The compute_factors = false simply tells the system to use the SPH factors directly from the equivalence library. The "CFEM-Diffusion" is defined in the [Transport Systems] block.
[Equivalence]
type = SPH
transport_system = CFEM-Diffusion
equivalence_library = 'sph/vtr_sph_verification.xml'
library_name = 'vtr_xs'
compute_factors = false
equivalence_grid_names = 'tfuel tcool cr'
equivalence_grid_variables = 'tfuel tcool cr'
[]AuxVariables and AuxKernels
There are three AuxVariables that are defined in this model; the fuel temperature, coolant temperature, and control rod fraction. These variables are defined to represent the state parameters that parameterize the cross section and equivalence libraries. The values for these variables are set with a parsed expression from the model parameters defined at the top of the input file.
[AuxVariables]
[tfuel]
initial_condition = ${initial_fuel_temperature}
[]
[tcool]
initial_condition = ${initial_salt_temperature}
[]
[cr]
family = MONOMIAL
order = CONSTANT
initial_condition = ${initial_cr_fraction}
[]
[]There are no AuxKernels that act on these variables since there are no derived quantities to obtain.
Initial Conditions and Functions
User specified functions are defined in this block. Here, we construct a function for the control rod movement. The type is set with the ConstantFunction parameter and the value equal to the all rods out case (1.92). This value can be changed to obtain ARO, ARI, and fractional rod cases.
[Functions]
[cr_withdrawal]
type = ConstantFunction
value = 1.92 # 1.92 = fully out, 0.88 = fully in
[]
[]Materials
Material cross sections are specified with the multi-group cross section library defined by library_file in the [GlobalParameters] block. In this model, we define three types of materials: fuel materials, materials for the reflector/shield/SR assemblies, and control rod materials. The first two materials are given the CoupledFeedbackMatIDNeutronicsMaterial type. This type uses neutronics properties based on mixed isotopes from the multi-group cross section library. It also allows the assignment of cross sections based on dynamic state variables (i.e. the AuxVariables defined above). In this example, we have a pseudo mixture and density that is defined in the [GlobalParameters] block. We must specify which blocks are fissile materials with the block card. The block ids (or names) can be identified on the mesh with an Exodus supported viewer (i.e. ParaView). For each block id, we also need to identify the densities and isotopes. For convenience, these parameters were put in the [GlobalParameters] block. The plus card indicates if absorption, fission, and kappa fission are to be used (for fissile materials).
The control rod materials are given by the CoupledFeedbackRoddedNeutronicsMaterial type. First, the segment_material_ids sets the materials of the control rod. The absorption region of the control rod is identified with the 116 id. The length of the rod segments are set with segment_material_ids. Only the active control rod segment lengths are declared. Next, a control rod function is specified for withdrawal as well as the withdrawal direction with front_position_function and rod_withdrawn_direction, respectively. The control rod function is defined in the [Functions] block. Lastly, the pseudo isotope mixture and density is declared for the control rod segments.
[Materials]
[fuel]
type = CoupledFeedbackMatIDNeutronicsMaterial
block = 1
plus = 1
[]
[refl_shield_sr]
# reflector, shield, SR assemblies
type = CoupledFeedbackMatIDNeutronicsMaterial
block = 2
plus = 0
[]
[control_rods]
# control rods modeled via RoddedNeutronicsMaterial
type = CoupledFeedbackRoddedNeutronicsMaterial
block = 3
segment_material_ids = '115 116 115' # 115 = non-abs, 116=abs
rod_segment_length = 1.0424
front_position_function = cr_withdrawal
rod_withdrawn_direction = y # Y-axis is the one vertical
isotopes = 'pseudo; pseudo; pseudo'
densities = '1.0 1.0 1.0'
plus = 0
[]
[]The power density is not calculated implicitly and must be defined in the [PowerDensity] block. This block provides the power density with a user specified power and the proper flux scaling factor in the postprocessor. The two required arguments are power and power_density_variable. The power is in units of Watts. The power_density_variable simply sets the variable name. Optional arguments include the post-processors to use. The integrated_power_postprocessor sets the desired name for the integrated power post-processor. Likewise, the power_scaling_postprocessor sets the desired name for the power scaling post-processor. Lastly, the family and order parameters are the polynomial representations of the power density corresponding to the underlying FEM.
[PowerDensity]
power = ${tpow}
power_density_variable = power_density
integrated_power_postprocessor = integrated_power
power_scaling_postprocessor = power_scaling
family = MONOMIAL
order = CONSTANT
[]Transport System
The transport system is the heart of the Griffin simulation. We define the transport particle and eigenvalue equation with particle and equation_type, respectively. The number of energy groups is set to "6" with the parameter G. Vacuum boundary conditions are imposed on side sets "1 2 3". The [CFEM-Diffusion] sub-block defines the method to solve the problem. We use Griffin's diffusion method built on the MOOSE continuous finite element method solver. Finite element parameters are set with the family and order, which define the family of functions used to approximate the solution and polynomial order. The last two options, assemble_scattering_jacobian and assemble_fission_jacobian are required to use the "PJFNKMO" method defined later in the [Executioner]. This tells the finite element solver to not update the material cross sections at each linear iteration.
[TransportSystems]
particle = neutron
equation_type = eigenvalue
G = 6
VacuumBoundary = '1 2 3'
[CFEM-Diffusion]
scheme = CFEM-Diffusion
n_delay_groups = 6
family = LAGRANGE
order = FIRST
assemble_scattering_jacobian = true
assemble_fission_jacobian = true
[]
[]Executioner
The [Executioner] block tells the solver what type of problem it needs to solve. Here, we select Eigenvalue as the executioner type which will solve the eigenvalue problem for the criticality and multi-group scalar flux (eigenvalue and eigenvectors). We also specify to solve with a Preconditioned Jacobian Free Newton Krylov Matrix Only method by setting solve_type equal to "PJFNKMO". The Matrix Only method forces the solver to not update material properties (i.e. cross sections) in the linear iterations of the solve. To use this solver, the constant_matrices parameter must be set to "true". There are a number of optional arguments that may also be included. For example, we define a few petsc options and non-linear solver tolerances.
The quadrature sub-block adds a custom quadrature rule. In this case, it is set to a fourth order.
[Executioner]
type = Eigenvalue
solve_type = 'PJFNKMO'
constant_matrices = true
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 100'
[Quadrature]
order = FOURTH
[]
[]Multi-apps and Transfers
This model is to run Griffin neutronics only. In the multi-physics model, we will declare sub-apps that include BISON and SAM.
Post-processors, Debug, and Outputs
The last blocks are for post-processors, debug options, and outputs. A post-processor can be thought of as an operator to compute a single scalar quantity of interest from the solution. For example, the total core power. This is defined with the ElementIntegralVariablePostprocessor type and integrates the power density variable we created earlier, power_density. The execute_on parameter declares at what point in the calculation to compute the total power.
[Postprocessors]
[power]
type = ElementIntegralVariablePostprocessor
variable = power_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[]Finally, the output block sets the output files from the simulation. Two of the most common options include the exodus and csv file. The Exodus file can be viewed with the same software as the mesh, but now will show the solution and solution derived quantities such as the multi-group scalar flux and power distribution. The csv file stores a summary of the solution that includes the k-effective. The perf_graph parameter is helpful to evaluate the computational run time.
[Outputs]
file_base = out_griffin_only
exodus = true
csv = true
perf_graph = true
[console]
type = Console
[]
[]How to run the model
The model can be ran with Griffin in serial or parallel. In serial, the solve takes approximately 1 min 40 sec on a Xeon Platinum 8268 processor. In parallel, with 48 processors, the solve takes approximately 6 seconds.
Run it via:
griffin-opt -i griffin_only.i
and in parallel,
mpirun -n 48 griffin-opt -i griffin_only.i
Multi-physics Model
The complete input file for the multi-physics model is shown below.
# ==============================================================================
# Model description
# Application : Griffin
# ------------------------------------------------------------------------------
# Idaho Falls, INL, June 23rd, 2022
# Author(s): Nicolas Martin
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================
# - VTR GRIFFIN neutronics input
# - Main Application
# ==============================================================================
# - The Model has been built based on [1-2].
# ------------------------------------------------------------------------------
# [1] Heidet, F. and Roglans-Ribas, J. Core design activities of the
# versatile test reactor -- conceptual phase. EPJ Web Conf. 247(2021)
# 01010. doi:10.1051/epjconf/202124701010.
# URL https://doi.org/10.1051/epjconf/202124701010.
# ==============================================================================
# MODEL PARAMETERS
# ==============================================================================
# ==============================================================================
# GLOBAL PARAMETERS
# ==============================================================================
[GlobalParams]
is_meter = true
grid_names = 'tfuel tcool cr'
grid_variables = 'tfuel tcool cr'
library_file = 'xs/vtr_xs.xml'
library_name = 'vtr_xs'
isotopes = 'pseudo'
densities = '1.0'
[]
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Mesh]
[fmg]
type = FileMeshGenerator
file = mesh/vtr_core.e # Y-axis = vertical direction for consistency with BISON solid mechanics RZ model
exodus_extra_element_integers = 'equivalence_id material_id'
[]
[eqvid]
type = ExtraElementIDCopyGenerator
input = fmg
source_extra_element_id = equivalence_id
target_extra_element_ids = 'equivalence_id'
[]
parallel_type = replicated
displacements = 'disp_x disp_y disp_z'
[]
[Equivalence]
type = SPH
transport_system = CFEM-Diffusion
equivalence_library = 'sph/vtr_sph.xml'
library_name = 'vtr_xs'
compute_factors = false
equivalence_grid_names = 'tfuel tcool cr'
equivalence_grid_variables = 'tfuel tcool cr'
[]
# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================
[Variables]
[]
[Kernels]
[]
# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
[tfuel]
initial_condition = 900.
[]
[tcool]
initial_condition = 700.
[]
[disp_x] # from grid plate expansion
initial_condition = 0
[]
[disp_z] # from grid plate expansion
initial_condition = 0
[]
[disp_y] # from fuel expansion
initial_condition = 0
[]
[cr]
family = MONOMIAL
order = CONSTANT
initial_condition = 0.
[]
[]
# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[ICs]
[]
[Functions]
[cr_withdrawal]
type = ConstantFunction
value = 1.9183 # 1.9183 = fully out, 0.876 = fully in
[]
[]
# ==============================================================================
# MATERIALS AND USER OBJECTS
# ==============================================================================
[Materials]
[fuel]
type = CoupledFeedbackMatIDNeutronicsMaterial
block = 1
plus = 1
displacements = 'disp_x disp_y disp_z'
[]
[refl_shield_sr] # reflector, shield, SR assemblies
type = CoupledFeedbackMatIDNeutronicsMaterial
block = 2
plus = 0
displacements = 'disp_x disp_y disp_z'
[]
[control_rods] # control rods modeled via RoddedNeutronicsMaterial
type = CoupledFeedbackRoddedNeutronicsMaterial
block = 3
segment_material_ids = '115 116 115' # 115 = non-abs, 116=abs
rod_segment_length = 1.0424
front_position_function = cr_withdrawal
rod_withdrawn_direction = y # Y-axis is the one vertical
isotopes = 'pseudo; pseudo; pseudo'
densities = '1.0 1.0 1.0'
plus = 0
displacements = 'disp_x disp_y disp_z'
[]
[]
[PowerDensity]
power = 300e6
power_density_variable = power_density # name of AuxVariable to be created
integrated_power_postprocessor = integrated_power
power_scaling_postprocessor = power_scaling
family = MONOMIAL
order = CONSTANT
[]
[UserObjects]
[power_density1_UO]
type = NearestPointLayeredAverage
block = 1
variable = power_density
direction = y
num_layers = 40
points_file = coordinates_orifice_1.txt
execute_on = TIMESTEP_END
[]
[power_density2_UO]
type = NearestPointLayeredAverage
block = 1
variable = power_density
direction = y
num_layers = 40
points_file = coordinates_orifice_2.txt
execute_on = TIMESTEP_END
[]
[power_density3_UO]
type = NearestPointLayeredAverage
block = 1
variable = power_density
direction = y
num_layers = 40
points_file = coordinates_orifice_3.txt
execute_on = TIMESTEP_END
[]
[]
# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[BCs]
[]
# ==============================================================================
# TRANSPORT SYSTEMS
# ==============================================================================
[TransportSystems]
particle = neutron
equation_type = eigenvalue
G = 6
VacuumBoundary = '1 2 3'
[CFEM-Diffusion]
scheme = CFEM-Diffusion
n_delay_groups = 6
family = LAGRANGE
order = FIRST
use_displaced_mesh = true
assemble_scattering_jacobian = true
assemble_fission_jacobian = true
[]
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Executioner]
automatic_scaling = True
type = Eigenvalue
solve_type = 'PJFNKMO'
constant_matrices = true
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 100'
free_power_iterations = 4
nl_abs_tol = 1e-8
nl_rel_tol = 1e-7
#fixed_point specific
fixed_point_force_norms = true
fixed_point_abs_tol = 1e-10
fixed_point_rel_tol = 1e-10
fixed_point_max_its = 3
accept_on_max_fixed_point_iteration = true
[]
# ==============================================================================
# MULTIAPPS AND TRANSFER
# ==============================================================================
[MultiApps]
#one BISON/SAM channel per orifice - currently 3 zones
[bison_1]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions_file = 'coordinates_orifice_1.txt'
input_files = bison_thermal_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
cli_args = "MultiApps/sam/cli_args='m_dot_in=29.8'"
[]
[bison_2]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions_file = 'coordinates_orifice_2.txt'
input_files = bison_thermal_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
cli_args = "MultiApps/sam/cli_args='m_dot_in=23.2'"
[]
[bison_3]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions_file = 'coordinates_orifice_3.txt'
input_files = bison_thermal_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
cli_args = "MultiApps/sam/cli_args='m_dot_in=15.9'"
[]
[core_support_plate]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions = '0 0 0'
input_files = core_support_plate_3d.i
execute_on = 'INITIAL' # no need for other calls - since tinlet is fixed only one update of the XS required
max_procs_per_app = 1
output_in_position = true
#no_restore = true
[]
[]
[Transfers]
# master -> sub-app
#-----------------------
# power_density to BISON
#-----------------------
[pdens_to_bison_1]
type = MultiAppProjectionTransfer
to_multi_app = bison_1
source_variable = power_density
variable = power_density
[]
[pdens_to_bison_2]
type = MultiAppProjectionTransfer
to_multi_app = bison_2
source_variable = power_density
variable = power_density
[]
[pdens_to_bison_3]
type = MultiAppProjectionTransfer
to_multi_app = bison_3
source_variable = power_density
variable = power_density
[]
# sub-app -> master
#------------------
# tfuel from BISON
#------------------
[tfuel_from_bison_1]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_1
source_variable = tfuel
variable = tfuel
[]
[tfuel_from_bison_2]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_2
source_variable = tfuel
variable = tfuel
[]
[tfuel_from_bison_3]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_3
source_variable = tfuel
variable = tfuel
[]
#------------------
# tcool from BISON/SAM
#------------------
[tcool_from_bison_1]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_1
source_variable = tcool
variable = tcool
[]
[tcool_from_bison_2]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_2
source_variable = tcool
variable = tcool
[]
[tcool_from_bison_3]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_3
source_variable = tcool
variable = tcool
[]
#------------------
# axial expansion = vertical displacements disp_y
#------------------
[disp_y_from_bison_1]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_1
source_variable = disp_y
variable = disp_y
[]
[disp_y_from_bison_2]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_2
source_variable = disp_y
variable = disp_y
[]
[disp_y_from_bison_3]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_3
source_variable = disp_y
variable = disp_y
[]
#------------------
# Radial expansion = horizontal displacements disp_x, disp_z from core_support_plate
#------------------
[disp_x_from_core_support_plate]
type = MultiAppNearestNodeTransfer
from_multi_app = core_support_plate
source_boundary = 'plateTop'
source_variable = disp_x
variable = disp_x
fixed_meshes = true
[]
[disp_z_from_core_support_plate]
type = MultiAppNearestNodeTransfer
from_multi_app = core_support_plate
source_boundary = 'plateTop'
source_variable = disp_z
variable = disp_z
fixed_meshes = true
[]
[bison_reporter_1]
type = MultiAppReporterTransfer
to_reporters = 'max_tfuel_r/value max_tclad_r/value max_tcool_r/value max_thcond_r/value'
from_reporters = 'max_tfuel/value max_tclad/value max_tcool/value max_thcond/value'
from_multi_app = bison_1
subapp_index = 0
[]
[]
# ==============================================================================
# RESTART
# ==============================================================================
[RestartVariables]
[]
# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
[avg_tfuel]
type = ElementAverageValue
variable = tfuel
block = 1
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tfuel]
type = ElementExtremeValue
variable = tfuel
block = 1
value_type = max
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tfuel]
type = ElementExtremeValue
variable = tfuel
block = 1
value_type = min
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tcool]
type = ElementAverageValue
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tcool]
type = ElementExtremeValue
variable = tcool
value_type = max
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tcool]
type = ElementExtremeValue
variable = tcool
value_type = min
execute_on = 'INITIAL TIMESTEP_END'
[]
[ptot]
type = ElementIntegralVariablePostprocessor
variable = power_density
execute_on = ' INITIAL TIMESTEP_END'
[]
[disp_x_max]
type = ElementExtremeValue
variable = disp_x
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[disp_y_max]
type = ElementExtremeValue
variable = disp_y
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[disp_z_max]
type = ElementExtremeValue
variable = disp_z
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[max_tfuel_r]
type = Receiver
[]
[max_tclad_r]
type = Receiver
[]
[max_tcool_r]
type = Receiver
[]
[max_thcond_r]
type = Receiver
[]
[]
[Debug]
[]
[Outputs]
exodus = true
csv = true
perf_graph = true
[]The neutronics input is similar to the standalone model that is described in the previous section. Still, we will highlight input additions that allow for multi-physics coupling (multi-apps and transfers). We'll also go into details on the sub-app input files (BISON, SAM, Tensor Mechanics).
MultiApps
In addition to the main-app neutronics input, there are now BISON inputs, SAM inputs, and a MOOSE Tensor Mechanics input. These sub-apps are defined in the [MultiApps] block. Note that only the sub-apps that are directly coupled to Griffin are defined here. In other words, the SAM model will be declared as a sub-app in the BISON model.
[MultiApps]
#one BISON/SAM channel per orifice - currently 3 zones
[bison_1]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions_file = 'coordinates_orifice_1.txt'
input_files = bison_thermal_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
cli_args = "MultiApps/sam/cli_args='m_dot_in=29.8'"
[]
[bison_2]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions_file = 'coordinates_orifice_2.txt'
input_files = bison_thermal_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
cli_args = "MultiApps/sam/cli_args='m_dot_in=23.2'"
[]
[bison_3]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions_file = 'coordinates_orifice_3.txt'
input_files = bison_thermal_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
cli_args = "MultiApps/sam/cli_args='m_dot_in=15.9'"
[]
[core_support_plate]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions = '0 0 0'
input_files = core_support_plate_3d.i
execute_on = 'INITIAL' # no need for other calls - since tinlet is fixed only one update of the XS required
max_procs_per_app = 1
output_in_position = true
#no_restore = true
[]
[]First, we start with the BISON sub-apps. There are three BISON MultiApps that correspond to the fuel rods in the three orifices of the VTR. For each of these models, we give the sub-app a name (i.e. bison_1, bison_2, bison_3). The type and app_type is set with FullSolveMultiApp and BlueCrabApp, respectively. A different command line argument is passed to override a single parameter (the mass flow rate) in the input files to form the three different simulations. The former performs a complete simulation every time the sub-app executes. The latter informs the code what type of app to build. We want BlueCrab for multi-physics simulations because it is a package of Griffin, BISON, and SAM. Next, the positions_file points to a file that contains spatial coordinates that define a point for the app location. The input_files simply points to the sub-app input file. We'll also set the sub-app to run its simulation at the end of the main-app solve (Griffin) with execute_on = 'TIMESTEP_END'. Next, we specify the sub-app to run on one processor with max_procs_per_app. Finally, the output_in_position allows the output to be 'moved' by its position vector.
The core support plate is modeled with the Tensor Mechanics module included in MOOSE. The sub-app is defined in a similar manner as the BISON apps. However, we only need to run this sub-app at the beginning of the simulation since the inlet temperature is fixed (the resulting radial displacements are also fixed). This is done by setting execute_on = 'INITIAL'.
Transfers
Variables are transfered internally to each app with MOOSE. We just need to specify which variables are transfered to each app and in which direction.
[Transfers]
# master -> sub-app
#-----------------------
# power_density to BISON
#-----------------------
[pdens_to_bison_1]
type = MultiAppProjectionTransfer
to_multi_app = bison_1
source_variable = power_density
variable = power_density
[]
[pdens_to_bison_2]
type = MultiAppProjectionTransfer
to_multi_app = bison_2
source_variable = power_density
variable = power_density
[]
[pdens_to_bison_3]
type = MultiAppProjectionTransfer
to_multi_app = bison_3
source_variable = power_density
variable = power_density
[]
# sub-app -> master
#------------------
# tfuel from BISON
#------------------
[tfuel_from_bison_1]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_1
source_variable = tfuel
variable = tfuel
[]
[tfuel_from_bison_2]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_2
source_variable = tfuel
variable = tfuel
[]
[tfuel_from_bison_3]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_3
source_variable = tfuel
variable = tfuel
[]
#------------------
# tcool from BISON/SAM
#------------------
[tcool_from_bison_1]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_1
source_variable = tcool
variable = tcool
[]
[tcool_from_bison_2]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_2
source_variable = tcool
variable = tcool
[]
[tcool_from_bison_3]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_3
source_variable = tcool
variable = tcool
[]
#------------------
# axial expansion = vertical displacements disp_y
#------------------
[disp_y_from_bison_1]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_1
source_variable = disp_y
variable = disp_y
[]
[disp_y_from_bison_2]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_2
source_variable = disp_y
variable = disp_y
[]
[disp_y_from_bison_3]
type = MultiAppGeometricInterpolationTransfer
from_multi_app = bison_3
source_variable = disp_y
variable = disp_y
[]
#------------------
# Radial expansion = horizontal displacements disp_x, disp_z from core_support_plate
#------------------
[disp_x_from_core_support_plate]
type = MultiAppNearestNodeTransfer
from_multi_app = core_support_plate
source_boundary = 'plateTop'
source_variable = disp_x
variable = disp_x
fixed_meshes = true
[]
[disp_z_from_core_support_plate]
type = MultiAppNearestNodeTransfer
from_multi_app = core_support_plate
source_boundary = 'plateTop'
source_variable = disp_z
variable = disp_z
fixed_meshes = true
[]
[bison_reporter_1]
type = MultiAppReporterTransfer
to_reporters = 'max_tfuel_r/value max_tclad_r/value max_tcool_r/value max_thcond_r/value'
from_reporters = 'max_tfuel/value max_tclad/value max_tcool/value max_thcond/value'
from_multi_app = bison_1
subapp_index = 0
[]
[]At a high level, the scheme works as follows. The power density is solved with Griffin and passed to BISON. BISON solves for the temperature distribution of the fuel and transfers the wall temperature to SAM. SAM calculates the coolant temperature and heat transfer coefficient, which then gets passed back to BISON. Following the BISON and SAM iterations, the fuel temperature gets passed back to Griffin and the cross sections are updated.
In a separate iteration between the thermal and mechanical BISON models, the mechanical model solves and returns the thermal expansion of the fuel. The thermal expansion of the core support plate is solved in a native MOOSE Tensor Mechanics module and the radial displacements are passed back to the Griffin main-app.
The transfers are organized in the following way: transfer type, direction, and variable. For example, the power density transfer from Griffin to BISON is specified with the MultiAppProjectionTransfer type and direction with the to_multi_app parameter. The source_variable and variable is the power_density.
How to run the model
The model can be ran by executing the BlueCrab app that consists of Griffin, BISON, and SAM. In parallel, with 48 processors, the solve takes approximately 52 seconds.
Run it via:
mpirun -n 48 blue_crab-opt -i griffin_multiphysics.i
BISON Thermal Model
The input file for the BISON thermal model of a VTR fuel rod is displayed below.
################################################################################
## VTR fuel rod, derived from BISON assessment case IFR1.i ##
## Bison Sub-Application ##
## Thermal conduction in fuel, closed gap and clad ##
## Convection on clad outer surface with SAM coupling ##
################################################################################
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# Units are in standard SI: J, K, kg, m, Pa, s.
# This input file does not include friction.
# radial/axial dimensions at hot irradiated condition to be consistent with griffin/SAM model
gap = 0. #0.5e-3
rod_outside_diameter = 6.282e-3 # 6.25e-3
clad_thickness = 0.503e-3 # 0.5e-3
slug_diameter = 5.277e-3 # 4.547e-3
fuel_height = 842.4e-3 # 800.0e-3
plenum_height = 782.2e-3 # 778.0e-3
# ==============================================================================
# GLOBAL PARAMETERS
# ==============================================================================
[GlobalParams]
# Parameters that are used in multiple blocks can be included here so that
# they only need to be specified one time.
order = FIRST
family = LAGRANGE
temperature = Temperature
#the following are needed in multiple UPuZr Materials
X_Zr = 0.225 # U-20Pu-10Zr
X_Pu = 0.171
density = 11120.0 # kg/m3 at hot operating condition (UPuZr fuel)
[]
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Mesh]
coord_type = RZ
# based on x447.i from examples
# rod specific parameters - dimensions are for fresh fuel at room temperature
[smeared_pellet_mesh]
type = FuelPinMeshGenerator
clad_thickness = ${clad_thickness}
pellet_outer_radius = '${fparse slug_diameter/2}'
pellet_height = ${fuel_height}
clad_top_gap_height = ${plenum_height} # fixme assumes no Na bond sodium
clad_gap_width = ${gap} # no gap for hot condition after irradiation (2 % fima)
clad_bot_gap_height = 0.
bottom_clad_height = ${clad_thickness}
top_clad_height = ${clad_thickness}
# meshing parameters
clad_mesh_density = customize
pellet_mesh_density = customize
nx_p = 6 # number of fuel elements in radial direction
ny_p = 300 # number of fuel elements in axial direction
nx_c = 3 # number of clad elements in radial direction
ny_c = 300 # number of clad elements in axial direction
ny_cu = 1 #?? number of cladding elements in upper plug in axial direction (default=1)
ny_cl = 1 #?? number of cladding elements in lower plug in axial direction (default=1)
pellet_quantity = 1
elem_type = QUAD4
[]
[add_side_clad]
type = SubdomainBoundingBoxGenerator
location = INSIDE
restricted_subdomains = clad
block_id = '4'
block_name = 'side_clad'
input = smeared_pellet_mesh
bottom_left = '${fparse rod_outside_diameter/2 - clad_thickness/3} 0. 0.'
top_right = '${fparse rod_outside_diameter/2} ${fparse fuel_height + plenum_height} 0.'
[]
[]
# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================
[Variables]
[Temperature]
initial_condition = 900
[]
[]
[Kernels]
active = 'heat_conduction heat_source_from_power'
[heat_time]
type = HeatConductionTimeDerivative
variable = Temperature
[]
[heat_conduction]
type = HeatConduction
variable = Temperature
[]
[heat_source_from_power]
type = CoupledForce
variable = Temperature
block = pellet
v = power_density_scaled
[]
[heat_source_from_fission]
type = FissionRateHeatSource
variable = Temperature
fission_rate = fission_rate
block = pellet
[]
[]
# ==============================================================================
# MODULES
# ==============================================================================
[ThermalContact]
#Action to control heat transfer in regions without meshes. Specifically
#the gap and the coolant.
[thermal_contact]
type = GapHeatTransfer
variable = Temperature
primary = clad_inside_right
secondary = pellet_outer_radial_surface
quadrature = true
gap_conductivity = 61.0 # [Fink and Leibowitz, 1995.] pg. 181 (840 K)
min_gap = 0.38e-03 # [Greenquist et al., 2020] pg. 4
[]
[]
# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
[power_density]
block = pellet
initial_condition = 4.2730e+08 #homogeneous power density - from python script get_vol.py
[]
[power_density_scaled]
block = pellet
[]
[tfuel]
block = pellet
[]
[twall]
block = 4 # side_clad
[]
[tcool]
initial_condition = 750
[]
[htc]
initial_condition = 1e5
[]
[disp_x]
initial_condition = 0
[]
[disp_y]
initial_condition = 0
[]
[]
[AuxKernels]
active = 'GetPowerDensity_from_griffin SetTwall SetTfuel' # replace _from_griffin by _from_func for standalone runs
[GetPowerDensity_from_griffin]
type = NormalizationAux
variable = power_density_scaled
source_variable = power_density
normal_factor = 2.66105 # Vhom/Vhet, from python script get_vol.py
block = pellet
[]
[GetPowerDensity_from_func]
type = FunctionAux
variable = power_density_scaled
function = power_density_shaped
block = pellet
[]
[SetTwall]
type = CoupledAux
block = 4
variable = twall
coupled = Temperature
[]
[SetTfuel]
type = CoupledAux
block = pellet
variable = tfuel
coupled = Temperature
[]
[]
# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[ICs]
[]
[Functions]
[lhgr] # W/m
type = ConstantFunction
value = 24.868e3
[]
# power density in W/m3 - normally from griffin
[power_density_hom]
type = ConstantFunction
value = 4.2725e8 # W/m3
[]
[power_density_het]
type = ConstantFunction
value = 11.369e8 # W/m3
[]
[]
# ==============================================================================
# MATERIALS AND USER OBJECTS
# ==============================================================================
[Materials]
#---------------
# thermal only
#---------------
active = 'fuel_thermal fuel_density clad_thermal clad_density'
#fuel
#thermal materials
[fuel_thermal]
type = UPuZrThermal
block = pellet
spheat_model = savage
thcond_model = lanl
porosity = 0
outputs = all
[]
[fuel_density]
type = Density
block = pellet
[]
#cladding
#thermal materials
[clad_thermal]
type = HT9Thermal
block = 'clad 4'
[]
[clad_density]
type = Density
block = 'clad 4'
density = 7800
[]
[]
[UserObjects]
[disp_y_UO]
type = LayeredAverage
variable = disp_y
direction = y
num_layers = 20
block = pellet
execute_on = TIMESTEP_END
use_displaced_mesh = true
[]
[]
# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[BCs]
#for coupling with SAM (provides both htc and tcool)
[convection_outer_clad]
type = CoupledConvectiveHeatFluxBC
boundary = 'clad_outside_bottom clad_outside_right clad_outside_top'
variable = Temperature
T_infinity = tcool
htc = htc
[]
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Preconditioning]
# Used to improve the solver performance
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Steady
automatic_scaling = true
solve_type = 'PJFNK'
nl_rel_tol = 1e-7
nl_abs_tol = 1e-8
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
petsc_options_value = 'hypre boomeramg 100'
[]
# ==============================================================================
# MULTIAPPS AND TRANSFER
# ==============================================================================
[MultiApps]
[sam]
type = FullSolveMultiApp
app_type = SamApp
positions = '0 0 0'
input_files = sam_channel.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
[]
[bison_mechanics]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions = '0 0 0'
input_files = bison_mecha_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
[]
[]
[Transfers]
[twall_to_sam]
type = MultiAppCoordSwitchNearestNodeTransfer
source_variable = twall
variable = T_wall_external # SAM variable
to_multi_app = sam
displaced_target_mesh = true
fixed_meshes = true
[]
[tcool_from_sam]
type = MultiAppCoordSwitchNearestNodeTransfer
source_variable = temperature
variable = tcool
from_multi_app = sam
displaced_source_mesh = true
fixed_meshes = true
[]
[htc_from_sam]
type = MultiAppCoordSwitchNearestNodeTransfer
source_variable = htc_external
variable = htc
from_multi_app = sam
displaced_source_mesh = true
fixed_meshes = true
[]
[temperature_to_mechanics]
type = MultiAppCopyTransfer
to_multi_app = bison_mechanics
source_variable = Temperature
variable = Temperature
[]
[disp_x_from_mechanics]
type = MultiAppCopyTransfer
from_multi_app = bison_mechanics
source_variable = disp_x
variable = disp_x
[]
[disp_y_from_mechanics]
type = MultiAppCopyTransfer
from_multi_app = bison_mechanics
source_variable = disp_y
variable = disp_y
[]
[sam_reporter]
type = MultiAppReporterTransfer
to_reporters = 'max_tcool_r/value'
from_reporters = 'max_tcool/value'
from_multi_app = sam
subapp_index = 0
[]
[]
# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
[ptot_bison]
type = ElementIntegralVariablePostprocessor
block = pellet
variable = power_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[ptot_scaled_bison]
type = ElementIntegralVariablePostprocessor
block = pellet
variable = power_density_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tfuel]
type = ElementAverageValue
variable = Temperature
block = pellet
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tfuel]
type = ElementExtremeValue
variable = Temperature
value_type = max
block = pellet
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tfuel]
type = ElementExtremeValue
variable = Temperature
value_type = min
block = pellet
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_twall]
type = ElementAverageValue
variable = twall
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_twall]
type = ElementExtremeValue
variable = Temperature
value_type = min
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_twall]
type = ElementExtremeValue
variable = Temperature
value_type = max
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tclad]
type = ElementExtremeValue
variable = Temperature
value_type = max
block = 'clad 4'
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tclad]
type = ElementExtremeValue
variable = Temperature
value_type = min
block = 'clad 4'
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tclad]
type = ElementAverageValue
variable = Temperature
block = 'clad 4'
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_htc]
type = ElementAverageValue
variable = htc
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tcool]
type = ElementAverageValue
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tcool]
type = ElementExtremeValue
value_type = max
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tcool]
type = ElementExtremeValue
value_type = min
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[disp_x_max]
type = ElementExtremeValue
variable = disp_x
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[disp_y_max]
type = ElementExtremeValue
variable = disp_y
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[max_thcond]
type = ElementExtremeValue
variable = thermal_conductivity
[]
[max_tcool_r]
type = Receiver
[]
[]
[Outputs]
# csv = true
# exodus = true
print_nonlinear_converged_reason = false
print_linear_converged_reason = false
[]The top of the input file houses the model parameters. For a BISON model, this would include the rod diameter, clad thickness, fuel height, etc.
Global Parameters
Global parameters are parameters that are used in more than one block. For example, the family of functions and polynomial order are defined for the finite element solver.
[GlobalParams]
# Parameters that are used in multiple blocks can be included here so that
# they only need to be specified one time.
order = FIRST
family = LAGRANGE
temperature = Temperature
#the following are needed in multiple UPuZr Materials
X_Zr = 0.225 # U-20Pu-10Zr
X_Pu = 0.171
density = 11120.0 # kg/m3 at hot operating condition (UPuZr fuel)
[]Geometry and Mesh
The pellet mesh is setup with the SmearedPelletMeshGenerator type. This object generates a smeared pellet mesh for use in simulating fuel rods with pellet fuel. In a smeared mesh, the fuel stack is treated as a single rectangle of fuel (i.e., dishes and chamfers are not included). There are a number of geometry and meshing options that are shown but will not be discussed. The clad is added around the pellet mesh with the SubdomainBoundingBoxGenerator type. More information about the mesh generation may be found in the BISON documentation for the SmearedPelletMeshGenerator and the SubdomainBoundingBoxGenerator.
[Mesh]
coord_type = RZ
# based on x447.i from examples
# rod specific parameters - dimensions are for fresh fuel at room temperature
[smeared_pellet_mesh]
type = FuelPinMeshGenerator
clad_thickness = ${clad_thickness}
pellet_outer_radius = '${fparse slug_diameter/2}'
pellet_height = ${fuel_height}
clad_top_gap_height = ${plenum_height} # fixme assumes no Na bond sodium
clad_gap_width = ${gap} # no gap for hot condition after irradiation (2 % fima)
clad_bot_gap_height = 0.
bottom_clad_height = ${clad_thickness}
top_clad_height = ${clad_thickness}
# meshing parameters
clad_mesh_density = customize
pellet_mesh_density = customize
nx_p = 6 # number of fuel elements in radial direction
ny_p = 300 # number of fuel elements in axial direction
nx_c = 3 # number of clad elements in radial direction
ny_c = 300 # number of clad elements in axial direction
ny_cu = 1 #?? number of cladding elements in upper plug in axial direction (default=1)
ny_cl = 1 #?? number of cladding elements in lower plug in axial direction (default=1)
pellet_quantity = 1
elem_type = QUAD4
[]
[add_side_clad]
type = SubdomainBoundingBoxGenerator
location = INSIDE
restricted_subdomains = clad
block_id = '4'
block_name = 'side_clad'
input = smeared_pellet_mesh
bottom_left = '${fparse rod_outside_diameter/2 - clad_thickness/3} 0. 0.'
top_right = '${fparse rod_outside_diameter/2} ${fparse fuel_height + plenum_height} 0.'
[]
[]Variables and Kernels
A variable is a desired quantity to be solved for in the finite element solution. In the BISON thermal model, we want the temperature distribution of the fuel element.
[Variables]
[Temperature]
initial_condition = 900
[]
[]Kernels are the inner product terms that make up a weak form of a differential equation. The kernels define the differential equation to solve. The active kernels that are used in the solve are identified with active. For example, a time derivative kernel is defined but not activated as we are interested in the steady state solution.
[Kernels]
active = 'heat_conduction heat_source_from_power'
[heat_time]
type = HeatConductionTimeDerivative
variable = Temperature
[]
[heat_conduction]
type = HeatConduction
variable = Temperature
[]
[heat_source_from_power]
type = CoupledForce
variable = Temperature
block = pellet
v = power_density_scaled
[]
[heat_source_from_fission]
type = FissionRateHeatSource
variable = Temperature
fission_rate = fission_rate
block = pellet
[]
[]Modules
ThermalContact is a MOOSE module to control the heat transfer in regions without meshes. In this model, we use [ThermalContact] to model the gap between the pellet and clad. We do this by defining the GapHeatTransfer type and corresponding parameters. Detailed information about the ThermalContact module may be found in the MOOSE documentation.
[ThermalContact]
#Action to control heat transfer in regions without meshes. Specifically
#the gap and the coolant.
[thermal_contact]
type = GapHeatTransfer
variable = Temperature
primary = clad_inside_right
secondary = pellet_outer_radial_surface
quadrature = true
gap_conductivity = 61.0 # [Fink and Leibowitz, 1995.] pg. 181 (840 K)
min_gap = 0.38e-03 # [Greenquist et al., 2020] pg. 4
[]
[]AuxVariables and AuxKernels
AuxVariables are variables that are derived from the solution variable (in this case, temperature). These can be quantities such as the fuel temperature, wall temperature, coolant temperature, etc. We also define the power_density, a variable transfered from the Griffin main-app.
[AuxVariables]
[power_density]
block = pellet
initial_condition = 4.2730e+08 #homogeneous power density - from python script get_vol.py
[]
[power_density_scaled]
block = pellet
[]
[tfuel]
block = pellet
[]
[twall]
block = 4 # side_clad
[]
[tcool]
initial_condition = 750
[]
[htc]
initial_condition = 1e5
[]
[disp_x]
initial_condition = 0
[]
[disp_y]
initial_condition = 0
[]
[]AuxKernels are kernels that act on a variable to derive an AuxVariable. Consider the fuel temperature and wall temperature. These AuxKernels are of type CoupledAux. These copy the values from the solution variable, Temperature, into the respective variables twall and tfuel on the desired blocks. AuxKernels are also defined to obtain (and scale) the power density from Griffin.
[AuxKernels]
active = 'GetPowerDensity_from_griffin SetTwall SetTfuel' # replace _from_griffin by _from_func for standalone runs
[GetPowerDensity_from_griffin]
type = NormalizationAux
variable = power_density_scaled
source_variable = power_density
normal_factor = 2.66105 # Vhom/Vhet, from python script get_vol.py
block = pellet
[]
[GetPowerDensity_from_func]
type = FunctionAux
variable = power_density_scaled
function = power_density_shaped
block = pellet
[]
[SetTwall]
type = CoupledAux
block = 4
variable = twall
coupled = Temperature
[]
[SetTfuel]
type = CoupledAux
block = pellet
variable = tfuel
coupled = Temperature
[]
[]Initial Conditions and Functions
Users may specify custom functions in the Functions block. For example, we have defined the linear heat generation rate, heterogenous and homogenous power densities. These functions would be used in a standalone BISON model if the power density is not provided by Griffin.
[Functions]
[lhgr]
# W/m
type = ConstantFunction
value = 24.868e3
[]
# power density in W/m3 - normally from griffin
[power_density_hom]
type = ConstantFunction
value = 4.2725e8 # W/m3
[]
[power_density_het]
type = ConstantFunction
value = 11.369e8 # W/m3
[]
[]Materials and UserObjects
Material characteristics are defined in the [Materials] block. They include the thermal properties and density of the fuel and clad.
[Materials]
#---------------
# thermal only
#---------------
active = 'fuel_thermal fuel_density clad_thermal clad_density'
#fuel
#thermal materials
[fuel_thermal]
type = UPuZrThermal
block = pellet
spheat_model = savage
thcond_model = lanl
porosity = 0
outputs = all
[]
[fuel_density]
type = Density
block = pellet
[]
#cladding
#thermal materials
[clad_thermal]
type = HT9Thermal
block = 'clad 4'
[]
[clad_density]
type = Density
block = 'clad 4'
density = 7800
[]
[]Boundary Conditions
The boundary conditions in the model are used to couple with SAM, where SAM will provide the coolant temperature and heat transfer coefficient. As such, the type is defined as a CoupledConvectiveHeatFluxBC and coupled to the Temperature solution variable.
[BCs]
#for coupling with SAM (provides both htc and tcool)
[convection_outer_clad]
type = CoupledConvectiveHeatFluxBC
boundary = 'clad_outside_bottom clad_outside_right clad_outside_top'
variable = Temperature
T_infinity = tcool
htc = htc
[]
[]Execution Parameters
Preconditioners transform the given problem into a form that is more suitable for the numerical solver, thus improving solver performance. The single matrix preconditioner (SMP) is defined here that builds a preconditioner using user defined off-diagonal parts of the Jacobian matrix.
[Preconditioning]
# Used to improve the solver performance
[SMP]
type = SMP
full = true
[]
[]The [Executioner] block tells the solver what type of problem it is to solve. Here, we select Steady as the executioner type to identify that the problem is steady state. We also specify to solve with a Preconditioned Jacobian Free Newton Krylov method by setting solve_type equal to "PJFNK".
[Executioner]
type = Steady
automatic_scaling = true
solve_type = 'PJFNK'
nl_rel_tol = 1e-7
nl_abs_tol = 1e-8
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
petsc_options_value = 'hypre boomeramg 100'
[]MultiApps and Transfers
There are two sub-apps that extend from this BISON sub-app input deck; a SAM thermal fluids model and a BISON mechanics model. Both are defined here in the same manner as discussed in MultiApps.
[MultiApps]
[sam]
type = FullSolveMultiApp
app_type = SamApp
positions = '0 0 0'
input_files = sam_channel.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
[]
[bison_mechanics]
type = FullSolveMultiApp
app_type = BlueCrabApp
positions = '0 0 0'
input_files = bison_mecha_only.i
execute_on = 'TIMESTEP_END'
max_procs_per_app = 1
output_in_position = true
#no_restore = true
[]
[]To accompany the multi-apps, variable transfers must be specified for coupling. We want to pass the outer wall temperature calculated in the BISON thermal model to SAM. SAM will then compute the coolant temperature and heat transfer coefficient to pass back to the BISON thermal model. There is also coupling between the thermal BISON model and mechanical model. A transfer is defined to pass the temperature distribution from the thermal model to the mechanical model. The mechanical model then solves for the thermal expansion and transfers the x and y displacements back to the thermal model.
[Transfers]
[twall_to_sam]
type = MultiAppCoordSwitchNearestNodeTransfer
source_variable = twall
variable = T_wall_external # SAM variable
to_multi_app = sam
displaced_target_mesh = true
fixed_meshes = true
[]
[tcool_from_sam]
type = MultiAppCoordSwitchNearestNodeTransfer
source_variable = temperature
variable = tcool
from_multi_app = sam
displaced_source_mesh = true
fixed_meshes = true
[]
[htc_from_sam]
type = MultiAppCoordSwitchNearestNodeTransfer
source_variable = htc_external
variable = htc
from_multi_app = sam
displaced_source_mesh = true
fixed_meshes = true
[]
[temperature_to_mechanics]
type = MultiAppCopyTransfer
to_multi_app = bison_mechanics
source_variable = Temperature
variable = Temperature
[]
[disp_x_from_mechanics]
type = MultiAppCopyTransfer
from_multi_app = bison_mechanics
source_variable = disp_x
variable = disp_x
[]
[disp_y_from_mechanics]
type = MultiAppCopyTransfer
from_multi_app = bison_mechanics
source_variable = disp_y
variable = disp_y
[]
[sam_reporter]
type = MultiAppReporterTransfer
to_reporters = 'max_tcool_r/value'
from_reporters = 'max_tcool/value'
from_multi_app = sam
subapp_index = 0
[]
[]Post-processors, Debug, and Outputs
Postprocessors are used to derive desired quantities from the solution variable(s). Some of the quantities that we may be interested in include the average, maximum, and minimum fuel, clad, and wall temperatures.
[Postprocessors]
[ptot_bison]
type = ElementIntegralVariablePostprocessor
block = pellet
variable = power_density
execute_on = 'INITIAL TIMESTEP_END'
[]
[ptot_scaled_bison]
type = ElementIntegralVariablePostprocessor
block = pellet
variable = power_density_scaled
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tfuel]
type = ElementAverageValue
variable = Temperature
block = pellet
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tfuel]
type = ElementExtremeValue
variable = Temperature
value_type = max
block = pellet
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tfuel]
type = ElementExtremeValue
variable = Temperature
value_type = min
block = pellet
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_twall]
type = ElementAverageValue
variable = twall
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_twall]
type = ElementExtremeValue
variable = Temperature
value_type = min
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_twall]
type = ElementExtremeValue
variable = Temperature
value_type = max
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tclad]
type = ElementExtremeValue
variable = Temperature
value_type = max
block = 'clad 4'
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tclad]
type = ElementExtremeValue
variable = Temperature
value_type = min
block = 'clad 4'
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tclad]
type = ElementAverageValue
variable = Temperature
block = 'clad 4'
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_htc]
type = ElementAverageValue
variable = htc
block = 4
execute_on = 'INITIAL TIMESTEP_END'
[]
[avg_tcool]
type = ElementAverageValue
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_tcool]
type = ElementExtremeValue
value_type = max
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[min_tcool]
type = ElementExtremeValue
value_type = min
variable = tcool
execute_on = 'INITIAL TIMESTEP_END'
[]
[disp_x_max]
type = ElementExtremeValue
variable = disp_x
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[disp_y_max]
type = ElementExtremeValue
variable = disp_y
value_type = max
execute_on = ' INITIAL TIMESTEP_END'
use_displaced_mesh = true
[]
[max_thcond]
type = ElementExtremeValue
variable = thermal_conductivity
[]
[max_tcool_r]
type = Receiver
[]
[]The outputs files are supressed but may be turned on if desired.
[Outputs]
# csv = true
# exodus = true
print_nonlinear_converged_reason = false
print_linear_converged_reason = false
[]SAM Thermal-hydraulic Model
The input file for the SAM thermal hydraulic model of a VTR assembly is displayed below.
################################################################################
## VTR assembly model ##
## SAM Sub-Application ##
## Steady state 1D thermal hydraulics ##
################################################################################
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
m_dot_in = 29.8 #kg/s
Nrods = 217 #number of rods
P_init = 1.0e5 # Pa
T_init = 623.15 # K (=350 degrees C)
rod_od = ${fparse 6.282e-3} # rod outside diameter, m
duct_f2f = 0.11154 # duct inside flat-to-flat, m
ww_pitch = 0.20419 # wire wrap pitch m
ww_od = 1.1e-3 # wire wrap outer diameter m
A_channel = ${fparse duct_f2f*duct_f2f*sqrt(3.)/2-Nrods*pi*rod_od*rod_od/4-Nrods*pi*ww_od*ww_od/4} #0.004048 channel flow area, m2
Ph_channel = ${fparse Nrods*pi*rod_od} # heated perimeter, m
Pw_channel = ${fparse Nrods*pi*rod_od+2*pi*sqrt(3)*duct_f2f} # wetted perimeter, m
D_channel = ${fparse 4* A_channel / Pw_channel} # hydraulic diameter, m
Dheat_channel = ${fparse 4*A_channel/Ph_channel} # heated diameter, m
#fuel height
fuel_h = 0.84235 # m
plenum_full_h = 0.23952 # m
plenum_empty_h = 0.54263 # m
fuel_height = ${fparse fuel_h+plenum_empty_h+plenum_full_h}
P_out = 1.0e5 # 0.1 MPa
rho_in = 868.137 # kg/m^3
V_init = ${fparse m_dot_in/rho_in/A_channel}
# ==============================================================================
# GLOBAL PARAMETERS
# ==============================================================================
[GlobalParams]
eos = eos # The equation-of-state name
scaling_factor_var = '1 1e-3 1e-6' # Scaling factors for fluid variables (p, v, T)
global_init_P = ${P_init}
global_init_T = ${T_init}
global_init_V = ${V_init}
[]
# ==============================================================================
# EQUATIONS OF STATE
# ==============================================================================
[EOS]
[eos] # EOS name
type = PBSodiumEquationOfState # Using the sodium equation-of-state
[]
[]
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Components]
[fuel_channel]
type = PBOneDFluidComponent
position = '0 0 0' # The origin position of this component
orientation = '0 1 0' # The orientation of the component
htc_ext = htc_external # name of HTC external function
Dh = ${D_channel} # Equivalent hydraulic diameter
length = ${fuel_height} # Length of the component
n_elems = 40 # Number of elements used in discretization
A =${A_channel} # Area of the One-D fluid component
Ph = ${Ph_channel} # heated channel perimeter
D_heated = ${Dheat_channel} # Equivalent heated diameter
PoD = 1.183 # rod pitch/diameter ratio
HoD = ${fparse ww_pitch/rod_od} # wire wrap pitch / rod outside diameter
WF_geometry_type = WireWrap # Wall friction geometry type
WF_user_option = ChengTodreas # wall friction correlation
HTC_geometry_type = Bundle # HTC geometry type (rod bunde or pipe)
[]
[inlet]
type = PBTDJ
input = 'fuel_channel(in)' # Name of the connected components and the end type
v_bc = ${V_init} # Velocity boundary condition
T_bc = ${T_init} # Temperature boundary condition
[]
[outlet]
type = PBTDV
input = 'fuel_channel(out) ' # Name of the connected components and the end type
p_bc = ${P_out} # Pressure boundary condition
[]
[heat_transfer_from_fuel]
type = HeatTransferWithExternalHeatStructure
flow_component = fuel_channel
initial_T_wall = 800
htc_name = htc_external
T_wall_name = T_wall_external
[]
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Preconditioning]
[SMP_PJFNK]
type = SMP # Single-Matrix Preconditioner
full = true # Using the full set of couplings among all variables
solve_type = 'PJFNK' # Using Preconditioned JFNK solution method
petsc_options_iname = '-pc_type' # PETSc option, using preconditiong
petsc_options_value = 'lu' # PETSc option, using ‘LU’ precondition type in Krylov solve
[]
[] # End preconditioning block
[Executioner]
type = Steady # This is a transient simulation
nl_rel_tol = 1e-8 # Relative nonlinear tolerance for each Newton solve
nl_abs_tol = 1e-9 # Relative nonlinear tolerance for each Newton solve
nl_max_its = 20 # Number of nonlinear iterations for each Newton solve
l_tol = 1e-4 # Relative linear tolerance for each Krylov solve
l_max_its = 100 # Number of linear iterations for each Krylov solve
petsc_options_iname = '-ksp_gmres_restart' # Additional PETSc settings, name list
petsc_options_value = '100' # Additional PETSc settings, value list
# transient specific
# start_time = 0.0
# end_time = 1
# steady_state_detection = true
# dtmin = 1e-4
# [TimeStepper]
# type = IterationAdaptiveDT
# dt = 0.01
# []
[Quadrature]
type = TRAP # Using trapezoid integration rule
order = FIRST # Order of the quadrature
[]
[] # close Executioner section
# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
[Tin]
type = ComponentBoundaryVariableValue
variable = temperature
input = fuel_channel(in)
execute_on = 'TIMESTEP_END'
[]
[Tout]
type = ComponentBoundaryVariableValue
variable = temperature
input = fuel_channel(out)
execute_on = 'TIMESTEP_END'
[]
[Tavg]
type = ElementAverageValue
variable = temperature
execute_on = 'TIMESTEP_END'
[]
[max_tcool]
type = ElementExtremeValue
value_type = max
variable = temperature
execute_on = 'TIMESTEP_END'
[]
[Pin]
type = ComponentBoundaryVariableValue
variable = pressure
input = fuel_channel(in)
execute_on = 'TIMESTEP_END'
[]
[Pout]
type = ComponentBoundaryVariableValue
variable = pressure
input = fuel_channel(out)
execute_on = 'TIMESTEP_END'
[]
[delta_t_core]
type = DifferencePostprocessor
value1 = Tout
value2 = Tin
execute_on = 'TIMESTEP_END'
[]
[delta_p_core]
type = DifferencePostprocessor
value1 = Pin
value2 = Pout
execute_on = 'TIMESTEP_END'
[]
[Twall_avg]
type = ElementAverageValue
variable = T_wall_external
execute_on = 'TIMESTEP_END'
[]
[Twall_max]
type = ElementExtremeValue
value_type = max
variable = T_wall_external
execute_on = 'TIMESTEP_END'
[]
[Twall_min]
type = ElementExtremeValue
value_type = min
variable = T_wall_external
execute_on = 'TIMESTEP_END'
[]
[htc_avg]
type = ElementAverageValue
variable = htc_external
execute_on = 'TIMESTEP_END'
[]
[htc_min]
type = ElementExtremeValue
value_type = min
variable = htc_external
execute_on = 'TIMESTEP_END'
[]
[htc_max]
type = ElementExtremeValue
value_type = max
variable = htc_external
execute_on = 'TIMESTEP_END'
[]
[]
[Outputs]
# csv = true
print_nonlinear_converged_reason = false
print_linear_converged_reason = false
#[out]
# type = Exodus
# use_displaced = true
# sequence = false
#[]
[]The top of input file houses the model parameters. For a SAM model, this could include the number of rods, rod diameter, channel dimensions, etc.
Global Parameters
Global parameters are parameters that are used in more than one block. This model includes the initial pressure, temperature, and velocity variables of the fluid.
[GlobalParams]
eos = eos # The equation-of-state name
scaling_factor_var = '1 1e-3 1e-6' # Scaling factors for fluid variables (p, v, T)
global_init_P = ${P_init}
global_init_T = ${T_init}
global_init_V = ${V_init}
[]Equations of State
The equations of state (EOS) block is unique to SAM. Because the VTR is sodium cooled, the EOS type is given by PBSodiumEquationofState.
[EOS]
[eos]
# EOS name
type = PBSodiumEquationOfState # Using the sodium equation-of-state
[]
[]Components
This block defines the components of the simulation. In SAM, a component is defined as the physics modeling (fluid flow and heat transfer) and mesh generation of a reactor component.
[Components]
[fuel_channel]
type = PBOneDFluidComponent
position = '0 0 0' # The origin position of this component
orientation = '0 1 0' # The orientation of the component
htc_ext = htc_external # name of HTC external function
Dh = ${D_channel} # Equivalent hydraulic diameter
length = ${fuel_height} # Length of the component
n_elems = 40 # Number of elements used in discretization
A = ${A_channel} # Area of the One-D fluid component
Ph = ${Ph_channel} # heated channel perimeter
D_heated = ${Dheat_channel} # Equivalent heated diameter
PoD = 1.183 # rod pitch/diameter ratio
HoD = '${fparse ww_pitch/rod_od}' # wire wrap pitch / rod outside diameter
WF_geometry_type = WireWrap # Wall friction geometry type
WF_user_option = ChengTodreas # wall friction correlation
HTC_geometry_type = Bundle # HTC geometry type (rod bunde or pipe)
[]
[inlet]
type = PBTDJ
input = 'fuel_channel(in)' # Name of the connected components and the end type
v_bc = ${V_init} # Velocity boundary condition
T_bc = ${T_init} # Temperature boundary condition
[]
[outlet]
type = PBTDV
input = 'fuel_channel(out) ' # Name of the connected components and the end type
p_bc = ${P_out} # Pressure boundary condition
[]
[heat_transfer_from_fuel]
type = HeatTransferWithExternalHeatStructure
flow_component = fuel_channel
initial_T_wall = 800
htc_name = htc_external
T_wall_name = T_wall_external
[]
[]The fuel channel is modelled with a PBOneDFluidComponent type that simulates 1-D fluid flow using the primitive variable based fluid model. The channel geometry is declared with a number of parameters. Additional information can be found in the SAM User's Manual.
The channel inlet and outlet are specified with the PBTDJ and PBTDV types, respectively. The PBTDJ component is an inlet boundary in which the flow velocity and temperature are provided by pre-defined functions. The PBTDV component is a boundary in which pressure and temperature conditions are provided by pre-defined functions.
The last block specifies heat transfer from the fuel with the The HeatTransferWithExternalHeatStructure type. Further information can be found in the SAM User's Manual.
Executioner
The same preconditioner is used as previous BISON input. See Execution Parameters.
[Preconditioning]
[SMP_PJFNK]
type = SMP # Single-Matrix Preconditioner
full = true # Using the full set of couplings among all variables
solve_type = 'PJFNK' # Using Preconditioned JFNK solution method
petsc_options_iname = '-pc_type' # PETSc option, using preconditiong
petsc_options_value = 'lu' # PETSc option, using ‘LU’ precondition type in Krylov solve # End preconditioning block
[]
[][Executioner]
type = Steady # This is a transient simulation
nl_rel_tol = 1e-8 # Relative nonlinear tolerance for each Newton solve
nl_abs_tol = 1e-9 # Relative nonlinear tolerance for each Newton solve
nl_max_its = 20 # Number of nonlinear iterations for each Newton solve
l_tol = 1e-4 # Relative linear tolerance for each Krylov solve
l_max_its = 100 # Number of linear iterations for each Krylov solve
petsc_options_iname = '-ksp_gmres_restart' # Additional PETSc settings, name list
petsc_options_value = '100' # Additional PETSc settings, value list
# transient specific
# start_time = 0.0
# end_time = 1
# steady_state_detection = true
# dtmin = 1e-4
# [TimeStepper]
# type = IterationAdaptiveDT
# dt = 0.01
# []
[Quadrature]
type = TRAP # Using trapezoid integration rule
order = FIRST # Order of the quadrature # close Executioner section
[]
[]Post-processors, Debug, and Outputs
Postprocessors are used to derive desired quantities from the solution variable(s). Some of the quantities that we may be interested in include the inlet and outlet temperatures, heat transfer coeficient, and temperature maximums and minimums.
[Postprocessors]
[Tin]
type = ComponentBoundaryVariableValue
variable = temperature
input = fuel_channel(in)
execute_on = 'TIMESTEP_END'
[]
[Tout]
type = ComponentBoundaryVariableValue
variable = temperature
input = fuel_channel(out)
execute_on = 'TIMESTEP_END'
[]
[Tavg]
type = ElementAverageValue
variable = temperature
execute_on = 'TIMESTEP_END'
[]
[max_tcool]
type = ElementExtremeValue
value_type = max
variable = temperature
execute_on = 'TIMESTEP_END'
[]
[Pin]
type = ComponentBoundaryVariableValue
variable = pressure
input = fuel_channel(in)
execute_on = 'TIMESTEP_END'
[]
[Pout]
type = ComponentBoundaryVariableValue
variable = pressure
input = fuel_channel(out)
execute_on = 'TIMESTEP_END'
[]
[delta_t_core]
type = DifferencePostprocessor
value1 = Tout
value2 = Tin
execute_on = 'TIMESTEP_END'
[]
[delta_p_core]
type = DifferencePostprocessor
value1 = Pin
value2 = Pout
execute_on = 'TIMESTEP_END'
[]
[Twall_avg]
type = ElementAverageValue
variable = T_wall_external
execute_on = 'TIMESTEP_END'
[]
[Twall_max]
type = ElementExtremeValue
value_type = max
variable = T_wall_external
execute_on = 'TIMESTEP_END'
[]
[Twall_min]
type = ElementExtremeValue
value_type = min
variable = T_wall_external
execute_on = 'TIMESTEP_END'
[]
[htc_avg]
type = ElementAverageValue
variable = htc_external
execute_on = 'TIMESTEP_END'
[]
[htc_min]
type = ElementExtremeValue
value_type = min
variable = htc_external
execute_on = 'TIMESTEP_END'
[]
[htc_max]
type = ElementExtremeValue
value_type = max
variable = htc_external
execute_on = 'TIMESTEP_END'
[]
[]BISON Mechanical Model
The BISON mechanical model is constructed in the same way as the thermal model, however, the solver focuses on the mechanical behavior. As such, we'll focus our discussion on input blocks that are different than that of the thermal model and refer the reader to the BISON Thermal Model for more information.
The input file for the BISON mechanical model of a VTR fuel rod is displayed below.
################################################################################
## VTR fuel rod, derived from BISON assessment case IFR1.i ##
## Bison lower level Sub-Application ##
## Fuel and clad mechanical deformation ##
################################################################################
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# Notes
# Units are in standard SI: J, K, kg, m, Pa, s.
# This input file does not include friction.
# radial/axial dimensions at cold condition, will be thermally expanded based on fuel temperature from bison_thermal_only.i
gap = 0. #0.5e-3
rod_outside_diameter = 6.282e-3 # 6.25e-3
clad_thickness = 0.503e-3 # 0.5e-3
slug_diameter = 5.277e-3 # 4.547e-3
fuel_height = 842.4e-3 # 800.0e-3
plenum_height = 782.2e-3 # 778.0e-3
# ==============================================================================
# GLOBAL PARAMETERS
# ==============================================================================
[GlobalParams]
# Parameters that are used in multiple blocks can be included here so that
# they only need to be specified one time.
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y'
temperature = Temperature
#the following are needed in multiple UPuZr Materials
X_Zr = 0.225 # U-20Pu-10Zr
X_Pu = 0.171
[]
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Mesh]
coord_type = RZ
# based on x447.i from examples
# rod specific parameters - dimensions are for fresh fuel at room temperature
[smeared_pellet_mesh]
type = FuelPinMeshGenerator
clad_thickness = ${clad_thickness}
pellet_outer_radius = ${fparse slug_diameter/2}
pellet_height = ${fuel_height}
clad_gap_width = ${gap}
bottom_clad_height = ${clad_thickness}
top_clad_height = ${clad_thickness}
clad_bot_gap_height =${clad_thickness}
clad_top_gap_height = ${plenum_height} # fixme assumes no Na bond sodium
# meshing parameters
clad_mesh_density = customize
pellet_mesh_density = customize
nx_p = 6 # number of fuel elements in radial direction
ny_p = 300 # number of fuel elements in axial direction
nx_c = 3 # number of clad elements in radial direction
ny_c = 300 # number of clad elements in axial direction
ny_cu = 1 #?? number of cladding elements in upper plug in axial direction (default=1)
ny_cl = 1 #?? number of cladding elements in lower plug in axial direction (default=1)
pellet_quantity = 1
elem_type = QUAD4
[]
[add_side_clad]
type = SubdomainBoundingBoxGenerator
location = INSIDE
restricted_subdomains = clad
block_id = '4'
block_name = 'side_clad'
input = smeared_pellet_mesh
bottom_left = '${fparse rod_outside_diameter/2 - clad_thickness/3} 0. 0.'
top_right = '${fparse rod_outside_diameter/2} ${fparse fuel_height + plenum_height} 0.'
[]
# mesh options
patch_size = 50
patch_update_strategy = iteration
partitioner = centroid
centroid_partitioner_direction = y
[]
# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================
[Variables]
[]
[Kernels]
[]
# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
# Only Temperature is specified here. disp_x and disp_y are generated by the
# [Physics/SolidMechanics] block.
[Temperature]
initial_condition = 900
[]
[]
[AuxKernels]
[]
# ==============================================================================
# MODULES
# ==============================================================================
[Physics]
# Modules are prepackaged actions to create objects in other blocks.
[SolidMechanics]
[QuasiStatic]
add_variables = true
strain = SMALL # changed from finite for enabling steady solve
generate_output = 'stress_xx stress_yy strain_xx strain_yy' # 'hoop_stress vonmises_stress hydrostatic_stress volumetric_strain'
[fuel_mechanics]
block = pellet
eigenstrain_names = fuel_thermal_strain # 'fuel_thermal_strain fuel_gaseous_strain fuel_solid_strain'
[]
[cladding_mechanics]
block = 'clad 4'
eigenstrain_names = clad_thermal_strain #'clad_thermal_strain clad_swelling_strain'
#additional_generate_output = 'hoop_creep_strain hoop_elastic_strain hoop_strain'
[]
[]
[]
[]
# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[ICs]
[]
[Functions]
# Any parameter that is a pre-specified function of space and time goes here.
[power_history]
type = ConstantFunction
value = 24.868e3
[]
[axial_peaking_factors]
type = ConstantFunction
value = 1.
[]
[]
# ==============================================================================
# MATERIALS AND USER OBJECTS
# ==============================================================================
[Materials]
#---------------
# mechanics only
#---------------
#fuel+clad
active = 'fuel_elasticity_tensor fuel_elastic_stress fuel_thermal_expansion
clad_elasticity_tensor clad_elastic_stress clad_thermal_expansion
fission_rate fuel_gaseous_swelling'
[fission_rate] # only needed with heat_source_standalone kernel
type = UPuZrFissionRate
rod_linear_power = power_history
axial_power_profile = axial_peaking_factors
#X_Pu_function = 0.0
# coeffs = '0.8952 -1.2801 '
pellet_radius = ${fparse slug_diameter/2}
block = pellet
[]
#fuel
#mechanics materials
[fuel_elasticity_tensor]
type = UPuZrElasticityTensor
block = pellet
porosity = porosity
output_properties = 'youngs_modulus poissons_ratio'
outputs = all
[]
[fuel_elastic_stress]
type = ComputeLinearElasticStress
block = pellet
[]
[fuel_thermal_expansion]
type = UPuZrThermalExpansionEigenstrain
block = pellet
eigenstrain_name = fuel_thermal_strain
stress_free_temperature = 900 # geometry already expanded at hot condition
[]
[fuel_gaseous_swelling] # needed to get porosity material property set
type = UPuZrGaseousEigenstrain
block = pellet
eigenstrain_name = fuel_gaseous_strain
anisotropic_factor = 0.5 # [Greenquist et al., pg. 8]
bubble_number_density = 2.09e18 # [Casagranda, 2020]
interconnection_initiating_porosity = 0.125 # [Casagranda, 2020]
interconnection_terminating_porosity = 0.2185 # [Casagranda, 2020]
fission_rate = fission_rate
output_properties = porosity
outputs = all
[]
#cladding
#mechanics materials
# [clad_elasticity_tensor]
# type = ComputeIsotropicElasticityTensor
# youngs_modulus = 1.645e11 # [Hofman et al., 1989] pg. E.1.1.6 (733 K)
# poissons_ratio = 0.35 # [Hofman et al., 1989] pg. E.1.1.6 (733 K)
# block = 'clad 4'
# []
[clad_elasticity_tensor]
type = HT9ElasticityTensor
block = 'clad 4'
[]
[clad_elastic_stress]
type = ComputeLinearElasticStress
block = 'clad 4'
[]
[clad_thermal_expansion]
type = HT9ThermalExpansionEigenstrain
block = 'clad 4'
eigenstrain_name = clad_thermal_strain
stress_free_temperature = 295 # [Greenquist et al., 2020] pg. 7
[]
[]
# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[BCs]
[no_x_all]
type = DirichletBC
variable = disp_x
boundary = centerline
value = 0.0
[]
[no_y_all]
type = DirichletBC
variable = disp_y
boundary = 'clad_outside_bottom bottom_of_bottom_pellet'
value = 0.0
[]
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Preconditioning]
# Used to improve the solver performance
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Steady
automatic_scaling = true
solve_type = 'PJFNK'
nl_rel_tol = 1e-7
nl_abs_tol = 1e-8
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
petsc_options_value = 'hypre boomeramg 100'
[]
# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
[disp_x_max]
type = NodalExtremeValue
variable = disp_x
value_type = max
use_displaced_mesh = true
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
use_displaced_mesh = true
[]
[disp_y_max]
type = NodalExtremeValue
variable = disp_y
value_type = max
use_displaced_mesh = true
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
use_displaced_mesh = true
[]
[]
[Outputs]
# csv = true
# exodus = true
print_nonlinear_converged_reason = false
print_linear_converged_reason = false
[]Global Parameters
The displacements disp_x and disp_y have been added to the global parameters list, this will make sure every object uses the displaced mesh instead of the original mesh.
[GlobalParams]
# Parameters that are used in multiple blocks can be included here so that
# they only need to be specified one time.
order = FIRST
family = LAGRANGE
displacements = 'disp_x disp_y'
temperature = Temperature
#the following are needed in multiple UPuZr Materials
X_Zr = 0.225 # U-20Pu-10Zr
X_Pu = 0.171
[]Modules
The mechanical model utilizes the solid mechanics module in MOOSE. The model is setup such that the mechanical behavior is modeled in both the fuel and clad. For a detailed explanation of the TensorMechanics model, the authors defer to the SolidMechanics section of the MOOSE documentation.
[Physics]
# Modules are prepackaged actions to create objects in other blocks.
[SolidMechanics]
[QuasiStatic]
add_variables = true
strain = SMALL # changed from finite for enabling steady solve
generate_output = 'stress_xx stress_yy strain_xx strain_yy' # 'hoop_stress vonmises_stress hydrostatic_stress volumetric_strain'
[fuel_mechanics]
block = pellet
eigenstrain_names = fuel_thermal_strain # 'fuel_thermal_strain fuel_gaseous_strain fuel_solid_strain'
[]
[cladding_mechanics]
block = 'clad 4'
eigenstrain_names = clad_thermal_strain #'clad_thermal_strain clad_swelling_strain'
#additional_generate_output = 'hoop_creep_strain hoop_elastic_strain hoop_strain'
[]
[]
[]
[]Materials and UserObjects
Material properties are defined in the [Materials] block. The properties include the elasticity tensor, elastic stress, thermal expansion, and density of the fuel and clad.
[Materials]
#---------------
# mechanics only
#---------------
#fuel+clad
active = 'fuel_elasticity_tensor fuel_elastic_stress fuel_thermal_expansion
clad_elasticity_tensor clad_elastic_stress clad_thermal_expansion
fission_rate fuel_gaseous_swelling'
[fission_rate]
# only needed with heat_source_standalone kernel
type = UPuZrFissionRate
rod_linear_power = power_history
axial_power_profile = axial_peaking_factors
#X_Pu_function = 0.0
# coeffs = '0.8952 -1.2801 '
pellet_radius = '${fparse slug_diameter/2}'
block = pellet
[]
#fuel
#mechanics materials
[fuel_elasticity_tensor]
type = UPuZrElasticityTensor
block = pellet
porosity = porosity
output_properties = 'youngs_modulus poissons_ratio'
outputs = all
[]
[fuel_elastic_stress]
type = ComputeLinearElasticStress
block = pellet
[]
[fuel_thermal_expansion]
type = UPuZrThermalExpansionEigenstrain
block = pellet
eigenstrain_name = fuel_thermal_strain
stress_free_temperature = 900 # geometry already expanded at hot condition
[]
[fuel_gaseous_swelling]
# needed to get porosity material property set
type = UPuZrGaseousEigenstrain
block = pellet
eigenstrain_name = fuel_gaseous_strain
anisotropic_factor = 0.5 # [Greenquist et al., pg. 8]
bubble_number_density = 2.09e18 # [Casagranda, 2020]
interconnection_initiating_porosity = 0.125 # [Casagranda, 2020]
interconnection_terminating_porosity = 0.2185 # [Casagranda, 2020]
fission_rate = fission_rate
output_properties = porosity
outputs = all
[]
#cladding
#mechanics materials
# [clad_elasticity_tensor]
# type = ComputeIsotropicElasticityTensor
# youngs_modulus = 1.645e11 # [Hofman et al., 1989] pg. E.1.1.6 (733 K)
# poissons_ratio = 0.35 # [Hofman et al., 1989] pg. E.1.1.6 (733 K)
# block = 'clad 4'
# []
[clad_elasticity_tensor]
type = HT9ElasticityTensor
block = 'clad 4'
[]
[clad_elastic_stress]
type = ComputeLinearElasticStress
block = 'clad 4'
[]
[clad_thermal_expansion]
type = HT9ThermalExpansionEigenstrain
block = 'clad 4'
eigenstrain_name = clad_thermal_strain
stress_free_temperature = 295 # [Greenquist et al., 2020] pg. 7
[]
[]Boundary Conditions
Dirichlet boundary conditions are specified at the centerline and boundary of the problem for each of the solution variables (disp_x and disp_y). Because the problem is defined in RZ geometry, disp_x refers to the radial displacement and disp_y, the axial displacement.
[BCs]
[no_x_all]
type = DirichletBC
variable = disp_x
boundary = centerline
value = 0.0
[]
[no_y_all]
type = DirichletBC
variable = disp_y
boundary = 'clad_outside_bottom bottom_of_bottom_pellet'
value = 0.0
[]
[]Post-processors, Debug, and Outputs
Postprocessors are used to derive desired quantities from the solution variable(s). Some of the quantities that we may be interested in include the maximum displacement in the radial and axial directions, and the strains.
[Postprocessors]
[disp_x_max]
type = NodalExtremeValue
variable = disp_x
value_type = max
use_displaced_mesh = true
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
use_displaced_mesh = true
[]
[disp_y_max]
type = NodalExtremeValue
variable = disp_y
value_type = max
use_displaced_mesh = true
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
use_displaced_mesh = true
[]
[]Core Support Plate
The input file for the core support plate solid mechanics model is displayed below. The objective of this simulation is to capture the radial thermal expansion of the support plate.
################################################################################
## 3D core support plate thermal expansion input ##
## Tensor Mechanics input model ##
################################################################################
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# given an inlet temperature (nominal = 350 degrees C = 623.15 K),
# computes the displacements along x and y
# transfer them to griffin for cross section adjustment (radial expansion feedback)
Tinlet = 623.15 # inlet temperature (nominal temperature for the )
Tref = 293.15 # reference temperature for the linear thermal expansion for SS316
#Tref = ${Tinlet}
#Tsf = 623.15 # stress-free temperature for the ComputeMeanThermalExpansionFunctionEigenstrain
# ==============================================================================
# GLOBAL PARAMETERS
# ==============================================================================
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
temperature = temp
[]
# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================
[Mesh]
[fmg]
type = FileMeshGenerator
file = mesh/cyl_plate_3d.e # Y - vertical axis
[]
[add_nodeset] # pick node at (0,0,0) and rename it "center"
type = BoundingBoxNodeSetGenerator
input = fmg
new_boundary = center
top_right = '1e-3 0 1e-3'
bottom_left = '-1e-3 0 -1e-3'
[]
[add_outer_nodeset]
type = BoundingBoxNodeSetGenerator
input = add_nodeset
new_boundary = outer_node
top_right = '1.443 0 1e-3'
bottom_left = '1.437 0 -1e-3'
[]
parallel_type = REPLICATED
[]
# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================
[Variables]
[]
[Kernels]
[]
# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
[temp] #core support plate temperature is set to inlet temperature
initial_condition = ${Tinlet}
[]
[]
[AuxKernels]
[]
# ==============================================================================
# MODULES
# ==============================================================================
[Physics/SolidMechanics/QuasiStatic]
[plate]
add_variables = true
strain = SMALL
eigenstrain_names = thermal_expansion
generate_output = 'stress_xx strain_xx stress_yy strain_yy stress_zz strain_zz'
[]
[]
# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[ICs]
[]
[Functions]
# from TEV 3749
[ss316_alphaMean_vtr]
type = ParsedFunction
expression = 'a+b*t+c*t*t'
symbol_names = 'a b c'
symbol_values = '1.789e-5 2.398e-9 3.269e-13'
[]
[]
# ==============================================================================
# MATERIALS AND USER OBJECTS
# ==============================================================================
[Materials]
[elasticity_tensor_ss316]
type = SS316ElasticityTensor
elastic_constants_model = legacy_ifr
[]
[stress]
type = ComputeLinearElasticStress
[]
[thermal_expansion_ss316]
type = SS316ThermalExpansionEigenstrain
stress_free_temperature = ${Tref}
eigenstrain_name = thermal_expansion
[]
[]
[UserObjects]
[]
# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[BCs]
[x_0] # center of support plate bottom surface
type = DirichletBC
variable = disp_x
boundary = center # node at (0,0,0)
value = 0
[]
[xz_0] # outer edge along x-axis (remove rotation motion)
type = DirichletBC
variable = disp_z
boundary = outer_node # node at (1.44,0,0)
value = 0
[]
[z_0] # center of support plate bottom surface
type = DirichletBC
variable = disp_z
boundary = center # node at (0,0,0)
value = 0
[]
[bottom_disp_y] # at bottom of support plate
type = DirichletBC
variable = disp_y
boundary = plateBottom
value = 0
[]
[]
# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Preconditioning]
active = 'SMP_PJFNK'
[SMP_PJFNK]
type = SMP
full = true
solve_type = 'PJFNK'
petsc_options_iname = '-ksp_gmres_restart -pc_type'
petsc_options_value = '100 lu'
[]
[]
[Executioner]
type = Steady
automatic_scaling = true
nl_forced_its = 3
#use default convergence criterion
[]
# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
[disp_x_max]
type = NodalExtremeValue
variable = disp_x
value_type = max
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[disp_y_max]
type = NodalExtremeValue
variable = disp_y
value_type = max
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
[]
[disp_z_max]
type = NodalExtremeValue
variable = disp_z
value_type = max
[]
[strain_zz]
type = ElementAverageValue
variable = strain_zz
[]
[]
[Debug]
[]
[Outputs]
print_nonlinear_converged_reason = false
print_linear_converged_reason = false
# exodus = true
# csv = true
# perf_graph = true
[]The top of input file houses the model parameters. Here, we specify the inlet temperature and reference temperature.
Global Parameters
The displacements that we wish to solve for (disp_x, disp_y, and disp_z) are specified as a global parameter. This is a 3D model of the core support plate and all coordinate directions will be solved for, however, we are only interested in the radial displacements (x and z).
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
temperature = temp
[]Geometry and Mesh
This block defines the computational mesh at which the solution is to be computed on. The mesh is supplied with an Exodus file, titled 'cyl_plate_3d.e' and read in using the FileMeshGenerator type.
Next, the mesh is centered at coordinate position (0,0,0) by defining a nodeset with the BoundingBoxNodeSetGenerator type.
[Mesh]
[fmg]
type = FileMeshGenerator
file = mesh/cyl_plate_3d.e # Y - vertical axis
[]
[add_nodeset]
# pick node at (0,0,0) and rename it "center"
type = BoundingBoxNodeSetGenerator
input = fmg
new_boundary = center
top_right = '1e-3 0 1e-3'
bottom_left = '-1e-3 0 -1e-3'
[]
[add_outer_nodeset]
type = BoundingBoxNodeSetGenerator
input = add_nodeset
new_boundary = outer_node
top_right = '1.443 0 1e-3'
bottom_left = '1.437 0 -1e-3'
[]
parallel_type = REPLICATED
[]AuxVariables and AuxKernels
There is one AuxVariable that is defined for this model; the core support plate temperature. The temperature is set with the inlet temperature.
[AuxVariables]
[temp]
#core support plate temperature is set to inlet temperature
initial_condition = ${Tinlet}
[]
[]There are no AuxKernels.
Modules
The core support plate model utilizes the solid mechanics module that is native to MOOSE. We are interested in the thermal expansion of the core support plate so we set the eigenstrain_names to thermal_expansion and generate a stress and strain output in each coordinate direction. For a detailed explanation of the TensorMechanics model, the authors defer to the SolidMechanics section of the MOOSE documentation.
[Physics]
[SolidMechanics]
[QuasiStatic]
[plate]
add_variables = true
strain = SMALL
eigenstrain_names = thermal_expansion
generate_output = 'stress_xx strain_xx stress_yy strain_yy stress_zz strain_zz'
[]
[]
[]
[]Initial Conditions and Functions
User specified functions are defined in this block. Here, we construct a function for the alpha mean of 316 stainless steel. The type is set with the ParsedFunction parameter and the value equal to a second order polynomial with specified coefficients.
[Functions]
# from TEV 3749
[ss316_alphaMean_vtr]
type = ParsedFunction
expression = 'a+b*t+c*t*t'
symbol_names = 'a b c'
symbol_values = '1.789e-5 2.398e-9 3.269e-13'
[]
[]Materials and UserObjects
Material characteristics are defined in the [Materials] block. Characteristics include the mechanical properties of 316 stainless steel such as the elasticity tensor, stress, and thermal expansion. For a detailed explanation of the material types, the authors defer the reader to the SolidMechanics section of the MOOSE documentation.
[Materials]
[elasticity_tensor_ss316]
type = SS316ElasticityTensor
elastic_constants_model = legacy_ifr
[]
[stress]
type = ComputeLinearElasticStress
[]
[thermal_expansion_ss316]
type = SS316ThermalExpansionEigenstrain
stress_free_temperature = ${Tref}
eigenstrain_name = thermal_expansion
[]
[]Boundary Conditions
Dirichlet boundary conditions are set at the problem boundary for each of the solution variables (disp_x, disp_y, and disp_z). The radial x and z directional boundary is the bottom surface of the core support plate center. The axial y directional boundary is specified at the bottom of the support plate.
[BCs]
[x_0]
# center of support plate bottom surface
type = DirichletBC
variable = disp_x
boundary = center # node at (0,0,0)
value = 0
[]
[xz_0]
# outer edge along x-axis (remove rotation motion)
type = DirichletBC
variable = disp_z
boundary = outer_node # node at (1.44,0,0)
value = 0
[]
[z_0]
# center of support plate bottom surface
type = DirichletBC
variable = disp_z
boundary = center # node at (0,0,0)
value = 0
[]
[bottom_disp_y]
# at bottom of support plate
type = DirichletBC
variable = disp_y
boundary = plateBottom
value = 0
[]
[]Execution Parameters
The same preconditioner is used as previous BISON inputs. See Execution Parameters. For proper convergence the nonlinear iterations are forced to three in the Executioner. Note that without nl_forced_its the default criteria would result in two nonlinear iterations and the solution convergence would not be satisfied.
[Executioner]
type = Steady
automatic_scaling = true
nl_forced_its = 3
#use default convergence criterion
[]Post-processors, Debug, and Outputs
Some of the quantities that we may be interested in include the maximum displacement in the coordinate directions as well as the strain.
[Postprocessors]
[disp_x_max]
type = NodalExtremeValue
variable = disp_x
value_type = max
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[disp_y_max]
type = NodalExtremeValue
variable = disp_y
value_type = max
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
[]
[disp_z_max]
type = NodalExtremeValue
variable = disp_z
value_type = max
[]
[strain_zz]
type = ElementAverageValue
variable = strain_zz
[]
[]