Dynamic Modeling of the Two-Phase Leakage Process of Natural Gas Liquid Storage Tanks

Abstract

The leakage process simulation of a Natural Gas Liquid (NGL) storage tank requires the simultaneous solution of the NGL’s pressure, temperature and phase state in the tank and across the leak hole. The methods available in the literature rarely consider the liquid/vapor phase transition of the NGL during such a process. This paper provides a comprehensive pressure-temperature-phase state method to solve this problem. With this method, the phase state of the NGL is predicted by a thermodynamic model based on the volume translated Peng-Robinson equation of state (VTPR EOS). The tank’s pressure and temperature are simulated according to the pressure-volume-temperature and isenthalpic expansion principles of the NGL. The pressure, temperature, leakage mass flow rate across the leak hole are calculated from an improved Homogeneous Non-Equilibrium Diener-Schmidt (HNE-DS) model and the isentropic expansion principle. In particular, the improved HNE-DS model removes the ideal gas assumption used in the original HNE-DS model by using a new compressibility factor developed from the VTPR EOS to replace the original one derived from the Clausius-Clayperon equation. Finally, a robust procedure of simultaneously solving the tank model and the leak hole model is proposed and the method is validated by experimental data. A variety of leakage cases demonstrates that this method is effective in simulating the dynamic leakage process of NGL tanks under critical and subcritical releasing conditions associated with vapor/liquid phase change.

1. Introduction

Natural Gas Liquids (NGLs) are flammable mixtures consisting of light hydrocarbon products (ethane, propane, isobutane, pentane and some heavier species). NGL has a variety of applications, being mainly used for heating appliances, cooking equipment, vehicles and in the chemical industry [1]. For safety and convenience, the tank plays an important role in the NGL storage and transportation processes. Although NGL tanks’ operators always take care during all operations, accidental releases of NGL still happen worldwide, and the subsequent dispersion followed by possible fires and explosions usually leads to fatal incidents and property damages. Over the past few years, NGL tank leakage accidents have killed hundreds of people [2,3]. Aiming at controlling of such catastrophic accidents, the leakage simulation technology of NGL storage tanks was developed to capture two essential parameters: (1) the release mass flow rate, to estimate the fire or explosion consequences; (2) the total time required to empty the tank, to evaluate the fire duration and time limit for rescue.

Accurate simulation of NGL release from storage tanks is a tough task because the leakage process of NGL is very complex [4,5]. Normally, the NGL in the tank is a saturated liquid [6]. Once the pressure drops down to the saturation vapor pressure of the NGL during the leakage process, the NGL would evaporate, causing the coexistence of liquid and vapor phases in the tank and the leak hole [7,8]. Additionally, the temperature would also change with the decreasing pressure due to the energy conservation law [9]. On the other hand, the phase state and thermodynamic properties (density, heat capacity, etc.) of the NGL in turn affect the leakage mass flow rate and related dynamic changes in the pressure and temperature. Therefore, a comprehensive leakage simulation model of NGL tanks should have four features: (1) be applicable to vapor/liquid single-phase and two-phase mixtures, (2) accurate prediction of the thermodynamic properties associated with the pressure and temperature changes, (3) accurate prediction of dynamic leakage mass flow rate and the pressure drop caused by mass reduction, and (4) simultaneous modeling and solving of points (2) and (3).

Over the past years, many achievements regarding the simulation of the tank leakage process have been published and can be classified into two main categories. The first category mainly aims at single-phase fluid tanks, especially natural gas tanks and oil tanks [10,11]. Obviously, those models are not applicable to the leakage simulation of NGL tanks associated with liquid-vapor phase transitions.

The second category mainly focuses on calculating the mass flow rate of a two-phase fluid flowing through specific equipment, e.g., safety valves, orifices, and nozzles. A notable method proposed by Henry and Fauske [12] accounts for the interphase heat, mass and momentum transfers of single-component fluid flowing through convergent nozzles. The Henry-Fauske model is suitable for single-phase and two-phase compressible fluids and thus has been extensively adopted by many commercial tools, e.g., Schlumberger OLGA. Another important model, named the Omega method, was developed by Leung [13,14,15], and is recommended to size the pressure relief valves for either flashing or non-flashing flow according to the API RP520 document [16]. However, in this method, the vapor-liquid equilibrium curve described by the Clausius-Clapeyron equation is not valid for multicomponent mixtures. In particular, the change in the mixed system volume with the pressure is described by an approximately linear relation [17]. Diener and Schmidt proposed a Homogeneous Non-Equilibrium-Diener/Schmidt (HNE-DS) model that adds a boiling delay coefficient to the compressibility factor of the original Omega method [18,19,20,21]. The HNE-DS method is adopted as a basis of the ISO 4126-10 International Standard [22] and covers the thermodynamic non-equilibrium of the boiling nucleation of a saturated liquid, and thus it is able to predict the evaporation process of the saturated NGL. Moreover, the HNE-DS model has similar accuracy to the Henry-Fauske model but has simpler formulas [23]. Unlike the previous methods, Raimondi [17] proposed a method of calculating the critical flow condition for the pressure safety device, which uses a specific technique to limit the maximum flow velocity of fluid flowing through the leak hole to the local sonic velocity and calculates the mass flow rate according to the product of the leak hole area, fluid flow velocity and density. Kanes et al. [24] also assumed that the maximum fluid flow velocity is equal to the local sonic velocity. In fact, the maximum flow velocity of the two-phase fluid flowing through a leak hole cannot always reach the local sonic velocity. We will discuss this problem in Section 3.

In summary, there are many limitations in existing methods of calculating the leakage mass flow rate through leak holes. Since the leakage mass flow rate is the basis for the simulation of a tank’s release process, a new method applicable to multi-component NGL should be developed. Moreover, a comprehensive model for the NGL tank’s leakage process simulation is needed to be built, which should be able to predict the NGL thermodynamic properties, the phase state, and dynamic changes in pressures and temperatures during the leakage process.

The purpose of this paper is to establish a dynamic leakage simulation model as well as its solution method. In what follows, NGL thermodynamic property models based on the volume translated Peng-Robinson Equation of State (VTPR EOS) are briefly introduced. Then, a comprehensive model to simulate the dynamic leakage process of the NGL tank is built, which covers the compressibility, liquid/vapor phase transition and non-equilibrium evaporation effect of the NGL. After that, a robust solution procedure of the model is studied and validated by experimental data. Finally, eight cases are demonstrated to show the features of the model.

2. NGL Thermodynamic Properties Model

2.1. Phase Behavior of NGLs

NGL is composed mainly of ethane, propane, butane, pentane or their mixtures. Compositions of four typical NGL mixtures are listed in

