Commit ba9e61c4 authored by Wysocki, Aaron's avatar Wysocki, Aaron
Browse files

Adding grid enhancement factors and dynamic gap model inputs

Description:
Allows setting grid heat transfer enhancement coefficients and
blockages. Turned on by default if grids exist in the model. Allows
enabling the dynamic gap model (disabled by default). Both models use
default VERAIn gap parameters if custom values not provided.

Gitlab ticket - 4116
parent f551ae6a
Loading
Loading
Loading
Loading
+4 −1
Original line number Diff line number Diff line
*.pyc
deck.inp
*.inp
!/tests/testDriver/deck.inp
!/tests/testDriver/deck_same.inp
!/tests/testDriver/deck_diff.inp
SubKit.egg-info
MCFR
.DS_Store
+108 −1
Original line number Diff line number Diff line
@@ -72,6 +72,15 @@ class Model:
        me.formLossChannels = []
        # Same size as formLosses.  Gives the momentum level in the channel of the loss.
        me.formLossLevels = []
        # Number of grid types, global across all assem types
        me.numGridTypes = 0
        # Dict of grid heat transfer enhancement values with the following structure:
        #     gridTypeNum (global grid type number across all assem types)
        #         'level' (float)
        #         'chans' (list of all chans associated with this grid type)
        #         'blockageRatio' (float)
        #         'enhancementCoeffs' (list of length 4)
        me.gridEnhancementDict = {}
        # Average linear heat rate and direct heating
        me.qprime = 0.0  # kW/m
        me.dhfrac = 0.0
@@ -811,7 +820,7 @@ class Model:
            assert(k > 0.0)
            kList = [k] * len(levelsList)
        channelsIn = []
        if not channels:
        if channels is None:
            for sec in me.sections.keys():
                for ch in me.sections[sec].channels.keys():
                    channelsIn.append(ch)
@@ -825,6 +834,7 @@ class Model:
                for ch in channels:
                    assert(isinstance(ch, int))
                channelsIn = channels

        for ch in channelsIn:
            for i, level in enumerate(levelsList):
                me.formLosses.append(kList[i])
