Numerical Study on Gaseous CO2 Leakage and Thermal Characteristics of Containers in a Transport Ship

Abstract

In general, gaseous CO2, compressed and liquefied, is stored in tanks for efficient transportation [8]. Because the tanks are used for transportation on the ocean [9], the possibility of gas leakage should be examined for confirming the safety of the transportation facilities [10]. [...]to provide a better understanding when designing CO2 cargo facilities, the present study aims to numerically investigate CO2 propagation characteristics inside a container and to examine thermal effects on the cargo walls on which thermal damage might occur. From the results, the initial leakage mass flow rates were estimated to be 2.175 kg/s, 8.698 kg/s, and 19.569 kg/s for crack sizes of 10 mm, 20 mm, and 30 mm, respectively. Because the phase change occurs after the leakage of liquefied CO2, we tried to estimate the required time for phase change theoretically based on the dispersion of flashing jets model [18]. [...]for the sake of simplicity, the present study introduced the effective area concept, as depicted in Figure 1b, for imposing the reasonable boundary condition with the same mass flow rate of gaseous CO2.

1. Introduction

Carbon dioxide (CO₂) is one of the most significant greenhouse gases; however, it is still emitted into the environment from large industrial facilities. Carbon capture and storage (CCS) technologies have been developed as an option to stabilize the concentration of CO₂ in the atmosphere [1]. Generally, the CCS chain consists of three stages: capture, transportation, and storage [2]. Ocean storage is becoming the most promising approach to sequestering CO₂ due to the large capacity of deep saline aquifers [3]. Although pipelines have recently been developed for transporting CO₂, there are major obstacles related to the initial installation, as well as maintenance costs [4]. Alternatively, CO₂ shipping vessels (or water carriers) are more attractive alternatives when considering maintenance costs [5]. To be specific, the capacity of a CO₂ transport ship is determined by the daily production of CO₂ and the sailing period of the ship [6]. Much effort to build CO₂ transport ships has been conducted in the Republic of Korea, because Korea is a peninsula surrounded by the sea [7]. Currently, CO₂ transport is still in the development stage but, because of its importance, a persistent effort is being made regarding the better design of transport ships.

In general, gaseous CO₂, compressed and liquefied, is stored in tanks for efficient transportation [8]. Because the tanks are used for transportation on the ocean [9], the possibility of gas leakage should be examined for confirming the safety of the transportation facilities [10]. The transport ships may be damaged if an accidental CO₂ leak occurs in the container during sailing [11]. When an accidental leak occurs, liquefied CO₂ undergoes rapid vaporization and expansion, causing a significant pressure rise [12]. Even if CO₂ is not classified as a dangerous substance, two important factors clearly affect safety, one of which is the occurrence of personal injury such as suffocation by workers exposed to an environment with a high CO₂ concentration [13]. The other is associated with thermal damage to the tank structure during leakage, because very cold CO₂ gas might impinge on the cargo walls, causing a local change in solid temperature and thermal deformation [14]. Thus, we need a good deal of information on thermal characteristics depending on different scenarios of gas leakage.

In spite of its importance, as noted above, it is very hard to directly measure temperature and CO₂ concentration in a real-scale cargo to obtain information on thermal behavior during leakage. As an alternative, the computational fluid dynamics (CFD) simulation is a promising way to get useful information on the variations in concentration and thermal characteristics during the leakage of CO₂ gas. Thus, to provide a better understanding when designing CO₂ cargo facilities, the present study aims to numerically investigate CO₂ propagation characteristics inside a container and to examine thermal effects on the cargo walls on which thermal damage might occur. We also suggest a new engineering model with an effective area concept to consider the phase change during leakage. Moreover, the influence of crack size on the transient evolution of CO₂ concentration and structure temperature is analyzed.

2. Modeling of Gas Leakage

In the literature [15,16], it is reported that there were serious gas leak accidents through various size of cracks near the pipe connections. Based on these reports, it is assumed that CO₂ gas escape caused by a leakage under the tank flows in a downward vertical direction at the corner of the container. There might be many factors related to leak scenarios, but the current work focuses mainly on the crack size effect on the leakage characteristics. Indeed, various crack sizes, depending on the defect types, have been reported [17]. Based on the reports, we used three leakage scenarios for a crack diameter ranging from 10 mm to 30 mm by assuming that the crack is caused by external damage, not by corrosion defects in the weld portion. The mass flow rate through the crack had to be determined because it must be used as input at the crack position in the CFD simulation. We used the commercial code of DNV PHAST (V.7.0), which includes the phase change model, for predicting the mass flows rate. From the results, the initial leakage mass flow rates were estimated to be 2.175 kg/s, 8.698 kg/s, and 19.569 kg/s for crack sizes of 10 mm, 20 mm, and 30 mm, respectively.

