Griffin-BISON Multiphysics Steady State Model

KRUSTY Mesh Model

The KRUSTY mesh employed in the simulations was generated by using MOOSE mesh generator tools including ParsedCurveGenerator, XYDelaunayGenerator, CircularBoundaryCorrectionGenerator, and PeripheralTriangleMeshGenerator, etc. Combinations of these tools facilitate the formation of sophisticated meshes while maintaining their volume integrity. Figure 1 (A) shows the comprehensive core mesh from MOOSE, while a detailed 2D mesh sample from the core's center is provided in Figure 1 (B). An input example of combining ParsedCurveGenerator and XYDelaunayGenerator to generate a 2D mesh of irregular shape (blue mesh in Figure 1 (C), which is part of the 2D mesh of the KRUSTY fuel region (Figure 1 (B))) is shown in Listing 1.

The KRUSTY mesh is 1/4 angular symmetric and could be divided into half or quarter meshes using PlaneDeletionGenerator or CartesianMeshTrimmer objects. Figure 2 shows the quarter mesh generated following these steps. The quarter core mesh was adopted in the simulations.

Figure 1: (A) KRUSTY model; (B) the 2D mesh of KRUSTY fuel region with heat pipes; (C) the enlarged framed area in (B).

[Mesh]
  [Region2]
    type = ParsedCurveGenerator
    x_formula = 't1:=t;
                 t2:=t-1;
                 t3:=t-2;
                 t4:=t-3;
                 x1:=r_hp_gap*cos(t1*(th1_end-th1_start)+th1_start)+hp_x;
                 x2:=r13*cos(t2*(th2_end-th2_start)+th2_start);
                 x3:=r_slot*cos(t3*(th3_end-th3_start)+th3_start)+hp_x;
                 x4:=r12*cos(t4*(th4_end-th4_start)+th4_start);
                 if(t<1,x1,if(t<2,x2,if(t<3,x3,x4)))'
    y_formula = 't1:=t;
                 t2:=t-1;
                 t3:=t-2;
                 t4:=t-3;
                 y1:=r_hp_gap*sin(t1*(th1_end-th1_start)+th1_start);
                 y2:=r13*sin(t2*(th2_end-th2_start)+th2_start);
                 y3:=r_slot*sin(t3*(th3_end-th3_start)+th3_start);
                 y4:=r12*sin(t4*(th4_end-th4_start)+th4_start);
                 if(t<1,y1,if(t<2,y2,if(t<3,y3,y4)))'
    section_bounding_t_values = '0 1 2 3 4'
    constant_names = 'pi              r12           r13            r_slot          r_hp_gap             hp_x       th1_start     th1_end        th2_start      th2_end           th3_start      th3_end        th4_start       th4_end'
    constant_expressions = '${fparse pi} ${fparse r12} ${fparse r13} ${fparse r_slot} ${fparse r_hp_gap} ${fparse hp_x}    ${fparse pi}  1.4983163888   0.12251793865   0.18118584459  1.5707963268     2.5324547478  0.12272974563       0'
    nums_segments = '4 2 4 4'
    is_closed_loop = true
    # show_info = true
  []
[]

Figure 2: Quarter-core KRUSTY mesh generated using MOOSE

KRUSTY Cross Section Generation

KRUSTY is a small fast spectrum core with a significant amount of neutrons leaking out of the fuel region. The energy spectrum shifts due to neutron leakages and thermal neutrons returning from the reflectors. Therefore, the KRUSTY fuel region was divided into seven radial regions and six axial zones as shown in Figure 1 for applying different cross sections for different fuel regions.

Both Serpent and MC2-3 were used to prepare cross sections for modeling KRUSTY. Sensitivity analysis showed that MC2-3 provides better anisotropic scattering cross sections in the BeO reflector regions, and Serpent provides better cross sections in the fuel region for accounting for the fuel Doppler reactivity feedback. Therefore, a hybrid multigroup cross section set that used MC2-3 cross sections in the reflector region and used Serpent cross sections in other regions was prepared for modeling KRUSTY.

