**1. Introduction**

An efficient trickle flow of high-temperature melts in coke beds is necessary in ironmaking blast furnaces, as it ensures a smooth continuous process, which is a prerequisite for high productivity in operations that use lower amounts of reducing agents. The trickle flow is driven by the gravitational force that acts on the liquid iron and molten slag (two immiscible liquids), and the dripping behavior is balanced by viscous and interfacial forces. When the volume of the high-strength coke, which acts as a spacer in the lower part of the blast furnace, is reduced, the resistance to the melt passage increases, which can cause instability. It is, therefore, essential to maintain a desired trickle flow through the coke bed to prevent clogging problems that may occur if the pressure drop in the bed becomes excessive or shows large variations. This approach emphasizes the need to understand the dripping behavior of the two liquids in the coke bed.

Over the years, liquid flow in blast furnaces has been investigated experimentally and mathematically. Several studies focused on the volume fraction of the liquid in the packed bed, called "liquid hold-up", and the effect of various in-furnace conditions, including the physical properties of the liquid and the wettability of coke. The relationship between hold-up and dimensionless number has been determined by applying an experimental system using room temperature media [1–3] or actual melt [4–6]. Since the hold-up of CaO-SiO2-Al2O3-based molten slag (>1773 K) is significantly affected by the wettability of coke [7], several experimental studies have been conducted on the static contact angle between molten slag and "carbonaceous material" [8–10]. The molten slag wetting behavior relatively depends on the carbonaceous material substrate, i.e., not only the equilibrium contact angles but the periods to reach a stable contact angle vary considerably. If basicity (CaO/SiO2 mass %) ≤ 1, the decrease in the contact angle with time as the basicity decreases is generally attributed to the formation of interface product SiC [11], but the change in the state of the Fe-Si-C system is still unclear because this change depends on the Si form (which may exist in a gaseous form). In investigating coke-slag wettability, the differences in small amounts of components (such as FeO and MgO) in slag and the differences in the specific properties (such as crystal structure and void structure) of carbonaceous substrates should be considered.

Computational fluid dynamics (CFD) is useful for understanding this phenomenon. The advantage of this approach is the precise tracking of time-varying fluid interfaces, which is difficult to observe experimentally. Most recent CFD studies demonstrate that, with regard to the wettability of the packed bed, the geometry of the bed has a significant effect on the interface dynamics, leading to the co-existence of complex interfaces; thus, the change in wettability leads to various interfacial movements that are unidentified under static conditions. Molten slag flow in a packed coke bed was investigated by using the three-dimensional combined discrete element method (DEM) and CFD [12–14]. A series of two-dimensional [15,16] and three-dimensional [17] high-resolution direct numerical simulations were also conducted with the aim of understanding the pore-scale fluid dynamics. Although the two-liquid flow in packed beds has received attention but usually focused on room temperature media [18,19], there is no guarantee that the information available from the water-oil-rock system can be applied to the slag-iron-coke system in blast furnaces. For example, the density difference between iron and slag (≈4000 kg/m3) is considerably higher than that of water and oil (≈200 kg/m3) [20]. It is essential to evaluate the interfacial force in an unsteady condition under the condition where the momentum in the gravity direction appears. CFD studies have been conducted for such cases. The phase field modeling for the behavior of molten slag confirmed that the increment of the interfacial tension between liquid iron and molten slag was the driving force of the molten slag separation from the melt [21]. Using the volume of fluid simulation, the driving force of the liquid iron penetration into the coke bed was proposed due to the total energy reduction by extending the area covered with the liquid iron [22]. On the other hand, the fully Lagrangian approach, which can track moving calculation points, has been adopted for solving multi-fluid physics problems during high-temperature metallurgical processes with complicated interfaces [23,24]. The smoothed particle hydrodynamics (SPH) method discretizes a continuous fluid phase by moving particles and is suitable for analyzing interfacial flow, even for numerous dispersed phases. This method can track the dripping flow [25], the movement of both the gas and the liquid phase [26], fluid flow in bed structures having different shapes (coupled with multi-sphere DEM) [27,28], and solid particle penetration into the liquid iron bath [29,30] directly. Recently, with the increase in computer capacity, the SPH method has been applied to very complicated interface problems [31]. For convenience, surface forces are converted into body forces in the SPH method, but recent expansions in computer capacity have alleviated this weakness. In this study, the trickle flow of a liquid iron and molten slag was simulated in the lower part of the blast furnace, and the effect of the wettability of molten slag and coke was investigated.

#### **2. Methods**

#### *2.1. Basic Formulation of SPH*

From the basic idea of SPH, it is evident that fluid motion is represented using a set of particle motion equations. A kernel function is introduced as an integral interpolation to solve a differential equation. In other words, the formulation is based on an interpolation scheme that allows the estimation of a vector or scalar function *f* at position **r** in terms of the values of the function at the discretization points.

