Loading SubKit/build/SquareLatticeLWR_Nodal.py +0 −2 Original line number Diff line number Diff line Loading @@ -1168,12 +1168,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): p = list(np.array(p) / np.average(p)) me.model.addAxialPowerShape('fuelrods', z, p) # me.model.addAxialPowerShape('guidetubes', z, p) radPowers = [] weights = [] for obj in me.model.solidObjects.values(): # obj.powID = 'allrods' # Need to weight the rod power when normalizing radial power shape assert(obj.symtype == 1) if radialPower: Loading tests/nodal_p7_uniformPow/buildDeck.py~deleted 100755 → 0 +0 −306 Original line number Diff line number Diff line #!/usr/bin/env python # Author: Bob Salko # Date: 1/23/19 # Description: # Makes a model of CASL problem 7 using 4 subchannels per assembly # # The following new features are being tested by this script: # - A custom name is used for the input file, ensuring that SubKit script mode can make a model with non # default name # - The ability to set percentTheoreticalDensity and the dynamic gap model is being tested # - The ability to set conductor options related to fuel material and gap models (Card 9.1 parameters) # - The ability to disable the main text output in Card Group 14 # - The ability to provide the total core power to be used to calculate the average LHR # - The ability to set global multipliers on wide and narrow gap widths and gap length # - The ability to set RTWFP and convergence criteria for transient input decks import sys sys.path.insert(0, '../../') import SubKit.build as cdb import SubKit.utils.UnitConversions as units from SubKit.build.SquareLatticeLWR_Nodal import SquareLatticeLWR_Nodal import copy import math import numpy as np def main(): dFuelOuter = 0.475e-2 * 2 # m dFuelInner = 0.418e-2 * 2 # m dFuelPellet = 0.4096e-2 * 2 # m dGuideTubeOuter = 0.602e-2 * 2 # m dGuideTubeInner = 0.561e-2 * 2 # m assemPitch = 21.5e-2 # m pitch = 1.26e-2 # m gapWidthMultWide = 0.1 gapWidthMultNarrow = 0.5 gapLengthMult = 2.0 baffleGap = 0.19e-2 # m inletMassFlux = 0.3483429E+04 # kg/m**2/s inletTemp = 291.85 # C outletPressure = 155.132039 # bar corePow = 3.411e6 # kW directHeat = 0.026 plenumPressure = 20.0 # bar plenumVolume = 8.783E-06 # m^3 gridLocations = list( np.array([13.884, 75.2, 127.4, 179.6, 231.8, 284.0, 336.2]) * units.t_cm_m) gridHeights = list( np.array([3.866, 3.810, 3.810, 3.810, 3.810, 3.810, 3.810]) * units.t_cm_m) axial_edit_bounds = list(np.array([11.951, 15.817, 24.028, 32.239, 40.45, 48.662, 56.873, 65.084, 73.295, 77.105, 85.17, 93.235, 101.3, 109.365, 117.43, 125.495, 129.305, 137.37, 145.435, 153.5, 161.565, 169.63, 177.695, 181.505, 189.57, 197.635, 205.7, 213.765, 221.83, 229.895, 233.705, 241.77, 249.835, 257.9, 265.965, 274.03, 282.095, 285.905, 293.97, 302.035, 310.1, 318.165, 326.23, 334.295, 338.105, 346.0262, 353.9474, 361.8686, 369.7898, 377.711]) * units.t_cm_m) formLoss = 0.9070 # Setup axial mesh gridLocations = list(np.array(gridLocations) - axial_edit_bounds[0]) axial_edit_bounds = list( np.array(axial_edit_bounds) - axial_edit_bounds[0]) dz = [] for j in range(1, len(axial_edit_bounds)): dz.append(axial_edit_bounds[j] - axial_edit_bounds[j - 1]) numLevels = len(dz) gridBottoms = list(np.array(gridLocations) - np.array(gridHeights) / 2.0) gridLevels = [] for grid in gridBottoms: found = False for level, z in enumerate(axial_edit_bounds): if np.isclose(grid, z): gridLevels.append(level + 1) found = True assert(found) assert(len(gridLevels) == len(gridLocations)) # The map of the core. Each 1 represents an assembly coreMap = [[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0], [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0], [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]] # Fuel rod is type 1, but mark it with an "f" to make it more visible # Guide tube is type 2 f = 1 pwr17x17 = [ [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f], [f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f], [f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f]] # The multiple gap connects are not working with CTF parallel for some reason rodDomains = [[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0], [0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], [0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], [0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], [3, 3, 3, 3, 3, 3, 3, 4, 2, 2, 2, 2, 2, 2, 2], [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4], [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4], [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4], [0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], [0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], [0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0], [0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], # [ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 0, 0], # [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 6, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6], # [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0], # [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0], # [ 0, 0, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 0, 0], # [ 0, 0, 0, 0, 4, 4, 5, 5, 5, 6, 6, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0], # [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], # [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], # [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], # [ 0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0], # [ 0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]] rodDomains = [[0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 3, 0, 0, 0, 0], [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 0, 0], [0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0], [0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], [1, 1, 1, 1, 7, 7, 7, 2, 2, 2, 3, 3, 3, 3, 3], [1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 3], [1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6], [4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6], [4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6], [4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], [0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], [0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], [0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0], [0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 1, 1, 1, 1, 1, 1, 5, 5, 5, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2], # [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2], # [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4], # [ 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4], # [ 3, 3, 3, 3, 3, 3, 5, 5, 5, 4, 4, 4, 4, 4, 4], # [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], # [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], # [ 0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0], # [ 0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]] print np.unique(np.array(rodDomains), return_counts=True) model = cdb.Model() model.setTitle("4-Loop nodal mesh loss of flow") solidGeoms = {1: cdb.FuelRod(d_outer=dFuelOuter, d_inner=dFuelInner, d_pellet=dFuelPellet, percentTheoreticalDensity=0.9450, gapModel='dynamic', # Fuel rod plenumPressure=plenumPressure, plenumVolume=plenumVolume), # Fuel rod 2: cdb.Tube(d_outer=dGuideTubeOuter, d_inner=dGuideTubeInner)} # Guide tube # Radial power shape radPow = [[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00], [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91, 1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00], [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18, 0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00], [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93, 1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00], [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22, 0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75], [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94, 1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76], [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12, 0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77], [0.78, 0.91, 1.18, 0.93, 1.22, 0.94, 1.12, 0.90, 1.12, 0.94, 1.22, 0.93, 1.18, 0.91, 0.78], [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12, 0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77], [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94, 1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76], [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22, 0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75], [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93, 1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00], [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18, 0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00], [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91, 1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00], [0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00]] builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch, baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels, rodDomains=rodDomains, modelGuideTubes=True, gapWidthMultWide=gapWidthMultWide, gapWidthMultNarrow=gapWidthMultNarrow, gapLengthMult=gapLengthMult) # Run to 1 second to let the solution reach steady state # Test the ability to provide rtwfp during transients (used for system coupled calcs) # model.addTransientGroup(tEnd=1.0, rtwfp=1e4) # Run a 10 second transient # model.addTransientGroup(tEnd=11.0, editInterval=1.0, rtwfp=1e4) # Write card group 19 conv criteria to the input file, even for this transient model model.setSteadyState(transPrintConv=True, energyBalance=1e-4, massBalance=1e-5, maxIterations=200000) # Ramp the power down over time time = [1.0, 4.0, 4.3, 4.6, 4.9, 5.2, 5.5, 5.8, 6.1, 6.4, 6.7, 7.0, 11.0] power = [1.0, 1.0, 0.99, 0.95, 0.90, 0.80, 0.68, 0.44, 0.27, 0.21, 0.19, 0.18, 0.16] # model.addTotalPowerRamp(time, power) # Set the core boundary conditions # Apply a flow ramp to represent the loss of flow time = [1.0, 1.1, 1.2, 1.4, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.8, 4.0, 4.6, 5.4, 6.5, 7.9, 9.2, 10.0] flow = [1.00, 1.00, 0.90, 0.83, 0.78, 0.75, 0.71, 0.70, 0.67, 0.66, 0.64, 0.62, 0.61, 0.58, 0.55, 0.51, 0.46, 0.43, 0.40] builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow]) builder.setCoreOutletBC(outletPressure, inletTemp) # Set the nominal power and direct heating to coolant model.setTotalCorePower(power=corePow, dhfrac=directHeat) # Turn on rod VTK and DNB edits. They are both off by default. builder.setEditOptions(dnbEdits=True, veracsHDF5=True, mainTextOutput=False) # Use default fuel modeling parameters model.setConductorModelOptions( irelf=0, iconf=0, imwr=0, ifswell=0, ifdens=0, ifreloc=0, ifcladcreep=0) coreHeight = np.sum(np.array(dz)) def axPow(z): return np.cos((z - coreHeight / 2) * np.pi / 2 / coreHeight) # Pass to the core builder builder.setCorePowerShape(radialPower=radPow, axialPower=axPow) model.generateModel('nodal_p7.inp') if __name__ == '__main__': main() Loading
SubKit/build/SquareLatticeLWR_Nodal.py +0 −2 Original line number Diff line number Diff line Loading @@ -1168,12 +1168,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): p = list(np.array(p) / np.average(p)) me.model.addAxialPowerShape('fuelrods', z, p) # me.model.addAxialPowerShape('guidetubes', z, p) radPowers = [] weights = [] for obj in me.model.solidObjects.values(): # obj.powID = 'allrods' # Need to weight the rod power when normalizing radial power shape assert(obj.symtype == 1) if radialPower: Loading
tests/nodal_p7_uniformPow/buildDeck.py~deleted 100755 → 0 +0 −306 Original line number Diff line number Diff line #!/usr/bin/env python # Author: Bob Salko # Date: 1/23/19 # Description: # Makes a model of CASL problem 7 using 4 subchannels per assembly # # The following new features are being tested by this script: # - A custom name is used for the input file, ensuring that SubKit script mode can make a model with non # default name # - The ability to set percentTheoreticalDensity and the dynamic gap model is being tested # - The ability to set conductor options related to fuel material and gap models (Card 9.1 parameters) # - The ability to disable the main text output in Card Group 14 # - The ability to provide the total core power to be used to calculate the average LHR # - The ability to set global multipliers on wide and narrow gap widths and gap length # - The ability to set RTWFP and convergence criteria for transient input decks import sys sys.path.insert(0, '../../') import SubKit.build as cdb import SubKit.utils.UnitConversions as units from SubKit.build.SquareLatticeLWR_Nodal import SquareLatticeLWR_Nodal import copy import math import numpy as np def main(): dFuelOuter = 0.475e-2 * 2 # m dFuelInner = 0.418e-2 * 2 # m dFuelPellet = 0.4096e-2 * 2 # m dGuideTubeOuter = 0.602e-2 * 2 # m dGuideTubeInner = 0.561e-2 * 2 # m assemPitch = 21.5e-2 # m pitch = 1.26e-2 # m gapWidthMultWide = 0.1 gapWidthMultNarrow = 0.5 gapLengthMult = 2.0 baffleGap = 0.19e-2 # m inletMassFlux = 0.3483429E+04 # kg/m**2/s inletTemp = 291.85 # C outletPressure = 155.132039 # bar corePow = 3.411e6 # kW directHeat = 0.026 plenumPressure = 20.0 # bar plenumVolume = 8.783E-06 # m^3 gridLocations = list( np.array([13.884, 75.2, 127.4, 179.6, 231.8, 284.0, 336.2]) * units.t_cm_m) gridHeights = list( np.array([3.866, 3.810, 3.810, 3.810, 3.810, 3.810, 3.810]) * units.t_cm_m) axial_edit_bounds = list(np.array([11.951, 15.817, 24.028, 32.239, 40.45, 48.662, 56.873, 65.084, 73.295, 77.105, 85.17, 93.235, 101.3, 109.365, 117.43, 125.495, 129.305, 137.37, 145.435, 153.5, 161.565, 169.63, 177.695, 181.505, 189.57, 197.635, 205.7, 213.765, 221.83, 229.895, 233.705, 241.77, 249.835, 257.9, 265.965, 274.03, 282.095, 285.905, 293.97, 302.035, 310.1, 318.165, 326.23, 334.295, 338.105, 346.0262, 353.9474, 361.8686, 369.7898, 377.711]) * units.t_cm_m) formLoss = 0.9070 # Setup axial mesh gridLocations = list(np.array(gridLocations) - axial_edit_bounds[0]) axial_edit_bounds = list( np.array(axial_edit_bounds) - axial_edit_bounds[0]) dz = [] for j in range(1, len(axial_edit_bounds)): dz.append(axial_edit_bounds[j] - axial_edit_bounds[j - 1]) numLevels = len(dz) gridBottoms = list(np.array(gridLocations) - np.array(gridHeights) / 2.0) gridLevels = [] for grid in gridBottoms: found = False for level, z in enumerate(axial_edit_bounds): if np.isclose(grid, z): gridLevels.append(level + 1) found = True assert(found) assert(len(gridLevels) == len(gridLocations)) # The map of the core. Each 1 represents an assembly coreMap = [[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0], [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0], [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0]] # Fuel rod is type 1, but mark it with an "f" to make it more visible # Guide tube is type 2 f = 1 pwr17x17 = [ [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f], [f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f, 2, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, 2, f, f, f, f, f, f, f, f, f, 2, f, f, f], [f, f, f, f, f, 2, f, f, 2, f, f, 2, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f], [f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f, f]] # The multiple gap connects are not working with CTF parallel for some reason rodDomains = [[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0], [0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], [0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], [0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], [3, 3, 3, 3, 3, 3, 3, 4, 2, 2, 2, 2, 2, 2, 2], [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4], [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4], [3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4], [0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], [0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], [0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0], [0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], # [ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 0, 0], # [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 1, 2, 2, 2, 6, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6], # [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0], # [ 0, 4, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 0], # [ 0, 0, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 0, 0], # [ 0, 0, 0, 0, 4, 4, 5, 5, 5, 6, 6, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 0, 0], # [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], # [ 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], # [ 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], # [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], # [ 0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], # [ 0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0], # [ 0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]] rodDomains = [[0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 3, 0, 0, 0, 0], [0, 0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 0, 0], [0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0], [0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 0], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3], [1, 1, 1, 1, 7, 7, 7, 2, 2, 2, 3, 3, 3, 3, 3], [1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 3], [1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6], [4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6], [4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6], [4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6], [0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], [0, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 0], [0, 0, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 0, 0], [0, 0, 0, 0, 4, 5, 5, 5, 5, 5, 6, 0, 0, 0, 0]] # rodDomains = [[ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0], # [ 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 0], # [ 1, 1, 1, 1, 1, 1, 5, 5, 5, 2, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2], # [ 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2], # [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 2, 2, 2, 2], # [ 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4], # [ 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4], # [ 3, 3, 3, 3, 3, 3, 5, 5, 5, 4, 4, 4, 4, 4, 4], # [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], # [ 0, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 0], # [ 0, 0, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 0, 0], # [ 0, 0, 0, 0, 3, 3, 3, 4, 4, 4, 4, 0, 0, 0, 0]] print np.unique(np.array(rodDomains), return_counts=True) model = cdb.Model() model.setTitle("4-Loop nodal mesh loss of flow") solidGeoms = {1: cdb.FuelRod(d_outer=dFuelOuter, d_inner=dFuelInner, d_pellet=dFuelPellet, percentTheoreticalDensity=0.9450, gapModel='dynamic', # Fuel rod plenumPressure=plenumPressure, plenumVolume=plenumVolume), # Fuel rod 2: cdb.Tube(d_outer=dGuideTubeOuter, d_inner=dGuideTubeInner)} # Guide tube # Radial power shape radPow = [[0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00], [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91, 1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00], [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18, 0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00], [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93, 1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00], [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22, 0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75], [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94, 1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76], [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12, 0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77], [0.78, 0.91, 1.18, 0.93, 1.22, 0.94, 1.12, 0.90, 1.12, 0.94, 1.22, 0.93, 1.18, 0.91, 0.78], [0.77, 1.10, 0.97, 1.19, 0.95, 1.19, 0.94, 1.12, 0.94, 1.19, 0.95, 1.19, 0.97, 1.10, 0.77], [0.76, 0.93, 1.17, 0.94, 1.19, 0.98, 1.19, 0.94, 1.19, 0.98, 1.19, 0.94, 1.17, 0.93, 0.76], [0.75, 1.10, 0.93, 1.12, 0.95, 1.19, 0.95, 1.22, 0.95, 1.19, 0.95, 1.12, 0.93, 1.10, 0.75], [0.00, 0.75, 1.12, 0.93, 1.12, 0.94, 1.19, 0.93, 1.19, 0.94, 1.12, 0.93, 1.12, 0.75, 0.00], [0.00, 0.00, 0.74, 1.12, 0.93, 1.17, 0.97, 1.18, 0.97, 1.17, 0.93, 1.12, 0.74, 0.00, 0.00], [0.00, 0.00, 0.00, 0.75, 1.10, 0.93, 1.10, 0.91, 1.10, 0.93, 1.10, 0.75, 0.00, 0.00, 0.00], [0.00, 0.00, 0.00, 0.00, 0.75, 0.76, 0.77, 0.78, 0.77, 0.76, 0.75, 0.00, 0.00, 0.00, 0.00]] builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch, baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels, rodDomains=rodDomains, modelGuideTubes=True, gapWidthMultWide=gapWidthMultWide, gapWidthMultNarrow=gapWidthMultNarrow, gapLengthMult=gapLengthMult) # Run to 1 second to let the solution reach steady state # Test the ability to provide rtwfp during transients (used for system coupled calcs) # model.addTransientGroup(tEnd=1.0, rtwfp=1e4) # Run a 10 second transient # model.addTransientGroup(tEnd=11.0, editInterval=1.0, rtwfp=1e4) # Write card group 19 conv criteria to the input file, even for this transient model model.setSteadyState(transPrintConv=True, energyBalance=1e-4, massBalance=1e-5, maxIterations=200000) # Ramp the power down over time time = [1.0, 4.0, 4.3, 4.6, 4.9, 5.2, 5.5, 5.8, 6.1, 6.4, 6.7, 7.0, 11.0] power = [1.0, 1.0, 0.99, 0.95, 0.90, 0.80, 0.68, 0.44, 0.27, 0.21, 0.19, 0.18, 0.16] # model.addTotalPowerRamp(time, power) # Set the core boundary conditions # Apply a flow ramp to represent the loss of flow time = [1.0, 1.1, 1.2, 1.4, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.8, 4.0, 4.6, 5.4, 6.5, 7.9, 9.2, 10.0] flow = [1.00, 1.00, 0.90, 0.83, 0.78, 0.75, 0.71, 0.70, 0.67, 0.66, 0.64, 0.62, 0.61, 0.58, 0.55, 0.51, 0.46, 0.43, 0.40] builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow]) builder.setCoreOutletBC(outletPressure, inletTemp) # Set the nominal power and direct heating to coolant model.setTotalCorePower(power=corePow, dhfrac=directHeat) # Turn on rod VTK and DNB edits. They are both off by default. builder.setEditOptions(dnbEdits=True, veracsHDF5=True, mainTextOutput=False) # Use default fuel modeling parameters model.setConductorModelOptions( irelf=0, iconf=0, imwr=0, ifswell=0, ifdens=0, ifreloc=0, ifcladcreep=0) coreHeight = np.sum(np.array(dz)) def axPow(z): return np.cos((z - coreHeight / 2) * np.pi / 2 / coreHeight) # Pass to the core builder builder.setCorePowerShape(radialPower=radPow, axialPower=axPow) model.generateModel('nodal_p7.inp') if __name__ == '__main__': main()