The MC2-3 calculation takes a two-step approach. The first step creates energy self-shielded cross sections and the number of energy-groups used is larger than 1000 groups. These cross sections are imported into an approximated RZ model of KRUSTY as shown in Figure 3. A neutron transport calculation was performed by TWODANT to solve the neutron fluxes in each material region in the RZ model. The neutron fluxes calculated by TWODANT have considered both the neutron leakage effects and spatial self-shielding effects and were imported into the MC2-3 second step. Multigroup cross sections were condensed based on the TWODANT fluxes in each material zone. Table 1 lists the 22-group structure used in the calculations. The MC2-3 input for the first step and the input for TWODANT are documented in mc_s1.in. The MC2-3 input for the second step calculation is similar to the input in the first step but has a different control block as shown in Listing 3.

Figure 3: KRUSTY RZ model for TWODANT calculation

\$control
     c_group_structure       =ANL1041
     i_number_region         =64
     c_geometry_type         =mixture
     i_scattering_order      =5
     l_twodant               =T
     c_twodant_group         =BG
     l_pendf                 =T
/

\$control
     c_group_structure       =USER
     i_number_region        =64
     c_geometry_type         =mixture
     c_externalspectrum_ufg  =rzmflx
/

Griffin Stand-Alone Neutronic Model

The standalone Griffin input for the neutronic model of KRUSTY is included. Listing 4 shows the mesh block of this input that was used to generate the unstructured file Griffin_mesh.e. The volume for each material region in the mesh file is preserved and made to be consistent with the KRUSTY Serpent reference model. All heat pipe locations are filled. The list in the subdomains block references block ids in the mesh file, and the extra_element_ids is the list of material ids filling the mesh block.

[Mesh]
 [fmg]
  type = FileMeshGenerator
  file = '../../gold/MESH/Griffin_mesh.e'
 []
 [rot]
  type = TransformGenerator
  input = fmg
  transform = ROTATE
  vector_value = '-90 0 0'
 []
 [id]
   type=SubdomainExtraElementIDGenerator
   input = rot
   subdomains = '1 2 3 4 5 6
                 64 71 72 73
                 8 9 10 11
                 12 1212
                 13 14 15 16 17 18 19
                 20 21 22 23 24
                 2081 2091 2101 2111 2121 2131 2141
                 2082 2092 2102 2112 2122 2132 2142
                 3081 3091 3101 3111 3121 3131 3141
                 3082 3092 3102 3112 3122 3132 3142
                 4081 4091 4101 4111 4121 4131 4141
                 4082 4092 4102 4112 4122 4132 4142'

   extra_element_id_names = 'material_id'
   extra_element_ids ='99 17 11 15 70 18
                       99 6 6 60
                       8 5 10 40
                       51 52
                       4 20 18 40 50 14 13
                       6 7 12 11 9
                       300  310  320  330  340  350  360
                       301  311  321  331  341  351  361
                       202  212  222  232  242  252  262
                       203  213  223  233  243  253  263
                       104  114  124  134  144  154  164
                       105  115  125  135  145  155  165'
  []

  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim= 3
    nx = 20   # about 2 cm may need to adjust
    ny = 20   # about 2 cm may need to adjust
    nz = 40   # about 5 cm may need to adjust to test convergence or performance
    xmin = -2
    xmax =  52
    ymin = -2
    ymax =  52
    zmin = -1
    zmax =  135
  []
   [assign_coarse_id]
     type = CoarseMeshExtraElementIDGenerator
     input = id
     coarse_mesh = coarse_mesh
     extra_element_id_name = coarse_element_id
   []
[]

[AuxVariables]
  [Tf]
    initial_condition = 300
    order = CONSTANT
    family = MONOMIAL
  []

  [Ts]
      initial_condition = 300
      order = CONSTANT
      family = MONOMIAL
  []
[]
[TransportSystems]
  particle = neutron
  equation_type = eigenvalue
  G =22
  VacuumBoundary = 'Core_top Core_bottom Core_outer_boundary'
  ReflectingBoundary = 'Mirror_X_surf Mirror_Y_surf'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 3
    NAzmthl = 5

    NA = 3
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering  = true
    n_delay_groups = 6
  []
[]

[PowerDensity]
  power = 1275.0
  power_density_variable = power_density
[]