Because the phase change occurs after the leakage of liquefied CO₂, we tried to estimate the required time for phase change theoretically based on the dispersion of flashing jets model [18]. In the case of the pressurized swirling nozzle at 6.8 bars (which was the same as the CO₂ storage pressure condition), the spray cone angle and the breakup length were estimated to be 82.13° and 9.34 mm, respectively, based on the atomization model with Re = 2.05 × 10⁶ and We = 3.99 × 10⁵ [19,20]. The predicted Sauter mean diameter (SMD) was 0.0204 μm, obtained from mechanical break-up criterion [21]. For the in-flight evaporation of water drops with a diameter smaller than 230 μm, the lifetime of a drop can be estimated by t_life = D₀²/(q₀ ΔT), where D₀ is the initial drop diameter [22]. Here, ΔT is the temperature difference between the ambient temperature and the drop temperature, and q₀ is an empirical constant [22]. Finally, the lifetime of CO₂ droplets after the break-up was estimated to be 2.40 × 10⁻²⁶ s. With the droplet velocity and lifetime, a distance required for phase change from liquid to gas was estimated to be 5.7 × 10⁻²⁴ m from the crack, which means that the liquefied CO₂ changes into the gas phase almost instantly, immediately after leaking. Thus, considering break-up length (L_bu), the total distance required for vaporization of CO₂ droplets is less than 0.5% of the total height (H) from the bottom surface, as shown in Figure 1a. It has also been reported that dry ice due to phase change might form and accumulate on the impingement surface according to the operating and environmental conditions. From simulations by the DNV PHAST program under the present conditions, it was found that there was no accumulation of dry ice on the impingement surface because the break-up length was much smaller than the distance between the crack and the bottom surface.

Next, the appropriate inlet boundary condition at the crack position where a rapid phase change occurs had to be determined. Because the phase change occurs rapidly, direct multiphase simulation near the crack becomes very difficult and complicated, resulting in high computational costs as well as causing high numerical uncertainty. As noted before, the total distance required for drop evaporation is much smaller than the impingement distance. Hence, for the sake of simplicity, the present study introduced the effective area concept, as depicted in Figure 1b, for imposing the reasonable boundary condition with the same mass flow rate of gaseous CO₂. Here, it is assumed that the liquefied CO₂ exposed from the crack evaporates spontaneously. Gas phase velocity, v_gas, which is predicted by the DNV PHAST program, considering substantial expansion and phase change, is used as an inlet boundary condition for CFD simulations. The DNV PHAST software assumes one-dimensional homogeneous flow along the expansion zone, which is in thermal equilibrium with zero air entrainment [21]. It also contains the expansion model and the orifice model, based on the conservation of energy and the isentropic process, and is able to estimate gas velocity and mass flow rate. Based on mass conservation, the effective area was calculated by considering the expansion and was used as the modified inlet area, with CO₂ gas velocity and mass flow rate obtained by the DNV PHAST code. The present study estimated the effective crack sizes to be 66 mm, 132 mm, and 198 mm for different crack sizes of 10 mm, 20 mm, and 30 mm, respectively.

It is important to make an appropriate and rapid response after CO₂ leakage. In this work, a transient simulation was conducted until a flow time of 300 seconds to analyze the gas propagation behavior in the early stage of gas leakage was achieved. More detailed numerical conditions are summarized in Table 1.

3. Mathematical Representations for CFD Simulations

For the simulations, a set of conservation equations, as expressed below, were solved simultaneously for mass, momentum, energy, and species.

∂ρ∂t+∂(ρ_ui)∂_xi=0

∂∂t(ρ_ui)+ρ_uj∂_ui∂_xj=−∂p∂_xi+∂_τji∂_xj+_SM

∂∂t(ρE)+∂∂_xi[_ui(ρE+p)]=∂∂_xj[_keff∂T∂_xj−∑k_{hk Jjk}]

∂∂t(ρ_Yk)+∂(ρ_{ui Yk})∂_xi=−∂_Jjk∂_xj

