Rupture calc

Manual.md

@@ -1000,6 +1000,65 @@ Both ASME F&D, DIN and 2:1 semielliptical are variants of a torispherical vessel
 Using the *fluids* library partial volumes, surface area (full and partial) and liquid level (from partial volume) can be calculated and used internally in *openthermo*.
 
 
+
+## Rupture evaluation 
+It is assumed that when the von Mises stress, $\sigma_e$ (MPa), exceeds the allowable tensile strength of the material, (ATS) (MPa), rupture will occur, i.e., when $\sigma_e$ > ATS.
+
+The ATS is calculated as [@scandpower]:
+
+$$ ATS = UTS  k_s k_y $$
+
+Where
+
+- $UTS$ is the material Ultimate Tensike Strength
+- $k_s$ is a general safety factor for a specific material with known material data. If typical material specific data is applied a factor of 0.85 is recommended.
+- $k_y$ is an additional factor used for material with missing or uncertain data. Normally this factor is 1.0.
+
+The von Mises stress is calculated by:
+
+$$ \sigma_e = \sqrt{3 \left(\frac{p D^2}{D^2 - d^2} \right)^2 +\sigma_a^2} $$
+
+Where 
+
+- $p$ is the pressure (MPa)
+- $D$ is the vessel/pipe external diameter
+- $d$ is the vessel/pipe internal diameter
+- $\sigma_a$ is the longitudinal stress due to the external force. Assumed to be 30 MPa [@scandpower]
+
+The evaluation of vessel rupture is performed as a post-calculation step following the actual depressurisation calculation. The depressurisation calculation is performed with the applicable back-ground heat load to generate the time dependent pressure profile of the vessel inventory. The background heat load is depending on the fire type as also summarised in [@Tbl:heatfluxes1]. During the depressurisation calaculation the internal heat flux is calculated. In the post-calculation step an energy balance is made for the vessel material, with the external heat is generated by the applicable peak heat load Stefan-Boltzmann formulation as also provided in [@Tbl:heatfluxes1], with the internal heat flux calculated using the back-gorund heat load. The heat balance for the vessel wall is used to solve for the vessel wall temperature as a function of time
+
+$$ \frac{dT}{dt} = \frac{q_{external}-q_{internal}}{C_p \rho dx} $$
+
+Where:
+
+- $T$ is the vessel wall temperature (K)
+- $q_{external}$ is time dependent peak fire heat flux (W/m$^2$ K)
+- $q_{internal}$ is the internal convective heat flux (W/m$^2$ K)
+- $C_p$ is the material temeprature dependent heat capacity (J/kg K)
+- $\rho$ is the materal density (assumed constant) (kg/m$^3$)
+- $dx$ is the vessel wall material thickness (m)
+
+This oridinary differential equation is solved using a simple explicit Euler scheme. 
+
+The temperature dependent material properties has been sourced from the Scandpower guideline [@scandpower] and visualised in [@Fig:heat_capacity] and [@Fig:UTS] for heat capacity and Ultimate Tensile Strength, respectively. The materials implemented in *openthermo* are summarised in [@Tbl:materials].  
+
+![Steel heat capacities as a finction of temperature for the materials implemented in *openthermo*. The values have been sourced from [@scandpower] and missing values extrapolated to cover the same temperature range.](docs/img/heat_capacity.png){#fig:heat_capacity}
+
+
+![Steel Ultimate Tensile Strengt as a finction of temperature for the materials implemented in *openthermo*. The values have been sourced from [@scandpower] and missing values extrapolated to cover the same temperature range.](docs/img/UTS.png){#fig:UTS}
+
+
+| Steel type   | Type / alloy | ASME | DIN | ASTM |
+|--------------|--------------|------|-----|------|
+| Carbon steel | 235LT        |      |     | A-333 / A-671 |  
+|              | 360LT        |      |     |      |
+| Duplex (SS)  | 2205         | SA-770 | 1.4462 | A-790 |
+| Austenitic (SS) | 316       | A-358 316 | 1.4401 | A-320 |
+| Super austenitic (SS) | 6Mo |      | 1.4529 | B-677 |
+
+: Steel materials implemented in *openthermo*. {#tbl:materials}
+
+
 ## Model implementation
 A simple (naive) explicit Euler scheme is implemented to integrate the mass balance over time, with the mass rate being calculated from an orifice/valve equation as described in [@Sec:flow]:
 

src/hyddown/examples/rupture.yml

@@ -0,0 +1,27 @@
+vessel:
+  length: 9
+  diameter: 3
+  orientation: "vertical"
+  type: "Flat-end"
+  heat_capacity: 500
+  density: 7700
+  thickness: 0.136 
+initial:
+  temperature: 298.15
+  pressure: 11500000
+  fluid: "CH4" 
+calculation:
+  type: "energybalance"
+  time_step: 1
+  end_time: 900.
+valve:
+  flow: "discharge"
+  type: "relief"
+  set_pressure: 13500000
+  back_pressure: 101300.
+heat_transfer:
+  type: "s-b"
+  fire: "scandpower_pool"
+rupture:
+  material: "CS_360LT" 
+  fire: "scandpower_jet_peak_large"

src/hyddown/fire.py

@@ -69,6 +69,45 @@ def jet_fire_scandpower(Tvessel):
     return stefan_boltzmann(alpha, e_flame, e_surface, h, Tflame, Tradiative, Tvessel)
 
 
+def jet_fire_peak_large_scandpower(Tvessel):
+    """
+    Incident heat flux of 350 kW/m2
+    """
+    alpha = 0.85
+    e_flame = 1
+    e_surface = 0.85
+    h = 100
+    Tflame = 1429.61
+    Tradiative = 1429.61
+    return stefan_boltzmann(alpha, e_flame, e_surface, h, Tflame, Tradiative, Tvessel)
+
+
+def jet_fire_peak_small_scandpower(Tvessel):
+    """
+    Incident heat flux of 250 kW/m2
+    """
+    alpha = 0.85
+    e_flame = 1
+    e_surface = 0.85
+    h = 100
+    Tflame = 1279.29
+    Tradiative = 1279.29
+    return stefan_boltzmann(alpha, e_flame, e_surface, h, Tflame, Tradiative, Tvessel)
+
+
+def pool_fire_peak_scandpower(Tvessel):
+    """
+    Incident heat flux of 150 kW/m2
+    """
+    alpha = 0.85
+    e_flame = 1
+    e_surface = 0.85
+    h = 30
+    Tflame = 1212.54
+    Tradiative = 1212.54
+    return stefan_boltzmann(alpha, e_flame, e_surface, h, Tflame, Tradiative, Tvessel)
+
+
 def sb_fire(T_vessel, fire_type):
     if fire_type == "api_jet":
         Q = jet_fire_api521(T_vessel)
@@ -78,6 +117,12 @@ def sb_fire(T_vessel, fire_type):
         Q = pool_fire_scandpower(T_vessel)
     elif fire_type == "scandpower_jet":
         Q = jet_fire_scandpower(T_vessel)
+    elif fire_type == "scandpower_jet_peak_large":
+        Q = jet_fire_peak_large_scandpower(T_vessel)
+    elif fire_type == "scandpower_jet_peak_small":
+        Q = jet_fire_peak_small_scandpower(T_vessel)
+    elif fire_type == "scandpower_pool_peak":
+        Q = pool_fire_peak_scandpower(T_vessel)
     else:
         raise ValueError("Unknown Stefan-Bolzmann fire heat load")
     return Q

src/hyddown/hdclass.py

@@ -191,6 +191,13 @@ def read_input(self):
         # - functional mass flow
         self.thickness = 0
 
+        if "rupture" in self.input:
+            self.rupture_material = self.input["rupture"]["material"]
+            if "fire" in self.input["rupture"]:
+                self.rupture_fire = self.input["rupture"]["fire"]
+            else:
+                self.rupture_fire = "api_jet"
+
         # Reading heat transfer related data/information
         if "heat_transfer" in self.input:
             self.heat_method = self.input["heat_transfer"]["type"]
@@ -1042,38 +1049,22 @@ def run(self, disable_pbar=True):
 
                     self.Q_outer[i] = (
                         fire.sb_fire(self.T_vessel[i - 1], self.fire_type)
-                        # (
-                        #     self.sb_heat_load(self.time_array[i])
-                        #     - 0.85 * 5.67e-8 * self.T_vessel[i - 1] ** 4
-                        # )
                         * (self.surf_area_inner - wetted_area)
                         * self.surf_area_outer
                         / self.surf_area_inner
                     )
 
                     self.q_outer[i] = fire.sb_fire(self.T_vessel[i - 1], self.fire_type)
-                    # (
-                    #     self.sb_heat_load(self.time_array[i])
-                    #     - 0.85 * 5.67e-8 * self.T_vessel[i - 1] ** 4
-                    # )
 
                     self.Q_outer_wetted[i] = (
                         fire.sb_fire(self.T_vessel_wetted[i - 1], self.fire_type)
-                        # (
-                        #     self.sb_heat_load(self.time_array[i])
-                        #     - 0.85 * 5.67e-8 * self.T_vessel_wetted[i - 1] ** 4
-                        # )
                         * wetted_area
                         * self.surf_area_outer
                         / self.surf_area_inner
                     )
                     self.q_outer_wetted[i] = fire.sb_fire(
                         self.T_vessel_wetted[i - 1], self.fire_type
                     )
-                    # (
-                    #     self.sb_heat_load(self.time_array[i])
-                    #     - 0.85 * 5.67e-8 * self.T_vessel_wetted[i - 1] ** 4
-                    # )
 
                     if np.isnan(self.Q_outer_wetted[i]):
                         self.Q_outer_wetted[i] = 0
@@ -1701,6 +1692,112 @@ def plot_tprofile(self, filename=None, verbose=True):
         if verbose:
             plt.show()
 
+    def analyze_rupture(self, filename=None):
+        from hyddown.materials import steel_Cp, ATS, von_mises
+        from hyddown import fire
+
+        pres = lambda x: np.interp(x, self.time_array, self.P)
+        q_unwetted = lambda x: np.interp(x, self.time_array, self.q_inner)
+        q_wetted = lambda x: np.interp(x, self.time_array, self.q_inner_wetted)
+
+        T0_unwetted = self.T_vessel[0]
+        T0_wetted = self.T_vessel_wetted[0]
+
+        thk = self.thickness
+        rho = self.vessel_density
+        inner_diameter = self.diameter
+
+        dt = 10
+        max_time = self.time_array[-1]
+        tsteps = int(max_time / dt)
+
+        T_wetted_wall = np.zeros(tsteps + 1)
+        T_unwetted_wall = np.zeros(tsteps + 1)
+        T_wetted_wall[0] = T0_wetted
+        T_unwetted_wall[0] = T0_unwetted
+        peak_times = np.zeros(tsteps + 1)
+        peak_times[0] = 0
+
+        for i in range(tsteps):
+            peak_times[i + 1] = peak_times[i] + dt
+            q_fire_wetted = fire.sb_fire(T_wetted_wall[i], self.rupture_fire)
+            q_fire_unwetted = fire.sb_fire(T_unwetted_wall[i], self.rupture_fire)
+            T_wetted_wall[i + 1] = T_wetted_wall[i] + (
+                q_fire_wetted - q_wetted(peak_times[i])
+            ) * dt / (thk * rho * steel_Cp(T_wetted_wall[i], self.rupture_material))
+            T_unwetted_wall[i + 1] = T_unwetted_wall[i] + (
+                q_fire_unwetted - q_unwetted(peak_times[i])
+            ) * dt / (thk * rho * steel_Cp(T_unwetted_wall[i], self.rupture_material))
+
+        ATS_wetted = np.array([ATS(T, self.rupture_material) for T in T_wetted_wall])
+        ATS_unwetted = np.array(
+            [ATS(T, self.rupture_material) for T in T_unwetted_wall]
+        )
+        von_mises_wetted = von_mises_unwetted = np.array(
+            [von_mises(pres(time), inner_diameter, thk) for time in peak_times]
+        )
+
+        self.peak_times = peak_times
+        self.von_mises = von_mises_unwetted
+        self.ATS_unwetted = ATS_unwetted
+        self.ATS_wetted = ATS_wetted
+        self.peak_T_wetted = T_wetted_wall
+        self.peak_T_unwetted = T_unwetted_wall
+
+        if np.all(ATS_unwetted > von_mises_unwetted) == True:
+            self.rupture_time = None
+            # print("No rupture")
+        elif np.all(ATS_unwetted < von_mises_unwetted) == True:
+            self.rupture_time = 0
+            # print("Rupture at time=0")
+        else:
+            rupture_idx = np.where(ATS_unwetted < von_mises_unwetted)[0][0]
+            self.rupture_time = (
+                peak_times[rupture_idx - 1] + peak_times[rupture_idx]
+            ) / 2
+            # print("Rupture time +/- 5 s:", self.rupture_time)
+            # print("Rupture pressure (bar)", pres(peak_times[rupture_idx - 1]))
+
+        from matplotlib import pyplot as plt
+
+        plt.figure(1)
+        plt.plot(peak_times, von_mises_wetted / 1e6, label="von Mises stress")
+
+        plt.plot(peak_times, ATS_unwetted / 1e6, label="ATS unwetted wall")
+        if self.liquid_level.all() != 0:
+            plt.plot(peak_times, ATS_wetted / 1e6, label="ATS wetted wall")
+        plt.xlabel("Time (s)")
+        plt.ylabel("Allowable Tensile Strength / von Mises Stress (MPa)")
+        plt.legend(loc="best")
+        if filename is not None:
+            plt.savefig(
+                filename + "_ATS_vonmises.png",
+            )
+        plt.figure(2)
+        if self.liquid_level.all() != 0:
+            plt.plot(peak_times, T_wetted_wall - 273.15, label="T wetted wall")
+        plt.plot(peak_times, T_unwetted_wall - 273.15, label="T unwetted wall")
+        plt.xlabel("Time (s)")
+        plt.ylabel("Wall temperature (C)")
+        plt.legend(loc="best")
+        if filename is not None:
+            plt.savefig(
+                filename + "_peak_wall_temp.png",
+            )
+
+        plt.figure(3)
+        plt.plot(
+            peak_times,
+            np.array([pres(time) for time in peak_times]) / 1e5,
+            label="Pressure",
+        )
+        plt.xlabel("Time (s)")
+        plt.ylabel("Pressure (bar)")
+        plt.legend(loc="best")
+
+        if filename is None:
+            plt.show()
+
     def __str__(self):
         return "HydDown vessel filling/depressurization class"