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

: This study investigates numerically gaseous CO 2 leakage characteristics inside the containers of a transport ship and examines thermal e ﬀ ects on the structural damage that might happen in the containers. First, with consideration of the phase change, the ejected mass ﬂow rate was estimated using the commercial code of DNV PHAST. Based on this estimated mass ﬂow rate, we introduced an e ﬀ ective area model for accounting for the fast evaporation of liqueﬁed CO 2 occurring in the vicinity of a crack hole. Using this leakage modeling, along with a concept of the e ﬀ ective area, the computational ﬂuid dynamics (CFD) simulations for analyzing transient three-dimensional characteristics of gas propagation in a conﬁned space with nine containers, as well as the thermal e ﬀ ect on the walls on which the leaking gas impinges, were conducted. The commercial code, ANSYS FLUENT V. 17.0, was used for all CFD simulations. It was found that there are substantial changes in the pressure and temperature of the gas mixture for di ﬀ erent crack sizes. The CO 2 concentration at human nasal height, a measure of clear height for safety, was also estimated to be higher than the safety threshold of 10% within 200 s. Moreover, very cold gas created by the evaporation of liqueﬁed CO 2 can cool the cargo walls rapidly, which might cause thermal damage.


Introduction
Carbon dioxide (CO 2 ) 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 2 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 2 due to the large capacity of deep saline aquifers [3]. Although pipelines have recently been developed for transporting CO 2 , there are major obstacles related to the initial installation, as well as maintenance costs [4]. Alternatively, CO 2 shipping vessels (or water carriers) are more attractive alternatives when considering maintenance costs [5]. To be specific, the capacity of a CO 2 transport ship is determined by the daily production of CO 2 and the sailing period of the ship [6]. Much effort to build CO 2 transport ships has been conducted in the Republic of Korea, because Korea is a peninsula surrounded by the sea [7]. Currently, CO 2 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 2 , 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 2 leak occurs in the container during sailing [11]. When an accidental leak occurs, liquefied CO 2 undergoes rapid vaporization and expansion, causing a significant pressure rise [12]. Even if CO 2 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 2 concentration [13]. The other is associated with thermal damage to the tank structure during leakage, because very cold CO 2 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 2 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 2 gas. Thus, to provide a better understanding when designing CO 2 cargo facilities, the present study aims to numerically investigate CO 2 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 2 concentration and structure temperature is analyzed.

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 2 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 2 , 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 2 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 6 and We = 3.99 × 10 5 [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 0 2 /(q 0 ∆T), where D 0 is the initial drop diameter [22].
Here, ∆T is the temperature difference between the ambient temperature and the drop temperature, and q 0 is an empirical constant [22]. Finally, the lifetime of CO 2 droplets after the break-up was estimated to be 2.40 × 10 −26 s. With the droplet velocity and lifetime, a distance required for phase change from liquid to gas was estimated to be 5.7 × 10 −24 m from the crack, which means that the liquefied CO 2 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 2 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.
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 CO2. Here, it is assumed that the liquefied CO2 exposed from the crack evaporates spontaneously. Gas phase velocity, vgas, 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 onedimensional 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 CO2 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 CO2 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. 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 2 . Here, it is assumed that the liquefied CO 2 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 2 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 2 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.

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.
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: where the turbulent viscosity µ t can be obtained by Appl. Sci. 2019, 9,2536 5 of 12 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 This simulation used a target CO 2 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 2 ship. It provides useful information about the frequency index, severity index, consequence, location, and hazard event [24]. where the turbulent viscosity μt can be obtained by where Gk represents the generation of turbulent kinetic energy due to mean velocity gradients. Gb is the generation of turbulent kinetic energy due to buoyancy. YM 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 1 2 1.44, 1.92, 0.09, This simulation used a target CO2 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 CO2 ship. It provides useful information about the frequency index, severity index, consequence, location, and hazard event [24].  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 CO2 gas was leaking at −78.37 °C, which is the 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 2 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.
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 pressurevelocity 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 −6 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.

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 CO2 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. 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 −6 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.

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 2 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 CO2 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 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.  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 CO2 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 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 2 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 Appl. Sci. 2019, 9, 2536 8 of 12 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.
Appl. Sci. 2018, 8, x FOR PEER REVIEW 8 of 13 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 nonhomogeneous 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.

Accumulation of CO2 Concentration
Since humans can lose consciousness within a few minutes when exposed to CO2 gas at a concentration of over 10% [27,28], a CO2 concentration of 10% can be regarded as a safety threshold. The predicted CO2 concentration distribution passing through the crack is provided in Figure 7. The results show that the maximum CO2 concentration occurs at the corner region near the crack location. A much higher CO2 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 CO2 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.

Accumulation of CO 2 Concentration
Since humans can lose consciousness within a few minutes when exposed to CO 2 gas at a concentration of over 10% [27,28], a CO 2 concentration of 10% can be regarded as a safety threshold. The predicted CO 2 concentration distribution passing through the crack is provided in Figure 7. The results show that the maximum CO 2 concentration occurs at the corner region near the crack location. A much higher CO 2 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 2 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 2 concentration being over 10%, to characterize three-dimensional leakage behavior. Physically, the leaked CO 2 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. From an engineering viewpoint, as shown in Figure 8, it would be useful to define a dangerous volume of CO2 concentration being over 10%, to characterize three-dimensional leakage behavior. Physically, the leaked CO2 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.

Evolution of Wall Temperature
The leaked CO 2 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.

Evolution of wall temperature
The leaked CO2 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 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.

Conclusions
The present study investigated CO 2 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 2 concentration changed rapidly depending on the crack size. The present study provided the evolution of a dangerous volume with a CO 2 concentration above 10 %, showing an increasing tendency as the crack size increased.

2.
Numerical results showed that the leaked CO 2 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.