Table 1. Compositions of some NGL mixtures (mol %).

Component	NGL1	NGL 2	NGL 3	NGL 4
CH4	0.00	0.00	0.00	0.00
C2H6	8.65	8.73	13.07	14.30
C3H8	47.68	47.23	53.33	53.59
iC4H10	19.26	18.99	16.85	15.45
nC4H10	24.06	24.10	14.66	14.06
iC5H12	0.33	0.88	1.29	1.56
nC5H12	0.01	0.07	0.74	0.98
C6+	0.00	0.00	0.06	0.06

. The Pressure-Temperature phase diagrams of NGL1 and NGL3 samples are depicted in Figure 1.

Figure 1 shows that the P-T diagram is divided into three different regions by the phase envelope [25,26]. The region on the left side of the phase envelope is the liquid phase region, and the portion on the right side is the vapor phase region. The portion covered by the phase envelope is the two-phase region. Normally, the NGL remains in the liquid phase region in the tank. However, if a leakage accident happens, the tank pressure and temperature may follow a T-P curve depicted in Figure 1 [17]. When the pressure in the tank drops down to the saturation vapor pressure of the NGL, part of the liquid NGL evaporates and the vapor phase appears. Finally, there might be no liquid phase in the tank due to the tank’s pressure will eventually reach the atmospheric pressure. In other words, a complete leakage process typically involves three stages: liquid leakage, liquid-vapor two-phase leakage, and the possible vapor leakage. Thus, the accurate calculation on the NGL phase behavior is a fundament of the leakage simulation. In this paper, the thermodynamic models based on the VTPR EOS are employed to calculate the phase change as well as the corresponding thermodynamic properties (e.g., density, entropy, and enthalpy) of the NGL. This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation as well as the experimental conclusions that can be drawn.

2.2. Thermodynamic Model

The VTPR EOS [27] is built based on the Peng-Robison (PR) EOS. The PR EOS [28] is given by Equations (1)–(3) as follows:

$$P=\frac{RT}{v_m-b}-\frac{a}{v_m\left(v_m+b\right)+b\left(v_m-b\right)}$$

(1)

$$a=0.45724\frac{R^2{T}_c^2}{P_c}\alpha$$

(2)

$$b=0.0778\frac{R{T}_c^2}{P_c}$$

(3)

where P and P_c refer to the fluid pressure and critical pressure, respectively, kPa; T and T_c refer to the temperature and critical temperature, respectively, K; v_m is the molar volume, m³/kmol; R is the gas constant, 8.314 kJ/(kmol·K). P_c, T_c and ω for each component are given in many references [29].

For hydrocarbon components, the following volume translated terms are applied to correct the molar volume and co-volume parameters [30]:

$$\tilde{b}=b+c$$

(4)

$${\tilde{v}}_m=v_m+c$$

(5)

$$c=c_1+c_2\left(T-288.15\right)$$

(6)

where ${\tilde{v}}_m$ and $\tilde{b}$ are specific molar volume and co-volume in the VTPR EOS, m³/kmol; c is the volume translation term, m³/kmol. c₁ and c₂ are component specified constants defined in the references [30]. When the VTPR EOS is applied to mixtures, the classical Van der Waals mixing rule is employed to calculate the parameters a, b, and c for mixtures.

Based on the VTPR EOS, the thermodynamic properties regarding the leakage simulation of NGL tanks are expressed as follows:

The density:

$$\rho =\frac{P}{ZRT}$$

(7)

The enthalpy:

$$h=h^0+\left(Z-1\right)RT+\frac{T\frac{da}{dT}}{2\sqrt{2}b}\ln \left[\frac{Z+\left(1+\sqrt{2}\right)bp/\left(RT\right)}{Z+\left(1-\sqrt{2}\right)bp/\left(RT\right)}\right]$$

(8)

The entropy:

$$s-s^0=-R\ln \frac{\rho RT}{101.325}-\frac{\ln \left(3-2\sqrt{2}\right)\frac{da}{dT}}{2\sqrt{2}b}-\frac{\frac{da}{dT}}{2\sqrt{2}b}\ln \left[\frac{1+\left(1+\sqrt{2}\right)b\rho }{1+\left(1-\sqrt{2}\right)b\rho }\right]$$

(9)

where Z is the compressibility factor; ρ is the density, kg/m³; h and h⁰ refer to the fluid enthalpy and ideal gas enthalpy, respectively, J/kg; s and s⁰ are the fluid entropy and ideal gas entropy, respectively, J/(kg·K); The methods of calculating h⁰ and s⁰ can be found in the API Technical Databook [31].

In addition to the thermophysical properties, the liquid/vapor phase change during the leakage process can also be captured by use of the flash algorithm based on the VTPR EOS. Once the phase state and the composition of each phase are obtained from the flash algorithm, corresponding thermodynamic properties can be calculated from Equations (1) to (9). The details with regard to the flash algorithm are referred to Michelsen [25,32] and Zhu and Okuno [33]. The thermodynamical model provides a basis for the compositional simulation of the tank’s release process.

3. The NGL Tank Leakage Simulation Model

3.1. Leakage Process of the NGL Tank

The leakage process of a pressure vessel is generally divided into two flow states, namely the critical flow and the subcritical flow [24,34]. The critical flow refers to the condition that the leakage mass flow rate keeps the maximum value and is independent of the pressure downstream of the leak hole when the upstream pressure is held constant. Conversely, if the leakage mass flow rate is dependent on the downstream pressure when the upstream pressure is held constant, the flow state is attributed to be a subcritical flow. The critical flow occurs when the tank pressure is relatively high. With the decreasing of the tank pressure, the critical flow gradually transfers to the subcritical flow [34,35]. A simulation model should be applicable to both of the critical and subcritical flow conditions.

When calculating the leakage mass flow rate, the tank’s pressure is usually treated as an a priori known parameter, and the pressure at the outlet of the leak hole could be regarded as the atmospheric pressure. This is true for the subcritical flow condition, so the leakage mass flow rate can be easily calculated according to the energy conservation law [20]. However, the outlet pressure under the critical flow condition might not be equal to the atmospheric pressure [17,34]. Here, the outlet pressure of the leak hole under the critical flow condition is also named as the critical flow pressure.

A number of methods based on the ideal gas assumption or on the specific iteration procedures are suggested to calculate the critical flow pressure. The API RP 581 document [36] uses Equation (10) to calculate the outlet pressure based on the ideal gas assumption as follows:

$$P_{\mathrm{out}}=P_{\mathrm{atm}}{\left(\frac{k_v+1}{2}\right)}^{\frac{k_v}{k_v-1}}$$

