Commit bd44addb authored by Gurecky, William's avatar Gurecky, William
Browse files

ensures radial power factors are set to 1.0 when a full 3d power distribution... 

ensures radial power factors are set to 1.0 when a full 3d power distribution is specified to LWR Nodal deck builder
parent 9dafc1e6

SubKit/build/SquareLatticeLWR_Nodal.py

+39 −7
Original line number Diff line number Diff line
@@ -895,9 +895,30 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 


   @property

   def z_centers(me):

      z_grid = np.cumsum(me.dz) - me.dz[0]

      z_grid[-1] = np.cumsum(me.dz)[-1] + 0.0001

      return z_grid

      """ Axial level z centers.  Has length = len(z_bounds) - 1 """

      z_center_grid = np.average(me.z_bounds_pairs, axis=1)

      return z_center_grid



   @property

   def z_power_grid(me):

      """ manually adjust z center endpoints to be consistent """

      z_pwr_grid = me.z_centers

      z_pwr_grid[0] = 0.0

      z_pwr_grid[-1] = np.cumsum(me.dz)[-1] + 0.0001

      return z_pwr_grid



   @property

   def z_bounds(me):

      """ Axial level z boundaries """

      z_bound = np.cumsum(me.dz)

      z_bound = np.insert(z_bound, 0, 0.0)

      return z_bound



   @property

   def z_bounds_pairs(me):

      """ Axial level z boundary pairs ((zstart_i, zend_i), ...) """

      z_pairs = np.asarray((me.z_bounds[1:], me.z_bounds[:-1])).T

      return z_pairs



   def setCorePowerShape3D(me, powerDist, wgts=None, coreStart=0.0):

      """
@@ -907,18 +928,24 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
      Args:

         powerDist (numpy.ndarray): 3 dim array. Contains power distribution with indixing

           (asm_i, asm_j, asm_axial)

         wgts (numpy.ndarray): Same shape as powerDist. Weights of powers (optional)

         coreStart (float): The axial location of the start of the core.  In CTF, the bottom of the

            model is zero.  If the core does not start at the bottom of the model, this must be provided

            so the power tables are correctly defined in CTF. [m]

      """

      assert len(powerDist.shape) == 3

      assert len(me.z_centers.shape) == 1

      assert len(me.z_power_grid.shape) == 1

      assert powerDist.shape[0] == me.chMapObj.getAssemDimLen()

      assert powerDist.shape[1] == me.chMapObj.getAssemDimLen()

      assert me.z_centers.size == powerDist.shape[2]

      assert me.z_power_grid.size == powerDist.shape[2]



      # expand powerDist to NodalMesh



      # normalize 3d powers

      if wgts is None:

          wgts = np.ones(powerDist.shape)

      assert wgts.shape == powerDist.shape

      powerDist *= wgts

      powerDist /= np.average(powerDist, weights=wgts)



      radPowers = []
@@ -932,5 +959,10 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core): 
                  rods = me.rodsInAssem[assemID]

                  for rod in rods:

                     powID = "pwr_" + str(rod)

                     me.model.addAxialPowerShape(powID, me.z_centers, powerDist[i, j, :])

                     me.model.addAxialPowerShape(powID, me.z_power_grid, powerDist[i, j, :])

                     obj.powID = powID

                     # radial power factors are implicit in the 3d power distribution

                     if powerDist[i, j, 0] == 0.0:

                         obj.power = 0.0

                     else:

                         obj.power = 1.0

tests/nodal_p7_pow3d/buildDeck.py

+252 −230
Original line number Diff line number Diff line
@@ -4,17 +4,17 @@ 
# Description:

# Makes a model of CASL problem 7 using 4 subchannels per assembly

#

import numpy as np

import math

import copy

from SubKit.build.SquareLatticeLWR_Nodal import SquareLatticeLWR_Nodal

import SubKit.utils.UnitConversions as units

import SubKit.build as cdb

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():

def main(power_3d_flag=True):



    dFuelOuter = 0.475e-2*2  # m

    dFuelInner = 0.418e-2*2  # m
@@ -29,8 +29,10 @@ def main(): 
    outletPressure = 155.132039  # bar

    linearHeatRate = 16.71977  # kW/m

    directHeat = 0.026

   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)

    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,
@@ -56,7 +58,6 @@ def main(): 
        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],
@@ -199,23 +200,36 @@ def main(): 


    # 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.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, \

    builder = SquareLatticeLWR_Nodal(model, coreMap, pwr17x17, pitch, dz, solidGeoms, assemPitch=assemPitch,

                                     baffleGap=baffleGap, assemLosses=formLoss, assemLossLevels=gridLevels, rodDomains=rodDomains,

                                     modelGuideTubes=True)



    # Run to 1 second to let the solution reach steady state
@@ -223,15 +237,14 @@ def main(): 
    # Run a 10 second transient

    model.addTransientGroup(tEnd=11.0, editInterval=1.0)



   # 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]

    # Ramp the power

    time = [1.0, 4.0]

    power = [1.0, 1.0]

    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]

    time = [1.0,  1.1]

    flow = [1.00, 1.00]

    builder.setCoreInletBC(inletMassFlux, inletTemp, flowRamp=[time, flow])

    builder.setCoreOutletBC(outletPressure, inletTemp)
@@ -242,20 +255,29 @@ def main(): 
    model.setEditOptions(dnbEdits=True)



    coreHeight = np.sum(np.array(dz))



    def axPow(z):

        return np.cos((z-coreHeight/2)*np.pi/2/coreHeight)



    if power_3d_flag:

        # set 3d power distribution

        axialMesh = builder.z_centers

        radPow_arr = np.array(radPow)

   radPow_arr_full = np.ones((radPow_arr.shape[0], radPow_arr.shape[1], len(axialMesh)))

        radPow_arr_full = np.ones(

            (radPow_arr.shape[0], radPow_arr.shape[1], len(axialMesh)))

        for z_i in range(len(axialMesh)):

       radPow_arr_full[:, :, z_i] = radPow_arr_full[:, :, z_i] * radPow_arr

            radPow_arr_full[:, :, z_i] = radPow_arr_full[:,

                                                         :, z_i] * radPow_arr

        axPow_arr = axPow(np.array(axialMesh))

        powerDist3d = radPow_arr_full * axPow_arr

        builder.setCorePowerShape3D(powerDist3d)

        model.generateModel('deck.inp')

    else:

        builder.setCorePowerShape(radialPower=radPow, axialPower=axPow)

        model.generateModel('orig_deck.inp')



   model.generateModel()



if __name__ == '__main__':

   main()
 
    # build 2 decks which should produce nearly the same CTF TH solution output

    main(power_3d_flag=True)

    main(power_3d_flag=False)