where ρ is the density; τ_ji is the stress; S_M is the external body force; k_eff is the effective conductivity; J_jk is the diffusion flux of species k; and Y_k is the mass fraction of species k. The standard k-ε turbulence model was used, and the standard wall function considering buoyancy effects was employed to describe the turbulent dissipation energy near the wall. The transport equations for the turbulent kinetic energy, k, and the rate of energy dissipation, ε, are given by [23] as follows:

∂∂t(ρk)+∂∂_xi(ρ_uik)=∂∂_xj[(μ+_μtσk)∂k∂_xj]+_Gk+_Gb−ρε+_YM

∂∂t(ρε)+∂∂_xi(ρ_uiε)=∂∂_xj[(μ+_μtσε)∂ε∂_xj]+_C1εk(_Gk+_Gb)−_C2ρ^ε2k

where the turbulent viscosity μ_t can be obtained by

_μt=_Cμρ^k2ε

where G_k represents the generation of turbulent kinetic energy due to mean velocity gradients. G_b is the generation of turbulent kinetic energy due to buoyancy. Y_M represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate to account for the compressibility effect. The empirical constants are given by

_C1=1.44,_C2=1.92,_Cμ=0.09,_σk=1.0,_σε=1.3.

This simulation used a target CO₂ ship that was designed using the hazard identification (HAZID) method, originally developed by the Korea Research Institute of Ships and Ocean Engineering (KRISO). The HAZID method is usually used for obtaining the basic design of a CO₂ ship. It provides useful information about the frequency index, severity index, consequence, location, and hazard event [24]. Figure 2 shows a three-dimensional computational model and the specifications of the CO₂ ship’s optimized container. The computational domain has a size of 30(L) m × 30(W) m × 21(H) m with a total cargo capacity of 18,000 m³. The container has nine tanks for storing liquefied CO₂ in a 3 × 3 square arrangement standing 8.4 m apart from each other. The height and diameter of the cylindrical tanks are 12.5 m and 7 m, respectively. The top and bottom of these tanks are shaped like a spherical cap to distribute high internal pressure up to 6.8 bars. Air pockets in the surrounding regions have an interval of 2 m between the boundary and inner walls, and they are divided by seven lateral and 12 horizontal plates to support the structure. The tween deck is placed horizontally at the height of 14.75 m from the bottom boundary to reinforce the structure of the container.

In the simulation, two separate volumes were constructed for considering conjugate heat transfer between the container and the air pocket. As shown in Figure 3, polyhedral meshes were employed for the entire domain to reduce computational costs. The number of grids selected was approximately 400,000, as determined by the results of the grid independence test. Finer meshes were generated in the vicinity of the crack position to capture rapid variations in temperature and flow motion. Boundary conditions for solving the governing equations were as follows: first, the inlet mass flow rates were estimated to be 2.175 kg/s, 8.698 kg/s, and 19.569 kg/s for different crack sizes of 10 mm, 20 mm, and 30 mm, respectively. The inlet CO₂ gas was leaking at −78.37 °C, which is the saturation temperature under atmospheric pressure. It was assumed that the container is filled with air at room temperature and atmospheric pressure.

The Pressure Implicit with Splitting of Operator (PISO) scheme was adopted for pressure-velocity coupling. Moreover, the first order upwind scheme was used in the discretization of all equations. The density of the mixture was calculated using the ideal gas law, and its specific heat was predicted by the mixture properties of temperature. The initial time step was fixed as 10⁻⁶ s, depending on the Courant number related to flow stability. The total computational time of each case was approximately two weeks for 300 s of flow time. The commercial CFD code (ANSYS FLUENT V.17.0) was used, and the simulation was carried out using a computer system with a 3.50 GHz processor (IntelR CoreTM i7-3770K CPU) and 8 GB RAM with four nodes.

4. Results and Discussion 4.1. Flow Characteristics

For a crack size of 10 mm, corresponding to the effective crack size of 66 mm as listed in Table 1, the velocity field at t = 300 s after leakage is shown in Figure 4a. As mentioned, the maximum velocity was predicted to be 237.7 m/s by using DNV PHAST, and it was used for the input boundary condition at the crack position. Rapid expansion occurs after leakage, as discussed above, and the gas phase with high velocity generated by the phase change and rapid volume expansion was considered in the CFD simulation. In Figure 4a, this CO₂ gas flow impinges on the bottom surface and spreads in a radial direction, similar to the impinging jet flow characteristics [25]. Owing to the flow diffusion, a rapid decrease in flow velocity was observed, e.g., the flow velocity decreased by up to 28 m/s approximately at y = 0.3 m. The velocity profiles at different locations from the bottom surface are presented in Figure 4b. A substantial change in the flow velocity is shown because of the rapid transport of momentum and energy through convection and diffusion in the free jet zone.

