Dynamic modeling of heat exchanger tube rupture

One fault that occurs with heat exchangers is a tube rupture, an overpressure scenario in which high pressure fluid flows into the low pressure region. It is a serious safety concern that may lead to significant damage. Accurate prediction of the pressure build-up after a rupture is critical to determine the appropriate size of a relief device and avoid exceeding allowable pressure limits. This paper describes a model-based step-by-step methodology to predict dynamic pressure profiles during tube rupture for liquid-liquid, vapor-liquid, and flashing liquid-liquid systems. The transient effects of the relief valve are considered. The effects of choked flow must also be considered for accurate maximum pressure predictions. Using a dimensionless analysis, the pressure ratio and density ratio are shown to significantly impact the severity of this incident. Results show that vapor-liquid systems result in the highest pressure surges.


Introduction
In process industries, pressure relief valves are commonly referred to as the last line of defense. Code requires that when a process fails, regardless of its built-in redundancies, there exists an independent protection system (e.g., pressure relief valve) powered by the medium it protects. [1] Therefore, correct pressure relief valve sizing is critical to ensure the protection of equipment, property, and life. Fortunately, for most overpressure scenarios, pressure relief standards set by API 521 help guide engineers on how to tackle various overpressure scenarios. For heat exchanger tube ruptures, an overpressure event that has been documented in industry, API 521 recommends a dynamic analysis be performed when there exist large differences in pressure. [2] However, a framework for dynamically modeling tube rupture scenarios is not given. A sketch of the tube rupture scenario is given in Fig. 1. This scenario involves a heat exchanger with high pressure on the tube-side and low pressure on the shell-side. Exchanger tubes over time can weaken (shown as a zigzag stress pattern in Fig. 1). If unnoticed or not replaced in time, this weak point may fail. High pressure tube-side *Correspondence: hasan@tamu.edu 1 Artie McFerrin Department of Chemical Engineering, Texas A&M University, TX 77843-3122 College Station, USA fluid enters the shell-side, increasing the system pressure. At the set point, a pressure relief valve opens, venting shell-side fluid to the flare. With heat exchangers in the offshore industry having tube-side pressures into the hundreds of bar, dynamic modeling is critical to selecting an appropriate relief device. An additional benefit of accurately sizing overpressure scenarios via dynamic analysis is the potential to reduce the weight of equipment, lowering capital costs. [3] This is especially important in the offshore industry, where topside weight reduction is an important design criterion.
Before modeling the tube rupture phenomena, it is helpful to understand its relevance to industry. Tube ruptures can occur in the chemical, oil, gas, offshore, and even the nuclear industry. Knowing where tube ruptures are more common allows one to scrutinize one set of systems over others. For reasons including the need to protect proprietary process information, reliability data on the failure of equipment in industry is hard to come by. However, there are efforts to counter this such as the Mary Kay O'Connor Process Safety Center's Instrument Reliability Network [4].
Because of this, it is often difficult to find statistics on the frequency of tube ruptures. Moreover, some failure rate databases may record a tube rupture failure under a Fig. 1 Sketch of heat exchanger tube rupture different name, making it even harder to estimate the likelihood of this event. Fortunately, there are some databases that give us a clue as to where tube ruptures occur. The United Kingdom's Hydrocarbon Release Database collected data on offshore incidents. [5] By filtering out exchanger related releases, Table 1 shows an estimate of which areas in a plant experienced a tube rupture. From the table, it is evident that most offshore tube ruptures occur in gas compression systems. This is expected since gas compression systems contain a series of intercoolers, with high pressure gases on the tube-side, and low pressure cooling water on the shell-side.
The model by Fowler et al. was one of the earliest attempts to describe the tube rupture transient process. [6] A shell-average pressure model was developed to capture the pressure-time relationship. This was followed by a stress analysis to determine if shell failure was likely to occur. Simpson developed a different shell-average pressure model based on the assumption that a tube failure will result in an expanding translating spherical explosion [7]. The models by Fowler et al. and Simpson both accounted for the inertia of the shell-side fluid and relief line. They did not however simulate the sudden pressure increase following a tube break (commonly referred to as water-hammer). Sumaria et al. modeled the transient pressure surge by using a lumped parameter model to incorporate the compressibility of the shell-side fluid and relief line. [8] This model was applied to a high pressure gas entering a low pressure liquid. Ennis et al. extended this work by considering the case of a tube-side Processing/gas/dehydration 3 5 8 Processing/gas/LPG/condensate 1 2 3 Processing/gas/sour gas 1 1 Processing/oil/oil treatment 1 1 Separation/gas production 1 1 Utilities/gas/fuel gas 4 4 Utilities/oil/heat transfer oil 1 1 Total 5 28 25 58 flashing liquid and accounted for the spatial and temporal aspects of the attached piping system. [9] Furthermore, Botros addressed the importance of accounting for piping systems and specifying appropriate boundary conditions [10]. Cassata et al used commercial transient analysis software to model a tube rupture by performing a series of simulations with the burst occurring at the U bend. [11] Pressure profiles at various points along the exchanger were generated. Ewan et al calculated pressure profiles when varying the location of the tube burst. [12] For relief devices close to the tube burst, pulse widths were higher than bursts far from the relief device. Nagpal evaluated several tube rupture cases using dynamic simulations and provides a chart for when dynamic simulations are recommended. [13] The simulation results show that pressure relief devices can be used to protect liquidfilled exchangers. Aside from simulating the severity of this event through consequence modeling, Acosta and Siu examined the risk of this scenario through the use of dynamic event trees. [14] Lastly, a Joint Industry Project conducted a full-scale heat exchanger tube rupture. [15] These experimental results were compared against blind simulations performed by consultants. The results support using one dimensional dynamic models to model the pressure profiles during tube ruptures.
One technique on estimating the relief rate during a tube rupture is to apply a steady state calculation. This is done by first assuming a complete tube rupture has occurred (referred to in the literature as a guillotine rupture). A standard hydraulic orifice calculation is then performed to determine what is the rate of fluid entering the shell-side. This number is then doubled due to the fact that fluid is entering from both ends of the ruptured tube. Finally, this rate is then inputted into a relief valve calculation in order to determine what size valve is required for adequate protection of the heat exchanger. There are some issues that exist with this method however. The first is the concern that the shell-side of the exchanger can fail before the relief valve can open. The reasoning behind this is that tube rupture scenarios result in fast-acting transient shockwaves which can reach peak pressure in milliseconds. [16] This is especially the case for high pressure differences between the shell and tube-side.
Another concern with this technique is that the orifice calculation method ignores the shell-side throughout the calculation. This is a problem because the tube rupture failure scenario is entirely focused on the shell-side failing since it is the low pressure side of the exchanger. Orificestyle calculations overlook how the shell-side pressure profile develops throughout the tube failure -the effects of mixing between the shell and tube-side fluids are ignored. Lastly, orifice calculations are more appropriate for tube rupture calculations that do not have a subcooled liquid on the tube-side. Subcooled liquids will flash once entering the low pressure shell-side, further complicating the calculation of the required relief rate. For heat exchanger tube ruptures, the goal is to prevent the shell-side of the exchanger from failing. Thus, a dynamic simulation gives a more complete picture of the overpressure scenario, increasing a plant's confidence in its ability to withstand a tube rupture.
The novel contributions of this paper include a step-bystep methodology on dynamically sizing a relief system for a tube rupture. The model will be useful in answering critical safety related questions such as what is the maximum pressure the shell-side reaches, how long after the tube rupture is maximum pressure reached, and how long is the shell design pressure exceeded. Examples on how to apply this methodology to a liquid-liquid, vapor-liquid, and flashing liquid-liquid system are given. This paper explains how to use process simulation software to generate the necessary thermophysical properties required for accurate pressure profile predictions. Finally, a dimensionless analysis of properties that affect tube ruptures is performed.
The article is organized as follows. "Background" section provides a background on this overpressure event, in addition to examining the common causes of heat exchanger tubes failing. "Background" section also provides a brief overview of ASME VIII pressure vessel code and its relation to the tube rupture scenario. This helps to specify an acceptable pressure limit during a tube rupture event. Lastly, "Background" section also describes the differences in relief devices used to protect exchangers and their associated response times. "Methods" section models the tube rupture. "Background" section first lists the assumptions in this model, followed by the modeling of each component that is fed into a tube rupture dynamic pressure profile algorithm. "Results & Discussion" section illustrates how the equations in "Methods" section are used in a practical setting, with examples on liquidliquid, vapor-liquid, and flashing liquid-liquid systems. A discussion on the differences between these systems is included. "Dimensionless Analysis of Tube Rupture" section describes a dimensionless analysis performed for a liquid-liquid tube rupture. "Mitigating Tube Rupture Scenarios" section includes strategies on mitigating a tube rupture. This is followed by this paper's conclusion, in "Conclusion" section.

