Updated 'Introduction to Bamboo'. Changed defaults 'h' models.

Now uses Bartz (with sigma) for gas side and Sieder-Tate for coolant side. Have modified the Sieder-Tate implementation to avoid a boil-off problem. Has not been checked thoroughly but seems to work as of now.

bamboo/cooling.py

@@ -49,8 +49,10 @@ def black_body(T):
     return SIGMA*T**4
 
 def h_gas_1(D, M, T, rho, gamma, R, mu, k, Pr):
-    """Get the convective heat transfer coefficient on the gas side, non-Bartz equation.
-    Uses Eqn (8-22) on page 312 or RPE 7th edition (Reference [2])
+    """Get the convective heat transfer coefficient on the gas side. Uses Eqn (8-22) on page 312 or RPE 7th edition (Reference [2]). I believe this is just a form of the Dittius-Boelter equation.
+    
+    Note:
+        Seems to give much lower wall temperatures than the Bartz equation, and is likely less accurate. h_gas_2 and h_gas_3 are likely more accurate.
 
     Args:
         D (float): Flow diameter (m)
@@ -75,9 +77,6 @@ def h_gas_2(D, cp_inf, mu_inf, Pr_inf, rho_inf, v_inf, rho_am, mu_am, mu0):
     """Bartz equation, 
     using Equation (8-23) from page 312 of RPE 7th edition (Reference [2]). 'am' refers to the gas being at the 'arithmetic mean' of the wall and freestream temperatures.
 
-    Note:
-        Seems to provide questionable results - may have been implemented incorrectly.
-
     Args:
         D (float): Gas flow diameter (m)
         cp_inf (float): Specific heat capacity at constant pressure for the gas, in the freestream
@@ -233,12 +232,13 @@ def __init__(self, E, sigma_y, poisson, alpha, k, **kwargs):
         self.perf_therm = (1 - self.poisson) * self.k / (self.alpha * self.E)   #Performance coefficient for thermal stress, higher is better
 
     def __repr__(self):
-        return f"""bamboo.cooling.Material Object \nYoung's modulus = {self.E/1e9} GPa 
-0.2% Yield Stress = {self.sigma_y/1e6} MPa 
-Poisson's ratio = {self.poisson}
-alpha = {self.alpha} strain/K
-Thermal conductivity = {self.k} W/m/K
-(may also have a specific heat capacity (self.c) and density (self.rho))"""
+        return f"""bamboo.cooling.Material Object
+        Young's modulus = {self.E/1e9} GPa 
+        0.2% Yield Stress = {self.sigma_y/1e6} MPa 
+        Poisson's ratio = {self.poisson}
+        alpha = {self.alpha} strain/K
+        Thermal conductivity = {self.k} W/m/K
+        (may also have a specific heat capacity (self.c) and density (self.rho))"""
 
     def relStrength(self, T, ignoreLowTemp = False, ignoreHighTemp = False):
         """Uses polynomial coefficients to determine the fraction of yield stress
@@ -341,7 +341,6 @@ def check_liquid(self, T, p):
                 return True
             else:
                 return False
-
 
     def k(self, T, p):
         """Thermal conductivity
@@ -361,7 +360,11 @@ def k(self, T, p):
             elif self.force_phase == 'g':
                 return self.thermo_object.kg
             else:
-                return self.thermo_object.k
+                #Manually check which phase we're in, and return the right conductivity (otherwise sometimes it seems to return odd results)
+                if self.thermo_object.phase == 'l':
+                    return self.thermo_object.kl
+                elif self.thermo_object.phase == 'g':
+                    return self.thermo_object.kg
 
         elif self.model == "CoolProp":
             return PropsSI("CONDUCTIVITY", "T", T, "P", p, self.coolprop_name)
@@ -400,6 +403,7 @@ def mu(self, T, p):
 
         elif self.model == "custom":
             return self.custom_mu
+
     def Pr(self, T, p):
         """Prandtl number
 
@@ -418,7 +422,11 @@ def Pr(self, T, p):
             elif self.force_phase == 'g':
                 return self.thermo_object.Prg
             else:
-                return self.thermo_object.Pr
+                #Manually check which phase we're in, and return the right property
+                if self.thermo_object.phase == 'l':
+                    return self.thermo_object.Prl
+                elif self.thermo_object.phase == 'g':
+                    return self.thermo_object.Prg
 
         elif self.model == "CoolProp":
             return PropsSI("PRANDTL", "T", T, "P", p, self.coolprop_name)