(10)

where k_v is the isentropic exponent; P_atm is the atmospheric pressure, kPa; P_out is the outlet pressure under the critical flow condition. Unlike the method based on the ideal gas assumption, another one method assumes that the fluid flow velocity is equal to the local sound speed. Thus, the outlet pressure can be calculated based on the energy and momentum conservation laws [17,34]. When this method is used, the momentum change across the leak hole can be calculated from Equation (11):

$$\Delta P=\frac{1}{2}{\rho }_{\mathrm{out}}{v}_{\mathrm{s}}^2$$

(11)

where, ρ_out is the fluid density at the outlet, kg/m³; v_s is the local sound speed at the critical flow.

These methods are not applicable to the two-phase NGL leakage process because the multi-component NGL is obviously inconsistent with an ideal gas hypothesis, and the NGL flow velocity cannot always reach the local sound speed as well [12]. A n-butane relief case given in Raimondi’s work [17] is taken as an example to express this phenomenon. In that paper, the pressure and temperature in the tank are 4000 kPa and 145 °C, and the critical flow pressure and temperature are 3400 kPa and 142.9 °C, respectively. The NGL density and the local sound speed at the critical flow condition are 285.51 kg/m³ and 209.1 m/s, respectively. Assuming the inlet flow velocity of the NGL is 0 m/s and the maximum liquid leakage flow velocity is equal to the local sound speed, the pressure difference calculated from Equation (11) is 6241 kPa, which indicates that the difference between the tank’s pressure and the critical flow pressure is not high enough to accelerate the NGL flow velocity to the local sound speed. Therefore, an efficient method is needed to predict the outlet pressure for further leakage simulations.

3.2. Simulation Model

3.2.1. The HNE-DS Model

This paper adopts the Homogeneous Non Equilibrium-Diener/Schmidt (HNE-DS) method [18,23] to calculate the outlet pressure, which accounts for the boiling delay effect of the fluid and can be regarded as an improved version of the widely used Omega method [14]. In the HNE-DS model, the outlet pressure ratio (η_b) and the non-equilibrium critical pressure ratio (η_cNE) are defined by Equations (12) and (13):

$${\eta }_b=\frac{P_{\mathrm{out}}}{P_{\mathrm{in}}}$$

(12)

$${\eta }_{cNE}=\frac{P_{\mathrm{cric}}}{P_{\mathrm{in}}}$$

(13)

where P_in is the inlet stagnation pressure [18]; P_out and P_cric refer to the outlet, and critical flow pressures, respectively.

The critical pressure ratio is related to the nonequilibrium compressibility factor (ω_NE) of the fluid, which is solved from Equation (14):

$${\eta }_{cNE}^2+\left({\omega }_{NE}^2-2{\omega }_{NE}\right){\left(1-{\eta }_{cNE}\right)}^2+2{\omega }_{NE}^2\ln {\eta }_{cNE}+2{\omega }_{NE}^2\ln \left(1-{\eta }_{cNE}\right)=0$$

(14)

Alternatively, η_cNE can be calculated from Equation (15) when ω_NE ≥ 2

$${\eta }_{cNE}^{}=0.55+0.217\ln {\omega }_{NE}-0.046{\left(\ln {\omega }_{NE}\right)}^2+0.004{\left(\ln {\omega }_{NE}\right)}^3$$

(15)

The parameter ω_NE is derived based on the Clausius-Clapeyron equation [14], given by Equations (16) and (17):

$${\omega }_{NE}=\frac{\frac{v}{v_{\mathrm{in}}}-1}{\frac{P_{\mathrm{in}}}{P}-1}=\frac{x_{\mathrm{in}}{v}_{\mathrm{g}.\mathrm{in}}}{v_{\mathrm{in}}}+\frac{C_{\mathrm{pl}.\mathrm{in}}{T}_{\mathrm{in}}{P}_{\mathrm{in}}}{v_{\mathrm{in}}}{\left(\frac{v_{\mathrm{g}.\mathrm{in}}-v_{\mathrm{l}.\mathrm{in}}}{\Delta h_{v.\mathrm{in}}}\right)}^{2}N$$

(16)

$$N={\left[x_{\mathrm{in}}+C_{\mathrm{pl}.\mathrm{in}}{T}_{\mathrm{in}}{P}_{\mathrm{in}}\left(\frac{v_{\mathrm{g}.\mathrm{in}}-v_{\mathrm{l}.\mathrm{in}}}{\Delta h_{v.\mathrm{in}}^2}\right)\ln \left(\frac{1}{{\eta }_c}\right)\right]}^{\tau }$$

(17)

where, x_in is the gas mass fraction at the inlet of the leak hole; T_in is the inlet temperature; v_g.in, v_l.in, v_in refer to the gas, liquid, and mixture specific volumes at the inlet, respectively, m³/kg; v and P refer to the outlet specific volume and the pressure; C_pl.in is the liquid specific heat capacity, J/(kg·K); Δh_v.in is the latent heat of vaporization at inlet, J/kg; τ is a coefficient, τ = 0.6 for the orifice, control valve, short nozzles; τ = 0.4 for the safety valve; N is a coefficient that indicates the nonequilibrium effect of boiling delay; N ≤ 1 and N = 1 refer to the non-equilibrium and equilibrium state, respectively. Moreover, for the frozen (non-flashing) flow, N = 0; η_c is the critical pressure ratio at the equilibrium condition, which is solved from Equation (14) or Equation (15) by use of the equilibrium compressibility factor (ω_N=1) calculated from Equation (16) [18,23].

In Equations (16) and (17), the homogeneous specific volume of the vapor-liquid two-phase mixture is defined as:

$$v_{\mathrm{in}}=x_{\mathrm{in}}{v}_{\mathrm{g}.\mathrm{in}}+\left(1-x_{\mathrm{in}}\right)v_{\mathrm{l}.\mathrm{in}}$$

(18)

Under the critical flow condition, the outlet pressure must be equal to or greater than the critical pressure, so the case η_b ≤ η_cNE represents the critical flow condition, and the outlet pressure should be set as P_out = η_cNEP_in when calculating the leakage mass flow rate. Conversely, the case η_b > η_cNE refers to a subcritical flow, and the outlet pressure is equal to the atmospheric pressure.

Based on the critical pressure calculated from the HNE-DS model, the leakage mass flow rate is then calculated according to isentropic frictionless flow equation as follows [23]:

$$m_{\mathrm{ds}}=\psi \varphi A_{\mathrm{ori}}\sqrt{\frac{2P_{\mathrm{in}}}{v_{\mathrm{in}}}}$$

(19)