Background
A tube in an exchanger can fail from fatigue, corrosion, increased temperature, and conflict with baffling. [17] A case study from SABIC Europe shows baffle hammering as having the potential to damage tubes. [18] Depending on the tolerance between the baffle and tubes, excessive movement by the tubes will result in continued striking against the baffle. Over time, this has the possibility to form a cavity. This creates a weak point in the tube, effectively reducing the design pressure of the tube. Similarly, Shahrani et al studied the failure mechanism of heat exchanger tubes in a plant. [19] It was concluded that fretting between the baffles and tubes, in combination with the corrosion due to the presence of sulfur, led to the tube failing. Vibration from equipment such as compressors can also cause tube failures via baffle hammering. One way to counter this is to reduce the velocity to minimize vibrations. For new exchangers, materials of construction and/or adding design features such as baffle spacing may also allow for increased vibration.
Metal erosion from high velocity particulates over time can also thin tubes and result in a rupture. [17] The erosion rate primarily depends on the velocity, diameter, pipe thickness, and mass of particulates. The impact angle also affects the severity of erosion. Because of the need for the fluid to change direction, in a piping system, pipe bends frequently exhibit higher rates of erosion. [20] Another case study by Khilnaney examines the failures of tube ruptures at a nuclear plant cooling water exchanger. [21] In this study, hydrogen sulfide was detected in the cooling tower, thus suspecting a tube leak had occurred. It was found that in 50% of these tubes, fluid elastic instability across the tube may have played a role in their failure. Lastly, corrosion from harsh chemicals such as sulfur and caustic chemicals can weaken tubes to the point of failure. ASME VIII code governs the performance of pressure vessels. It specifically applies to pressure vessels operating at or exceeding 15 psig. [22] Most exchangers in a plant have operating pressures that fall under this criterion and are consequently governed by ASME VIII code. The code also defines a maximum allowable working pressure (MAWP) based on an equipment's material, thickness, and other mechanical properties. During overpressure scenarios, pressure vessels must vent excess pressure to a flare, recycled to the process, or to the environment (depending on the fluid being vented). For the most part, pressure relief valves are set at or below the MAWP. There are exceptions to this when dealing with multiple relief valves. Note that the MAWP and the design pressure often refer to the same pressure and are used interchangeably. This paper makes no distinction between the two terms.
For a single pressure relief valve protecting a system, a 10% accumulation in pressure beyond the MAWP is allowed [22]. For multiple relief valves providing protection, the accumulation is allowed to increase to 16% above the MAWP [22]. The last common overpressure scenario is for dealing with external fire, in which case, accumulation can increase up to 21% beyond the MAWP [22]. For pressure vessels, the hydrotest pressure is tested anywhere between 1.3 to 1.5 times the MAWP [23]. This work assumes a hydrotest pressure of 1.5 times the design pressure. This work can easily be reformulated to apply to a hydrotest pressure of 1.3 times the design pressure, however. Lastly, because of the low probability of tube ruptures occurring in comparison with other more common relief events, and because a tube rupture occurs over a short period of time, the shell-side hydrotest pressure is set as the upper bound. Some plants may not even consider an overpressure scenario up to the hydrotest pressure to be a layer of protection analysis (LOPA) scenario. [24] Figure 2 shows an overview of where the different pressures are located in relation to one another. The range of pressures for a tube rupture scenario can theoretically vary from the shell-side operating pressure to the tubeside hydrotest pressure. A normal operating window for the shell-side is between its operating pressure and design pressure. Somewhere between these pressures, a pressure relief device is set to open. Similarly, the tube-side also has an operating window between its operating pressure and design pressure, as well as a relief device set between these pressures. The next threshold after the shell-side design pressure is its hydrotest pressure. When modeling the dynamic pressure profile, any pressure between the operating pressure and hydrotest pressure will be considered acceptable. Beyond this pressure, the severity of a tube rupture and the likelihood of the shell-side failing increase with pressure (shown as the region in yellow-to-red region in Fig. 2).
For most overpressure scenarios, relief valves are able to open in time for pressure increases. With heat exchanger tube ruptures however, there is some concern that tube ruptures (which can generate peak pressures in milliseconds) can overpressure the shell-side of an exchanger before a pressure relief valve can take effect. A 2002 study prepared by Pipeline Simulation and Integrity Ltd for the United Kingdom's Health and Safety Executive tested the respons times of pressure relief valves and rupture disks.  [16] The study found that fast-acting relief valves do exist. Note that if a vendor includes a datasheet that lists the response time, that should be used in the model. An ifstatement can be included that doesn't allow any fluid to vent through the relief valve from the initial rupture until the response time has been met.

