Numerical Simulation of Propagation Characteristics of Hazardous Noxious Substances Spilled from Transport Ships

Abstract

This study numerically investigates the propagation characteristics of hazardous noxious substances (HNSs) spilled from transport ships and suggests the metal model for predicting the HNS propagation velocity varied with the current velocity and HNS density. The commercial computational fluid dynamics (CFD) code ANSYS FLUENT (V. 17.2) was used for two-dimensional simulation based on the Reynolds-averaged Navier–Stokes (RANS) equation together with the standard k–ε model. The scalar transport equation was also solved to estimate the spatial and transient behaviors of HNS. The main parameters to analyze the near-field propagation characteristics of HNSs spilled from the ship were layer thickness, HNS concentration, and propagation velocity. It was found that advection becomes more dominant in propagating an HNS layer that becomes thinner as the current velocity increases. When the current velocity increased beyond a certain level (~0.75 m/s), the mixing effect made the HNS layer less dense but thicker. Consequently, lower-density HNS causes increased HNS concentrations at sea level. As the current velocity increased, the concentration distribution became homogeneous regardless of HNS density. In particular, the second-order response surface model provided for three variables on the basis of the numerical results for 15 cases with the use of the general least-squares regression method, showing a good fit. This model would be useful in estimating the propagation velocity of HNS spilled from a ship.

1. Introduction

Hazardous and noxious substances (HNSs) are widely known as liquid and mixed-liquid substances that can have harmful effects on the marine environment [1]. Large quantities of HNSs are transported by sea. If an HNS-related accident occurs in a port or on a ship, harmful substances are diffused in the vicinity of the ships in different ways depending on the type of spill material (e.g., sinking, floating, evaporating, dissolving) [2]. Indeed, it is difficult to cope with HNS accidents because of the high risk associated with various HNSs with bioaccumulative, biodegradable, toxic, and explosive characteristics [3,4,5,6,7,8]. To date, many government authorities and environmental organizations have concentrated their effort on establishing emergency planning against HNS disasters. There is still a lack of information on HNS characteristics under HNS accident scenarios.

Indeed, it is difficult to perform experiments directly on HNS propagation in real circumstances because of the toxicity of HNSs with various chemical properties. Nevertheless, Fuhrer et al. [2] dissolved a small amount of styrene, an HNS, in the sea and monitored the residual concentration over time. They reported that the spilled styrene was above 50% of the total HNS amount only for 1 h in the early stage of the accident. This report can be useful in making subtle plans for accident response. Because they conducted the experiment for directly releasing HNS in the actual sea, however, they could cause severe safety problems and destroy the marine ecosystem in the experiment. Considering the safety and environmental issues, direct experiments need to be strongly restricted. This is why the available experimental data related to HNS propagation is hardly available.

The computational fluid dynamics (CFD) approach has been a promising way to investigate the diffusion characteristics of HNSs near ships in the sea [9,10,11,12]. Elhakeem et al. [9] developed a three-dimensional prediction model involving the hydrodynamic model and the oil-spill model. They compared the real photos taken for oil-spill propagation in the Arabian Gulf region with the oil-spill model. Their model also made an accurate prediction for a credible analysis of oil spills. Berry et al. [10] developed the OILTRANS model that consists of two inter-linked models, namely, the particle-transport model, based on the Lagrangian equation, and the oil-fates model. This model was applied to a spill accident that had occurred during a ship-to-ship fuel transfer in the Celtic Sea in 2009. The model was compared with aerial observations of the movement of the oil layer and showed good agreement in estimating the oil trajectory within ten days. Zhong et al. [11] focused mainly on the validation of their own three-dimensional and high-resolution oil-spill model. They constructed a high-resolution grid with a size ranging from 440 m × 440 m to 500 m × 500 m. They simulated the trajectory of the spilled oil and compared the results with the historical information from the motor vessel (M/V) Marathassa oil spill. Janeiro et al. [12] also suggested a prediction model of oil-spill trajectory, which was coupled with hydrodynamics, oil transport, and a weathering model. They focused mainly on the validation of the developed model by comparison with the observed data from oil spills in the Tuscany Archipelago region.

Most of the previous works mentioned above mainly reported oil-spill accidents, not general HNS accidents [13]. Their simulation results were obtained in domains with tens to hundreds of kilometers. This approach is called the far-field analysis, by which the local propagation of a substance spilled from the ship cannot be predicted. The main purpose of this far-field analysis is to show the bulk motion of a spilled substance based on acquired data of winds and sea currents. This approach fails to intrinsically examine the near-field propagation characteristics in the vicinity of the transport ship because it treats such a finite region near the ship as a point source, numerically. When establishing a safety strategy, it is essential to understand the HNS spill and propagation characteristics within one kilometer or less, especially in the early stages of the incident. Thus, the present study attempts to obtain useful information on HNS propagation near the spill location by using CFD, which is called “near-field analysis” in this study.