$$f(\mathbf{r}) \cong \int f(\mathbf{r'}) \mathcal{W}(\mathbf{r}, h)dV \tag{1}$$

In Equation (1), **r** denotes arbitrary coordinates, **r** denotes particle positions, *V* denotes the volume, *W* is the smoothing kernel function, and *h* is the radius of influence. The summation of function *f* can be replaced with a summation over particles only within the distance *h* from **r***<sup>i</sup>* owing to the compactness; thus, *W* **r***ij*, *h* = 0 when **r***ij* <sup>&</sup>gt; *<sup>h</sup>*. The kernel must possess a form symmetrical to **r***ij* <sup>=</sup> 0. Here, *i* is the particle index, and *j* is the index of the neighboring particle around *i*. The kernel has at least a continuous first derivative and must satisfy the normalization condition as % *W* **r***ij*, *h d***r** = 1. Within the *h* → 0 limit, the kernel is required to reduce to a Dirac delta function δ **r***ij* . Wendland's kernel can be applied to prevent having various kernel artifacts in a multiphase system [32]:

$$\mathcal{W}\{\mathbf{r}\_{\rm ij},h\} = \frac{21}{16\pi h^3} \begin{cases} \left(1 - \frac{q}{2}\right)^4 (2q+1), & q < 2\\ 0, & q \ge 2 \end{cases} \tag{2}$$

where **r***ij* <sup>0</sup> is the interparticle distance corresponding to the initial conditions, *<sup>q</sup>* <sup>=</sup> **r***ij* /*h*, and it is assumed that *h* = 1.05 **r***ij* <sup>0</sup> [27]. The gradient form of Equation (1) can be represented by using the divergence theorem as follows:

$$
\nabla f(\mathbf{r}) \cong -\int f\left(\mathbf{r'}\right) \nabla W(\mathbf{r}, h)dV\tag{3}
$$

When applying this approximation to dispersed fluids, the discontinuities in the density distribution of the fluids become significant, which increases numerical errors. Nonuniform distribution or insufficiency of particles in the regions near the interface results in a significant fluctuation of pressure. The moving least squares (MLS) method is useful for approximating the function to solve this problem [33], which is related to the numerical fluctuations in the pressure at the nearby interface. The pressure is a function of the local density, and thus, the smooth density field of a bulk phase leads to a continuous pressure distribution. The MLS method improves the mass-area-density consistency and filters small-scale pressure oscillations, as described briefly in the following section.

#### *2.2. Density Approximation*

The SPH formulation can be transformed into a particle-based format to express the mass-area-density consistency process. The density of the particles was expressed in terms of the sum of the kernel functions of the N particles present within the radius of influence as follows:

$$\rho\_{\bar{i}} = \sum\_{j=1}^{N} m\_{j} \mathcal{W}(\mathbf{r}\_{ij\prime}, h) \tag{4}$$

where the subscripts *i* and *j* denote the particle indices, *m* is the mass, and ρ is the local density around the particle. Further, *mj* is the mass of particle *j*. Therefore, the kernel function around particle *i* can be discretized using the following equation, which is derived from Equations (1)–(3):

$$f\_i(\mathbf{r}) = \sum\_{j=1}^{N} \frac{m\_j}{\rho\_j} f(\mathbf{r}\_j) \mathcal{W}(\mathbf{r}\_{ij}, h) \tag{5}$$

The gradient of *fi* can be represented as expressed in Equation (6):

$$\nabla f\_i(\mathbf{r}) = -\sum\_{j=1}^{N} \frac{m\_j}{\rho\_j} f(\mathbf{r}\_j) \nabla \mathcal{W}(\mathbf{r}\_{ij}, h) \tag{6}$$

MLS involves a first-order consistent gradient approximation, which allows pressure smoothing, and its first derivative values are obtained using the method for the homogenous bulk phase mentioned above. The method of the least square interplant with constraint condition (CLS) represents an improved scheme, leading to a more accurate approximation than the MLS method around the sampling points [33,34]. In the CLS method, the particles can be made to directly represent a physical quantity by extending the MLS method in the one-dimensional error space for multiple dimensions. In the three-dimensional space, the CLS method approximates the values of various physical parameters around the particles based on the following equation: \* ρ + *<sup>i</sup>* = *a*<sup>0</sup> + *a*1(*x* − *xi*) + *a*2(*y* − *yi*) + *a*3(*z* − *zi*), where, *x*, *y*, and *z* are the coordinates of the sampling points and *a*0, *a*1, *a*2, and *a*<sup>3</sup> are the undetermined coefficients. If *x* = *xi*, *y* = *yi*, and *z* = *zi*, then \* ρ*i* + = *a*<sup>0</sup> at particle *i*. Refer to a previous report for parameters determining the procedure [26].