Methods
This paper examines the three most common tube rupture scenarios: the liquid-liquid system, the vapor-liquid system, and the flashing liquid-liquid system. For all shell and tube exchanger configurations, tube ruptures can be divided into four distinct phases [25]: • Phase I: Sudden rupture • Phase II: Development of transient shell-side pressure wave • Phase III: Shell-side liquid discharges through pressure relief system • Phase IV: Shell-side vapor discharges through pressure relief system The first phase generates a shock wave from the rupture of the tube. This is, however, considered to be nearly instantaneous and is not modeled. With an open path, the high pressure side fluid enters the low pressure side. Because the tube-side fluid acts against the shell-side fluid, the shell-side fluid is pushed into the relief system. Once this is complete, the tube-side fluid (either vapor or liquid) has a path to exit the exchanger. Phases II, III, and IV are modeled with the following assumptions: • Only one tube is assumed to rupture.
• The rupture is assumed to be a full-bore rupture.
• A tube rupture exposes both ends of the severed tube, allowing fluid to enter the shell-side at twice the cross-sectional area of one tube. • The tube-side fluid enters the shell-side as an isentropic nozzle flow. • The effects of temperature are ignored (e.g., increased heating resulting in the buildup of pressure). • An infinite reservoir of tube-side fluid exists. This is a conservative assumption. Without this assumption, one would need to determine what is the quantity of tube-side fluid available (from upstream and downstream units). One can also determine what is the rate of tube-side pressure loss (and decrease in velocity) as the fluid is being lost to the shell-side. • Pressure relief systems have an instantaneous response time (i.e., zero milliseconds). • No credit is taken for outflow via inlet or outlet piping. The phrase "taking credit for outflow" allows one to reduce the relieving requirements for overpressure scenarios by accounting for fluid leaving the system. This however comes with a set of requirements to ensure that an open path exists (and is adequately sized) for fluid to escape. For this model, all excess pressure must be vented through available relief systems.
• All areas of the exchanger's shell-side are in equilibrium. The effects of localized pressure are ignored. • The exchanger shell-side can safely handle pressures up to the hydrotest pressure.
What follows describes how one can model and predict the liquid and vapor flow rates, flashing liquid flow rates, vapor densities, mass fluxes from the tube to shell-side and from the shell-side to outside through a pressure relief valve, and the dynamic pressure profiles as a consequence of tube rupture in a shell and tube heat exchanger.

Liquid Flow Rate
The volumetric energy balance for isentropic nozzle flow is given by Eq. 1. [26] Where G is the tube-side mass flux entering the shell (kg/s-m 2 ), v is the specific volume of the tube-side fluid (m 3 /kg), P 1 and P 2 are the pressure (Pa) the tube-side fluid experiences at the inlet and outlet, respectively, and t is the fluid condition at the throat of the nozzle where the crosssectional area is minimized. Assuming a constant pressure step size, the integration portion in Eq. 1 can be expressed via Eq. 2.
Where v is the specific volume of the tube-side fluid (m 3 /kg), h is the constant pressure step size chosen for summation purposes (Pa), n is the index for fluid conditions at the assumed throat pressure (i.e., the assumed endpoint pressure for the integral), and j is the increment counter used for summation purposes. Combining Equations 1 and 2, we are able to obtain the mass flux equation for liquid flow through a tube (Eq. 3).
Once G is determined for different pressure intervals, the mass flux should be fitted against pressure by using a quadratic fit in the following format (Eq. 4), where A1, A2, and A3 are the coefficients that give the best fit (e.g., method of least squares). Smaller pressure intervals will increase the accuracy of G. The end result is a continuous function that yields a mass flux for all pressures between the shell and tube-side pressure.

Vapor Flow Rate
The methodology for the calculation of vapor flow through the ruptured tube is very similar to that of the liquid flow rate. A series of isentropic flashes from high pressure to low pressure conditions are performed.
Using the relationship between pressure and specific volume, Equations 1 -3 are solved to yield a mass flux as a function of pressure. Up to this point, the solution method for the liquid flow rate and the vapor flow rate are identical. However, for most vapor-liquid tube ruptures, a choked flow condition will be present. The mass flux function, G, when calculated, will show whether or not a choked condition exists. If that is the case, the mass flux versus pressure must be adjusted to reflect the choked condition. This is fitted using a cubic fit in the following format (Eq. 5), where α1, α2, α3, and α4 are the coefficients that give the best fit (e.g., method of least squares). An example of how this is performed is presented later.