[Executioner]
  type = SweepUpdate

  richardson_abs_tol = 1e-6
  richardson_rel_tol = 1e-20

  richardson_max_its = 1000
  inner_solve_type = GMRes # SI #GMRes

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  max_diffusion_coefficient = 10.0
  diffusion_eigen_solver_type = newton
  diffusion_prec_type = lu
  prolongation_type = multiplicative

[]

[Materials]
 [all]
  type=CoupledFeedbackMatIDNeutronicsMaterial
  block='1 2 3 4 5 6
         64 71 72 73
         8 9 10 11
         12 1212
         13 14 15 16 17 18 19
         20 21 22 23 24
         2081 2091 2101 2111 2121 2131 2141
         2082 2092 2102 2112 2122 2132 2142
         3081 3091 3101 3111 3121 3131 3141
         3082 3092 3102 3112 3122 3132 3142
         4081 4091 4101 4111 4121 4131 4141
         4082 4092 4102 4112 4122 4132 4142'
  library_file = '../Serp_hbrid_reflector_updated.xml'
  library_name ='krusty_serpent_ANL_endf70_g22'
  isotopes = 'pseudo'
  densities = '1.0'
  grid_names = 'Tfuel Tsteel'
  grid_variables = 'Tf Ts'
  plus = 1
  is_meter = true
  dbgmat = false
 []
[]

[Outputs]
  csv = true
  exodus = true
  [pgraph]
    type = PerfGraphOutput
    level = 2
  []
[]

The KRUSTY mesh file has modeled 1/4th of the core with reflecting boundaries applied at the center-cut planes. Vacuum boundary conditions are applied on the KRUSTY top, bottom surface, as well as on its outermost radial surface as shown in Listing 6.

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue
  G =22
  VacuumBoundary = 'Core_top Core_bottom Core_outer_boundary'
  ReflectingBoundary = 'Mirror_X_surf Mirror_Y_surf'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 3
    NAzmthl = 5

    NA = 3
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering  = true
    n_delay_groups = 6
  []
[]

The Griffin neutronic model solves the and the flux distributions at critical state. Listing 7 shows the transport system block. It used the DFEM-SN solver. Spatial variables are discretized by the DFEM method with shape functions from the first order MONOMIAL family. Angular variables are discretized by the SN method using the Gauss-Chebyshev quadratures. Sensitivity analysis were performed and results were included in Table 2 with different orders of anisotropic scattering terms and different number of angles used in discretizing the angular space. It showed that the calculated was converged after using NPolar=3 angles per octant in the polar direction and NAzmthl=5 angles per octant in the azimuthal direction. With scattering term NA=3, Griffin neutronic model provides better agreement to the Serpent reference among these cases. Including higher scattering order terms may further improve the accuracy of the neutronic results but will require a lot more computational resources. In the Multiphysics model, the neutronic model used NA=3, NPolar=2 and NAzmthl=3 in order to further reduce computational time. As shown in Listing 7, six groups of delayed neutron group families were considered in the numerical calculations. Six delayed neutron data were obtained using Serpent and the values for each group were included in the ISOXML cross section file.

Listing 8 shows the executioner block in the Griffin input where the maximum number of allowed iterations and convergence criteria are specified. The SweepUpdate executioner unique for DFEM-SN was used. The Coarse-Mesh Finite Difference (CMFD) method was deployed to accelerate the DFEM-SN solver. The coarse mesh for the low-order diffusion calculation was a regular mesh that was generated in the mesh block as shown in Listing 5. The maximum diffusion coefficient was set to 10.0, as a limit value facilitating convergence in problems that contain void regions. It can be adjusted and will not significantly affect the accuracy of the results. Newton’s method was used to solve the low-order diffusion equation. Another option to use is the krylovshur method. It was found that it is important to use preconditioners for solving the diffusion equation in order to converge the DFEM-SN solver with CMFD acceleration. In this case, the LU decomposition was used and was found to be effective. The multiplicative prolongation was chosen to update all transport flux moments after solving the low-order diffusion calculation.

