Hexagonal Duct Bowing - Linear Thermal Gradient (IAEA VP1)

High Level Summary of Model

This model examines the free (unrestrained) thermal bowed deformation of a single, 3D hexagonal assembly (duct) which is subject to thermal gradients along the radial and axial directions. The duct is mechanically fixed at the bottom. The model specification is taken directly from Verification Problem 1 (VP1) in the series of benchmark problems published in IAEA (1990) to support the verification and validation of Liquid Metal Fast Breeder Reactor (LMFBR) analysis codes organized by The International Atomic Energy Agency (IAEA)’s International Working Group on Fast Reactors (IWGFR). Eleven participating agencies in nine different countries contributed numerical solutions to the set of benchmarks. This model employs the MOOSE Tensor Mechanics Module to model the thermo-mechanical deformation response due to the thermal gradients in the duct. This type of deformation is present in liquid-metal cooled fast reactors and serves as an important reactivity feedback effect. The results obtained from this MOOSE-based model demonstrate proper set up of a single assembly thermal duct bowing problem which can be used to build larger and more complex models.

Computational Model Description

In this example, known as IAEA Verification Problem 1 (VP1), a known linearly varying thermal gradient is applied to the hexagonal duct along the axial height in the active core region, as shown in Figure 1. Below the active core, the duct is at 400°C which varies linearly from the core inlet to 550°C at the left corner midwall and 500°C at the far right corner midwall. Along the cross-section, the temperature follows a linear gradient at each axial location. Above the active core, the temperatures are held constant from the core outlet to the top of the duct. The cross-sectional thermal gradient produces a differential thermal expansion where the hotter side of the duct expands more than the cooler side, causing a bowed displacement from hot side towards cold side.

Figure 1: Representation of thermal gradient along duct axial location, from 400°C below the core extending to 500°C and 550°C on opposite corners of the duct.

Figure 2 shows the geometry of the duct in both elevation and cross-section view, with the primary geometric parameters labeled. Table 1 summarizes the values for the geometric parameters. The duct is 4000 mm long, with the active core region extending from 1500 mm along the axial direction to 2500 mm. There is an above core load pad (ACLP) located at 3000 mm axial height. The duct is a sharp cornered hexagon with a distance between the mid-walls of the corners of 150 mm and a wall thickness of 3 mm at room temperature, 20°C. The ACLP has a thickness on top of the duct wall of 2.75 mm and an arbitrary height of 100 mm. The duct and load pad are made of a fictitious steel material with constant material properties (not dependent on temperature) which are given in Table 2.

Figure 2: Schematic showing duct model geometric parameters.

The geometry and phenomenon in this model are relevant to core-bowing phenomena present in liquid-metal cooled fast reactor design and analysis. The core bowing present in reactor assemblies induces reactivity feedback effects due to sensitivity of fast neutron reactors to leakage pathways. The MOOSE Reactor Module is used to generate the duct mesh, and the MOOSE Tensor Mechanics module is used to calculate the thermo-mechanical response of the duct due to differential thermal expansion. This work was originally performed as part of an assessment of the MOOSE Tensor Mechanics module for calculating radial core expansion for fast reactors Wozniak et al. (2021) funded under the NEAMS program. The temperature along the cold () corner of the duct is defined by Eq. (1):

(1)

The temperature along the hot () corner of the duct is defined by Eq. (2):

(2)

The bending moment along the duct axial height can be represented by the moment in a beam induced by thermal gradients given by the following equation:

(3)

Where is the Young’s modulus, is the area moment of inertia, is the thermal expansion coefficient per , and are the temperatures evaluated at the hot and cold corners of the cross section, is the depth of the beam along the cross section in the direction of the bending deflection, and is the axial height above the fixed boundary condition (located ). It should be noted that the moment is also defined as a piecewise continuous function along the same values as the temperature. The deflection of the duct, , can then be found as the deflection of a cantilever beam. Based on beam theory Beer et al. (2012), the second derivative of the deflection multiplied by and is equal to bending moment of the beam, which is defined in Eq. (4):

(4)

The deflection of the beam is then found through double integration of the moment and applying continuity boundary equations to the piecewise continuous sections:

(5)

Mesh Block

The mesh was created with the MOOSE Reactor Module. First, a 2D duct cross-section is created with PolygonConcentricCircleMeshGenerator, followed by removal of the duct interior by BlockDeletionGenerator to create a thin-walled hexagonal duct. The 2D duct is extruded to 3D using MeshExtruderGenerator. We note that while the presence of load pads is not important physically for this particular problem as no contact occurs, load pads are modeled so that this mesh can be re-used in more complex multi-assembly models. Therefore, the same steps are taken to create the load pad, which is moved axially along the duct to the correct position using TransformGenerator and stitched to the duct with StitchedMeshGenerator. Sidesets are renamed using RenameBoundaryGenerator, which are useful to name the duct and load pad faces for creating contact sets or for post-processing. The sideset on the bottom of the duct is renamed to “fixed” so the boundary conditions can be applied by easily identifying the sideset by name. The mesh of the top-down view and the duct at the load pad is shown in Figure 3.