Flashing Liquid Flow Rate
The flow rate of a flashing liquid entering the shell-side incorporates techniques used in calculating the liquid flow rate and vapor flow rate. First, a series of isentropic flashes from the tube-side pressure to the shell-side pressure are performed. Equations 1 -3 are solved to yield a mass flux as a function of pressure. If a choked condition does result, then the mass flux should be fitted as a cubic equation (Eq. 5). These steps are identical to the ones performed in the "Vapor Flow Rate" section. The main difference, however, is that the downstream phase in a flashing liquid may take the form of a vapor, twophase, or liquid, depending on its pressure. Thus, the vapor phase fraction at each pressure interval needs to be calculated in addition to the total mass flux. The vapor phase fraction, y t , is expressed as a piecewise function as follows: Where C1 and C2 are the coefficients that give the best fit, and P tbp is the tube-side fluid's bubble point pressure. The vapor phase fraction can then be multiplied with the total mass flux, G, to yield the tube-side vapor mass flux, G v , and one minus the vapor phase fraction multiplied by G yields the tube-side liquid mass flux, G l :

Density of Tube-Side Vapor, ρ tv
With liquid-liquid systems, the density is assumed to be constant. However, with vapor-liquid systems, density must be changed as a function of pressure. If the vapor density is treated as an ideal gas, it would be inaccurate at high pressures and low temperatures. If density versus pressure is fitted using experimental data, the results will be accurate. The problem with this, however, is that the number of gases available to deploy this method are limited. A general way is to use an equation of state to determine the density versus pressure relationship. This can be accomplished by using widely available process simulation software (e.g., Aspen HYSYS). Density and pressure are expressed via a linear relationship as follows: Where B1 and B2 are the coefficients that give the best fit.

Rate of Tube-Side Fluid Entering Shell-Side
With the tube-side mass flux, G, determined, the area of a single tube and the tube-side mass flow rate can be obtained by Equations 10, 11, and 12. Note that the mass flow rate is listed as two times the flux multiplied by the area. This is because the flow entering the shell-side will be entering from upstream of one ruptured tube and downstream of the other ruptured tube.

Rate of Shell-Side Fluid Exiting Relief Device
Assuming a pressure relief valve can be modeled as flow through an orifice, the following gives the mass flow rate exiting the relief device: Note that the equations used for flow through a pressure relief valve can be modified to model spring loaded, balanced bellows, and pilot valves. An example of this performed is by Singh who developed a one-dimensional dynamic model for a spring loaded pressure safety valve (PSV). [27] Moreover, any equations of flow through a PSV provided by a vendor can also be substituted here. The expression can also be modified to reflect the reseat pressure of a relief valve.

Bulk Modulus of Elasticity
Interpolation should be used to determine the bulk modulus for a fluid, B, shown in Eq. 14. While using the same temperature as that of the system, the high and low side pressures (P 1 and P 0 ) and their respective volumes (V 1 and V 0 ) are chosen to calculate the bulk modulus. Although the bulk modulus does vary, it is reasonable to assume it is constant throughout the entire range of pressures during the overpressure event.

Tube-Side Vapor and Liquid Volumes Present in Shell, V tv and V tl
V tv and V tl represent the tube-side vapor and tube-side liquid volumes that are present in the exchanger shell. At the beginning of a tube rupture, no tube-side fluid is present in the shell. As fluid starts entering the shell-side, V tv and V tl , are calculated by the following Equations: Whereṁ tv and ρ tv are the tube-side vapor mass flow rate (kg/s) and density (kg/m 3 ), respectively. Andṁ tl and ρ tl are the tube-side liquid mass flow rate (kg/s) and density (kg/m 3 ), respectively.

Shell Pressurization
A conservation of volume balance on the exchanger shell yields Eq. 17. [9] This allows one to determine the change in pressure over time during a tube rupture. For liquidliquid systems, the vapor associated termsṁ tv , ρ tv , V tv , and c tv0 are removed from the equation. Similarly, for vapor-liquid systems, the termsṁ tl , ρ tl , V tl , and B tl are removed from the equation. In the case of a flashing tubeside liquid with liquid on the shell-side, all terms in Eq. 17 would be present.
Whereṁ tv and ρ tv are the tube-side vapor mass flow rate (kg/s) and density (kg/m 3 ), respectively.ṁ tl and ρ tl are the tube-side liquid mass flow rate (kg/s) and density (kg/m 3 ), respectively.ṁ psv and ρ sl are the PSV mass flow rate (kg/s) and shell-side liquid density (kg/m 3 ), respectively. represents the tube-side volume and bulk modulus terms. V sl and B sl are the volume of shell-side liquid (m 3 ) remaining in the shell-side and the shell-side liquid's bulk modulus (Pa), respectively. V shell and B shell are the shell volume (m 3 ) and the bulk modulus (Pa) of the shell material of construction (e.g., carbon steel), respectively. The three cases for are given by Eq. 18. V tv and c tv0 are the volume of tube-side vapor (m 3 ) that has entered the shell-side and the tube-side vapor's speed of sound (m/s), respectively. V tl and B tl are the volume of tube-side liquid (m 3 ) that has entered the shell-side and the tube-side liquid's bulk modulus (Pa), respectively. This differential equation is the main mathematical model that will be solved and is valid for liquid-liquid, vaporliquid, and flashing liquid-liquid systems. The previous variables, described in earlier subsections all feed into this model.
Since the initial shell-side pressure in Eq. 17 is known, this differential model can be viewed as an initial value problem. Furthermore, this equation is a first order differential equation. A simple explicit method can be used to solve the model. Euler's method is used for this paper. Using Matlab, a custom algorithm based on Euler's method is developed [28]. The accuracy of this method increases with a smaller step size. Because tube ruptures can exhibit peak pressures in milliseconds, a step size of no larger than 0.1 millisecond is appropriate. Figure 3 lists the major steps in generating a tube rupture pressure profile. The first step in modeling a tube rupture is collecting the appropriate data about the system. These include, but are not limited to, equipment diagrams, process simulation parameters (pressure, temperature, density, etc.), and relief valve datasheets. The second step is to obtain the necessary thermodynamic properties. One way to achieve this is by creating a table of isentropic flashes from the tube-side conditions to the shell-side conditions. Table 2 lists the isentropic flashes for the liquid-liquid system solved in "Liquid-Liquid System" section. The integral term in Table 2 refers to the solution of Eq. 2. These properties are used to determine the tube-side mass flux which can be used to calculate the mass flow rate entering the shell-side. The bulk modulus is also calculated at this stage. Once the above information are obtained, the equations relating the mass flux and pressurization of the shell are used to calculate the shell-side pressure. An appropriate time-step is selected to solve Eq. 17 numerically. If the time-step is too large, then the resolution will not be detailed enough. If this occurs, the dynamic pressure profile obtained may not display the true peak of the tube rupture phenomena. Too small a time-step, on the other hand, increases the number of iterations needed to achieve a solution.