where, m_ds is the discharge mass flow rate, kg/s; A_ori is the leak hole area, m²; ψ is the expansion coefficient given by Equation (21); φ is the two-phase slip correction given by Equation (21); φ = 1 for the single-phase flow:

$$\psi =\frac{\sqrt{{\omega }_{NE}\ln \left(\frac{1}{{\eta }_{cNE}}\right)-\left({\omega }_{NE}-1\right)\left(1-{\eta }_{cNE}\right)}}{\left[\omega \left(\frac{1}{{\eta }_{cNE}}-1\right)+1\right]}$$

(20)

$$\varphi =\sqrt{\frac{v_{\mathrm{in}}}{v_{\mathrm{l}.\mathrm{in}}}}{\left\{1+x_{\mathrm{in}}\left[{\left(\frac{v_{\mathrm{g}.\mathrm{in}}}{v_{\mathrm{l}.\mathrm{in}}}\right)}^{1/6}-1\right]\left[1+x_{\mathrm{in}}\left({\left(\frac{v_{\mathrm{g}.\mathrm{in}}}{v_{\mathrm{l}.\mathrm{in}}}\right)}^{5/6}-1\right)\right]\right\}}^{-1/2}$$

(21)

The HNE-DS model is designed to size the devices using only the tank’s pressure, temperature and fluid properties. Hence, the temperature change across the leak hole has no effects on calculation results. However, studies show that great temperature drop happens across the leak hole, and the outlet temperature is essential to predict the formation of the liquid pool outside of the hole [37]. In this paper, we predict the change in temperature by use of an isentropic expansion model and simultaneously solve the temperature with the pressure and leakage flow rate.

Assuming the process of the NGL flowing through the leak hole is frictionless and adiabatic, the corresponding change in temperature can be calculated from the isentropic expansion correlation [4], as expressed by Equation (22):

$$s_{\mathrm{in}}\left(T_{\mathrm{in}},P_{\mathrm{in}}\right)=s_{\mathrm{out}}\left(T_{\mathrm{out}},P_{\mathrm{out}}\right)$$

(22)

where s_in and s_out refer to the NGL’s specific entropy at the inlet and outlet side of the leak hole, respectively, J/(kg·K); P_out is calculated by the HNE-DS model; T_in, P_in can be obtained from the tank simulation method. Therefore, the only one unknown variable in Equation (22) is T_out [38].

3.2.2. Improvement of the HNE-DS Model

The HNE-DS model uses Equation (14) or Equation (15) to calculate the critical pressure ratio. The results are largely dependent on ω_NE solved from Equation (16). However, Equation (16) is built based on two assumptions: (1) the vapor phase behaves like the ideal gas [14]; (2) the single-component Clausius-Clapeyron equation which assumes that the mixture’s specific volume linearly changes with the pressure like the ideal gas [23]. These assumptions are not sufficient enough for the multicomponent NGL mixtures. In this paper, a rigorous ω_NE derived from the VTPR EOS is adopted to replace Equation (16), yielding an improved HNE-DS model.

Following the derivation of the ω_NE in the HNE-DS model [18], the new compressibility factor is given by Equation (23). The related derivation process of Equation (23) is shown in Appendix A. Since the VTPR EOS is suitable for single-phase and two-phase hydrocarbons, Equation (23) can be applied to the pure liquid/vapor and two-phase NGL mixtures:

$${\omega }_{NEnew}=-{\eta }_{cNEnew}\frac{P_{\mathrm{in}}}{v_{\mathrm{in}}}\left[x_{\mathrm{in}}\frac{\mathrm{d}{v}_{\mathrm{g}.\mathrm{in}}}{\mathrm{d}{P}_{\mathrm{in}}}+\left(1-x_{\mathrm{in}}\right)\frac{\mathrm{d}{v}_{\mathrm{l}.\mathrm{in}}}{\mathrm{d}{P}_{\mathrm{in}}}-\frac{C_{\mathrm{pl}.\mathrm{in}}}{\Delta h_{v.\mathrm{in}}}\left(v_{\mathrm{g}.\mathrm{in}}-v_{\mathrm{l}.\mathrm{in}}\right)N\frac{\mathrm{d}{T}_{\mathrm{in}}}{\mathrm{d}{P}_{\mathrm{in}}}\right]$$

(23)

$$\frac{\mathrm{d}{v}_{mi}}{\mathrm{d}{P}_{}}=-{\left[\frac{R{T}_i}{{\left(v_{mi}-b\right)}^2}-\frac{2a\left(v_{mi}+b\right)}{{\left(v_{mi}^2+2b{v}_{mi}-b^2\right)}^2}\right]}^{-1}$$

(24)

$$\frac{\mathrm{d}T}{\mathrm{d}P}={\left[\frac{R}{v_m-b}+\frac{\mathrm{d}a}{\mathrm{d}T}\frac{1}{v_m\left(v_m+b\right)+b\left(v_m-b\right)}\right]}^{-1}$$

(25)

where the subscript i = l or g; the subscript ‘new’ refers to parameters calculated by the improved model; the parameters a, b, v_m, are the same as those used in Equation (1). The specific volume (v) and molar volume (v_m) has the following relationship: v = v_m/M_w, where M_w is the molecular weight, kg/kmol.

In the original HNE-DS model [18], once ω_NE is calculated from Equation (16), η_cNE can then be easily calculated from the explicit Equation (15) when ω_NE ≥ 2. However, in the improved model, Equation (23) shows that ω_NEnew is a function of η_cNEnew, thus ω_NEnew cannot be obtained in priority of solving Equation (14) or (15) for η_cNEnew. In other words, substituting Equation (23) into Equation (14) or (15) results in a nonlinear equation in terms of η_cNEnew, which should be solved by use of the bisection method or the Newton method. According to the resulting η_cNE, ω_NE is then directly calculated from Equation (23). That is the major difference in solving the original and improved HNE-DS models.

Figure 2 shows the comparison of original ω_NE, η_cNE and the improved ω_NEnew, η_cNEnew. The composition of the NGL is set as C₂H₆ 50 mol % + C₃H₈ 50 mol %. The NGL’s bubble point and dew point pressure at 300 K are 1730 kPa and 2427 kPa, respectively. In order to observe the effect of the NGL’s phase state on ω_NE and η_cNE, in this case, the pressure is designed to change from 1600 to 2500 kPa while keeping the temperature at 300 K so that the NGL is able to transfer from the pure vapor phase to the pure liquid phase. Figure 2 demonstrates that the original and improved models give the similar tendency regarding ω_NE and ω_NEnew that both of them firstly increase with increasing vapor mass fraction, and then decreases when x_in is close to the unity. This tendency is in accordance with the fact that the compressibility factor has a low value for the single-phase fluid [14,15]. It should be noted that ω_NE in the original HNE-DS model [18] ignores the compressibility of the liquid phase, thus it yields ω_NE = 0 when x_in = 0. On the contrary, ω_NEnew is calculated to be 0.0066 by the improved model when x_in = 0.