As mentioned, the present study used the gas leakage model for appropriate inlet conditions at the prescribed crack position. The well-predicted velocity profile becomes very important because it significantly affects the flow behavior near the bottom surface. For validation, the present study conducted additional simulations under the same experimental conditions [26]. We compared the predicted velocity profile and the radius at half-maximum velocity with the experimental results [26], as shown in Figure 5. The prediction of the mean velocity profile shows very good agreement with the experimental measurements, showing that the current gas leakage model is acceptable for determining the inlet conditions.

After a leak, the decrease of the hydraulic level in the tank leads to a gas pressure increase in the container due to the accumulated mass of leaked CO₂ gas. The present simulation considered the change in the pressure difference between the tank pressure and the environmental pressure, resulting in a decrease in the mass flow rate of the leaked gas. We made a controllable user-defined function embedded in the commercial code of FLUENT to control the mass flow rate varying with time.

Figure 6 presents the predicted volume-averaged static pressure. The static pressure increased linearly with a gradient proportional to the leakage flow rates, depending on the effective crack area. The estimated pressures at t = 300 s were reported to be 1.04 bars, 1.14 bars, and 1.32 bars for crack sizes of 10 mm, 20 mm, and 30 mm, respectively. According to the design guidelines of the transport ship, such ventilation or decompression systems should be operated when the static pressure exceeds 1.2 bars. Thus, the estimated pressure would be useful information to have when assessing the needs of ventilation system operation. Additionally, we considered that the local pressure near the crack would exceed the average pressure value because the CFD simulation was considering non-homogeneous conditions. However, the dynamic pressure of leaking gas on the bottom surface was less than 0.01 bars, indicating that the dynamic pressure effect was negligible, compared to the order of static pressure.

4.2. Accumulation of CO₂ Concentration

Since humans can lose consciousness within a few minutes when exposed to CO₂ gas at a concentration of over 10% [27,28], a CO₂ concentration of 10% can be regarded as a safety threshold. The predicted CO₂ concentration distribution passing through the crack is provided in Figure 7. The results show that the maximum CO₂ concentration occurs at the corner region near the crack location. A much higher CO₂ concentration in comparison to other regions is also predicted around the corner. For efficient evacuation, it is important to analyze a critical time at which the CO₂ concentration reaches 10% at human nasal height. The first critical time is observed at the corner and estimated to be 195 s, 48 s, and 25 s for crack sizes of 10 mm, 20 mm, and 30 mm, respectively. With the increase in the crack size, the critical time is inversely proportional to the square of the crack size depending on the leakage mass flow rate.

From an engineering viewpoint, as shown in Figure 8, it would be useful to define a dangerous volume of CO₂ concentration being over 10%, to characterize three-dimensional leakage behavior. Physically, the leaked CO₂ gas is rarely propagated toward the upper tween deck region because of the buoyancy effect. In the case of the 30 mm crack size, it was reported that the whole region below the tween deck could become a dangerous volume within 300 s. To investigate the propagation tendency, the dangerous volume fraction, defined as the ratio of the dangerous volume to the entire cargo volume, is presented in Figure 9. We see that a substantial difference in dangerous volume depending on crack size is observed at the initial stage of leakage. As noted previously, the initial leakage mass flow rates were estimated to be 2.175 kg/s for a crack size of 10 mm, which is about nine times lower than a crack size of 30 mm. Thus, the diffusion effect is dominant in the early stage of leakage owing to a lower mass flow rate, showing a lower increasing rate of the dangerous volume, compared to the other two cases. When the crack size increases, the gradient at the initial stage becomes higher than the ratio of the mass flow rate. From the results, it was noted that the gas sensor systems should be installed near the bottom corners of the container for the rapid detection of leaked gas.

4.3. Evolution of Wall Temperature

