Seismic analysis of a base-isolated nuclear power plant building

This model was adopted from the list of examples on the MASTODON website. The inputs can be found here.

commentnote:Units of this model

GN, GPa, m, and sec

Part 2: NPP building analysis

Now that the seismic response with just the basemat is shown to be reasonable, the modeling and response of the seismically isolated building is presented in this section.

Modeling

The finite-element mesh of the building, developed in Cubit is presented in Figure 3 in Part 1. The building is founded on the basemat and the isolation system presented in Figure 2 in Part 1. The figures below show the plan view and an isometric view of the internal components of the building. The building is a one-story shear wall structure with buttresses on all four sides and the roof. It houses four steam generators (shown in purple and modeled as solid cylinders for simplicity) supported by individual concrete bases (in yellow), and head-supported reactor vessel is suspended by an internal concrete slab. The reactor vessel is contained in a concrete cylindrical housing that can be seen in red the figures below. The reactor vessel contains molten salt with properties that are assumed to be the same as water. The reactor vessel, as seen from the bottom without its housing or the internal walls, is shown in Figure 1 below.

Figure 2: Isometric view of the internal components of the NPP building.

Figure 3: Plan view of the internal components of the NPP building.

Figure 1: Head-supported reactor vessel as seen from the bottom.

All the materials in the building are modeled using a linear elastic material and 3D 8-noded HEX elements, except for the FP isolators, which are modeled using two noded link elements with the FP isolator material model see theory manual and user manual for more information. The reactor vessel is modeled as a thin cylindrical vessel using continuum elements. The fluid inside the reactor vessel is modeled using a linear elastic material with a very small shear modulus to reproduce fluid-like behavior. A more comprehensive, fluid-structure interaction (FSI) approach is also possible in MASTODON to model the fluid. More information on FSI modeling in MOOSE and MASTODON is provided here. Further information on the building model is also provided in Bolisetti et al. (2020).

In this demonstration, the building is subjected to ground motions in X, Y, and Z directions at the base of the isolation system. These ground motions are presented in Figure 1, and their spectral accelerations are shown in Figure 6, Figure 8, and Figure 10, all from Part 1. Pressure, temperature, and velocity dependencies of the isolators are switched on, as described in method 2 of Part 1 of this model. The first three timesteps of the analysis involve a static gravity analysis as described in the note above.

The input file for the simulation of Part 2 is listed below.

[Mesh]
  [mesh_gen]
    type = FileMeshGenerator
    file = mesh/full_structure_with_isolators_new.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [rot_x]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_y]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_z]
    block = 'isolator_elems upper_rigid_elems'
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [vel_y]
  []
  [accel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [rot_vel_x]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_vel_y]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_vel_z]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_accel_x]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_accel_y]
    block = 'isolator_elems upper_rigid_elems'
  []
  [rot_accel_z]
    block = 'isolator_elems upper_rigid_elems'
  []
  [Fb_x]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
  [Fb_y]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
  [Fb_z]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
  [Defb_x]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
  [Velb_x]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
  [Defb_y]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
  [Defb_z]
    block = 'isolator_elems'
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/LineElement/QuasiStatic]
  displacements = 'disp_x disp_y disp_z'
  rotations = 'rot_x rot_y rot_z'

  velocities = 'vel_x vel_y vel_z'
  accelerations = 'accel_x accel_y accel_z'
  rotational_velocities = 'rot_vel_x rot_vel_y rot_vel_z'
  rotational_accelerations = 'rot_accel_x rot_accel_y rot_accel_z'

  beta = 0.275625
  gamma = 0.55
  alpha = -0.05

  [rigid_beams]
    block = 'upper_rigid_elems'
    area = 130.06
    Iy = 24166.729
    Iz = 24166.729
    y_orientation = '0.0 0.0 1.0'
  []
[]

[Physics/SolidMechanics/Dynamic]
  displacements = 'disp_x disp_y disp_z'
  [all]
    strain = FINITE
    displacements = 'disp_x disp_y disp_z'
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases SGs int_wall int_slab RV_housing RV small_walls upper_basemat fluid_material RV_slab'
    hht_alpha = -0.05
    static_initialization = true
    stiffness_damping_coefficient = 0.0019
  []
[]