On the other hand, it is well known that the original method generally gives conservative sizing results for a safety valve or a nozzle for saturated two-phase flow and conversely, it overestimates the size when the fluid at the inlet is a sub-cooled liquid or a boiling liquid with low x_in [19]. This problem can be expressed by the under-estimation of η_cNE for the saturated two-phase flow and to the over-estimation of η_cNE for the sub-cooled liquid or boiling liquid with low x_in. It can be observed in Figure 2 that improved η_cNEnew is slightly higher than the original η_cNE when x_in is in the range of 0.41 to 0.95, whereas the improved η_cNEnew is lower than the original η_cNE when x_in is lower than 0.41. In particular, the improved η_cNEnew is significantly lower than the original η_cNE when x_in is approaching zero. Thus, η_cNEnew obtained from the improved model is probably better than η_cNE obtained from the original HNE-DS model_. In order to clarify the difference of calculating ω_NE and η_cNE by use of the original and improved models, a detailed example is shown in Appendix B.

Additionally, in the original HNE-DS model, Equations (14) and (15) are not valid to calculate η_cNE for the compressible single liquid phase fluid due to ω_NE is calculated as zero for the pure liquid phase. Consequently, Schmidt [19] used η_cNE = η_b for the single liquid phase flow, yielding a low η_cNE value for the single liquid-phase flow. The improved HNE-DS model accounts for the compressibility of the liquid phase, thus it is able to give a non-zero ω_NEnew for single liquid-phase flow, finally yielding a low η_cNEnew under the single liquid-phase flow condition. For example, if P_out = 101.325 kPa, P_in = 2500 kPa and T_in = 300 K, η_cNE obtained from Schmidt [19] is 0.04053, while η_cNEnew calculated from the improved model is 0.1035.

3.2.3. The Tank Simulation Model

In this subsection, the tank simulation model is built to simulate the change of the pressure and temperature in the tank over the leakage duration time. Here, following two assumption are adopted (1) the NGL in the tank is either a homogeneous liquid-vapor mixture or a pure liquid/vapor fluid; (2) there is no heat transfer to or from surroundings; (3) the whole dynamic leakage process can be discrete into a number of time steps, and the tank’s model is to be applied at each time step.

Based on these assumptions, the expansion process caused by the mass reduction in the tank can be treated as an isenthalpic process [38], which indicates that the specific enthalpy of the NGL in the tank is held a constant. In other words, if the tank pressure is given, the tank temperature can be calculated from the isenthalpic equation as expressed by Equation (26):

$$h\left(T_{\mathrm{in}},P_{\mathrm{in}}\right)=\mathrm{constant}$$

(26)

Additionally, one more equation is needed to calculate the tank’s pressure. This paper simulates the change in pressure by use of the PVT relationship of the NGL. At each time step, the released number of NGL moles is calculated by:

$$\frac{\mathrm{d}n}{\mathrm{d}t}=-\frac{m_{\mathrm{d}s}}{M}$$

(27)

where n is the number of NGL moles in the tank, kmol; m_ds is the leakage mass flow rate, kg/s; M is the NGL molar weight, kg/kmol; t is the time, s. Assuming the leakage mass flow rate at each time step (e.g., the k-th time step) is a constant [4], integrating Equation (27) with respect to the time yields Equation (28):

$$\Delta n^k=-\frac{m_{\mathrm{ds}}^k}{M}\Delta t^k$$

(28)

where the superscript k refers to the time step k; Δn represents the released number of NGL moles in a specified time step; Δt is the length of the time step k. Thus, the residual number of NGL moles in the tank at the beginning of the k + 1 time step can be calculated from Equation (29):

$$n_{}^{k+1}=n^k-\frac{m_{\mathrm{ds}}^k}{M}\Delta t^k$$

(29)

where n^k+1 is the number of NGL moles in the tank at beginning of the k + 1 time step; n^k is calculated from Equation (30):

$$n^k=\frac{{\rho }^{k}V}{M}$$

(30)

where, ρ^k is the NGL density in the tank at the k-th time step, kg/m³; V is the physical volume of the tank, m³. In Equation (29), the time step Δt is important for solving the model because it influences both of the computation speed and accuracy. Generally, both of the computation accuracy and required total computation time increase with decreasing the time interval. A smaller time step leads to a higher computation accuracy and a longer total computation time. Thus, a proper time step should be determined to balance the computation time and accuracy. Here, based on the ratio of the residual number of NGL moles in the tank to the current leakage flow rate, we proposed a method to calculate the time step, as expressed by Equation (31):

$$\Delta t^k=\frac{n^k}{Q\Delta n^k}$$

(31)

where Q is an adjustable parameter to control the time step. Figure 3 shows curves of the first-derivative of Equation (31) with respect to Q, where the range of n^k/Δn^k is from 100 to 100,000 which covers most of the leakage cases in industry. It is demonstrated that, even for the case n^k/Δn^k = 100,000, the curve approaches a value close to zero when Q is larger than 200. That means, increasing Q will not significantly decrease the time step if Q is larger than 200. Thus, we set Q = 200 in this paper. As a result, Equation (31) offers a way to dynamically adjust the time step according to the residual NGL in the tank and the current leakage mass flow rate.

In Equations (29) and (30), $m_{ds}^k$ and ρ^k are determined by the priority known $P_{\mathrm{in}}^k$, and $T_{\mathrm{in}}^k$, therefore n^k+1 can be calculated from Equation (29). The task of solving the tank model is to find a specified $P_{\mathrm{in}}^{k+1}$ and $T_{\mathrm{in}}^{k+1}$ at which the number of NGL moles in the tank is equal to n^k+1. According to the general equation of state, if the NGL is in the single liquid- or vapor-phase, the relationship between the tank’s pressure and temperature the is described as follows:

$$P_{\mathrm{in}}^{k+1}V=n^{k+1}ZR{T}_{\mathrm{in}}^{k+1}$$

(32)

For the two-phase mixture, Equation (32) is formulated as:

$$P_{\mathrm{in}}^{k+1}V=n^{k+1}\left[x_{\mathrm{in}}{Z}_{\mathrm{g}}+\left(1-x_{\mathrm{in}}\right)Z_{\mathrm{l}}\right]R{T}_{\mathrm{in}}^{k+1}$$

(33)