Liquid-Liquid System
Consider the following example. A carbon steel shell and tube exchanger has ethylene glycol in the tube-side and cooling water in the shell-side. The exchanger contains 150 tubes each 2 meters long with an inner diameter of 15 mm. The shell's volume is 7.5 m 3 . The high and low side operating pressures are 10 bar and 1 bar, respectively. The high and low side operating temperatures are 100°C and 20°C, respectively. The shell-side has a PSV set pressure of 1.2 bar (equal to its design pressure) and a hydrotest pressure of 1.8 bar. The following steps are used to determine the required PSV size in order to adequately protect this system.
Step 1: Obtain the physical properties of the liquidliquid system.
Using the Peng-Robinson equation of state, at 1 bar and 20°C, the density of water is 1011 kg/m 3 . Similarly, at 10 bar and 100°C, the density of ethylene glycol is 1055 kg/m 3 . We then determine the bulk modulus for each fluid. It is determined that water has a bulk modulus of 3.44931 ×10 9 Pa, and ethylene glycol has a bulk modulus of 0.89769 ×10 9 Pa. The bulk modulus of common metals can be found in engineering reference manuals [29]. Carbon steel has a bulk modulus (B shell ) of 159 ×10 9 Pa.
Step 2: Use process simulation software to calculate thermodynamic properties.
Using Aspen HYSYS [30], a property table is created for the tube-side fluid. Isentropic flash calculations (in increments of 1 bar) are performed for ethylene glycol from the tube-side pressure (10 bar) to the shell-side pressure (1 bar). The results are listed in Table 2.
Step 3: Calculate the mass flow rate through the tube, entering the shell-side.
Using the mass flux equation and property table, the mass flux can be calculated (shown in Table 2). Once the mass flux for each pressure is determined, the data is then fitted to a quadratic equation in order to obtain Fig. 3 Steps needed to generate dynamic pressure profile during tube rupture the mass flux as a function of pressure. For the ethyleneglycol water system described, Eq. 4 is generated with coefficients A1, A2, A3 equal to -434.4, 526.4, and 41854.5, respectively. Note that the pressure in Eq. 4 refers to the shell-side pressure. At a pressure of 1 bar and 10 bar, this equation yields the maximum and minimum tubeside flow rates. This is reasonable since the largest driving force exists when the shell-side pressure is 1 bar. Similarly, when the shell-side is at 10 bar, no driving force exists.
Step 4: Using an appropriate time-step, determine the pressure at time t.
At time is equal to zero, the shell-side pressure is set to 1 bar. Using Matlab [28], Equations 4, 9-13, and 15-16 are then simultaneously solved using a time-step of 1 millisecond. The required data are shown in Table 3. Note the calculations presented shall assume a zero second response time for relief valves.
The pressure profiles for the liquid-liquid system are shown in Fig. 4. The pressure (bar) is plotted against time (milliseconds) for various standard PSV orifice sizes. This figure uses a time axis of 500 milliseconds, enough to capture the pseudo steady state pressures for all orifice sizes. A dotted line representing the hydrotest pressure (1.8 bar) is shown. As stated earlier, we do not wish to exceed the hydrotest pressure for this scenario. At time is equal to 0 milliseconds, the tube rupture begins. This period in time has the largest difference in pressure between the two exchanger sides.For orifice sizes D, E, F, and G, the hydrotest pressure is exceeded at 25 milliseconds. An 10 additional miliseconds shows orifice size H exceeding the hydrotest pressure. These orifice sizes reach a pseudo steady state condition at approximately 300 miliseconds, with shell-side pressures ranging from 3 bar to above 9 bar. Thus, for orifice sizes D, E, F, G, and H, the shell-side failing is almost certain due to the hydrotest pressure being significantly exceeded. Orifice sizes J and K are the only PSV sizes that adequately protect against this overpressure event. While orifice size J does exceed the design pressure, it is comfortably below the hydrotest pressure. Orifice size K is even more conservative with shell-side pressures not exceeding the design pressure during this tube rupture. Due to the lower cost of the J orifice PSV (resulting from its smaller size), this PSV would most likely be preferred in a plant.