[Kernels]
  [inertia_x]
    type = InertialForce
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases SGs int_wall int_slab RV_housing RV small_walls upper_basemat fluid_material RV_slab'
    variable = disp_x
    eta = 0.038
    alpha = -0.05
  []
  [inertia_y]
    type = InertialForce
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases SGs int_wall int_slab RV_housing RV small_walls upper_basemat fluid_material RV_slab'
    variable = disp_y
    eta = 0.038
    alpha = -0.05
  []
  [inertia_z]
    type = InertialForce
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases SGs int_wall int_slab RV_housing RV small_walls upper_basemat fluid_material RV_slab'
    variable = disp_z
    eta = 0.038
    alpha = -0.05
  []
  [lr_disp_x]
    type = StressDivergenceIsolator
    block = 'isolator_elems'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    component = 0
    variable = disp_x
    static_initialization = true
    zeta = 0.0019
    alpha = -0.05
  []
  [lr_disp_y]
    type = StressDivergenceIsolator
    block = 'isolator_elems'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    component = 1
    variable = disp_y
    static_initialization = true
    zeta = 0.0019
    alpha = -0.05
  []
  [lr_disp_z]
    type = StressDivergenceIsolator
    block = 'isolator_elems'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    component = 2
    variable = disp_z
    static_initialization = true
    zeta = 0.0019
    alpha = -0.05
  []
  [lr_rot_x]
    type = StressDivergenceIsolator
    block = 'isolator_elems'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    component = 3
    variable = rot_x
    static_initialization = true
    zeta = 0.0019
    alpha = -0.05
  []
  [lr_rot_y]
    type = StressDivergenceIsolator
    block = 'isolator_elems'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    component = 4
    variable = rot_y
    static_initialization = true
    zeta = 0.0019
    alpha = -0.05
  []
  [lr_rot_z]
    type = StressDivergenceIsolator
    block = 'isolator_elems'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    component = 5
    variable = rot_z
    static_initialization = true
    zeta = 0.0019
    alpha = -0.05
  []
  [gravity]
    type = Gravity
    variable = disp_z
    value = -9.81
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases SGs int_wall int_slab RV_housing RV small_walls upper_basemat fluid_material RV_slab'
    alpha = -0.05
  []
[]

[AuxKernels]
  [Fb_x]
    type = MaterialRealCMMAux
    property = basic_forces
    row = 0
    column = 0
    variable = Fb_x
    block = 'isolator_elems'
  []
  [Fb_y]
    type = MaterialRealCMMAux
    property = basic_forces
    row = 1
    column = 0
    variable = Fb_y
    block = 'isolator_elems'
  []
  [Fb_z]
    type = MaterialRealCMMAux
    property = basic_forces
    row = 2
    column = 0
    variable = Fb_z
    block = 'isolator_elems'
  []
  [Defb_x]
    type = MaterialRealCMMAux
    property = deformations
    row = 0
    column = 0
    variable = Defb_x
    block = 'isolator_elems'
  []
  [Velb_x]
    type = MaterialRealCMMAux
    property = deformation_rates
    row = 0
    column = 0
    variable = Velb_x
    block = 'isolator_elems'
  []
  [Defb_y]
    type = MaterialRealCMMAux
    property = deformations
    row = 1
    column = 0
    variable = Defb_y
    block = 'isolator_elems'
  []
  [Defb_z]
    type = MaterialRealCMMAux
    property = deformations
    row = 2
    column = 0
    variable = Defb_z
    block = 'isolator_elems'
  []
  [accel_x]
    type = TestNewmarkTI
    displacement = disp_x
    variable = accel_x
    first = false
  []
  [vel_x]
    type = TestNewmarkTI
    displacement = disp_x
    variable = vel_x
  []
  [accel_y]
    type = TestNewmarkTI
    displacement = disp_y
    variable = accel_y
    first = false
  []
  [vel_y]
    type = TestNewmarkTI
    displacement = disp_y
    variable = vel_y
  []
  [accel_z]
    type = TestNewmarkTI
    displacement = disp_z
    variable = accel_z
    first = false
  []
  [vel_z]
    type = TestNewmarkTI
    displacement = disp_z
    variable = vel_z
  []
  [rot_accel_x]
    block = 'isolator_elems upper_rigid_elems'
    type = TestNewmarkTI
    displacement = rot_x
    variable = rot_accel_x
    first = false
  []
  [rot_vel_x]
    block = 'isolator_elems upper_rigid_elems'
    type = TestNewmarkTI
    displacement = rot_x
    variable = rot_vel_x
  []
  [rot_accel_y]
    block = 'isolator_elems upper_rigid_elems'
    type = TestNewmarkTI
    displacement = rot_y
    variable = rot_accel_y
    first = false
  []
  [rot_vel_y]
    block = 'isolator_elems upper_rigid_elems'
    type = TestNewmarkTI
    displacement = rot_y
    variable = rot_vel_y
  []
  [rot_accel_z]
    block = 'isolator_elems upper_rigid_elems'
    type = TestNewmarkTI
    displacement = rot_z
    variable = rot_accel_z
    first = false
  []
  [rot_vel_z]
    block = 'isolator_elems upper_rigid_elems'
    type = TestNewmarkTI
    displacement = rot_z
    variable = rot_vel_z
  []
