v0.6.3: fixed bug in low-fidelity model region energy balance with mesh disagreement and in gap subchannel adjacency map; added tests to demonstrate core-wide ebal with low-fidelity models

Author: Atz

dassh/__init__.py

@@ -14,7 +14,7 @@
 # permissions and limitations under the License.
 ########################################################################
 """
-date: 2021-05-25
+date: 2021-05-26
 author: Milos Atz, Micheal Smith
 """
 ########################################################################
@@ -44,4 +44,4 @@
 np.set_printoptions(threshold=sys.maxsize, linewidth=500)
 
 
-__version__ = '0.6.2-dev'
+__version__ = '0.6.3-dev'

dassh/core.py

@@ -14,7 +14,7 @@
 # permissions and limitations under the License.
 ########################################################################
 """
-date: 2021-05-03
+date: 2021-05-26
 author: matz
 Methods to describe the layout of assemblies in the reactor core and
 the coolant in the gap between them
@@ -189,6 +189,7 @@ def load(self, asms):
         self.gap_params['asm wp'] = self._calculate_asm_sc_wp()
         self.gap_params['area'] = self._calculate_sc_area()
         self.gap_params['L'] = self._calculate_dist_between_sc()
+        assert not np.any(0 in np.sum(self.gap_params['L'], axis=1))
         self.gap_params['de'] = self._calculate_sc_de()
 
         # Core-total parameters
@@ -318,6 +319,10 @@ def _map_asm_gap_adjacency(self):
             # Add the temporary arrays to main arrays
             asm_adj_sc.append(tmp_asm_adj_sc)
 
+        # Cast all value as int
+        for i in range(len(asm_adj_sc)):
+            for j in range(len(asm_adj_sc[i])):
+                asm_adj_sc[i][j] = [int(x) for x in asm_adj_sc[i][j]]
         return asm_adj_sc
 
     def _calculate_gap_xbnds(self):
@@ -911,7 +916,8 @@ def _find_adjacent_sc(self, asm_sc_adj):
         for ai in range(len(asm_sc_adj)):
             asm_sc = asm_sc_adj[ai]
             for side in range(len(asm_sc)):
-                for sci in range(len(asm_sc[side]) - 1):
+                # for sci in range(len(asm_sc[side]) - 1):
+                for sci in range(len(asm_sc[side])):
                     # Look to trailing corner on previous side
                     if sci == 0:
                         # Fill trailing corner's value into active index
@@ -931,15 +937,16 @@ def _find_adjacent_sc(self, asm_sc_adj):
                                 idx = np.where(sc_adj[sc - 1] == 0)[0][0]
                             sc_adj[sc - 1, idx] = asm_sc[side][sci]
                     # For the sc in current index: map the sc in next index
-                    sc = asm_sc[side][sci]
-                    if asm_sc[side][sci + 1] not in sc_adj[sc - 1]:
-                        idx = np.where(sc_adj[sc - 1] == 0)[0][0]
-                        sc_adj[sc - 1, idx] = asm_sc[side][sci + 1]
-                    # For the sc in next index; map the sc in current index
-                    sc = asm_sc[side][sci + 1]
-                    if asm_sc[side][sci] not in sc_adj[sc - 1]:
-                        idx = np.where(sc_adj[sc - 1] == 0)[0][0]
-                        sc_adj[sc - 1, idx] = asm_sc[side][sci]
+                    if sci < len(asm_sc[side]) - 1:
+                        sc = asm_sc[side][sci]
+                        if asm_sc[side][sci + 1] not in sc_adj[sc - 1]:
+                            idx = np.where(sc_adj[sc - 1] == 0)[0][0]
+                            sc_adj[sc - 1, idx] = asm_sc[side][sci + 1]
+                        # For the sc in next index; map the sc in current index
+                        sc = asm_sc[side][sci + 1]
+                        if asm_sc[side][sci] not in sc_adj[sc - 1]:
+                            idx = np.where(sc_adj[sc - 1] == 0)[0][0]
+                            sc_adj[sc - 1, idx] = asm_sc[side][sci]
         return sc_adj
 
     def _calculate_sc_wp(self):

dassh/reactor.py

@@ -1280,7 +1280,7 @@ def _calculate_asm_temperatures(self, asm, i, z, dz, dump_step):
         if self.core.model is None:
             gap_temp = np.ones(asm.duct_outer_surf_temp.shape[0])
             gap_htc = np.ones(asm.duct_outer_surf_temp.shape[0])
-        elif asm.active_region.is_rodded:
+        # elif asm.active_region.is_rodded:
             # gap_temp = dassh.mesh_functions.interpolate_quad(
             #     self.core.x_pts,
             #     self.core.adjacent_coolant_gap_temp(i),
@@ -1291,6 +1291,7 @@ def _calculate_asm_temperatures(self, asm, i, z, dz, dump_step):
             #     self.core.x_pts,
             #     self.core.adjacent_coolant_gap_temp(i),
             #     asm.x_pts)
+        else:
             gap_htc = dassh.mesh_functions.map_across_gap(
                 self.core.adjacent_coolant_gap_htc(i),
                 asm.active_region._map['gap2duct'])
@@ -1299,20 +1300,20 @@ def _calculate_asm_temperatures(self, asm, i, z, dz, dump_step):
                  * self.core.adjacent_coolant_gap_htc(i)),
                 asm.active_region._map['gap2duct'])
             gap_temp = gap_temp / gap_htc
-        else:
-            gap_htc = dassh.mesh_functions.map_across_gap(
-                self.core.adjacent_coolant_gap_htc(i),
-                asm.active_region._map['gap2duct'])
-            gap_temp = dassh.mesh_functions.map_across_gap(
-                self.core.adjacent_coolant_gap_temp(i),
-                asm.active_region._map['gap2duct']) / gap_htc
-            # poop = self.core.adjacent_coolant_gap_temp(i).reshape(6, -1)
-            # gap_temp = gap_temp.reshape(6, -1)
-            # gap_temp[:, -1] = poop[:, -1]
-            # gap_temp = gap_temp.flatten()
-            # gap_adj_temps = map_across_gap(
-            #     self.core.adjacent_coolant_gap_temp(i),
-            #     asm.active_region._map['gap2duct'])
+        # else:
+        #     gap_htc = dassh.mesh_functions.map_across_gap(
+        #         self.core.adjacent_coolant_gap_htc(i),
+        #         asm.active_region._map['gap2duct'])
+        #     gap_temp = dassh.mesh_functions.map_across_gap(
+        #         self.core.adjacent_coolant_gap_temp(i),
+        #         asm.active_region._map['gap2duct']) / gap_htc
+        #     # poop = self.core.adjacent_coolant_gap_temp(i).reshape(6, -1)
+        #     # gap_temp = gap_temp.reshape(6, -1)
+        #     # gap_temp[:, -1] = poop[:, -1]
+        #     # gap_temp = gap_temp.flatten()
+        #     # gap_adj_temps = map_across_gap(
+        #     #     self.core.adjacent_coolant_gap_temp(i),
+        #     #     asm.active_region._map['gap2duct'])
         asm.calculate(z, dz, gap_temp, gap_htc,
                       self._is_adiabatic, self._options['ebal'])
         # Write the results

dassh/region_unrodded.py

@@ -14,7 +14,7 @@
 # permissions and limitations under the License.
 ########################################################################
 """