Vapor-Liquid System
Consider the following example. A carbon steel shell and tube exchanger has methane in the tube-side and cooling water in the shell-side. The exchanger contains 100 tubes each 3 meters long with an inner diameter of 10 mm. The shell's volume is 7.5 m 3 . The high and low side operating pressures are 5 bar and 1 bar, respectively. The high and low side operating temperatures are 100°C and 20°C, respectively. The shell-side has a PSV set pressure of 1.2 bar (equal to its design pressure) and a hydrotest pressure of 1.8 bar. We would like to determine the required PSV size in order to adequately protect this system following the steps: Step 1: Obtain the physical properties of the vaporliquid system.
Using the Peng-Robinson equation of state, at 1 bar and 20°C, the density of water is 1011 kg/m 3 . Between 5 bar and 1 bar, the density of methane in kg/m 3 is given by Eq. 9 with coefficients B1 and B2 equal to 0.4747 and 0.58, respectively. As previously stated, water has a bulk modulus of 3.44931 ×10 9 Pa. Similarly, carbon steel has a bulk modulus (B shell ) of 159 ×10 9 Pa.
Step 2: Use process simulation software to calculate thermodynamic properties.
Using Aspen HYSYS [30], a property table is created for the tube-side fluid. Isentropic flash calculations (in increments of 0.4 bar) are performed for methane from the  The results are listed in Table 4.
Step 3: Calculate the mass flow rate through the tube, entering the shell-side.
Using the mass flux equation and property table, the mass flux can be calculated and is given in Table 4. Note that the methane will exhibit choked flow, and the mass flux must reflect this (see Fig. 5). Thus, the mass flux versus pressure data is fitted by a cubic polynomial (Eq. 5) with the coefficients α1, α2, α3, and α4 equal to -34.219, 219.62, -439.53, and 997.29, respectively. Between 1 bar and 5 bar, this equation yields the maximum and minimum tube-side flow rates.
Step 4: Using an appropriate time-step, determine the pressure at time t.
At time is equal to zero, the shell-side pressure is set to 1 bar. Using Matlab [28], Equations 5,[9][10][11][12][13][14], and 16 are then simultaneously solved using a time-step of 1 millisecond. The required data are shown in Table 3. Note the calculations presented shall assume a zero second response time for relief valves. Figure 5 shows the mass flux profile for methane at various downstream pressures. Two mass flux profiles for methane are included, one which takes into account choked flow (listed as red circles) and one that excludes the effects of choked flow (listed as blue triangles). The choked flow region is roughly between downstream pressures of 1 bar and 2.6 bar. In other words, from an initial pressure and temperature of 5 bar and 100°C, the mass flow rate of methane transported through an isentropic nozzle (as is the case in a tube rupture) increases with decreasing downstream pressure until 2.6 bar. Between downstream pressures of 1 bar and 2.6 bar, the mass flow rate of methane is constant and at its maximum (731 kg/s-m 2 ). Neglecting the effects of choked flow can result in a misleading mass flux, as can be seen by the negative parabolic shaped trendline in Fig. 5. Accepting this erroneous mass flux may lead to an undersized relief system. The pressure profiles for the methane water system are shown in Fig. 6. By plotting the pressure (bar) versus time (milliseconds), the effects of selecting different PSV orifice sizes (D, J, N, P, Q, R, and T) to protect this system are shown. A dotted line representing the hydrotest pressure (1.8 bar) is included. For this tube rupture, orifice size Q appears to provide adequate overpressure protection, without being conservative. One notable difference between this system and the liquid-liquid system is the speed at which pseudo steady state pressure is reached.   Fig. 4. This highlights the severity of a tube rupture in vapor-liquid systems.

Comparison of Liquid-Liquid and Vapor-Liquid Systems
Assuming the same tube and shell-side pressures, for an equivalent mass entering the shell-side, a vapor-liquid system reaches peak pressure significantly faster (and as a result is more severe) than that of a liquid-liquid system. This is true even though the fluid flow (on a mass basis) that enters the shell-side for vapors is significantly less than that of liquids. The vapor-liquid simulations indicate that the peak pressures are reached in less than 50 ms as opposed to the liquid-liquid cases which took between 100 and 300 ms, depending on the size of the PSV orifice. Shell-side transient peak pressures are able to significantly exceed their hydrotest pressures. For systems left with inadequate overpressure protection, these simulations demonstrate the consequences of tube ruptures.
The T and R orifices, as shown in Fig. 6, resulted in the pressure oscillating above and below the set pressure.
These pressure oscillations can be broken into a series of steps. First, the system pressure builds up due to an overpressure event, causing the PSV to open. Because the PSV is sufficiently large to handle the influx of fluid entering, the PSV vents all of the excess fluid along with additional fluid. At this point, the system pressure decreases, eventually dropping below the set pressure. Since the system pressure is below the set pressure, the PSV reseats to its original position. This cycle repeats for the duration of the overpressure event. More sophisticated modeling can include the effects of the PSV reseating. The rapid opening and closing of the relief valve is known as chattering and has been documented as affecting the reliability of a PSV. [31] Chattering can occur due to oversizing a relief valve, which appears to be the case for orifices T and R. For more common overpressure scenarios such as blocked outlet, control valve failures, and fires, plants may desire to select a relief valve that avoids chattering. However, because of the low frequency of tube ruptures in plants, chattering is not considered to be a primary concern.
Another distinction between vapor-liquid and liquidliquid systems is the presence of choked flow. For vaporliquid systems, ignoring choked flow can result in a relief system being undersized. Figure 6 shows the tube-side mass flux when including and excluding choked flow. For downstream pressures that are close to the upstream pressures, it is possible the fluid may not enter the choked region. This would yield the same mass flow rate whether choked flow is included or excluded. This is seen in the pressure interval between 5 bar and 2.6 bar in both Table 4 and Fig. 5. This should not be left to chance however. Testing for choked flow should always be performed.