The present study provides the two-dimensional CFD simulation results and makes a meta-model platform on the basis of the numerical predictions with the use of the response surface-regression method. The current simulation was conducted for a relatively high current velocity ranging from 0.5 to 1.5 m/s that was taken among various scenarios. This is because the highest risk scenario is generally considered in handling HNS accidents in the sea. Under this limitation, the two-dimensional analysis seems to be suitable for analyzing the propagation characteristics because of the very high Peclet number under which the convection effect becomes dominant. In fact, there may be many environmental factors in the case of a real HNS accident. Among them, the present study focuses on the influence of the current velocity on HNS propagation behavior using the commercial program, ANSYS FLUENT (V. 17.2). The HNS layer thickness, averaged HNS concentration at sea level, and average propagation velocity were predicted for different current velocities and HNS densities.

2. Numerical Methods

2.1. Governing Equations and Details for Simulation

The present study solves a set of conservation equations for mass, momentum, and species that are expressed as follows.

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u_i)}{\partial x_i}=0,$$

(1)

$$\frac{\partial }{\partial t}(\rho u_i)+\rho u_j\frac{\partial u_i}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ji}}{\partial x_j}+\rho g,$$

(2)

$$\frac{\partial }{\partial t}(\rho Y_m)+\frac{\partial (\rho u_i{Y}_m)}{\partial x_i}=-\frac{\partial J_{jm}}{\partial x_j},$$

(3)

where ρ, u, x, and t are the density, velocity, position, and time, respectively. p is the pressure and τ_ji is the stress, J_jm is the diffusion flux of species m, and Y_m is the mass fraction of species m. The standard k–ε turbulence model was used, and the standard wall function was employed to describe the turbulent dissipation energy near the wall. The equations for the turbulent kinetic energy k and the rate of energy dissipation ε are as follows.