In Equations (32) and (33), the parameters x_in, Z_g, and Z_l can be calculated based on the given $P_{\mathrm{in}}^{k+1}$ and $T_{\mathrm{in}}^{k+1}$, and the VTPR EOS. Thus, $P_{\mathrm{in}}^{k+1}$ and $T_{\mathrm{in}}^{k+1}$ are the only two independent variables in Equations (26) and (33). It seems that $P_{\mathrm{in}}^{k+1}$ and $T_{\mathrm{in}}^{k+1}$ could be simultaneously solved by sequentially applying a direct substitution algorithm [25] to Equations (26) and (33). This approach could firstly start from an initial estimate of $P_{\mathrm{in}}^{k+1}$, and then updates $T_{\mathrm{in}}^{k+1}$ according to Equation (26). Then, the updated $T_{\mathrm{in}}^{k+1}$is substituted into Equation (33) as a priority known parameter, and a new $P_{\mathrm{in}}^{k+1}$ is solved from the resulting equation. If the difference between the initial and the new $P_{\mathrm{in}}^{k+1}$ is less than a tolerance value, the algorithm converges; otherwise, the new $P_{\mathrm{in}}^{k+1}$ is set as a new initial estimate to repeat the previous iteration process until the algorithm converges. Nevertheless, it comes to a surprise that this direct substitution algorithm cannot stably converge due to x_in used in Equation (33) is extremely sensitive to the pressure $P_{\mathrm{in}}^{k+1}$ and temperature $T_{\mathrm{in}}^{k+1}$, and Z_g is the dominant factor in Equation (33) when vapor phase emerges in the tank.

Alternatively, one better option is using the number of moles as the convergence criterion to solve the model. In order to start the algorithm, the initial estimates of $T_{\mathrm{in}}^{k+1}$ and $P_{\mathrm{in}}^{k+1}$ are set as:

$$T_{\mathrm{in}}^{k+1}=T_{\mathrm{in}}^k$$

(34)

$$P_{\mathrm{in}}^{k+1}=\frac{n_{}^{k+1}}{n^k}{P}_{\mathrm{in}}^k$$

(35)

By using the estimated $T_{\mathrm{in}}^{k+1}$ and $P_{\mathrm{in}}^{k+1}$, the NGL density ρ^k+1 is then obtained from the thermodynamic models presented by Equations (1) to (7). Thus, the number of NGL moles at $T_{\mathrm{in}}^{k+1}$ and $P_{\mathrm{in}}^{k+1}$ is given by Equation (36):

$$n_{\mathrm{cal}}^{k+1}=\frac{{\rho }_{}^{k+1}V}{M}$$

(36)

The distance between the actual (n^k+1) and calculated ($n_{\mathrm{cal}}^{k+1}$) number of NGL moles in the tank is used to modify the $P_{\mathrm{in}}^{k+1}$, given by Equation (37). It is shown that $P_{\mathrm{in}}^{k+1}$ increases if $n_{\mathrm{cal}}^{k+1}<n^{k+1}$ and conversely, $P_{\mathrm{in}}^{k+1}$ decreases if $n_{\mathrm{cal}}^{k+1}>n^{k+1}$. Also, $P_{\mathrm{in}}^{k+1}$ keeps a constant when $n_{\mathrm{cal}}^{k+1}=n^{k+1}$. Once the new $P_{\mathrm{in}}^{k+1}$ is determined, $T_{\mathrm{in}}^{k+1}$ can be solved from Equation (26):

$$P_{\mathrm{in}}^{k+1}=P_{\mathrm{in}}^{k+1}\left[1-\kappa \frac{\left(n_{\mathrm{cal}}^{k+1}-n_{}^{k+1}\right)}{n_{}^{k+1}}\right]$$

(37)

where κ is a relaxing factor from 0 to 1.

The difference between n^k+1 and $n_{\mathrm{cal}}^{k+1}$ defined by Equation (38) is adopted here as the convergence criterion of the algorithm. Since n^k+1 in Equation (38) is a fixed value at each time step, the proposed algorithm is more stable than the method that uses the modification on $P_{\mathrm{in}}^{k+1}$ as the convergence criterion:

$$\left|n^{k+1}-n_{\mathrm{cal}}^{k+1}\right|<\epsilon$$

(38)

where ε is the convergence tolerance, we set ε = 0.001 kmol.

4. Solution Method and Model Validation

4.1. Solution Procedures

• In Section 3, the improved HNE-DS model and the NGL tank model have been described. The former one is used to predict the leakage mass flow rate and the temperature at the outlet of the leak hole. The latter one is used to simulate the pressure and temperature in the tank. Since the tank’s pressure and temperature are necessary parameters for solving the leak hole model, the tank model and the improved HNE-DS model should be sequentially solved at each time step until the tank pressure reaches the atmospheric pressure. In other words, at time step k, the tank’s pressure is updated based on the leakage mass flow rate calculated from the improved HNE-DS model. The new pressure and temperature in the tank are then set as the inlet parameters of the leak hole model in the next time step k + 1. Hence, the HNE-DS model and the tank model should be simultaneously solved. The solution procedure is described as follows:
• Step 1: Set the iteration time k = 0; input the initial tank pressure $P_{\mathrm{in}}^k$, tank temperature $T_{\mathrm{in}}^k$.
• Step 2: Set the simulation time interval Δt according to Equation (31).
• Step 3: Perform the flash calculation and calculate the molar weight, compressibility factor Z, specific volume, entropy, enthalpy and specific heat capacity of the NGL using the VTPR EOS and Equations (1)–(9).
• Step 4: Calculate η_cNEnew and ω_cNEnew from Equations (14) and (23).
• Step 5: Calculate η_b from Equation (12). If η_b ≤ η_cNEnew, set P_out = η_cNEnewP_in; otherwise, set P_out = P_atm.
• Step 6: Calculate the leakage mass flow rate from Equation (19).
• Step 7: Calculate the temperature at the outlet of the leak hole from Equation (22).
• Step 8: Calculate n^k+1 from Equations (29) and (30).
• Step 9: Calculate estimated tank parameters $P_{\mathrm{in}}^{k+1}$ and $T_{\mathrm{in}}^{k+1}$ using Equations (34) and (35).
• Step 10: Calculate $n_{\mathrm{cal}}^{k+1}$ from Equation (36), and check the convergence criterion Equation (38). If Equation (38) is satisfied, go to Step 12. Otherwise, go to Step 11.
• Step 11. Update $P_{\mathrm{in}}^{k+1}$ and $T_{\mathrm{in}}^{k+1}$ by subsequently use of Equations (37) and (26), and go back to step 10.
• Step 12: If $P_{\mathrm{in}}^{k+1}<P_{\mathrm{atm}}$, stop. Otherwise, set t = t + Δt, k = k + 1 and go back to step 2.