@@ -844,6 +854,103 @@ class Model:
                    raise RuntimeError(
                        "Channel " + str(ch) + " Level " + str(level) + " has multiple form losses defined")

    def addGridEnhancement(me, levels, gridEnhancementCoeffs=None, gridBlockageRatios=None, channels=None):
        """ Specify parameters in the Yao-Hochreiter-Leech formula for spacer grid convective heat transfer
               enhancement. The formula is:
                  Nu/Nu0 = [1+a*eps^2*exp(b*x/D)]*[1+A^2*tan^2(PHI)*exp(c*X/D)]]^d
            gridEnhancementCoeffs specifies parameters a, b, c, and d in the formula. gridBlockageRatios
            specifies parameter eps in the formula, which represents the blocked flow area divided by the
            bare flow area in the absence of spacer grids.

        Args:
           gridEnhancementCoeffs (list or list of lists): Coefficients a, b, c, and d in the Yao-Hochreiter-Leech
              model. A single list (len=4) may be provided (e.g. [0.1, 0.0, 0.3, -0.2], to be used for each
              location in levels in the specified channels. Alternately, a separate list (len=4) may be
              provided for each entry in levels, such that gridEnhancementCoeffs is a list of lists.  For
              example, if levels=[5,10,15], a valid gridEnhancementCoeffs would be [[0.1, 0.0, 0.3, -0.2],
              [2.0, 0.4, 0.3, 0.0],[-2.0, 1.3, 0.0, 0.0]]. If gridEnhancementCoeffs=None, the default
              parameter values will be applied for all levels and channels.
           gridBlockageRatios (float or list): eps in the Yao-Hochreiter-Leech model. Enter as a float to be
              applied to all locations in levels. Alternately, enter as a list of floats specifying the value
              at each location in levels. If None, a default value will be applied.
           levels (int or list): Can be one level or a list of levels where the form loss is applied. A unique
              grid type will be assigned to each level in this list. 
           channels (int or list): Can be one channel or a list of channels where the form loss is applied.
              Use channel IDs that were used when adding the channels to the model.  If None, the form loss
              is applied uniformly in all channels at the passed levels.
              

        """
        levelsList = []
        if isinstance(levels, int):
            levelsList = [levels]
        elif isinstance(levels, list):
            for j in levels:
                assert(isinstance(j, int))
            levelsList = levels
        else:
            raise RuntimeError(
                "Unexpected entry for levels in addGridEnhancement " + str(levels))
        if gridEnhancementCoeffs is not None:
            if isinstance(gridEnhancementCoeffs, list):
                if all([isinstance(entry,list) for entry in gridEnhancementCoeffs]):
                    assert(len(gridEnhancementCoeffs) == len(levelsList))
                    for coeffs in gridEnhancementCoeffs:
                        assert(all([isinstance(coeff,float) for coeff in coeffs]))
                        assert(len(coeffs)==4)
                    gridEnhancementCoeffsList = gridEnhancementCoeffs
                elif all([isinstance(entry,float) for entry in gridEnhancementCoeffs]):
                    assert(len(gridEnhancementCoeffs)==4)
                    gridEnhancementCoeffsList = [gridEnhancementCoeffs] * len(levelsList)
                else:
                    raise RuntimeError(
                        "gridEnhancementCoeffs must be a list of lists or a list of floats " + str(gridEnhancementCoeffs))
            else:
                raise RuntimeError(
                    "Unexpected entry for gridEnhancementCoeffs in addGridEnhancement " + str(gridEnhancementCoeffs))

        if gridBlockageRatios is not None:
            if isinstance(gridBlockageRatios, list):
                assert(all(gridBlockageRatios) >= 0.0)
                assert(len(gridBlockageRatios) == len(levelsList))
                gridBlockageRatiosList = gridBlockageRatios
            elif isinstance(gridBlockageRatios, float):
                assert(gridBlockageRatios > 0.0)
                gridBlockageRatiosList = [gridBlockageRatios] * len(levelsList)

        channelsIn = []
        if channels is None:
            for sec in me.sections.keys():
                for ch in me.sections[sec].channels.keys():
                    channelsIn.append(ch)
            if len(channelsIn) == 0:
                warnings.warn(
                    "No channels in model when calling gridEnhancementCoeffs.  Ensure channels were added first or call gridEnhancementCoeffs with channel indices specified.", Warning)
        else:
            if isinstance(channels, int):
                channelsIn = [channels]
            elif isinstance(channels, list):
                for ch in channels:
                    assert(isinstance(ch, int))
                channelsIn = channels

        # Put the channel IDs in order, since that's how xml2ctf prints them
        channelsIn.sort()

        startingGridNum = len(me.gridEnhancementDict)
        for (grid,level) in enumerate(levelsList):
            gridNum = startingGridNum + grid
            me.gridEnhancementDict[gridNum] = {}
            if gridBlockageRatios is not None:
                me.gridEnhancementDict[gridNum]['blockageRatio'] = gridBlockageRatiosList[grid]
            if gridEnhancementCoeffs is not None:
                me.gridEnhancementDict[gridNum]['enhancementCoeffs'] = gridEnhancementCoeffsList[grid]
            me.gridEnhancementDict[gridNum]['level'] = level
            me.gridEnhancementDict[gridNum]['chans'] = channelsIn

        # Sets the total number of grid types in the model
        me.numGridTypes = len(me.gridEnhancementDict)

    def addBoundaryCondition(me, bc, chID, lev):
        """ Add a boundary condition object to the model.

+17 −1
Original line number Diff line number Diff line
@@ -79,14 +79,30 @@ class FuelRod(SolidGeo):
       num_rings (int): Number of rings to model in the pellet region.
       percentTheoreticalDensity (float): The percent of theoretical density of the pellet material in the pellet
          region.
       gapModel (string): 'constant' or 'dynamic'; the model used to determine gap conductance
       plenumPressure (float): Cold pressure of the fuel rod gap and plenum. Required if gapModel='dynamic'. Not
          used if gapModel='constant'.
       plenumVolume (float): Gas plenum volume at cold conditions. Required if gapModel='dynamic'. Not
          used if gapModel='constant'.

    """
    def __init__(me, d_outer, d_inner, d_pellet, pellet_material=None, clad_material=None, num_rings=10, percentTheoreticalDensity=1.0):
    def __init__(me, d_outer, d_inner, d_pellet, pellet_material=None, clad_material=None, num_rings=10, percentTheoreticalDensity=1.0,
                 gapModel=None, plenumPressure=None, plenumVolume=None):
        me.d_outer = d_outer
        me.d_inner = d_inner
        me.d_pellet = d_pellet
        me.num_rings = num_rings
        me.percentTheoreticalDensity = percentTheoreticalDensity
        me.gapModel = gapModel
        if gapModel == 'dynamic':
            assert(plenumPressure is not None)
            assert(plenumVolume is not None)
            me.plenumPressure = plenumPressure
            me.plenumVolume = plenumVolume
            
            # Gap thickness is equal to: (cladding inner radius) - (pellet outer radius) = 0.5 * (d_inner - d_pellet)
            # Only used if the dynamic gap model is enabled.
            me.gapThickness = 0.5 * (d_inner - d_pellet)