@@ -445,7 +453,11 @@ def cp(self, T, p):
             elif self.force_phase == 'g':
                 return self.thermo_object.Cpg
             else:
-                return self.thermo_object.Cp
+                #Manually check which phase we're in, and return the right property
+                if self.thermo_object.phase == 'l':
+                    return self.thermo_object.Cpl
+                elif self.thermo_object.phase == 'g':
+                    return self.thermo_object.Cpg
 
         elif self.model == "CoolProp":
             return PropsSI("CPMASS", "T", T, "P", p, self.coolprop_name)
@@ -469,8 +481,12 @@ def rho(self, T, p):
             elif self.force_phase == 'g':
                 return self.thermo_object.rhog
             else:
-                return self.thermo_object.rho
-
+                #Manually check which phase we're in, and return the right property
+                if self.thermo_object.phase == 'l':
+                    return self.thermo_object.rhol
+                elif self.thermo_object.phase == 'g':
+                    return self.thermo_object.rhog
+
         elif self.model == "CoolProp":
             return PropsSI("DMASS", "T", T, "P", p, self.coolprop_name)
 

bamboo/main.py

@@ -1417,17 +1417,17 @@ def regen_ablative_thermal_circuit(self, r, h_gas, h_coolant, wall_material, inn
 
         return q_dot, R_gas, R_ablative, R_wall, R_coolant,
 
-    def steady_heating_analysis(self, number_of_points=1000, h_gas_model = "1", h_coolant_model = "1", to_json = "heating_output.json"):
+    def steady_heating_analysis(self, number_of_points=1000, h_gas_model = "3", h_coolant_model = "2", to_json = "heating_output.json"):
         """Steady state heating analysis. Can be used for regenarative cooling, or combined regenerative and ablative cooling.
 
         Args:
             number_of_points (int, optional): Number of discrete points to divide the engine into. Defaults to 1000.
-            h_gas_model (str, optional): Equation to use for the gas side convective heat transfer coefficients. Options are '1', '2' and '3'. Defaults to "1".
-            h_coolant_model (str, optional): Equation to use for the coolant side convective heat transfer coefficients. Options are '1', '2' and '3'. Defaults to "1".
+            h_gas_model (str, optional): Equation to use for the gas side convective heat transfer coefficients. Options are '1', '2' and '3'. Defaults to "3".
+            h_coolant_model (str, optional): Equation to use for the coolant side convective heat transfer coefficients. Options are '1', '2' and '3'. Defaults to "2".
             to_json (str or bool, optional): Directory to export a .JSON file to, containing simulation results. If False, no .JSON file is saved. Defaults to 'heating_output.json'.
 
         Note:
-            h_gas_model = '2' seems to provide questionable results (if it works at all) - use it with caution. h_coolant_model = '2' can raise errors if using the 'force_phase' setting with your coolant TransportProperties object. See the functions h_gas_1(), h_gas_2(), h_coolant_1(), etc.. in the documentation for details on each model.
+            See the bamboo.cooling module for explanations of each h_gas and h_coolant option. Defaults are Bartz (using sigma correlation) for gas side, and Sieder-Tate for coolant side. These are believed to be the most accurate.
 
         Returns:
             dict: Results of the simulation. Contains the following dictionary keys: 
@@ -1526,18 +1526,34 @@ def steady_heating_analysis(self, number_of_points=1000, h_gas_model = "1", h_co
                                             Pr_gas[i])
 
                 elif h_gas_model == "2":
-                    #We need the previous wall temperature to use h_gas_3. If we're on the first step, then just use h_gas_1()
+                    #We need the previous wall temperature to use h_gas_2. If we're on the first step, assume wall temperature = freestream temperature.
                     if i == 0:
-                        h_gas[i] = cool.h_gas_1(2*self.y(x),
-                                                self.M(x),
-                                                T_gas[i],
-                                                self.rho(x),
-                                                self.perfect_gas.gamma,
-                                                self.perfect_gas.R,
-                                                mu_gas[i],
-                                                k_gas[i],
-                                                Pr_gas[i])
-                    #Use h_gas_2() for all subsequent steps                            
+                        gamma = self.perfect_gas.gamma
+                        R = self.perfect_gas.R
+                        D = 2*self.y(x)            #Flow diameter
+
+                        #Freestream properties
+                        p_inf = self.p(x)
+                        T_inf = T_gas[i]
+                        rho_inf = self.rho(x)
+                        M_inf = self.M(x)
+                        v_inf = M_inf * (gamma*R*T_inf)**0.5    #Gas velocity
+                        mu_inf = mu_gas[i]
+                        Pr_inf = Pr_gas[i]
+                        cp_inf = self.perfect_gas.cp
+
+                        ##Properties at arithmetic mean of T_wall and T_inf. Assume wall temperature = freestream temperature for the first step.
+                        T_am = T_inf
+                        mu_am = self.exhaust_transport.mu(T = T_am, p = p_inf)
+                        rho_am = p_inf/(R*T_am)                                 #p = rho R T - pressure is roughly uniform across the boundary layer so p_inf ~= p_wall
+
+                        #Stagnation properties
+                        p0 = self.chamber_conditions.p0
+                        T0 = self.chamber_conditions.T0
+                        mu0 = self.exhaust_transport.mu(T =  T0, p = p0)
+
+                        h_gas[i] = cool.h_gas_2(D, cp_inf, mu_inf, Pr_inf, rho_inf, v_inf, rho_am, mu_am, mu0)
+
                     else:
                         gamma = self.perfect_gas.gamma
                         R = self.perfect_gas.R
@@ -1566,16 +1582,18 @@ def steady_heating_analysis(self, number_of_points=1000, h_gas_model = "1", h_co
                         h_gas[i] = cool.h_gas_2(D, cp_inf, mu_inf, Pr_inf, rho_inf, v_inf, rho_am, mu_am, mu0)
 
                 elif h_gas_model == "3":
-                    #We need the previous wall temperature to use h_gas_3. If we're on the first step, then just use h_gas_1()
+                    #We need the previous wall temperature to use h_gas_3. If we're on the first step, assume wall temperature = freestream temperature.
                     if i == 0:
-                        h_gas[i] = cool.h_gas_1(2*self.y(x),
-                                                self.M(x),
-                                                T_gas[i],
-                                                self.rho(x),
-                                                self.perfect_gas.gamma,
-                                                self.perfect_gas.R,
-                                                mu_gas[i],
-                                                k_gas[i],
+                        h_gas[i] = cool.h_gas_3(self.c_star,
+                                                self.nozzle.At, 
+                                                self.A(x), 
+                                                self.chamber_conditions.p0, 
+                                                self.chamber_conditions.T0, 
+                                                self.M(x), 
+                                                T_gas[i], 
+                                                mu_gas[i], 
+                                                self.perfect_gas.cp, 
+                                                self.perfect_gas.gamma, 
                                                 Pr_gas[i])
 
                     #Use h_gas_3() for all subsequent steps
@@ -1636,22 +1654,33 @@ def steady_heating_analysis(self, number_of_points=1000, h_gas_model = "1", h_co
                                                     rho_coolant[i])
 
                 elif h_coolant_model == "2":
-                    #Use '3' for the first step (model '2' relies on the wall temperature, which hasn't yet been calculated for the first step).
+                    #This model requires the cooling channel wall temperature, which hasn't been calculated in the first step. 
+                    #Assume the wall temperature = coolant temperature for the first step.
                     if i == 0:
-                        h_coolant[i] = cool.h_coolant_3(rho_coolant[i], 
+                        h_coolant[i] = cool.h_coolant_2(rho_coolant[i], 
                                                         v_coolant[i], 
                                                         self.cooling_jacket.D(x=x, y=self.y(x=x, up_to = "wall in")), 
                                                         mu_coolant[i], 
+                                                        self.cooling_jacket.coolant_transport.mu(T = T_coolant[i], p = p_coolant[i]), 
                                                         Pr_coolant[i], 
                                                         k_coolant[i])
 
-                    #Use model '2' for every other step.
                     else:
+                        #If the wall temperature is above the boiling temperature of the fluid, cap the wall temperature used in Sieder-Tate to the boiling temperature of the liquid.
+                        #Currently doesn't do anything for CoolProp models.
+                        if self.cooling_jacket.coolant_transport.model == "thermo":
+                            self.cooling_jacket.coolant_transport.thermo_object.calculate(P = p_coolant[i]) 
+
+                            if self.cooling_jacket.coolant_transport.thermo_object.Tb < T_wall_outer[i-1]:
+                                liquid_wall_temp = self.cooling_jacket.coolant_transport.thermo_object.Tb - 0.001 #Make it a bit smaller just to avoid thermo using the gas phase.
+                            else:
+                                liquid_wall_temp = T_wall_outer[i-1]
+
                         h_coolant[i] = cool.h_coolant_2(rho_coolant[i], 
                                                         v_coolant[i], 
                                                         self.cooling_jacket.D(x=x, y=self.y(x=x, up_to = "wall in")), 
                                                         mu_coolant[i], 
-                                                        self.cooling_jacket.coolant_transport.mu(T = T_wall_outer[i-1], p = p_coolant[i]), 
+                                                        self.cooling_jacket.coolant_transport.mu(T = liquid_wall_temp, p = p_coolant[i]), 
                                                         Pr_coolant[i], 
                                                         k_coolant[i])
 
@@ -1707,7 +1736,7 @@ def steady_heating_analysis(self, number_of_points=1000, h_gas_model = "1", h_co
 
         #Not sure if the Sieder-Tate equation is valid with your coolant above the boiling temperature at the wall.
         if h_coolant_model == '2' and self.cooling_jacket.coolant_transport.check_liquid(T = np.amax(T_wall_outer[i]), p = p_coolant[-1]) == False:
-            print("Coolant temperature at the wall was above its boiling point when using the Sieder-Tate equation (h_coolant_model = '2') - results should be used with caution.")
+            print("Wall tempreature was above coolant boiling point when using the Sieder-Tate equation (h_coolant_model = '2') - coolant boiling temperature was used instead of wall temperature.")
 
         #Dictionary containing results
         output_dict = {"x" : list(discretised_x),

docs/_autosummary/bamboo.cooling.html

@@ -199,7 +199,7 @@
 <td><p>Dittus-Boelter equation for convective heat transfer coefficient.</p></td>
 </tr>
 <tr class="row-odd"><td><p><a class="reference internal" href="../index.html#bamboo.cooling.h_gas_1" title="bamboo.cooling.h_gas_1"><code class="xref py py-obj docutils literal notranslate"><span class="pre">h_gas_1</span></code></a>(D, M, T, rho, gamma, R, mu, k, Pr)</p></td>
-<td><p>Get the convective heat transfer coefficient on the gas side, non-Bartz equation.</p></td>
+<td><p>Get the convective heat transfer coefficient on the gas side.</p></td>
 </tr>
 <tr class="row-even"><td><p><a class="reference internal" href="../index.html#bamboo.cooling.h_gas_2" title="bamboo.cooling.h_gas_2"><code class="xref py py-obj docutils literal notranslate"><span class="pre">h_gas_2</span></code></a>(D, cp_inf, mu_inf, Pr_inf, rho_inf, …)</p></td>
 <td><p>Bartz equation, using Equation (8-23) from page 312 of RPE 7th edition (Reference [2]).</p></td>
@@ -219,7 +219,7 @@
 <tr class="row-odd"><td><p><a class="reference internal" href="../index.html#bamboo.cooling.Ablative" title="bamboo.cooling.Ablative"><code class="xref py py-obj docutils literal notranslate"><span class="pre">Ablative</span></code></a>(ablative_material, wall_material[, …])</p></td>
 <td><p>Container for refractory or ablative properties.</p></td>
 </tr>
-<tr class="row-even"><td><p><a class="reference internal" href="../index.html#bamboo.cooling.CoolingJacket" title="bamboo.cooling.CoolingJacket"><code class="xref py py-obj docutils literal notranslate"><span class="pre">CoolingJacket</span></code></a>(geometry, inner_wall, …[, …])</p></td>
+<tr class="row-even"><td><p><a class="reference internal" href="../index.html#bamboo.cooling.CoolingJacket" title="bamboo.cooling.CoolingJacket"><code class="xref py py-obj docutils literal notranslate"><span class="pre">CoolingJacket</span></code></a>(geometry, inner_wall, inlet_T, …)</p></td>
 <td><p>Container for cooling jacket information - e.g.</p></td>
 </tr>
 <tr class="row-odd"><td><p><a class="reference internal" href="../index.html#bamboo.cooling.Material" title="bamboo.cooling.Material"><code class="xref py py-obj docutils literal notranslate"><span class="pre">Material</span></code></a>(E, sigma_y, poisson, alpha, k, **kwargs)</p></td>