Flashing Liquid-Liquid System
Consider the following example. A carbon steel shell and tube exchanger has propane in the tube-side and cooling water in the shell-side. The exchanger contains 100 tubes each 4 meters long with an inner diameter of 10 mm. The shell's volume is 7.5 m 3 . The high and low side operating pressures are 30 bar and 6 bar, respectively. The high and low side operating temperatures are 60°C and 6°C, respectively. The shell-side has a PSV set pressure of 7.2 bar (equal to its design pressure) and a hydrotest pressure of 10.8 bar. We would like to determine the required PSV size in order to adequately protect this system following the steps: Step 1: Obtain the physical properties of the flashing liquid-liquid system.
Using the Peng-Robinson equation of state, the density of 60°C liquid propane between 30 bar and 6 bar is on average 446 kg/m 3 . The density of vapor propane between 21 bar and 6 bar is given by Eq. 9 with coefficients B1 and B2 equal to 2.32 and -1.5468, respectively. As previously stated, water has a density of 1011 kg/m 3 and a bulk modulus of 3.44931 ×10 9 Pa. Liquid propane has a bulk modulus of 0.1536 ×10 9 Pa. Carbon steel has a bulk modulus (B shell ) of 159 ×10 9 Pa.
Step 2: Use process simulation software to calculate thermodynamic properties.
Using Aspen HYSYS [30], a property table is created for the tube-side fluid. Isentropic flash calculations (in increments of 1.5 bar) are performed for propane from the tube-side pressure (30 bar) to the shell-side pressure (6 bar). Because vapor is generated, the vapor fraction is included in the property table. The results are listed in Table 5.
Step 3: Calculate the mass flow rate through the tube, entering the shell-side.
Using the mass flux equation and property table, the mass flux can be calculated and is given in Table 5. Methane exhibits choked flow with a maximum mass flux of 27170 kg/s-m 2 . The corrected mass flux versus pressure is fitted against a cubic polynomial (Eq. 5) with the coefficients α1, α2, α3, and α4 equal to -8.131, 323.33, -3295.7, and 27649, respectively. The vapor fraction is expressed as Eq. 6 with coefficients C1 and C2 equal to -0.025 and Step 4: Using an appropriate time-step, determine the pressure at time t.
At time is equal to zero, the shell-side pressure is set to 6 bar. Using Matlab [28], Equations 6-16, and 16 are then simultaneously solved using a time-step of 1 millisecond. The required data are shown in Table 3. Note the calculations presented shall assume a zero second response time for relief valves. Figure 7 describes the vapor and liquid mass fluxes for propane with changes in pressure. The x-axis is the downstream pressure in bar. This downstream pressure is equivalent to the shell-side pressure. Thus, Fig. 7 can be interpreted as beginning from the left of the graph (at a downstream pressure of 6 bar), and shifting to the right as the tube rupture develops. If no PSV is present, this would continue until the downstream pressure is in equilibrium with the upstream pressure (30 bar). In addition to the x-axis, two y-axis are shown, one being the mass flux (kg/s-m 2 ) and the other being the vapor fraction. Between downstream pressures of 21 bar and 30 bar, the vapor fraction is zero. Between pressures 6 bar and 21 bar, the vapor fraction increases with decreasing downstream pressure. The vapor mass flux is generated by multiplying the total mass flux with the vapor fraction. Therefore, no vapor mass flux is present between the downstream pressures 21 bar and 30 bar.
The pressure profiles for the propane water system are shown in Fig. 8. Compared to the previous two cases, this system requires roughly between 100 to 200 milliseconds before reaching pseudo steady state. A dotted line line at 10.8 bar represents the hydrotest pressure. The smallest PSV orifice size that does not exceed this pressure is K. One interesting feature of this system is that it functions as both a liquid-liquid and vapor-liquid tube rupture. In other words, as the tube rupture begins, the high pressure liquid flashes when exposed to the shell-side. This phase of the tube rupture is closely associated with a vapor-liquid rupture. Over time, pressure builds up and the amount of flashing decreases. Eventually, the liquid exiting the tube-side will fully remain a liquid in the shell-side, exhibiting characteristics similar to that of a liquid-liquid rupture. Orifice sizes D and E are both piece-wise functions since around 20 bar, the shell-side is entirely a liquid due to the high pressure.

Dimensionless Analysis of Tube Rupture
The severity of a tube rupture can be affected by a myriad of properties. These include the fluid's density, phase, pressure, and temperature. Exchanger properties such as the shell volume and tube diameter affect the flux and the rate of pressure increase. Lastly, pressure relief properties (e.g. valve size and opening times) affect the rate of fluid exiting the shell. To better understand the different parameters that can affect a tube rupture, properties were grouped together to form dimensionless numbers. The following numbers are proposed: P ratio = P tube P shell (20) P max = P(t) shell max P shell design (21) ρ ratio = ρ tube ρ shell (22) τ in Eq. 19 evaluates the size of the tube in relation to the exchanger shell. A larger tube diameter increases the severity of a tube rupture, while a larger shell volume slows the accumulation of pressure. For Eq. 19, characteristic lengths of 1e-4 m and 1 m 3 for D c and V c were selected. The pressure ratio (Eq. 20) is investigated because it is most likely the primary driver in the severity of a tube rupture. In fact, API 521 recommends a dynamic simulation when pressure differences exceed 1000 psi. The dimensionless number P max in Eq. 21 is the maximum pressure the shell-side experiences during a tube rupture divided by the shell design pressure. Lastly, the effects of the density ratio's (Eq. 22) impact on tube rupture will be studied. A series of liquid-liquid tube rupture simulations were performed. In addition to the dimensionless numbers proposed, these simulations also measure the time it takes to reach the maximum pressure. The simulations contained ethylene-glycol and water on the tube and shellside, respectively. In Fig. 9, for a constant pressure ratio of 10, the maximum pressure was plotted against the density ratio (which was adjusted from 0.1 to 10). The dimensionless number, τ , was set at 13.3, 20, and 40. These values for τ are obtained from a 7.5 m 3 shell with 1 cm, 1.5 cm, and 3 cm diameter tubes, respectively. Shell-side peak pressure vs. density ratio for varying τ with a pressure ratio of 10 From Fig. 9, it is seen that the density ratio can significantly affect the ability of an exchanger to pressurize. A high density ratio reduces the ability of a tube rupture to result in overpressure. There exists for each system a density ratio threshold, in which no increase in pressure is present for this overpressure scenario. This relationship can be explained by thinking of an infinitely large density ratio. This would mean that the tube-side fluid is infinitely small and cannot pressurize the system. In mathematical terms, an infinite tube-side density forces the dP dt term in Eq. 17 to zero. Alternatively, low density ratios were able to reach 10 bar regardless of the tube diameter and shell volume. When looking at the effects of τ , for the maximum pressure a system experiences, a low value for τ results in a lower peak pressure. This is due to the fact that a low τ value results in a small tube diameter in relation to the shell volume. Thus, the cross sectional area of flow will be less, allowing a pressure relief device to mitigate the tube rupture's impact.
A second set of simulations were performed with the aim of studying the time it takes to reach the system peak pressure, P max . This parameter is important to accurately capturing the severity of a tube rupture. A scenario that reaches peak pressure quickly may result in the shell-side failing due to the PSV not opening in time. This is in contrast to a longer time to reach P max , which allows the PSV to open, limiting our concern to the maximum pressure the system reaches. The same pressure ratio (held constant at 10), density ratios, and values for τ as Fig. 9 were used. The results are shown in Fig. 10. All else being equal, it is preferred to have a tube rupture with a longer time to reach maximum pressure. For an equivalent density ratio, a faster time to reach P max is observed for higher τ values. For a fixed τ , a higher density ratio increases the time it takes for a system to reach its maximum pressure. In addition, for all values of τ , the maximum pressure is achieved in less than 350 milliseconds.
The third set of dimensionless analysis simulations (Fig. 11) study the impact the pressure ratio has on the severity of a tube rupture. The pressure ratio was varied, in increments of 2, from 2 bar to 10 bar. Two values of τ were used, 20 and 40. Similarly, three values for the density ratio were used, 0.49, 1.09, and 5.04. From Fig. 11, a higher pressure ratio leads to a higher maximum pressure. This is not surprising. What is interesting to note, however, is how strong the density ratio can affect this peak pressure. As an example, a pressure ratio of 10 with a τ of 20, can yield a maximum pressure of 1.2 bar, 3 bar, or 7 bar. This wide variation in maximum pressure is due to the systems having density ratios of 5.04, 1.09, and 0.49, respectively. An increase in τ also shows an increase in a tube rupture's maximum pressure, although not as severe as the density ratio. A pressure ratio of 8, with a density ratio of 0.49 yields a maximum pressure of 5 bar (τ equal to 20) and 8 bar (τ equal to 40).