$$\frac{\partial }{\partial t}(\rho k)+\frac{\partial }{\partial x_i}\left(\rho u_{i}k\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho \epsilon ,$$

(4)

$$\frac{\partial }{\partial t}(\rho \epsilon )+\frac{\partial }{\partial x_i}\left(\rho u_i\epsilon \right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_1\frac{\epsilon }{k}(G_k+G_b)-C_2\rho \frac{{\epsilon }^2}{k},$$

(5)

where G_k is the production of turbulent kinetic energy due to mean velocity gradients and G_b is the production of turbulence kinetic energy due to buoyancy. C₁ and C₂ are model constants and σ_k and σ_ε are the turbulent Prandtl numbers of k and ε, respectively. The model constants are as follows.

$$\begin{array}{l}{C}_1=1.44,C_2=1.92,C_{\mu }=0.09\\ {\sigma }_k=1.0,{\sigma }_{\epsilon }=1.3\end{array}$$

(6)

In addition, μ_t denotes the turbulent viscosity obtained by

$${\mu }_t=C_{\mu }\rho \frac{k^2}{\epsilon }.$$

(7)

Figure 1 shows the schematic of the two-dimensional computational domain of the HNS accident ship. The shape of the ship was provided by Korea Research Institute of Ship and Ocean Engineering (KRISO). It is assumed that the ship moves parallel to the direction of the current. The current scenario considers the specific crack position and leakage flow rate, as illustrated in Figure 1, and assumes that the spilled HNS is mainly affected by the sea current. In spite of the two-dimensional analysis, the present study is expected to provide a useful analysis of HNS propagation affected by the current velocity and substance density. The domain size was 1000 m in length and 100 m in depth (y). The length of the accident ship (L_ship) was 80 m and its height (H_draft) was 8 m. The crack area was 0.1 m², located 30 m from the top of the ship. The number of grids in the computation domain was taken as 355,496, determined from grid-independent tests. In general, the sea current consists of three different components: tide-induced, density-driven, and wind-driven currents. Among them, the density-driven current and tide-induced current have almost uniform distributions in the depth direction [14]. Furthermore, the wind-driven current velocity is smaller than the other two. Thus, the ocean current is assumed to be a simplified uniform velocity profile for calculation.

The sea-current velocity was varied in the range from 0. 5 to 1.5 m/s for the present simulation. The Neumann condition was applied to the outlet boundary surface and the slip condition was used for the top surface of the computational domain. The floating HNS refers to a substance whose density is less than seawater density (1024 kg/m³), and the HNS density was conventionally in the range from 700 to 900 kg/m³ [2]. The present simulation varied the HNS density in this range, resulting in floating behavior. The viscosity and diffusion coefficient of the HNS were set to 0.00076 kg/m·s and 8 × 10⁻⁸ m²/s, based on styrene. It was assumed that the HNS spilled with a mass flow rate of 45 kg/s·m², estimated from a possible accident scenario. A first-order upwind scheme was used, and the density of the mixture was calculated by using the volume-weighted mixing law, as follows.

$$\rho =\frac{1}{{\displaystyle \sum _i^k\frac{Y_i}{{\rho }_i}}}$$

(8)

The time step was fixed as 0.25 s, and its averaged Courant number during the calculation was 0.183. The total flow time was taken as 1000 s for all calculations.

2.2. Response Surface Model

The meta-model can be suggested for providing a surrogate mathematical model of the original CFD simulation [15,16,17]. The present study used the response surface model (RSM) [17,18,19]. First, the surface function of interest was set as

$$y\left(x\right)=f\left(x\right)+e,$$

(9)

where y(x) is the surface function of interest, f(x) is a polynomial function of x, and e is the error or noise observed in y(x). The polynomial function f(x) used to approximate y(x) is assumed to be of typical quadratic form, expressed by

$$y\prime ={\beta }_0+{\displaystyle \sum _{i=1}^k{\beta }_i{x}_i}+{\displaystyle \sum _{i=1}^k{\beta }_{ii}{x}_i^2}+{\displaystyle \sum _{i<j}^{}{\beta }_{ij}{x}_i^{}{x}_j+e}.$$

(10)

The regression coefficients, β₀, β_i, β_ii, and β_ij, were determined by the least-squares method that minimizes the sum of the squares of the residual between $y\prime$ and y(x). The regression coefficients can be determined by

$$\beta ={\left[X^{\prime }X\right]}^{-1}{X}^{\prime }y,$$

(11)

where y is an n × 1 column vector of the response values, X is a k × n matrix of the sample points, and X′ is its transpose.

2.3. Validation of CFD Simulation

Unfortunately, it was not possible to find any experimental or numerical results for the above-mentioned scenario. As an alternative, additional CFD simulations were conducted to evaluate the flow-rate estimation made in this study by comparing the results from a previous report [20]. This is because the flow rate model becomes an important boundary condition that strongly affects the HNS propagation. Indeed, Tavakoli et al. [20] conducted a lab-scale experiment for the oil-spill behavior with the damaged tank below the waterline. According to Reference [20], the present study constructed the computational domain, as shown in Figure 2. The size of the room was 12 m × 5 m × 5 m and the depth of water was 3 m. An oil tank with a size of 1 m × 0.5 m × 1 was placed at the center of the water bath. The draft of the tank refers to the depth of the submerged part of the tank, and its value was taken as 0.6 m in this study. A puncture with a diameter of 2.2 cm was drilled 0.1 m from the bottom of the oil tank. The number of grids in the computational domain was taken to be approximately 370,000, determined from grid-independent tests. The oil tank was fully filled with oil, initially. The top surface was set as an atmospheric-pressure outlet condition.

The oil in the tank began to flow out through the puncture initially, and the volume flow rate was estimated by the change in the oil head in the tank. After the spill began, the oil level and volume flow rate decreased owing to the change in pressure difference. The predicted volume flow rate of spilled oil was compared with the experimental result, as shown in Figure 3. It shows reasonable agreement between the prediction and measurement, showing a difference of less than 8%.

3. Results and Discussion

3.1. HNS Propagation Characteristics

Figure 4 shows the distributions of HNS mass fractions at t = 1000 s after leakage. After 1000 s, the flow field reached the steady state, and the shape of the HNS layer was almost constant over time. The mass fraction of HNS in the wake region was measured in the range from 0.02 to 0.1. The dissolved HNS floated near the sea surface because it had a lower density than seawater. Here, the HNS layer formed behind the wake was mainly analyzed. The maximum mass fraction appeared at sea level, and it showed decreasing patterns in the depth direction. For quantitative analysis, three main parameters, namely, the HNS layer thickness (δ_layer), HNS propagation velocity (v_p), and averaged HNS mass fraction at sea level ($\varphi$_s−l) were introduced in the present study.

Figure 5a represents the change in the HNS layer thickness for different current velocities and HNS densities at t = 1000 s. The HNS layer thickness is defined as the depth at which the integrated HNS mass reaches 90% of the total amount of spilled HNS by integrating the mass fraction of HNS ($\varphi$_hns) from sea level at t = 1000 s, as follows.

$${\frac{{\displaystyle {\int }_0^{{\delta }_{layer}}{\displaystyle {\int }_0^{x_{\max }}\varphi \left(x,y\right)\rho dxdy}}}{{\displaystyle {\int }_0^{y_{\max }}{\displaystyle {\int }_0^{x_{\max }}\varphi \left(x,y\right)\rho dxdy}}}|}_{t=1000}=0.9,$$

(12)

where x_max and y_max represent the length of the computational domain in the x and y directions, respectively. The denominator represents the total amount of spilled HNS. The HNS layer thickness indicates the upper limit for the y-axis of the integral in the numerator, as shown in Figure 5b. From the results, the lower the HNS density, the thinner the layer that was formed under the same current velocity conditions, because of the buoyancy effect. When the density of HNS decreases, the dissolved HNS floats rapidly to the sea surface and forms a thinner layer. The change in HNS layer thickness is presented in Figure 5a. The HNS layer thickness decreases with the increase in current velocity, but it increases after a certain value of current velocity because of the competitive role between advection and diffusion. The inflection point appears at approximately v_c = 0.75 m/s. When v_c is 0.5 m/s, the corresponding δ_layer values are 4.56, 4.21, and 4.01 m for ρ_hns values of 700, 800, and 900 kg/m³, respectively. As v_c increases, the advection becomes dominant to overcome the diffusion effect and acts as the main driving force in HNS propagation. However, when v_c increases above a certain value (~0.75 m/s) at which an inflection point appears, the mixing of HNS becomes significant and makes the HNS layer less dense but thicker.

Figure 6a shows the transient evolution of the HNS interface for different current velocities in the case that ρ_hns = 700 kg/m³. The HNS interface is defined as the location at which the HNS mass fraction is 1.0 × 10⁻⁴. Knowledge of the maximum propagation length, defined as the distance from the crack position to the interface in the mainstream direction, is important in predicting the damage region that becomes wider with time. In the case of v_c = 1.5 m/s, the interface moves in the down-stream direction by approximately 40 m every 25 s. The second main parameter, HNS propagation velocity (v_p), is defined as the time-averaged velocity of the HNS interface in the downstream direction:

$$v_p=\frac{1}{T}{\displaystyle \underset{t=0}{\overset{T}{\int }}{v}_x\left(t\right)dt},$$

(13)

where t and v_x are the time and velocity of the HNS interface in the downstream direction. T denotes the period of integration, which is adopted as the whole computation time in the present study. Figure 6b shows the maximum propagation length of the HNS interface in the downstream direction. The maximum propagation length linearly increases with time and is significantly affected by the current velocity. As shown in Figure 6c, however, very similar tendencies are found for different densities because the advection and diffusion effects become more dominant than the buoyancy effect. The propagation velocity v_p is slightly greater than v_c. For example, in the case of v_c = 1.5 m/s, v_p is approximately 1.6 m/s. In addition, v_p varies with HNS density. When v_c is 0.5 m/s, v_p is 0.56, 0.60, and 0.65 m/s with respect to the HNS density of 700, 800, and 900 kg/m³, respectively. The reason is as follows. The dissolved HNS floats with vertical-directional momentum because of the buoyancy-driven force. The vertical-directional momentum of dissolved HNS will change to the streamwise-directional momentum at the sea surface region. Therefore, floating HNS with a density of 700 kg/m³ propagates faster than that with a density of 900 kg/m³. This effect is valid even in the case of v_c = 1.5 m/s but is not significant at higher-velocity conditions. This is because the vertical direction momentum component of the dissolved HNS decreases as v_c increases.

Figure 7 represents the change in the averaged HNS mass fraction at sea level for different current velocity and HNS density. The averaged HNS mass fraction at the sea level, $\varphi$_s−l, is defined as the area-averaged mass fraction of the HNS at sea level (y = 0) as follows:

$${\varphi }_{s-l}={\frac{1}{A}{\displaystyle \int {\displaystyle \int \varphi \left(\overrightarrow{x}\right)dA}}|}_{t=1000s}$$

(14)

where $\varphi$(x) represents the HNS mass fraction at sea-level and A is the area at the sea surface. As v_c increases, $\varphi$_s−l decreases proportionally, since spilled HNS is not stagnant but is actively propagated by advection. From 1 to 1.5 m/s v_c, the averaged HNS mass fraction at sea level increased for less dense HNS because buoyancy is a significant driving force for floating HNS at sea level. However, when v_c is 0.5 m/s, $\varphi$_s−l is estimated to be about 0.021 regardless of HNS density. This is because the spilled HNS with a density of 700 kg/m³ propagates faster than that with a density 900 kg/m³ due to the buoyancy force, and it makes the HNS layer less dense. Therefore, the $\varphi$_s−l values become similar regardless of HNS density at a v_c of 0.5 m/s.

3.2. Response Surface Model for HNS Spill Predictions

The CFD simulation provides deterministic solutions that are bounded by the boundary and operating conditions, which vary with different scenarios. Thus, obvious limitations exist in taking the CFD approach whenever accident scenarios are changed. To reduce the computational cost, the statistical approach has been widely used in many engineering applications and is also used in the present study. As in the meta-modeling, the present study adopted the RSM to provide the surrogate mathematical model from the CFD simulation results. To our knowledge, this is the first study of the near-field propagation characteristics. Of course, it is not easy to create a meta-model with high accuracy because there are insufficient simulation results. However, a more accurate model can be made by accumulating the simulation data and updating the model continuously. The main objective of this work is, therefore, to construct the numerical procedure required for making the meta-model of HNS propagation. The three main parameters, HNS layer thickness (δ_layer), HNS propagation velocity (v_p), and averaged HNS mass fraction at sea level ($\varphi$_s−l) are the dependent variables used in the mathematical model. The second-order RSM for the three main parameters is fitted to 15 cases using a general least-squares regression method. Each input/dependent variable is obtained from the CFD simulation, as discussed above. The general form of second-order RSM is expressed by

$$y\prime =a_o+a_1{x}_1+a_2{x}_2+a_3{x}_1{}^2+a_4{x}_2{}^2+a_5{x}_1{x}_2,$$

(15)

where, y′ is the main parameter, x₁ = v_c, and x₂ = ρ_hns; a₀ to a₅ are model constants. The obtained model constants are listed in

Table 1. Model constant of second-order response surface model (RSM).

Run No.	Model Constants
	a0	a1	a2	a3	a4	a5
δlayer	3.6602	−1.9772	−4.5864	0.5883	2.3693	0.9463
vp	0.9969	0.3154	−1.1469	0.1170	0.3343	0.3788
$\varphi$s−l	0.1889	−1.0484	4.0449	0.4659	−1.8737	−0.7809

. Three-dimensional response surface contours of each dependent variable are shown in Figure 8. The predicted values are compared with the corresponding CFD simulation results (black dots), and the goodness of fit is measured. The variables and scatter of the data points around the fitted regression curve are always evaluated between 0 and 1. As listed in

Table 2. Error testing of response surface models.

Measure	δlayer	vp	$\mathit{\varphi }$s−l
$R_{}^2$	0.9446	0.9988	0.9938
$R_{adj}^2$	0.9354	0.9986	0.9927

, both values are greater than 0.9, which means that the prediction model provides a good fit.

4. Conclusions

The present study conducted a numerical simulation of an HNS spill and developed a prediction model. The HNS spill and propagation behaviors were evaluated by numerical analysis according to the water-current velocity conditions and HNS density. The following three main parameters were defined: HNS layer thickness, HNS propagation velocity, and averaged HNS mass fraction at sea level. In addition, corresponding mathematical prediction models were developed using the response surface-model method. The results are as follows, and such data would be useful in HNS prediction technology and understanding the detailed physics behind the near-field characteristics of HNS spills and propagation.

(1)

Lower-density HNS forms thinner layers in all water-current velocity ranges. As current velocity increases, the advection effect becomes a significant driving force of HNS propagation. In this case, the spilled HNS was not stagnant, but was actively propagated by advection, and the layer thickness decreased. However, when the current velocity increased above a certain level (~0.75 m/s), the mixing effect became significant and made the HNS layer less dense, but the layer thickness increased.

(2)

The current velocity strongly affected the propagation velocity, which linearly increased in proportion to current velocity. The averaged HNS mass fraction at sea level increased for less dense HNS because buoyancy was a significant driving force for floating HNS at sea level. Moreover, the HNS mass fraction decreased as the current velocity increased.

(3)

Statistical approximations were used in the present study to construct a metamodel, providing a surrogate mathematical model of the original CFD simulation. The mathematical model provides meaningful results, in seconds, for the numerical conditions that have not been calculated.