-date: 2021-05-20
+date: 2021-05-26
 author: matz
 Methods for unrodded axial regions; to be used within Assembly objects
 """
@@ -147,7 +147,8 @@ def __init__(self, z_lo, z_hi, duct_ftf, vf_cool, flow_rate,
                     range(0, len(duct_ftf), 2)]
         self.duct_ftf = duct_ftf[-1]
         self.duct_thickness = (self.duct_ftf[1] - self.duct_ftf[0]) / 2
-        self.duct_perim = self.duct_ftf[0] * 6 / _sqrt3
+        # self.duct_perim = self.duct_ftf[0] * 6 / _sqrt3
+        self.duct_perim = self.duct_ftf[1] * 6 / _sqrt3
         # One sixth of the duct perimeter for the coolant temp calc
         # (self.duct_perimeter defined in SingleNodeHomogeneous)
         self.duct_perim_over_6 = self.duct_perim / 6
@@ -357,11 +358,18 @@ def calculate(self, dz, power, t_gap, htc_gap,
         None
 
         """
+        # Interior coolant temperatures: calculate using coolant
+        # properties from previous axial step
         self.temp['coolant_int'] += self._calc_coolant_temp(
             dz, power, adiabatic_duct, ebal)
+
+        # Update coolant properties for the duct wall calculation
+        self._update_coolant_params(self.temp['coolant_int'][0])
+
+        # Duct temperatures: calculate with new coolant properties
         self._calc_duct_temp(t_gap, htc_gap, adiabatic_duct)
 
