Fixing Component density to agree with 3D density

armi/materials/zr.py

@@ -130,7 +130,7 @@ def linearExpansionPercent(self, Tk=None, Tc=None):
         Tk = getTk(Tc, Tk)
         self.checkPropertyTempRange("linear expansion percent", Tk)
 
-        if Tk >= 291.62 and Tk < 1137:
+        if Tk >= 293 and Tk < 1137:
             return (
                 -0.111 + (2.325e-4 * Tk) + (5.595e-7 * Tk ** 2) - (1.768e-10 * Tk ** 3)
             )

armi/reactor/blueprints/tests/test_blockBlueprints.py

@@ -321,9 +321,6 @@ def test_explicitFlags(self):
 
     # TODO: This test passes, but shouldn't.
     def test_densityConsistentWithComponentConstructor(self):
-        # when comparing to 3D density, the comparison is not quite correct.
-        # We need a bigger delta, this will be investigated/fixed in another PR
-        biggerDelta = 0.001  # g/cc
         a1 = self.blueprints.assemDesigns.bySpecifier["IC"].construct(
             self.cs, self.blueprints
         )
@@ -336,7 +333,6 @@ def test_densityConsistentWithComponentConstructor(self):
         self.assertAlmostEqual(
             clad.getMassDensity(),
             clad.material.density3(Tc=clad.temperatureInC),
-            delta=biggerDelta,
         )
 
         self.assertAlmostEqual(

armi/reactor/components/component.py

@@ -232,27 +232,11 @@ def __init__(
         self.temperatureInC = Thot
         self.material = None
         self.setProperties(material)
-        self.tInputWarning(Tinput)
         self.applyMaterialMassFracsToNumberDensities()  # not necessary when duplicating...
         self.setType(name)
         self.p.mergeWith = mergeWith
         self.p.customIsotopicsName = isotopics
 
-    def tInputWarning(self, Tinput):
-        """
-        Check whether thermal expansion factor is 0.0% exactly at T=Tinput
-        """
-        expansionFactor = (
-            self.material.linearExpansionPercent(Tc=self.inputTemperatureInC) / 100.0
-        )
-        if not (abs(expansionFactor) < 1.0e-6):
-            runLog.warning(
-                f"Thermal expansion for {self.material} at Tinput = {Tinput} is non-zero "
-                f"({expansionFactor}). The modeled density for this material will be off "
-                f"by a factor of {(1 + expansionFactor) ** 2}.",
-                single=True,
-            )
-
     @property
     def temperatureInC(self):
         """Return the hot temperature in Celsius."""
@@ -353,42 +337,50 @@ def applyMaterialMassFracsToNumberDensities(self):
           due to the difference in self.inputTemperatureInC and self.temperatureInC
         - After the expansion, the density of the component should reflect the 3d
           density of the material
-
-        See Also
-        --------
-        self.applyHotHeightDensityReduction
         """
         # note, that this is not the actual material density, but rather 2D expanded
         # `density3` is 3D density
+        # call getProperty to cache and improve speed
         density = self.material.getProperty("density", Tc=self.temperatureInC)
 
         self.p.numberDensities = densityTools.getNDensFromMasses(
             density, self.material.p.massFrac
         )
-        self.applyHotHeightDensityReduction()
 
-    def applyHotHeightDensityReduction(self):
+        # material needs to be expanded from the material's cold temp to hot,
+        # not components cold temp, so we don't use mat.linearExpansionFactor or
+        # component.getThermalExpansionFactor.
+        # Materials don't typically define the temperature for which their references
+        # density is defined so linearExpansionPercent must be called
+        coldMatAxialExpansionFactor = (
+            1.0 + self.material.linearExpansionPercent(Tc=self.temperatureInC) / 100
+        )
+        self.changeNDensByFactor(1.0 / coldMatAxialExpansionFactor)
+
+    def adjustDensityForHeightExpansion(self, newHot):
         """
-        Adjust number densities to account for hot block heights (axial expansion)
-        (crucial for preserving 3D density).
+        Change the densities in cases where height of the block/component is changing with expansion.
 
         Notes
         -----
-        - We apply this hot height density reduction to account for pre-expanded
-          block heights in blueprints.
-        - This is called when inputHeightsConsideredHot: True.
+        Call before setTemperature since we need old hot temp.
+        This works well if there is only 1 solid component.
+        If there are multiple components expanding at different rates during thermal
+        expansion this becomes more complicated and, and axial expansion should be used.
+        Multiple expansion rates cannot trivially be accommodated.
+        See AxialExpansionChanger.
+        """
+        self.changeNDensByFactor(1.0 / self.getHeightFactor(newHot))
 
-        See Also
-        --------
-        self.applyMaterialMassFracsToNumberDensities
+    def getHeightFactor(self, newHot):
         """
-        # this is the same as getThermalExpansionFactor but doesn't fail
-        # on non-fluid materials that have 0 or undefined thermal expansion
-        # (we don't want materials to fail on  __init__ which calls this)
-        axialExpansionFactor = 1.0 + self.material.linearExpansionFactor(
-            self.temperatureInC, self.inputTemperatureInC
-        )
-        self.changeNDensByFactor(1.0 / axialExpansionFactor)
+        Return the factor by which height would change by if we did 3D expansion.
+
+        Notes
+        -----
+        Call before setTemperature since we need old hot temp.
+        """
+        return self.getThermalExpansionFactor(Tc=newHot, T0=self.temperatureInC)
 
     def getProperties(self):
         """Return the active Material object defining thermo-mechanical properties."""

armi/reactor/converters/axialExpansionChanger.py

@@ -152,8 +152,8 @@ def applyColdHeightMassIncrease(self):
         -----
         A cold 1 cm tall component will have more mass that a component with the
         same mass/length as a component with a hot height of 1 cm. This should be
-        called when the setting `inputHeightsConsideredHot` is used. This basically
-        undoes component.applyHotHeightDensityReduction
+        called when the setting `inputHeightsConsideredHot` is used. This adjusts
+        the expansion factor applied during applyMaterialMassFracsToNumberDensities.
         """
         for c in self.linked.a.getComponents():
             axialExpansionFactor = 1.0 + c.material.linearExpansionFactor(

armi/reactor/tests/test_blocks.py

@@ -2274,11 +2274,16 @@ def test_coldMass(self):
         # set ref (input/cold) temperature.
         Thot = fuel.temperatureInC
         Tcold = fuel.inputTemperatureInC
+
+        # change temp to cold
         fuel.setTemperature(Tcold)
         massCold = fuel.getMass()
         fuelArea = fuel.getArea()
+        # we are at cold temp so cold and hot area are equal
+        self.assertAlmostEqual(fuel.getArea(cold=True), fuel.getArea())
         height = self.b.getHeight()  # hot height.
-        rho = fuel.getProperties().density(Tc=Tcold)
+        rho = fuel.getProperties().density3(Tc=Tcold)
+        # can't use getThermalExpansionFactor since hot=cold so it would be 0
         dllHot = fuel.getProperties().linearExpansionFactor(Tc=Thot, T0=Tcold)
         coldHeight = height / (1 + dllHot)
         theoreticalMass = fuelArea * coldHeight * rho
@@ -2287,7 +2292,7 @@ def test_coldMass(self):
             massCold,
             theoreticalMass,
             7,
-            "Cold mass of fuel ({0}) != theoretical mass {1}.  "
+            msg="Cold mass of fuel ({0}) != theoretical mass {1}.  "
             "Check calculation of cold mass".format(massCold, theoreticalMass),
         )
 

armi/reactor/tests/test_components.py

@@ -499,9 +499,22 @@ def test_changeNumberDensities(self):
 class TestComponentExpansion(unittest.TestCase):
     # when comparing to 3D density, the comparison is not quite correct.
     # We need a bigger delta, this will be investigated/fixed in another PR
-    biggerDelta = 0.01  # g/cc
     tWarm = 50
     tHot = 500
+    coldOuterDiameter = 1.0
+
+    def test_ComponentMassIndependentOfInputTemp(self):
+        coldT = 20
+        circle1 = Circle("circle", "HT9", coldT, self.tHot, self.coldOuterDiameter)
+        coldT += 200
+        # pick the input dimension to get the same hot component
+        hotterDim = self.coldOuterDiameter * (
+            1 + circle1.material.linearExpansionFactor(coldT, 20)
+        )
+        circle2 = Circle("circle", "HT9", coldT, self.tHot, hotterDim)
+        self.assertAlmostEqual(circle1.getDimension("od"), circle2.getDimension("od"))
+        self.assertAlmostEqual(circle1.getArea(), circle2.getArea())
+        self.assertAlmostEqual(circle1.getMassDensity(), circle2.getMassDensity())
 
     def test_ExpansionConservationHotHeightDefined(self):
         """
@@ -513,9 +526,9 @@ def test_ExpansionConservationHotHeightDefined(self):
           inputHeightsConsideredHot = True (the default)
         """
         hotHeight = 1.0
-        coldOuterDiameter = 1.0
-        circle1 = Circle("circle", "HT9", 20, self.tWarm, coldOuterDiameter)
-        circle2 = Circle("circle", "HT9", 20, self.tHot, coldOuterDiameter)
+
+        circle1 = Circle("circle", "HT9", 20, self.tWarm, self.coldOuterDiameter)
+        circle2 = Circle("circle", "HT9", 20, self.tHot, self.coldOuterDiameter)
 
         # mass density is proportional to Fe number density and derived from
         # all the number densities and atomic masses
@@ -538,49 +551,50 @@ def test_ExpansionConservationHotHeightDefined(self):
 
         # material.density is the 2D density of a material
         # material.density3 is true density and not equal in this case
-        # density must be density by calling applyHotHeightDensityReduction
-        # or other methods (see rest of test).
         for circle in [circle1, circle2]:
-            # 2D density is not equal after application of applyHotHeightDensityReduction
+            # 2D density is not equal after application of coldMatAxialExpansionFactor
             # which happens during construction
             self.assertNotAlmostEqual(
                 circle.getMassDensity(),
                 circle.material.density(Tc=circle.temperatureInC),
             )
-            # 2D density is off by the thermal exp factor
+            # 2D density is off by the material thermal exp factor
+            percent = circle.material.linearExpansionPercent(Tc=circle.temperatureInC)
+            thermalExpansionFactorFromColdMatTemp = 1 + percent / 100
             self.assertAlmostEqual(
-                circle.getMassDensity() * circle.getThermalExpansionFactor(),
+                circle.getMassDensity() * thermalExpansionFactorFromColdMatTemp,
                 circle.material.density(Tc=circle.temperatureInC),
             )
             self.assertAlmostEqual(
                 circle.getMassDensity(),
                 circle.material.density3(Tc=circle.temperatureInC),
-                delta=self.biggerDelta,
             )
-        # Change temp forward and backward and show equal
+
+        # brief 2D expansion with set temp to show mass is conserved
+        # hot height would come from block value
+        warmMass = circle1.getMassDensity() * circle1.getArea() * hotHeight
+        circle1.setTemperature(self.tHot)
+        hotMass = circle1.getMassDensity() * circle1.getArea() * hotHeight
+        self.assertAlmostEqual(warmMass, hotMass)
+        circle1.setTemperature(self.tWarm)
+
+        # Change temp forward and backward with axial expansion and show equal
         oldArea = circle1.getArea()
         initialDens = circle1.getMassDensity()
-        # this math is done in applyHotHeightDensityReduction
-        applyHotHeightDensityReductionFactor = (
-            1.0
-            + circle1.material.linearExpansionFactor(
-                circle1.temperatureInC, circle1.inputTemperatureInC
-            )
+        # this math is done in applyMaterialMassFracsToNumberDensities
+        coldMatAxialExpansionFactor = 1.0 + circle1.material.linearExpansionFactor(
+            circle1.temperatureInC, circle1.inputTemperatureInC
         )
         factorToUndoHotHeight = circle1.getThermalExpansionFactor()
-        self.assertAlmostEqual(
-            applyHotHeightDensityReductionFactor,
-            factorToUndoHotHeight,
-        )
         # when block.setHeight is called (which effectively changes component height)
         # component.setNumberDensity is called (for solid isotopes) to adjust the number
         # density so that now the 2D expansion will be approximated/expanded around
         # the hot temp which is akin to these adjustments
 
-        # undo the old applyHotHeightDensityReduction
-        circle1.changeNDensByFactor(factorToUndoHotHeight)
+        # change to self.THot
+        heightFactor = circle1.getHeightFactor(self.tHot)
+        circle1.adjustDensityForHeightExpansion(self.tHot)  # apply temp at new height
         circle1.setTemperature(self.tHot)
-        circle1.applyHotHeightDensityReduction()  # apply at new temp
 
         # now its density is same as hot component
         self.assertAlmostEqual(
@@ -589,27 +603,18 @@ def test_ExpansionConservationHotHeightDefined(self):
         )
 
         # show that mass is conserved after expansion
-        circle1NewHotHeight = (
-            hotHeight * circle1.getThermalExpansionFactor() / factorToUndoHotHeight
-        )
+        circle1NewHotHeight = hotHeight * heightFactor
         self.assertAlmostEqual(
             mass1, circle1.getMassDensity() * circle1.getArea() * circle1NewHotHeight
         )
-        # you can calculate the height exp factor directly this way
-        self.assertAlmostEqual(
-            circle1.getThermalExpansionFactor() / factorToUndoHotHeight,
-            circle1.getThermalExpansionFactor(Tc=circle1.temperatureInC, T0=self.tWarm),
-        )
 
         self.assertAlmostEqual(
             circle1.getMassDensity(),
-            circle1.material.density3(Tc=circle2.temperatureInC),
-            delta=self.biggerDelta,
+            circle1.material.density3(Tc=circle1.temperatureInC),
         )
         # change back to old temp
-        circle1.changeNDensByFactor(circle1.getThermalExpansionFactor())
+        circle1.adjustDensityForHeightExpansion(self.tWarm)
         circle1.setTemperature(self.tWarm)
-        circle1.applyHotHeightDensityReduction()
 
         # check for consistency
         self.assertAlmostEqual(initialDens, circle1.getMassDensity())
@@ -628,17 +633,14 @@ def test_ExpansionConservationColdHeightDefined(self):
           inputHeightsConsideredHot = False
         """
         coldHeight = 1.0
-        circle1 = Circle("circle", "HT9", 20, self.tWarm, 1.0)
-        circle2 = Circle("circle", "HT9", 20, self.tHot, 1.0)
+        circle1 = Circle("circle", "HT9", 20, self.tWarm, self.coldOuterDiameter)
+        circle2 = Circle("circle", "HT9", 20, self.tHot, self.coldOuterDiameter)
         # same as 1 but we will make like 2
         circle1AdjustTo2 = Circle("circle", "HT9", 20, self.tWarm, 1.0)
 
         # make it hot like 2
-        circle1AdjustTo2.changeNDensByFactor(
-            circle1AdjustTo2.getThermalExpansionFactor()
-        )
+        circle1AdjustTo2.adjustDensityForHeightExpansion(self.tHot)
         circle1AdjustTo2.setTemperature(self.tHot)
-        circle1AdjustTo2.applyHotHeightDensityReduction()
         # check that its like 2
         self.assertAlmostEqual(
             circle2.getMassDensity(), circle1AdjustTo2.getMassDensity()
@@ -650,7 +652,6 @@ def test_ExpansionConservationColdHeightDefined(self):
             self.assertAlmostEqual(
                 circle.getMassDensity(),
                 circle.material.density3(Tc=circle.temperatureInC),
-                delta=self.biggerDelta,
             )
             # total mass consistent between hot and cold
             # Hot height will be taller
@@ -660,7 +661,6 @@ def test_ExpansionConservationColdHeightDefined(self):
                 * circle.getArea(cold=True)
                 * circle.material.density3(Tc=circle.inputTemperatureInC),
                 hotHeight * circle.getArea() * circle.getMassDensity(),
-                delta=self.biggerDelta,
             )