The corresponding solution flow chart is depicted in Figure 4.

4.2. Model Validation

Any experiments regarding NGL release from a tank are very risky because of the potential for explosions. To the best of our knowledge, rare experimental results have been reported in published papers. Botterill et al. [39] carried out an experiment on n-butane release from a short pipe which is assumed to be a tank in that paper. The pipe was 5 m in length and 12 mm in internal diameter. The initial pressure was 1.0 MPa, and the ambient temperature was 7 °C. The diameter of the leakage hole was set as 5 mm. Comparisons of the predicted tank pressure and temperature during the whole leakage process against the experimental data are shown in Figure 5 and Figure 6.

Based on the phase state of n-butane, each figure can be divided into two regions, namely the liquid phase region and the vapor phase region. This feature is in accordance with the pure n-butane’s P-T phase behavior, in which the saturation vapor pressure curve divides the whole P-T region into two sub-regions: the liquid phase region and the vapor phase [25,32]. When the pressure and temperature lie below the saturation vapor pressure curve on the P-T phase diagram, n-butane is in the vapor phase, otherwise, n-butane is in the vapor phase [40]. Since the liquid and vapor phases are different in thermophysical properties, the related leakage parameters suddenly change at the phase transition boundary. The inflection points on Figure 5 and Figure 6 show the phase transition boundary between the liquid phase and vapor phase. It is demonstrated that n-butane is in the liquid-phase state in the first 0.67 s of the leakage process and after that, n-butane evaporates.

Figure 5 presents that the pressures predicted by the original and improved models match the experimental data very well. It is shown that the pressure rapidly drops in the first 0.67 s due to the low compressibility of the pure liquid n-butane. Once n-butane evaporates, the pressure slowly drops because of the high compressibility of the vapor phase. It should be noted that the improved model gives slightly different results in the vapor phase region. Such differences are mainly caused by the different compressibility factors given by Equations (22) and (29). In particular, Equation (22) is developed based on the ideal gas assumption [14], whereas Equation (29) uses the VTPR EOS to represent the compressibility factor. These comparisons prove that the improved HNE-DS model, the tank simulation model, and the dynamic simulation procedures are capable of simulating the dynamic leakage process of NGL tanks associated with the liquid/vapor phase transition.

Figure 6 shows the calculated tank temperatures based on the original and the improved HNE-DS model. It should be noticed that both the original and improved HNE-DS model do not provide a temperature calculation method. That means, the temperatures given in Figure 6 are all calculated by Equation (26) used in the tank simulation model. Since Equation (26) has only two variables, namely the pressure and the temperature, the differences between the two curves in Figure 6 are mainly caused by the pressure differences. Also, the isenthalpic expansion law implied by Equation (26) indicates that the temperature change during the isenthalpic expansion process is positively related to the pressure drop [38]. Figure 5 shows that the improved model gives higher pressures in comparison to the original model, so the temperatures based on the improved model are naturally higher than temperatures based on the original model.

Figure 6 also shows that the calculated temperatures firstly increases in the liquid phase region, and then decreases in the vapor phase region. The leading reason is that the Joule-Thomson coefficient of the n-butane in the tank gradually transfers from a negative value to a positive value during the leakage process, which is true for the isenthalpic expansion process [41].

5. Results and Discussion

5.1. The Natural Gas Tank

In this subsection, we start from the simplest single-phase leakage process of the natural gas tank to discuss the simulation results. The natural gas mixtures are mainly composed of CH₄, and always keep as a single vapor phase during the leakage process. Due to the absence of the phase change, the leakage process of the natural gas tank is simpler than that of the NGL tank. The natural gas composition is given in

Table 2. Natural gas composition.

Component	N2	CO2	CH4	C2H6
mole fraction,%	1.0	2.0	95.0	2.0

. The initial tank pressure and temperature are 3000 kPa and 290 K, respectively. The tank diameter is 5 m, and the height is 3 m. The leak hole is set to be a circle with a diameter of 40 mm. The simulation results are depicted from Figure 7, Figure 8, Figure 9 and Figure 10.

Figure 7 depicts changes in the residual natural gas mass in the tank and the leakage mass flow rate over the leakage time. It is shown that the leakage process lasts 18.7 min. The residual mass in the tank and the leakage mass flow rate decrease over time. Since there is no phase change and the natural gas composition is constant, the dominant factor of influencing the leakage mass flow rate is the pressure difference between the inlet and outlet of the leak hole. Figure 5 demonstrates that such a pressure difference also decreases over the leakage duration time. The leakage process stops when the tank pressure drops down to the atmospheric pressure.

In the improved HND-DS model, the pressure difference between the inlet and outlet of the leak hole is represented by the parameter η_cNE in Equation (13), and the corresponding curve is shown in Figure 9. According to Step 5 of the solution procedure, the case η_b ≤ η_cNE refers to a critical flow, and the case η_b > η_cNE refers to a subcritical flow. If the case is a critical flow, the outlet pressure is calculated as η_cNEP_in, otherwise, the outlet pressure is set to the atmospheric pressure. In Figure 8, the intersection point of η_cNE and η_b is the point that transfers from the critical flow to subcritical flow.

Figure 10 presents the temperatures in the tank and at the outlet of the leak hole. The natural gas in the tank follows the iso-enthalpy expansion law, thus the tank temperature is dependent on the tank pressure. The results show that the tank temperature drops from 17 °C to 4.3 °C over time which is accordance with the tendency of the tank pressure. Unlike the temperature in the tank, the temperature at the outlet of the leak hole is dependent not only on the pressure difference but also on the tank temperature. Generally, the outlet temperature decreases with decreasing tank temperature, and increases with decreasing pressure difference. Under the critical flow condition, although the pressure difference keeps decreasing, the rapid decreasing tank temperature is a dominant factor of influencing the outlet temperature, which makes outlet temperature drops over time. On the contrary, under the subcritical flow condition, the tank temperature almost keeps a constant while the decreasing pressure difference becomes a dominant factor in influencing the outlet temperature, resulting in the increase of the outlet temperature. Finally, the tank temperature and the outlet temperature reach the same value since the pressure difference is equal to zero when the leakage process stops.

These results demonstrate that even for the simplest single-phase tank, the coupled behaviors between the temperature, pressure, leakage mass flow rate, are important in influencing the leakage process. In the following subsection, the model is to be applied to NGL tanks associated with the liquid/vapor phase transition.

5.2. The NGL Tank