[]

[Materials]
  [elasticity_concrete]
    type = ComputeIsotropicElasticityTensor
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases int_wall int_slab RV_housing small_walls RV_slab'
    youngs_modulus = 24.8 #GPa
    poissons_ratio = 0.2
  []
  [elasticity_rigid_concrete]
    type = ComputeIsotropicElasticityTensor
    block = 'upper_basemat'
    youngs_modulus = 99.2 #GPa # 4 x concrete for rigid basemat
    poissons_ratio = 0.2
  []
  [elasticity_steel_316]
    type = ComputeIsotropicElasticityTensor
    block = 'SGs RV'
    youngs_modulus = 170 #GPa
    poissons_ratio = 0.3
  []
  [elasticity_fluid]
    type = ComputeIsotropicElasticityTensor
    block = 'fluid_material'
    bulk_modulus = 2 #GPa #water
    poissons_ratio = 0.45 #water
  []
  [stress_1]
    type = ComputeFiniteStrainElasticStress
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases SGs int_wall int_slab RV_housing RV small_walls upper_basemat fluid_material RV_slab'
  []
  [concrete_density]
    type = GenericConstantMaterial
    block = 'roof ext_buttresses ext_walls int_buttresses SG_bases int_wall int_slab RV_housing small_walls upper_basemat RV_slab'
    prop_names = density
    prop_values = 2.4e-6 #e9kg/m3
  []
  [steel_density]
    type = GenericConstantMaterial
    block = 'SGs RV'
    prop_names = density
    prop_values = 7.85e-6 #e9kg/m3
  []
  [fluid_density]
    type = GenericConstantMaterial
    block = 'fluid_material'
    prop_names = density
    prop_values = 1.0e-6 #e9kg/m3 #water
  []

  [isolator_deformation]
    type = ComputeIsolatorDeformation
    sd_ratio = 0.5
    y_orientation = '1.0 0.0 0.0'
    displacements = 'disp_x disp_y disp_z'
    rotations = 'rot_x rot_y rot_z'
    velocities = 'vel_x vel_y vel_z'
    block = 'isolator_elems'
  []
  [elasticity]
    type = ComputeFPIsolatorElasticity
    mu_ref = 0.06
    p_ref = 0.006 # GPa
    block = 'isolator_elems'
    diffusivity = 4.4e-6
    conductivity = 18
    a = 100
    r_eff = 1.0 # meters. 2sec sliding period
    r_contact = 0.2
    uy = 0.001
    unit = 4
    beta = 0.275625
    gamma = 0.55
    pressure_dependent = false
    temperature_dependent = false
    velocity_dependent = false
    k_x = 78.53 # 7.853e10 N
    k_xx = 0.62282 # 622820743.6 N
    k_yy = 0.3114 # 311410371.8 N
    k_zz = 0.3114 # 311410371.8 N
  []

  [elasticity_beam_rigid]
    type = ComputeElasticityBeam
    youngs_modulus = 2e4
    poissons_ratio = 0.27
    shear_coefficient = 0.85
    block = 'upper_rigid_elems'
  []
  [stress_beam_rigid]
    type = ComputeBeamResultants
    block = 'upper_rigid_elems'
  []
[]

[Functions]
  [input_motion_x]
    type = PiecewiseLinear
    data_file = 'case2_scaled.csv'
    format = columns
    scale_factor = 9.81
    y_index_in_file = 1
    xy_in_file_only = false
  []
  [input_motion_y]
    type = PiecewiseLinear
    data_file = 'case2_scaled.csv'
    format = columns
    scale_factor = 9.81
    y_index_in_file = 2
    xy_in_file_only = false
  []
  [input_motion_z]
    type = PiecewiseLinear
    data_file = 'case2_scaled.csv'
    format = columns
    scale_factor = 9.81
    y_index_in_file = 3
    xy_in_file_only = false
  []