-        # Update pressure drop
+        # Update pressure drop (now that correlations are updated)
         self._pressure_drop += self.calculate_pressure_drop(dz)
 
     def activate2(self, previous_reg, t_gap, h_gap, adiabatic):
@@ -412,8 +420,7 @@ def _calc_coolant_temp(self, dz, power, adiabatic=False, ebal=False):
 
         """
         # Update coolant parameters and material properties
-        Tj = self.avg_coolant_int_temp  # <-- dont wanna calc everytim
-        self._update_coolant_params(Tj)
+        # Tj = self.avg_coolant_int_temp  # <-- dont wanna calc everytim
 
         # Calculate change in temperature from heat generation
         dT = power['refl']
@@ -425,9 +432,11 @@ def _calc_coolant_temp(self, dz, power, adiabatic=False, ebal=False):
                      / self.duct.thermal_conductivity)
                 R += 1 / self.coolant_params['htc']
                 dT_duct = (self.duct_perim_over_6 / R
-                           * (self.temp['duct_mw'][0] - Tj))
+                           * (self.temp['duct_mw'][0]
+                              - self.temp['coolant_int'][0]))
             else:
-                dT_duct = ((self.temp['duct_surf'][0, 0] - Tj)
+                dT_duct = ((self.temp['duct_surf'][0, 0]
+                            - self.temp['coolant_int'][0])
                            * self.coolant_params['htc']
                            * self.duct_perim_over_6)
 
@@ -478,17 +487,24 @@ def _calc_duct_temp(self, temp_gap, htc_gap, adiabatic=False):
 
         else:
             # Update coolant parameters and material properties
-            self._update_coolant_params(self.avg_coolant_int_temp)
+            # self._update_coolant_params(self.avg_coolant_int_temp)
             self.duct.update(self.avg_duct_mw_temp[0])
             L_over_2 = self.duct_thickness / 2
+            # c1 = (self.coolant_params['htc']
+            #       * (temp_gap - self.temp['coolant_int'])
+            #       / (self.coolant_params['htc'] * self.duct_thickness
+            #          + (self.duct.thermal_conductivity
+            #             * (1 + self.coolant_params['htc'] / htc_gap[1]))))
+            # c2 = (temp_gap
+            #       - c1 * (L_over_2
+            #               + self.duct.thermal_conductivity / htc_gap[1]))
             c1 = (self.coolant_params['htc']
                   * (temp_gap - self.temp['coolant_int'])
                   / (self.coolant_params['htc'] * self.duct_thickness
                      + (self.duct.thermal_conductivity
-                        * (1 + self.coolant_params['htc'] / htc_gap[1]))))
+                        * (1 + self.coolant_params['htc'] / htc_gap))))
             c2 = (temp_gap
-                  - c1 * (L_over_2
-                          + self.duct.thermal_conductivity / htc_gap[1]))
+                  - c1 * (L_over_2 + self.duct.thermal_conductivity / htc_gap))
             self.temp['duct_mw'][0] = c2
             self.temp['duct_surf'][0, 0] = c1 * -L_over_2 + c2
             self.temp['duct_surf'][0, 1] = c1 * L_over_2 + c2

tests/test_assembly.py

@@ -14,7 +14,7 @@
 # permissions and limitations under the License.
 ########################################################################
 """