[Executioner]
  type = SweepUpdate

  richardson_abs_tol = 1e-6
  richardson_rel_tol = 1e-20

  richardson_max_its = 1000
  inner_solve_type = GMRes # SI #GMRes

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  max_diffusion_coefficient = 10.0
  diffusion_eigen_solver_type = newton
  diffusion_prec_type = lu
  prolongation_type = multiplicative
[]

[Materials]
  [all]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '1 2 3 4 5 6
         64 71 72 73
         8 9 10 11
         12 1212
         13 14 15 16 17 18 19
         20 21 22 23 24
         2081 2091 2101 2111 2121 2131 2141
         2082 2092 2102 2112 2122 2132 2142
         3081 3091 3101 3111 3121 3131 3141
         3082 3092 3102 3112 3122 3132 3142
         4081 4091 4101 4111 4121 4131 4141
         4082 4092 4102 4112 4122 4132 4142'
    library_file = '../Serp_hbrid_reflector_updated.xml'
    library_name = 'krusty_serpent_ANL_endf70_g22'
    isotopes = 'pseudo'
    densities = '1.0'
    grid_names = 'Tfuel Tsteel'
    grid_variables = 'Tf Ts'
    plus = 1
    is_meter = true
    dbgmat = false
  []
[]

[AuxVariables]
  [Tf]
    initial_condition = 300
    order = CONSTANT
    family = MONOMIAL
  []

  [Ts]
    initial_condition = 300
    order = CONSTANT
    family = MONOMIAL
  []
[]

Listing 9 shows the material block in the Griffin input. The CoupledFeedbackMatIDNeutronicsMaterial, which interpolates multigroup cross sections on-the-fly among a tabulated cross section set, was used for multiphysics simulation. The material ID was read from the mesh block. The hybrid cross section set that was generated from Serpent and MC2-3 was used. The cross sections within the file are zone-averaged macroscopic cross sections. Therefore, the isotopes were assigned to be ‘pseudo’ and its density was 1.0.

The cross section set included a 2-D cross section table. The indices to lookup the 2-D table were the fuel temperature “Tf” and the structure temperature “Ts” that covers all other nonfuel regions. The state temperature points include fuel at 300 K, 600 K, 800 K, 1000 K and 1200 K and the structure material at 300 K, 600 K, 800 K and 1000 K, respectively. Listing 10 shows how different fuel and structure temperatures can be specified in the stand-alone neutronic calculations.

Finally, Listing 11 shows the steps to tally the axial power distribution in the KRUSTY innermost (“r1”) layer. To integrate the fission power within each fuel layer (seven radial layers), an AuxVariable in MOOSE that can be processed by the VectorPostprocessor was first defined to hold the layer-integrated total fission power for the axial regions within that layer as shown in the first block of Listing 11. A UserObject as shown in the second block of Listing 11 was necessary to perform the actual integration among the different axial layers. The integrated values were then placed into the AuxVariable as shown in the third block of Listing 11 and were printed out by the VectorPostprocessor as shown in the fourth block of Listing 11.

[AuxVariables]
  [axial_r1_z]
    order = FIRST
    family = MONOMIAL
  []
[]
[UserObjects]
  [axial_r1_z]
    type = LayeredIntegral
    variable = power_density
    direction = z
    bounds = '35 36 37 38 39 40.083435 41.16687 42.250305 43.33374 43.97763
              44.62152 44.85012 45.11512 45.38012 46.19006 47 48.09006 49.18012
              49.71012 49.95142 50.80945 51.66748 52.50061 53.33374 54.16687
              55 55.8416779 56.6833558 57.5250337 58.3667116 59.1839658 60.00122'
    block = '2081 2082 3081 3082 4081 4082'
  []
[]
[AuxKernels]
  [axial_r1_z]
    type = SpatialUserObjectAux
    variable = axial_r1_z
    execute_on = timestep_end
    user_object = axial_r1_z
  []
[]
[VectorPostprocessors]
  [axial_r1_z_elm]
    type = ElementValueSampler
    variable = axial_r1_z
    sort_by = z
    execute_on = timestep_end
  []
[]

BISON Thermal-Mechanical Model