In this subsection, the proposed model is applied to an NGL tank. The major difference between the leakage process of the NGL tank and that of the natural gas tank lies in the phase state of the fluid. The leakage process of the NGL tank usually couples with the liquid/vapor phase transition, whereas the natural gas keeps on vapor phase all the time. In this case, the diameter, height, initial pressure and temperature, and the leak hole of the NGL tank are set the same values as those of the natural gas tank used in the previous case. The NGL composition is set as the NGL1 sample given in Table 1, and corresponding density curve at 290 K is depicted in Figure 11. The simulation results regarding the leakage mass flow rate, residual NGL mass in the tank, the vapor phase fractions, pressures and temperatures in the tank and at the outlet of the leak hole are depicted from Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16.

Figure 11 demonstrates that the NGL density rapidly increases with increasing pressure in both of the vapor phase and two-phase states. However, in the liquid phase region corresponding to the pressure higher than 822 kPa at 290 K, the density changes very slowly with the pressure, which proves the low compressibility of the liquid NGL. Figure 12 shows the volume fractions of the vapor phase in the tank and at the outlet of the leak hole. It is shown that the NGL is in a single liquid phase state in the first one minute and after that, the vapor phase emerges in the tank. The vapor fraction in the tank also gradually increases over time due to the decrease of the tank pressure that we will discuss later.

Figure 13 shows that the whole leakage process lasts for 56.3 min. In the first one minute, the leakage mass flow rate drops from 48.0 kg/s to 18.0 kg/s, and the accumulated released NGL mass is 1856 kg, which takes 6.31% of the total NGL mass in the tank. In the following 55.3 min, the leakage mass flow rate gradually drops to the value of 0 kg/s and the NGL in the tank is totally released. These results demonstrate that the vapor phase fraction is an essential factor of influencing the NGL leakage mass flow rate for the two-phase leakage process due to the liquid density is one order of magnitude higher than the vapor density.

Figure 14 demonstrates that the tank pressure abruptly drops from 3000 kPa to 758 kPa in the first one minute of the leakage process, and then slowly drops to the atmospheric pressure. This particular phenomenon is not observed in the leakage process of the natural gas tank and can be attributed to the low compressibility of the liquid phase and the liquid/vapor phase change of the NGL. When the NGL is in the single liquid phase, the small releasing amount can greatly reduce the tank pressure due to the low compressibility of the liquid phase. Once the pressure drops down to the saturation vapor pressure of the NGL, changes in the tank temperature and pressure become slower because the volume previously taken by the released NGL can be continually filled by the evaporated NGL.

Figure 12 and Figure 16 demonstrate that the NGL keeps the two-phase state at the end of the leakage process. That is because the tank temperature is low enough to keep the NGL stay in the two-phase region as shown in Figure 14, although the tank pressure reaches the atmospheric pressure. Similarly, the NGL at the outlet of the leak hole is always in the two-phase region, which implies the formation a possible liquid pool outside of the tank [6]. Traditional methods typically ignore the change in the tank temperature, thus they are not able to predict the residual liquid NGL at the end of the leakage process.

5.3. The Analysis on Influencing Factors

5.3.1. The Tank Pressure

We set three different tank pressures at 500 kPa, 1000 kPa and 3000 kPa to research the effect of the tank pressure on the leakage process. The initial tank temperature is 290 K. The other parameters of the tank and the NGL composition are set the same values as those in the previous section. The simulated results regarding the leakage mass flow rate and the residual NGL mass in the tank are depicted in Figure 17 and Figure 18.

The results show that except for the portion in the first one minute, the results are almost identical when the initial tank pressures are equal to 1000 kPa and 3000 kPa. Actually, the NGL densities at 1000 kPa and 3000 kPa are equal to 494.6 kg/m³ and 501.4 kg/m³, respectively. Thus, under these two conditions, the difference between the total NGL mass in the tank is 317 kg, which takes only 1.43% of the total mass. Such a small difference in terms of the NGL mass in the tank doesn’t have a significant influence on the leakage mass flow rate and the total leakage duration time.

However, it should be noted here that the case with the initial pressure 500 kPa is significantly different from other two cases. That is because the NGL is in the two-phase state at 500 kPa and 290 K. The corresponding total NGL mass in the tank is 1286 kg. Hence, this case gives a lower leakage mass flow rate and the shorter leakage duration time. From these results, we can draw the conclusion that initial phase state of the NGL in the tank has a great effect on the leakage process. However, as long as the initial phase state of the NGL is in the liquid phase, the pressure is not a dominant factor in influencing the leakage process anymore.

5.3.2. The Tank Pressure

The leak hole diameter is also an important factor in influencing the leakage process. We set three diameters 20 mm, 30 mm and 40 mm for the leak hole. The simulated results are depicted in Figure 19 and Figure 20. It is easy to observe that the total leakage duration time decreases with increasing the hole diameter. The total duration time of these cases is equal to 225 min, 102 min, and 56 min, respectively. However, the results do not give a quantitative correlation between the hole diameter and the total duration time of the whole leakage process.

6. Conclusions

This paper proposes a comprehensive method to simulate the dynamic leakage process of NGL tanks. The advantages of this method lie in four aspects: (1) it uses a method based on VTPR EOS to calculate the physical properties of the NGL, which is a rigorous way covering wide pressure, temperature and composition ranges; (2) the compressibility factor derived from the Clausius-Clayperon equation in the original HNE-DS model is replaced by a new one developed based on the VTPR EOS. The improved HNE-DS model is capable of predicting the leakage mass flow rate of the pure liquid/vapor phase and two-phase multicomponent mixtures; (3) the temperature changes in the tank and at outlet of the leak hole are considered in the model; and (4) the improved HNE-DS model and the tank model are solved simultaneously.

The proposed method is validated by experimental data and applied to a series of leakage cases of natural gas and NGL tanks. The results reveal details regarding the changes in the leakage mass flow rate, residual NGL mass in the tank, pressures, and temperatures with the leakage duration time. Our findings demonstrated that the liquid/vapor phase transition of the NGL has a great effect on the leakage mass flow rate and total leakage duration time. Under the single-liquid phase releasing condition, the pressure drop rate in the tank and the leakage mass flow rate are much higher than those under the two-phase or single vapor phase releasing conditions. Moreover, the initial phase state of the NGL in the tank is a dominant factor of influencing the leakage process. However, if the NGL is in the single liquid phase at the beginning of the leakage process, the initial pressure in the tank will not obviously affect the leakage process due to the low compressibility of the liquid NGL.

The proposed method may provide more exact values regarding the released mass of NGL in accidents. Also, the simultaneous solution of the pressure, temperature and phase state offers a practical way to predict the potential risks of forming a leakage pool outside of the tank, giving a higher reliability to the calculation of atmospheric dispersion of NGL and related effects.