**2. Materials and Methods**

#### *2.1. Numerical Model*

The homogeneous Eulerian multiphase model [24,25] was adopted to model the phase change process after LNG leaked to the ground. The realizable *k*-*ε* model had higher accuracy in concentration distribution than the standard *k*-*ε* model by simulating Thorney's heavy gas diffusion (Freon-12) field test [26]. Therefore, the realizable *k*-*ε* model was selected for gas diffusion turbulence.

At present, the *k*-*ε* model is the most widely used turbulence model for the turbulence simulations of wind fields in large structures such as storage tanks. The standard *k*-*ε* model proposed by Launder and Spalding greatly improves the zero-equation model and one equation model, so it is widely used in engineering flow field calculation and has been well verified in practice. However, the applicability of the standard *k*-*ε* model for each component of Reynolds stress is not strong. For example, it is assumed that the turbulent viscosity coefficient is isotropic, while the turbulence is anisotropic in the case of curved wall flow, curved streamline flow, or strong swirling flow. Therefore, it is not recommended to use the standard *k*-*ε* model to calculate the wind field of a storage tank with curved wall flow; otherwise, it will produce a certain degree of distortion in calculation.

The equation of turbulent kinetic energy *k* and dissipation rate *ε* of the standard *k*-*ε* model is described as follows.

$$\rho \frac{Dk}{Dt} = \frac{\partial}{\partial \mathbf{x}\_i} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_k} \right) \frac{\partial k}{\partial \mathbf{x}\_i} \right] + \mathbf{G}\_k + \mathbf{G}\_b - \rho \varepsilon - \mathbf{Y}\_M \tag{1}$$

$$
\rho \frac{D\varepsilon}{Dt} = \frac{\partial}{\partial \mathbf{x}\_i} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_k} \right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_i} \right] + \mathbf{C}\_{1\varepsilon} \frac{\varepsilon}{k} (\mathbf{G}\_k + \mathbf{C}\_{3\varepsilon} \mathbf{G}\_b) - \mathbf{C}\_{2\varepsilon} \rho \frac{\varepsilon^2}{k} \tag{2}
$$

In the above equation, *Gk* represents the turbulent kinetic energy generated due to the average velocity gradient; *Gb* refers to the turbulent kinetic energy caused by buoyancy; *YM* refers to the effect of compressible turbulent pulsating expansion on the total dissipation rate.

The coefficient of turbulence viscosity is:

$$
\mu\_t = \rho \mathbb{C}\_{\mu} \frac{k^2}{\varepsilon}
$$

In FLUENT, these three parameters *C*1*ε*, *C*2*ε*, *C<sup>μ</sup>* are the default constant, *C*1*<sup>ε</sup>* = 1.4, *C*2*<sup>ε</sup>* = 1.92, and *C<sup>μ</sup>* = 0.09. The turbulent Prandt numbers of turbulent kinetic energy *k* and dissipation rate ε are respectively 1.0 and 1.3.

The transport equation of turbulent kinetic energy *k* and dissipation rate ε is described as follows.

The equation *k*:

$$
\rho \frac{Dk}{Dt} = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_k} \right) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + \mathbf{G}\_k + \mathbf{G}\_b - \rho \varepsilon - \mathbf{Y}\_M \tag{3}
$$

The equation *ε*:

$$\rho \frac{D\varepsilon}{Dt} = \frac{\partial}{\partial x\_j} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_t} \right) \frac{\partial \varepsilon}{\partial x\_j} \right] + \rho C\_1 S \varepsilon - \rho C\_2 \frac{\varepsilon^2}{k + \sqrt{\upsilon \varepsilon}} + C\_{1\varepsilon} \frac{\varepsilon}{k} C\_{3\varepsilon} G\_b \tag{4}$$

where,

$$C\_1 = \max\left[0.43, \frac{\eta}{\eta + 5}\right] \\ \eta = Sk/\varepsilon = \left(2E\_{i\bar{j}} \cdot E\_{i\bar{j}}\right)^{1/2} \\ \frac{k}{\varepsilon}E\_{i\bar{j}} = \frac{1}{2} \left(\frac{\partial u\_i}{\partial x\_{\bar{j}}} + \frac{\partial u\_{\bar{j}}}{\partial x\_{\bar{i}}}\right)$$

In the above equation, *k* is the kinetic energy of turbulence pulsation; ε is the dissipation rate of the turbulent pulsation kinetic energy; *Gk* is the turbulent kinetic energy caused by the average velocity gradient; *Gb* is the turbulent kinetic energy caused by buoyancy; *YM* is the effect of compressible turbulent pulsating expansion on the total dissipation rate. *C*<sup>2</sup> and *C*1*<sup>ε</sup>* are constant; *σ<sup>k</sup>* and *σε* is turbulent Prandt numbers of turbulent kinetic energy and dissipation rate, respectively.

In Fluent,

 $\mathbb{C}\_{1\varepsilon} = 1.44$ ,  $\mathbb{C}\_{2\varepsilon} = 1.9$ ,  $\mathbb{C}\_{3\varepsilon} = 0.09$ ,  $\mathbb{C}\_{2} = 1.9$ ,  $\sigma\_{k} = 1.0$ ,  $\sigma\_{\varepsilon} = 1.2$ .

#### *2.2. Parameter Setting*

A 16 × 104 m<sup>3</sup> large cylindrical LNG storage tank was chosen for numerical simulation, and its structural dimension are shown in Figure 1. The outer diameter of the tank was 82 m and the height was 50 m. The normal operating pressure of the storage tank was 25 kPa, and the maximum liquid level in the tank was 34.6 m. The origin of the computational domain was located at the center of the bottom of the tank. The coordinate of the leakage hole center point was (41, 10, 0), which was located on the leeward side of the tank. The leakage hole sizes were, respectively, 0.1 m × 0.1 m, 0.13 m × 0.13 m, 0.15 m × 0.15 m, 0.18 m × 0.18 m and 0.2 m × 0.2 m. Considering the calculation accuracy, the computational domain was determined to be 1000 m × 250 m × 500 m in the x, y, and z directions, and the tank that has a blocking rate of 2.78% was placed at a distance of 200 m downwind. The whole computational domain was discretized by the structured grid, and the specific grid division is shown in Figure 2a. In order to adapt to the change of flow field and ensure the accuracy of the solution, the grid around the leakage hole was encrypted by the block method. The independence of the grid and time step had been verified. The total number of cells in the calculation domain was finally determined to be 1,865,345, and the simulation time step was set to 0.1 s.

**Figure 1.** The settings of the large-scale LNG storage tank. (**a**) Geometric schematic of the tank; (**b**) the boundary settings of the tank.

**Figure 2.** Mesh division in this study. (**a**) The meshing of the computational watershed; (**b**) the meshing of LNG leakage diffusion experiment.

In order to represent the node coordinates more accurately and ensure the convergence of calculation, a double-precision solver and implicit method were used in the calculation. Figure 2b shows the meshing of the LNG leakage diffusion experiment. The calculation domain was established with a size of 900 m × 500 m × 50 m on the *x*-axis, *y*-axis, and *z*-axis, respectively. The *x*–*z* plane was placed on the ground, and the *y*-direction was the vertical height. Furthermore, the wind direction remained unchanged throughout the calculation domain. The boundary conditions on the left and right sides of the calculation

domain were the velocity–inlet and the pressure–outlet, respectively. Hexahedral mesh units were used for mesh generation, while the area around the pond was divided into fine meshes. A total of 803,287 cells were used for subsequent simulations.