Mitigating Tube Rupture Scenarios
Depending on how far along in the lifecycle of an exchanger a plant is, mitigating the severity of tube rupture incidents can be done through one of several ways. First, one should flag exchangers that will have large differences in pressure. These exchangers are most susceptible to failure in the event of a tube rupture. For exchangers that are still in the design phase, increasing the strength of the exchanger is one way of mitigating the effects of a tube rupture. This design change does increase the investment cost of the exchanger, however. If the exchanger is upgraded to comply with the two-thirds rule, then a pressure relief valve would not be required for that system. Another strategy to mitigate tube ruptures is preventing exchanger tubes from failing in the first place. Nowadays, nondestructive technology such as eddy current testing can help estimate the remaining life of exchanger tubes. [32] This technology can be used to more frequently monitor exchangers that are more susceptible to erosion and corrosion. Proactive prevention can also focus on detecting future operational hazards, such as by incorporating alarms within a process model. [33] One can also select an appropriately sized relief valve to handle a tube rupture scenario. For the ethylene-glycol water example in "Liquid-Liquid System" section, without a pressure relief valve, the shell-side would be pressurized from 1 bar to 9 bar in 300 ms. However, by installing either relief valve sizes J or K, the shell-side would exhibit peak pressures of 1.43 and 1.27 bar, respectively. Moreover, while relief size K exhibits chattering during this tube rupture, relief size J does not. Furthermore, the results from this paper indicate vapor-liquid and flashing liquidliquid ruptures to be more severe than liquid-liquid tube ruptures. Thus, more resources (e.g., increased design Fig. 11 Shell-side peak pressure vs. pressure ratio relationship for varying τ and density ratios pressure and larger PSV orifice sizes) can be dedicated to these higher risk exchangers. Lastly, it was found that the density ratio and pressure ratio can significantly affect the ability of a shell to pressurize. As a result, systems with a high density ratio and low pressure ratio are preferred for this overpressure scenario.

Conclusion
Shell and tube heat exchangers are commonly used in the oil, gas, chemical, and nuclear industries. Although rare, tube rupture overpressure events may compromise the mechanical integrity of an exchanger and can lead to the equipment's failure. This has the potential to result in catastrophic failures and should be modeled with rigorous sizing methods. This paper points out the challenges in modeling tube ruptures. The importance of accurately modeling tube rupture scenarios increases with large differences in pressure between the shell and tube. In addition, low tube-to-shell density ratios increase the severity of this event. By switching from an orifice-style calculation into a more rigorous dynamic simulation, one can better predict the effects of a tube rupture. This paper covered the main steps in generating a pressure profile, and a step-by-step calculation was performed for an ethylene-glycol water, methane water, and propane water system. A comparison of the PSV sizes reveals that vapor-liquid and liquid-liquid systems are the most and least severe cases, respectively. Flashing liquid-liquid tube rupture scenarios performed in between these two cases. The examples covered in this paper can serve as a basis for approaching liquid-liquid, vapor-liquid, and flashing liquid-liquid systems.

Abbreviations
A psv : Relief valve orifice area; B shell : Shell material of construction bulk modulus; B sl : Shell-side liquid bulk modulus; B tl : Tube-side liquid bulk modulus; C 0 : Coefficient of discharge; D c : Characteristic tube diameter; D tube : Tube diameter; P(t): Transient shell pressure; P max,norm : Normalized maximum transient shell-side pressure; P ratio : Dimensionless number relating tube and shell pressures; r tube : Tube radius; V c : Characteristic shell volume; V shell : Shell volume; V tl : Volume of tube-side liquid present in shell-side; V tv : Volume of tube-side vapor present in shell-side;ṁ psv : Relief valve mass flow rate;ṁ tl : Mass flow rate of tube-side liquid entering shell-side;ṁ tv : Mass flow rate of tube-side vapor entering shell-side; ρ ratio : Dimensionless number relating tube and shell densities; ρ sl : Shell-side liquid density; ρ tl : Tube-side liquid density; ρ tv : Tube-side vapor density; τ : Dimensionless number relating tube size to shell volume; MAWP: Maximum allowable working pressure; PSV: Pressure safety valve or pressure relief valve; G: Mass flux of tube-side fluid entering shell-side; h: Pressure step size; set: Shell-side set pressure