BISON was utilized to conduct thermo-mechanical simulations. In KRUSTY, the fuel disk and its radial shield hung down from the top stainless-steel plates. The bottom and radial reflectors sat on a movable platen below the fuel disk and were pushed up into the radial shield region to reach criticality. With very good insulation of the fuel disk, the top and bottom plates were about room temperature. Therefore, the reactor core was modeled between the fixed top and bottom surfaces with no total displacement in the axial direction but allowing radial thermal expansion or contraction. The top and bottom parts of KRUSTY were connected by “soft materials” to represent the vacuum and air region surrounding the fuel disk. In the KRUSTY thermo-mechanical model, the radial reflector expands upward into the “soft material” region when the fission power increases, and the fuel expansion is downward resulting in more part of fuel covered by the reflectors. Convective heat flux boundary conditions (h=5 W/m²·K) were applied to the top, bottom, and side surfaces of the core.

For a thermo-mechanical model, thermophysical properties (i.e., thermal conductivity, density, specific heat, and thermal expansion coefficient) and mechanical properties (i.e., elastic modulus) of all the involved components are needed. Table 3 lists the BISON/MOOSE models used for thermo-mechanical simulations. BISON doesn't include models for KRUSTY’s U-7.65Mo fuel material, so the thermal conductivity of U-10Mo and specific heat of U-8Mo are currently applied. Other thermal and mechanical attributes for U-7.65Mo and B4C were sourced from the open literature or based on certain assumptions. The displacements in the “soft materials” regions are not solved using MOOSE’s solid mechanics solver as the solid components. Instead, simple diffusion is used for these regions to simply propagate the deformation of the solid components.

As mentioned above, in the KRUSTY design, a thin insulation layer is used to thermally isolate the fuel from other components to reduce the loss of energy through heat transfer from fuel to reactor external surface. In order to simulate the thermal behavior of this thin insulation layer, MOOSE’s interface kernel is adopted to model the thermal resistance of the insulation layer without the need to explicitly mesh the thin structure. The thermal resistance of the insulation layer is tuned to match the reported behavior of KRUSTY (i.e., 30 W power leading to 200ºC fuel temperature).

MOOSE Multiphysics Coupling

In KRUSTY warm critical experiments, the power excursion transient was initiated by shifting the radial reflector of KRUSTY upward to insert positive reactivity when the reactor was in a cold/critical state (Poston et al., 2020). This VTB model includes the multiphysics model prepared for modeling this reactivity insertion test. A two-level MOOSE MultiApps hierarchy which tightly couples the aforementioned Griffin neutronics model and BISON thermo-mechanical model has been developed to simulate the KRUSTY warm critical tests. Griffin is set as the parent (main) application and uses the DFEM-SN(2,3) solver with CMFD acceleration and NA=3 (anisotropic scattering order), while BISON is set as the child application.

The power density profile, initially calculated by Griffin, is then transferred to BISON. In BISON, thermo-mechanical computations take place, allowing for the determination of temperature distribution within all solid components. This subsequently leads to the calculation of thermal expansion. Both the fuel temperature profile and displacement field are then sent back to Griffin as feedback for the neutronics simulation. The coupling of these two applications occurs through fixed point iteration. Notably, in this calculation, all heat is passively removed through the external boundary of the reactor, as no heat pipes are involved.

The movement of the axial reflector to insert the reactivity was modeled by forcing a Dirichlet boundary condition on the bottom of the solid assembly that includes the axial reflector to allow BISON to calculate the consequent displacement field through solving solid mechanics equations. Two steady-state eigenvalue calculations were included in the VTB model. These calculations were performed to confirm that an adequate amount of external reactivity was introduced for initiating the transient. The first one corresponds to the initial steady state for future transient simulation, and the second steady state corresponds to the first step of the transient after the reflector is shifted upward. In this calculation, the bias in the axial direction was set in the BISON input as shown in Listing 12 to generate the displacement in the mesh. The displacements were passed to Griffin neutronic for calculating the total reactivity inserted into the reactor using the MultiApp coupling framework.

reflector_disp = 0.0 # use 1.48e-3 for reactivity insertion
[BCs]
  [BottomSSFixZ]
    type = DirichletBC
    variable = disp_z
    boundary = 'ss_bot'
    value = ${reflector_disp}
  []
[]

References