[]

[BCs]
  [x_motion]
    type = PresetAcceleration
    acceleration = accel_x
    velocity = vel_x
    variable = disp_x
    beta = 0.2725625
    boundary = 'bottom_isolators'
    function = 'input_motion_x'
  []
  [y_motion]
    type = PresetAcceleration
    acceleration = accel_y
    velocity = vel_y
    variable = disp_y
    beta = 0.2725625
    boundary = 'bottom_isolators'
    function = 'input_motion_y'
  []
  [z_motion]
    type = PresetAcceleration
    acceleration = accel_z
    velocity = vel_z
    variable = disp_z
    beta = 0.2725625
    boundary = 'bottom_isolators'
    function = 'input_motion_z'
  []
  [fixrxbot]
    type = DirichletBC
    variable = rot_x
    boundary = 'bottom_isolators'
    value = 0.0
  []
  [fixrybot]
    type = DirichletBC
    variable = rot_y
    boundary = 'bottom_isolators'
    value = 0.0
  []
  [fixrzbot]
    type = DirichletBC
    variable = rot_z
    boundary = 'bottom_isolators'
    value = 0.0
  []
  [fixrxcon]
    type = DirichletBC
    variable = rot_x
    boundary = 'connections'
    value = 0.0
  []
  [fixrycon]
    type = DirichletBC
    variable = rot_y
    boundary = 'connections'
    value = 0.0
  []
  [fixrzcon]
    type = DirichletBC
    variable = rot_z
    boundary = 'connections'
    value = 0.0
  []
[]

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

[Executioner]
  type = Transient
  petsc_options = '-ksp_snes_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  solve_type = 'NEWTON'
  nl_rel_tol = 1e-7
  nl_abs_tol = 1e-15
  dt = 0.01
  end_time = 28
  timestep_tolerance = 1e-6
  automatic_scaling = true
  [TimeIntegrator]
    type = NewmarkBeta
    beta = 0.275625
    gamma = 0.5
    inactive_tsteps = 2
  []
[]

[Controls]
  [inertia_switch]
    type = TimePeriod
    start_time = 0.0
    end_time = 0.03
    disable_objects = '*/inertia_x */inertia_y */inertia_z
                       */vel_x */vel_y */vel_z
                       */accel_x */accel_y */accel_z
                       */rot_vel_x */rot_vel_y */rot_vel_z
                       */rot_accel_x */rot_accel_y */rot_accel_z'
    set_sync_times = true
    execute_on = 'timestep_begin timestep_end'
  []
[]

[Postprocessors]
  [inp_accel_x]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = accel_x
  []
  [inp_accel_y]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = accel_y
  []
  [inp_accel_z]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = accel_z
  []
  [basemat_accel_x]
    type = PointValue
    point = '0.0 0.0 0.0'
    variable = accel_x
  []
  [basemat_accel_y]
    type = PointValue
    point = '0.0 0.0 0.0'
    variable = accel_y
  []
  [basemat_accel_z]
    type = PointValue
    point = '0.0 0.0 0.0'
    variable = accel_z
  []
  [iso1_fb_axial]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = Fb_x
  []
  [iso1_defb_axial]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = Defb_x
  []
  [iso1_fb_shear1]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = Fb_y
  []
  [iso1_defb_shear1]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = Defb_y
  []
  [iso1_fb_shear2]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = Fb_z
  []
  [iso1_defb_shear2]
    type = PointValue
    point = '5.0 0.0 -1.3'
    variable = Defb_z
  []
[]