class _Material:
+169 −35
Original line number Diff line number Diff line
@@ -7,7 +7,6 @@ import Section
import utils.utils as utils
from scipy.interpolate import interp1d


class NodalChannelGeom(object):
    """ Calculates geometry of a nodal channel in the model

@@ -379,6 +378,18 @@ class AssemblyChannelRelations(object):
            ids.remove(0)
        return ids

    def getChansInAssemType(me,asmTypeID):
        """ Returns a list of all channels in all assemblies of the specified assembly type """

        # Loop over the core assembly map and for each match of this assemID, get the channels in that
        # assembly and add to the list.
        chansInAssemType = []
        for i, row in enumerate(me.assemTypeIdx):
            for j, assemTypeIdx in enumerate(row):
                if assemTypeIdx == asmTypeID:
                    chansInAssemType += me.getChansInAssem(i, j)
        return chansInAssemType

    def getDimLen(me):
        return len(me.chIdx)

@@ -621,14 +632,26 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
       sectionID (int): Index for the core section that will be added to the model.
       symOpt (int): Enter either 1 for full core geometry or 4 for quarter symmetry geometry.
       radialMesh (string): May be "nodal" only at this time.
       assemLosses (dict): A dictionary that takes the assembly ID and returns a list of loss coefficients in
          that assembly.  Must also provide the assemLossLevels dict if this is provided.  Alternatively,
       assemLosses (dict): A dictionary that takes the assembly type ID and returns a list of loss coefficients in
          that assembly type.  Must also provide the assemLossLevels dict if this is provided.  Alternatively,
          can also provide a single form loss coefficient if using the same value core-wide or a list of
          form loss coefficients (matching the size of assemLossLevels) if there is only one assembly type in
          the model.
       assemLossLevels (dict): A dictionary that takes the assembly ID and returns a list of levels where form
          form loss coefficients (matching the size of assemLossLevels) to apply this list to all
          assembly types.
       assemLossLevels (dict): A dictionary that takes the assembly type ID and returns a list of levels where form
          losses are provided.  Must also provide the assemLosses dict if this is provided.  Can also provide
          a list of levels (instead of a dict) if all assemblies will have losses at the same levels.
          a list of levels (instead of a dict) if all assembly types will have losses at the same levels.
       assemGridBlockageRatios (dict): A dictionary that takes the assembly type ID and returns a list of blockage
          ratios in that assembly type. Only the grid strap blockage ratio (ABLOC in Card 7.1.5) is considered at
          this time. Must also provide the assemLossLevels dict if this is provided, as the same levels will be
          used for the blockages.  Alternatively, can also provide a single blockage ratio if using the same
          value core-wide or a list of blockage ratios (matching the size of assemLossLevels) to apply this list
          to all assembly types.
       assemGridEnhancementCoeffs (dict): A dictionary that takes the assembly ID and returns a list of lists
          with Yao-Hochreiter-Leech grid heat transfer enhancement coefficients (len(assemLosses) x 4) in that
          assembly. Must also provide the assemLossLevels dict if this is provided. Alternatively, can also
          provide a single list of the 4 Yao-Hochreiter-Leech coefficients if using the same values core-wide,
          or a list of lists of coefficients (matching the size of assemLossLevels) to apply these to all assembly
          types.
       rodDomains (int): A list of lists, same dimensions as coreMap.  Gives an integer greater than or equal
          to zero.  The zeros are blank spots and must match the blank spots in coreMap.  The positive integers
          specify the domain the rod exists in.  This is only needed if modeling a parallel model.  The domains
@@ -640,8 +663,8 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
       sym_opt (int): Enter either 1 for full core geometry or 4 for quarter symmetry geometry.
    """
    def __init__(me, model, coreMap, assemMaps, pitch, dz, solidGeoms, startChanIndex=1, assemPitch=None,
                 baffleGap=0.0, sectionID=1, symOpt=1, assemLosses=None, assemLossLevels=None, rodDomains=None,
                 modelGuideTubes=False, sym_opt=1):
                 baffleGap=0.0, sectionID=1, symOpt=1, assemLosses=None, assemLossLevels=None,
                 assemGridBlockageRatios=None, assemGridEnhancementCoeffs=None, rodDomains=None, modelGuideTubes=False, sym_opt=1):
        assert(isinstance(model, Model.Model))
        assert(len(coreMap) > 0)
        assert(utils.isArraySquare(coreMap))