The leaked CO₂ temperature might drop rapidly to −78.37 °C under atmospheric pressure according to the phase diagram [9]. AH steel, which is a higher tensile steel of grade “A”, is typically used as a material for hull structures [14]. In the case of AH steel in the manufacturing industry, the design temperature was determined as −10 °C when considering the thickness of hull structure during thermal deformation, based on the Charpy V-notch test results. Based on the IGC code, this AH steel was determined as the material for hull construction, without a secondary barrier in container [29]. Figure 10 shows the temperature profile with respect to the vertical distance from the crack to the bottom surface at 300 s. The inlet gas temperature is −78.37 °C, and the gas impinges on the wall with an initial temperature of 20 °C in the CFD simulation. A steep gradient of gas temperature is shown in the vicinity of the crack due to intensive convective heat transfer. The temperature gradient increases with the crack size because of the increased thermal mass, proportional to the leakage flow rate. Variations in temperature affect the air pocket region where conjugated heat transfer occurs through the wall.

According to the regulations regarding ship classification, the design temperature of AH Steel for a hull structure is known to be −10 °C [14]. To avoid thermal damage, it is necessary to evaluate temperature evolution based on the given design temperature. The minimum solid temperature near the impingement surface was estimated to be −2.54 °C, −15.77 °C, and −28.17 °C for crack sizes of 10 mm, 20 mm, and 30 mm, respectively. For crack sizes larger than 20 mm, the impingement surface temperature decreased below −10 °C, which could cause structural damage induced by thermal stress. This result supports the fact that material selection and structural design should be carefully made by considering the threshold temperature to avoid thermal damage.

5. Conclusions

The present study investigated CO₂ gas propagation and thermal characteristics in a tank of a transport ship numerically by using two commercial codes. The local variations of temperature and gas concentration behavior were examined for different crack sizes. Moreover, an effective crack area was introduced to determine the boundary conditions for CFD simulations. The following conclusions were drawn.

1. Considering the phase change, the gas velocity exhausted from the crack was estimated by considering rapid gas expansion. According to the mass conservation principle, an effective crack area was calculated for making appropriate boundary conditions for CFD simulations. The ejected gas flow was rapidly diffused and impinged on solid surfaces. In particular, the CO₂ concentration changed rapidly depending on the crack size. The present study provided the evolution of a dangerous volume with a CO₂ concentration above 10 %, showing an increasing tendency as the crack size increased.

2. Numerical results showed that the leaked CO₂ gas rarely propagated toward the upper tween deck region because of the buoyancy effect. For a crack size of 30 mm, it was predicted that the entire region below the tween deck could become a dangerous volume within 300 s. Moreover, a substantial difference in dangerous volume with differing crack sizes was observed at the initial stage of leakage. The present CFD result supported the fact that, for quick detection before propagation, the gas sensors should be mounted near the bottom corners of the container of a transport ship.

3. The temperature profile from the crack to the air pocket exhibited a decreasing tendency as the crack size increased. The solid temperature of the impingement surface decreased rapidly below the design temperature when the crack size was larger than 20 mm. However, the sudden change in the wall temperature was delayed owing to the conjugate heat transfer from the air pocket.

Description (unit)	Detail	Note
Crack size (mm)	10/20/30
Mass flow rate (kg/s)	2.175/8.698/19.569	estimation by PHAST
Effective crack size (mm)	66/132/198
Density (kg/m3)	air@20°C = 1.225 CO2@20°C = 1.789	ideal gas law
Heat capacity (J/kg∙K)	air@20°C = 1006.43 CO2@20°C = 847.07	piecewise-polynomial profile
Thermal conductivity (W/m∙K)	air@20°C = 0.0242 CO2@20°C = 0.0145 steel@20°C = 16.27	piecewise-polynomial profile
Inlet conditions	TCO2 = −78.37 °C YCO2 = 1	saturation temperature of CO2 at 1 bar [9]
Initial environmental conditions	Tair = 20 °C Yair = 1	room temperature
Boundary conditions	Tamb = 20 °C	ambient temperature
Design temperature	Tdesign = −10 °C	damage to AH steel [14]
CO2 concentration limit	YCO2 = 0.1	injury to humans
Flow time (s)	300	sufficient time to analyze the gas propagation

Author Contributions

Conceptualization, S.H.L. and M.-S.S.; investigation, D.Y.K., C.H.J., B.J.P., and M.S.K.; methodology, D.Y.K. and C.H.J.; supervision, S.H.L.; writing of the original draft, D.Y.K.

Acknowledgments

This research was supported by the development of human resources for design specialists in applications of functional materials (Brain Korea (BK) 21).

Conflicts of Interest

The authors declare no conflict of interest.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).