[VectorPostprocessors]
  [accel_hist_x]
    type = ResponseHistoryBuilder
    variables = 'accel_x'
    nodes = '5252 2767 59044 24207 44503 41781 59152 38767 59100'
    # locations:
    # 5252-roof_edge
    # 2767-roof_center
    # 59044-RV_slab_center
    # 24207-SG_base
    # 44503-basemat_center-(0.35,-0.75,-1)(approx)
    # 41781-center_isolator_top-(5,0,-1)
    # 59152-center_isolator_bottom-(5,0,-1.3)
    # 38767-edge_isolator_top-(-45,30,-1)
    # 59100-edge_isolator_bottom(-45,30,-1.3)
    outputs = out1
  []
  [accel_spec_x]
    type = ResponseSpectraCalculator
    vectorpostprocessor = accel_hist_x
    regularize_dt = 0.01
    damping_ratio = 0.05
    start_frequency = 0.1
    end_frequency = 50
    outputs = out1
  []

  [accel_hist_y]
    type = ResponseHistoryBuilder
    variables = 'accel_y'
    nodes = '5252 2767 59044 24207 44503 41781 59152 38767 59100'
    # locations:
    # 5252-roof_edge
    # 2767-roof_center
    # 59044-RV_slab_center
    # 24207-SG_base
    # 44503-basemat_center-(0.35,-0.75,-1)(approx)
    # 41781-center_isolator_top-(5,0,-1)
    # 59152-center_isolator_bottom-(5,0,-1.3)
    # 38767-edge_isolator_top-(-45,30,-1)
    # 59100-edge_isolator_bottom(-45,30,-1.3)
    outputs = out1
  []
  [accel_spec_y]
    type = ResponseSpectraCalculator
    vectorpostprocessor = accel_hist_y
    regularize_dt = 0.01
    damping_ratio = 0.05
    start_frequency = 0.1
    end_frequency = 50
    outputs = out1
  []

  [accel_hist_z]
    type = ResponseHistoryBuilder
    variables = 'accel_z'
    nodes = '5252 2767 59044 24207 44503 41781 59152 38767 59100'
    # locations:
    # 5252-roof_edge
    # 2767-roof_center
    # 59044-RV_slab_center
    # 24207-SG_base
    # 44503-basemat_center-(0.35,-0.75,-1)(approx)
    # 41781-center_isolator_top-(5,0,-1)
    # 59152-center_isolator_bottom-(5,0,-1.3)
    # 38767-edge_isolator_top-(-45,30,-1)
    # 59100-edge_isolator_bottom(-45,30,-1.3)
    outputs = out1
  []
  [accel_spec_z]
    type = ResponseSpectraCalculator
    vectorpostprocessor = accel_hist_z
    regularize_dt = 0.01
    damping_ratio = 0.05
    start_frequency = 0.1
    outputs = out1
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
  csv = true
  [out1]
    type = CSV
    execute_on = 'final'
  []
[]
(msr/generic_msr/seismic_analysis/building_basemat_with_isolators_new.i)

Outputs

The output locations and responses for Part 2 include the same as those of Part 1 (an isolator at the center of the isolation system and the center of the basemat). In addition to these responses, the acceleration response at the roof, at the center of the reactor head, and the base of one of the steam generators is shown below.

The figures below present the isolator shear responses in XZ and YZ directions calculated in the same manner described in Part 1.

Figure 4: Shear force-displacement history in the XZ direction for an isolator at the center of the isolation system.

Figure 5: Shear force-displacement history in the YZ direction for an isolator at the center of the isolation system.

The figure below present the building responses for the same input ground motions shown in Part 1. The first row presents the spectral accelerations in all three directions at a node located at the center of the basemat of the building (same node as Part 1). The second, third, and fourth rows present these responses calculated at the center of the roof of the building, center of the reactor head, and at the base of one of the steam generators, respectively.

Figure 6: 5% damped acceleration response spectra in X direction (basemat center)

Figure 7: 5% damped acceleration response spectra in Y direction (basemat center)

Figure 8: 5% damped acceleration response spectra in Z direction (basemat center)

Figure 9: 5% damped acceleration response spectra in X direction (center of the building roof)

Figure 10: 5% damped acceleration response spectra in Y direction (center of the building roof)

Figure 11: 5% damped acceleration response spectra in Z direction (center of the building roof)

Figure 12: 5% damped acceleration response spectra in X direction (center of reactor vessel head)

Figure 13: 5% damped acceleration response spectra in Y direction (center of reactor vessel head)

Figure 14: 5% damped acceleration response spectra in Z direction (center of reactor vessel head)

Figure 15: 5% damped acceleration response spectra in X direction (base of one of the steam generators)

Figure 16: 5% damped acceleration response spectra in Y direction (base of one of the steam generators)

Figure 17: 5% damped acceleration response spectra in Z direction (base of one of the steam generators)

References

  1. C. Bolisetti, W. Hoffman, J. L. Coleman, S. S. Parsi, K. Lal, A. S. Whittaker, M. Cohen, K. Kramer, P. Kirchman, H. Bowers, and J. Redd. Seismic isolation of major advanced reactor systems for economic improvement and safety assurance. Technical Report INL/EXT-20-59608, Idaho National Laboratory, 2020.[BibTeX]