-date: 2021-04-02
+date: 2021-05-26
 author: matz
 Test the DASSH Assembly object
 """
@@ -73,7 +73,8 @@ def test_temperature_unrodded_asm(c_shield_asm):
     # Try calculating temperature; shouldn't fail
     dz = 0.005
     tgap = np.ones(6) * 623.15
-    htc_gap = [1e5, 1e5]
+    htc_gap = np.ones(6) * 1e5
+    c_shield_asm.active_region._update_coolant_params(623.15)
     c_shield_asm.calculate(dz, dz, tgap, htc_gap)
     assert c_shield_asm.avg_coolant_temp > 623.15
 

tests/test_inputs/VARPOW.out

@@ -0,0 +1,91 @@
+########################################################################
+# author: Milos Atz
+# date: 2021-05-26
+# comment: DASSH input example (2 rings, but only 2 assemblies)
+########################################################################
+# SETUP: Problem setup, user options, etc
+
+[Setup]
+    [[Units]]
+        temperature    = Celsius
+        length         = cm
+        mass_flow_rate = kg/s
+    [[Options]]
+        calc_energy_balance = True
+    [[Dump]]
+        coolant  = True
+        gap      = True
+        gap_fine = True
+        duct     = True
+
+
+########################################################################
+ # NEUTRONICS AND GEOMETRY FILES
+
+[ARC]
+    pmatrx = ../test_data/seven_asm_vac/PMATRX
+    geodst = ../test_data/seven_asm_vac/GEODST
+    ndxsrf = ../test_data/seven_asm_vac/NDXSRF
+    znatdn = ../test_data/seven_asm_vac/ZNATDN
+    labels = ../test_data/seven_asm_vac/LABELS
+    nhflux = ../test_data/seven_asm_vac/NHFLUX
+    ghflux = ../test_data/seven_asm_vac/GHFLUX
+
+
+########################################################################
+# CORE DETAILS: Describe characteristics of the reactor core
+
+[Core]
+    coolant_inlet_temp   = 350.0
+    fuel_material        = metal
+    fuel_alloy           = zr
+    coolant_heating      = sodium
+    coolant_material     = sodium_se2anl_425
+    length               = 375.0
+    assembly_pitch       = 12.0
+    bypass_fraction      = 0.1
+    gap_model            = flow
+    power_model          = pin_only
+    total_power          = 0.0
+
+
+########################################################################
+# ASSEMBLY DETAILS: Describe a group of assemblies
+# NAME identifies a group of assemblies (e.g. "INNER_DRIVER")
+# Can be repeated for as many assemblies as required by the user.
+
+[Assembly]
+
+    [[fuel]]
+        num_rings      = 3
+        pin_pitch      = 2.4200
+        pin_diameter   = 2.0000
+        clad_thickness = 0.0500
+        wire_pitch     = 20.000
+        wire_diameter  = 0.4100
+        duct_ftf       = 11.25, 11.75
+        duct_material  = HT9_se2anl_425
+
+    [[other]]
+        num_rings      = 5
+        pin_pitch      = 1.37
+        pin_diameter   = 1.13
+        clad_thickness = 0.05
+        wire_pitch     = 20.0
+        wire_diameter  = 0.24
+        duct_ftf       = 11.25, 11.75
+        duct_material  = HT9_se2anl_425
+        use_low_fidelity_model = True
+        convection_factor      = 1.0
+
+
+########################################################################
+# ASSIGN assemblies to positions in the core; assign fixed
+# flow rates or temperature limits to one or multiple assemblies.
+# Indicate which assemblies must be grouped together in an orificing
+# calculation.
+
+[Assignment]
+    [[ByPosition]]
+        fuel  = 1,  1,  1, FLOWRATE = 1.0
+        other = 2,  1,  1, FLOWRATE = 1.0