[Mesh]
  [cladding_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '8 8 8 8 8 8'
    background_intervals = 2
    background_block_ids = '10 10'
    polygon_size = 0.066451905
    polygon_size_style = 'apothem'
    duct_sizes_style = 'apothem'
    duct_sizes = '0.063451905'
    duct_intervals = '1'
    duct_block_ids = '1500'
    preserve_volumes = on
    quad_center_elements = true
    external_boundary_id = 100
    interface_boundary_id_shift = 100
    # To get each side of the pin
    generate_side_specific_boundaries = true
  []

  [cladding_center_removal]
    type = BlockDeletionGenerator
    block = '10'
    input = cladding_hex
  []

  [rename_cladding_1]
    # to avoid conflicts when generating load pad hex frame
    type = RenameBoundaryGenerator
    input = cladding_center_removal
    old_boundary = '10000 10001 15001 10002 15002 10003 15003 10004 15004 10005 15005 10006 15006'
    new_boundary = '1000 1001 1001 1002 1002 1003 1003 1004 1004 1005 1005 1006 1006'
  []

  [cladding_extrude]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 4.0'
    num_layers = 400
    bottom_sideset = '1201' # will be fixed boundary
    top_sideset = '1202'
    input = rename_cladding_1
  []

  ##################
  [load_pad_hex]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '8 8 8 8 8 8'
    background_intervals = 2
    background_block_ids = '10 10'
    polygon_size = 0.069201905
    polygon_size_style = 'apothem'
    duct_sizes_style = 'apothem'
    duct_sizes = '0.066451905'
    duct_intervals = '1'
    duct_block_ids = '1600'
    preserve_volumes = on
    quad_center_elements = true
    external_boundary_id = 200
    interface_boundary_id_shift = 200
    generate_side_specific_boundaries = true
  []

  [load_pad_center_removal]
    type = BlockDeletionGenerator
    block = '10'
    input = load_pad_hex
    new_boundary = '2200' # load pad inner surface to stitch
  []

  [load_pad_extrude]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 0.1'
    num_layers = 10
    input = load_pad_center_removal
  []

  [load_pad_translate]
    type = TransformGenerator
    input = load_pad_extrude
    transform = translate
    vector_value = '0 0 2.95'
  []

  [Cladding_load_pad_stitching]
    type = StitchedMeshGenerator
    inputs = 'cladding_extrude  load_pad_translate'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = '1000 2200'
  []

  [sidesets_rename_all]
    type = RenameBoundaryGenerator
    input = Cladding_load_pad_stitching
    old_boundary = '1201 1001 1002 1003 1004 1005 1006 10001 15001 10002 15002 10003 15003 10004 15004 10005 15005 10006 15006'
    new_boundary = 'fixed face3 face2 face1 face6 face5 face4 face3_aclp face3_aclp face2_aclp face2_aclp face1_aclp face1_aclp face6_aclp face6_aclp face5_aclp face5_aclp face4_aclp face4_aclp'
  []
  [block_rename]
    type = RenameBlockGenerator
    input = sidesets_rename_all
    old_block = '1500 1600'
    new_block = '1 1'
  []

  # Mesh parallel partitioning parameters
  patch_update_strategy = auto
  patch_size = 20
  partitioner = centroid
  centroid_partitioner_direction = z
[]

Figure 3: Single duct mesh top down and load pad location.

Variables

The TensorMechanics module solves for the displacement variables in , , and directions. These represent the components of the vector displacement at each node in the duct.

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
[]

Thermal Gradient Function

The applied temperature along the duct is defined with a MOOSE function, temp_func, of type ParsedFunction which specifies linear variation of temperature along the active core region in both the radial and axial directions as in Eq. (1) and Eq. (2).

[Functions]
  # The duct temperatures are defined at the corners and linearly vary in the axial direction
  # and along the face of the duct.
  [temp_func]
    type = ParsedFunction
    #At center of wall, y=+-0.075m
    #T varies across the cross-section from 500C to 550C, ramps up to that from 400C at z=1.5m to 2.5m
    expression = '400+if(z>2.5,t*(125-25/.075*y*t),if(z>1.5,t*(z-1.5)/1.0*(125-25/.075*y*t),0))'
  []
[]

The MOOSE AuxVariable, temp, is defined as the placeholder for temperature with the initial condition of 400°C, which is the temperature of the duct below the active core region. The Tensor Mechanics module will use this variable to calculate the duct thermal expansion.

[AuxVariables]
  [temp]
    initial_condition = 400
    order = FIRST
    family = LAGRANGE
  []
[]

Since the duct temperature is space-dependent rather than a constant value, the MOOSE AuxKernel of type FunctionAux instructs MOOSE to retrieve the temperature from the temp_func Function and stores it in the AuxVariable, temp.

[AuxKernels]
  [tfunc]
    type = FunctionAux
    variable = temp
    function = temp_func
    block = '1'
  []
[]

Boundary Conditions

The bottom surface of the duct is fixed from translation in all directions (, , and displacements) to simulate a cantilever beam bending. This is achieved through the assignment of Dirichlet boundary conditions in the BCs block. The variable parameter represents the variable on which this BC operates, and the value parameter is the value of the variable at the specified boundary. For this model, a “fixed” boundary sideset was named to easily reference the bottom of the duct to assign the boundary condition.

[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = fixed
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = fixed
    value = 0.0
  []
  [no_z]
    type = DirichletBC
    variable = disp_z
    boundary = fixed
    value = 0.0
  []
[]

Tensor Mechanics Module

The Tensor Mechanics Module is a library of code modules which models the mechanical deformation and stress of solids. The basic set of calculations required for stress calculation problems include computing the strain, the elasticity/constitutive relationship, and the stress tensors. In this model, finite strain formulation is used to account for large relative deformations in the thin wall cross-section. For thermal expansion calculations, an eigenstrain, which is a mechanical deformation not caused by external stresses, is defined as thermal_expansion.

[Modules]
  [TensorMechanics]
    [Master]
      [all]
        strain = FINITE
        volumetric_locking_correction = true
        add_variables = true
        eigenstrain_names = thermal_expansion
        decomposition_method = EigenSolution
        generate_output = 'vonmises_stress'
        temperature = temp
        use_finite_deform_jacobian = true
        extra_vector_tags = 'ref'
      []
    []
  []
[]

To compute the stresses and relate them with the strains, the elasticity (constitutive relationship) tensor is required. The constitutive relationship between the stress and the strain is shown in Eq. (6):

(6)

Where is the elasticity tensor, and is the incremental strain. Because the material is defined as steel, which is considered an isotropic material, the elasticity tensor is computed with ComputeIsotropicElasticityTensor material property, which requires the Young’s Modulus (modulus of elasticity) and Poisson’s ratio. The stress tensor is also required, which is calculated with ComputeFiniteStrainElasticStress, because the finite strain formulation is defined in the TensorMechanics block. To compute the thermal expansion due to the temperature, temp, the ComputeThermalExpansionEigenstrain property is added to compute the eigenstrain which is defined in the TensorMechanics block. Then the total strain in the duct is the summation of the elastic strain and the eigenstrain.

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.7e11
    poissons_ratio = 0.3
    block = 1
  []
  [small_stress]
    type = ComputeFiniteStrainElasticStress
    block = 1
  []
  [thermal_expansion_strain]
    type = ComputeThermalExpansionEigenstrain
    stress_free_temperature = 400.0
    thermal_expansion_coeff = 18.0e-6
    temperature = temp
    eigenstrain_name = thermal_expansion
    block = '1'
  []
[]

Results

As a debugging check, the temperature calculated with the defined function is plotted in Figure 4 to ensure the temperatures have been correctly mapped to the duct. We can visually verify that the temperature below the active core region is 400°C, the temperature increases through the active core region, and above the active core the temperature on each side of the duct is held constant. The deformation output from the MOOSE simulation is also plotted in Figure 4, magnified to show the deformed shape. The duct visibly bows from the high temperature side (left) towards the cooler temperature side (right).

Figure 4: Temperature gradient mapping of the duct (left) and duct bowing (right, magnified).

A closed form solution is available for the expected deformation of the duct centerline as shown in Eq. (5). Since MOOSE models only the duct walls, the centerline deformation of the duct is estimated by averaging the deformation of the flat faces at each axial height. The deformation from the closed-form solution and the MOOSE output are plotted in Figure 5 and agree well. The IAEA report IAEA (1990) compares the different participant codes and averages their results to compare with the closed form solution at the core top, the ACLP, and the top of the duct. Those values are presented in Table 3, comparing the closed form, the IAEA average, and the MOOSE results. The MOOSE simulation results agree well with both the code averages and the closed-form solution, with a deformation at the core top equal to 1.00 mm, deformation at the ACLP equal to 3.25 mm, and deformation at the top of the duct equal to 12.26 mm.

Figure 5: Duct bowing comparison between MOOSE and the closed-form solution.

Run Command

This model requires the openly available Reactor and Tensor Mechanics Modules which can be independently compiled from a MOOSE framework installation or leveraged from a MOOSE application which links to a version of MOOSE containing both these modules.

mpirun -np 20 /path/to/APP  -i hex_duct_linear_gradient.i

References