@@ -760,9 +783,15 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
                    # Connect it ot the model
                    model.addSolidOuterChan(solidID, chIdx, 1.0)

        if assemLosses:
        if assemLosses is not None:
            assert(assemLossLevels)
        if assemGridBlockageRatios is not None:
            # Grid enhancement parameters use assemLossLevels to define the enhancement locations, so assemLossLevels must be provided.
            assert(assemLossLevels)
        if assemLossLevels:
        if assemGridEnhancementCoeffs is not None:
            # Grid enhancement parameters use assemLossLevels to define the enhancement locations, so assemLossLevels must be provided.
            assert(assemLossLevels)
        if assemLossLevels is not None:
            assert(assemLosses)

            # Convert the passed assemLossLevels to a dictionary that takes assembly type ID and returns
@@ -779,38 +808,53 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
                    assert(assemTypeIdx in assemLossLevels)
                assemLossLevelsDict = assemLossLevels

            # Convert the passed assemLosses to a dictionary that takes assembly type ID and returns the
            # losses in the assembly.
            assert(isinstance(assemLosses, float) or isinstance(
                assemLosses, list) or isinstance(assemLosses, dict))
            if isinstance(assemLosses, float):
                assemLossesDict = {}
                for assemTypeIdx in chMapObj.getAssemTypeIds():
                    assemLossesDict[assemTypeIdx] = [assemLosses for x in range(
                        len(assemLossLevelsDict[assemTypeIdx]))]
            elif isinstance(assemLosses, list):
                assemLossesDict = {}
            assemLossesDict = me.convertAssemInputToDict(assemLosses,assemLossLevelsDict,chMapObj)
            # Grid blockage and enhancement are optional if loss levels were provided.
            if assemGridBlockageRatios is not None:
                assemGridBlockageRatiosDict = me.convertAssemInputToDict(assemGridBlockageRatios,assemLossLevelsDict,chMapObj)
            else:
                assemGridBlockageRatiosDict = {}
                for assemTypeIdx in chMapObj.getAssemTypeIds():
                    assert(len(assemLosses) == len(
                        assemLossLevelsDict[assemTypeIdx]))
                    assemLossesDict[assemTypeIdx] = assemLosses
            elif isinstance(assemLosses, dict):
                    assemGridBlockageRatiosDict[assemTypeIdx] = None
            if assemGridEnhancementCoeffs is not None:
                assemGridEnhancementCoeffsDict = me.convertAssemListInputToDict(assemGridEnhancementCoeffs,assemLossLevelsDict,chMapObj)
            else:
                assemGridEnhancementCoeffsDict ={}
                for assemTypeIdx in chMapObj.getAssemTypeIds():
                    assert(assemTypeIdx in assemLosses)
                    assert(len(assemLosses[assemTypeIdx]) == len(
                        assemLossLevelsDict[assemTypeIdx]))
                assemLossesDict = assemLosses

            for assemID in assemLossesDict.keys():
                # Loop over the core assembly map and for each match of this assemID, get the channels in that
                    assemGridEnhancementCoeffsDict[assemTypeIdx] = None

            # A unique grid type will be assumed for each loss level
            # in each assembly type. The total number of grid types
            # will be stored so that default grid enhancement
            # parameters will be written to Card group 7 if these
            # parameters were not explicitly provided.
            numGridTypes = 0
            for assemTypeID in assemLossLevelsDict.keys():
                numGridTypes += len(assemLossLevelsDict[assemTypeID])
            model.numGridTypes = numGridTypes

            for assemTypeID in assemLossesDict.keys():
                # Loop over the core assembly map and for each match of this assemTypeID, get the channels in that
                # assembly and apply the form losses.
                for i, row in enumerate(chMapObj.assemTypeIdx):
                    for j, assemTypeIdx in enumerate(row):
                        if assemTypeIdx == assemID:
                        if assemTypeIdx == assemTypeID:
                            chans = chMapObj.getChansInAssem(i, j)
                            for chIdx in chans:
                                model.addFormLosses(
                                    assemLossesDict[assemID], assemLossLevelsDict[assemID], chIdx)
                                    assemLossesDict[assemTypeID], assemLossLevelsDict[assemTypeID], chIdx)

            # Grid blockages and enhancement coeffs must be defined
            # for the same assembly type IDs.
            assert(assemGridBlockageRatiosDict.keys()==assemGridEnhancementCoeffsDict.keys())
            for assemTypeID in assemGridBlockageRatiosDict.keys():
                # Find all chans associated with this assembly type
                # and assign the grid enhancement parameters for them.
                chans = chMapObj.getChansInAssemType(assemTypeID)
                model.addGridEnhancement(assemLossLevelsDict[assemTypeID],
                                    gridBlockageRatios = assemGridBlockageRatiosDict[assemTypeID],
                                    gridEnhancementCoeffs = assemGridEnhancementCoeffsDict[assemTypeID],
                                    channels=chans)

        if rodDomains:
            assert(utils.isArraySquare(rodDomains))
@@ -849,6 +893,96 @@ class SquareLatticeLWR_Nodal(CoreBuilder.Core):
        me.mdotTotal = None
        me.tempInitial = None

    def convertAssemInputToDict(me,vals,levelsDict,chMapObj):
        # Convert the passed dict, list, or float of assembly values
        # (e.g. assemLosses, assemGridBlockageRatios) to a dictionary
        # that takes assembly type ID and returns the values in the
        # assembly.
        assert(isinstance(vals, float) or isinstance(
            vals, list) or isinstance(vals, dict))
        if isinstance(vals, float):
            valsDict = {}
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                valsDict[assemTypeIdx] = [vals for x in range(
                    len(levelsDict[assemTypeIdx]))]
        elif isinstance(vals, list):
            valsDict = {}
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                assert(len(vals) == len(
                    levelsDict[assemTypeIdx]))
                valsDict[assemTypeIdx] = vals
        elif isinstance(vals, dict):
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                assert(assemTypeIdx in vals)
                assert(len(vals[assemTypeIdx]) == len(
                    levelsDict[assemTypeIdx]))
            valsDict = vals
        return valsDict

    def convertAssemListInputToDict(me,vals,levelsDict,chMapObj):
        # Convert the passed dict, list, or list of lists containing
        # assembly values (e.g. assemGridEnhancementCoeffs) to a
        # dictionary that takes assembly type ID and returns the
        # values in the assembly. Similar to
        # convertAssemListInputToDict, but used when the basic
        # parameter being defined is a list (e.g. list of
        # coefficients) instead of a single float (e.g. loss factor)
        assert(isinstance(
            vals, list) or isinstance(vals, dict))

        if isinstance(vals, list):
            if all([isinstance(entry,float) for entry in vals]):
                # A single list was provided which contains a set of
                # values to be applied at all levels for all assembly
                # types.
                valsDict = {}
                for assemTypeIdx in chMapObj.getAssemTypeIds():
                    valsDict[assemTypeIdx] = [vals for x in range(
                        len(levelsDict[assemTypeIdx]))]
            elif all([isinstance(entry,list) for entry in vals]):
                # A list of lists was provided which contains a
                # separate set of values for each level. This will be
                # applied to all assembly types.
                valsDict = {}
                for assemTypeIdx in chMapObj.getAssemTypeIds():
                    assert(len(vals) == len(
                        levelsDict[assemTypeIdx]))
                    valsDict[assemTypeIdx] = vals
        elif isinstance(vals, dict):
            # A dict was provided which contains separate entries for
            # each assembly type.
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                assert(assemTypeIdx in vals)
                assert(len(vals[assemTypeIdx]) == len(
                    levelsDict[assemTypeIdx]))
            valsDict = vals
        return valsDict

    def convertAssemInputToDict(me,vals,levelsDict,chMapObj):
        # Convert the passed dict of assembly values (e.g. assemLosses, assemGridBlockageRatios)
        # to a dictionary that takes assembly type ID and returns the
        # values in the assembly.
        assert(isinstance(vals, float) or isinstance(
            vals, list) or isinstance(vals, dict))
        if isinstance(vals, float):
            valsDict = {}
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                valsDict[assemTypeIdx] = [vals for x in range(
                    len(levelsDict[assemTypeIdx]))]
        elif isinstance(vals, list):
            valsDict = {}
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                assert(len(vals) == len(
                    levelsDict[assemTypeIdx]))
                valsDict[assemTypeIdx] = vals
        elif isinstance(vals, dict):
            for assemTypeIdx in chMapObj.getAssemTypeIds():
                assert(assemTypeIdx in vals)
                assert(len(vals[assemTypeIdx]) == len(
                    levelsDict[assemTypeIdx]))
            valsDict = vals
        return valsDict

    def setEditOptions(me, chanVTK=None, rodVTK=None, ctfHDF5=None, veracsHDF5=None, chanASCII=None,
                       rodEdits=None, dnbEdits=None):
        """
+66 −4

File changed.

Preview size limit exceeded, changes collapsed.

Loading