The present study numerically investigated the liquid-gas-particles mixture flow during dam break using the computational fluid dynamics (CFD) and discrete element method (DEM). The adopted numerical methods were validated through comparison with previous researches of computation and experiments. The focus is on the effect of the particle friction coefficient (PFC) and the ratio of a particle size to a baseline size (PS) on behavior of the particles and the free surface. The PFC ranges from 0.1 to 0.9 and the PS varies from 50% to 100%. The liquid front position and the formation of the free-surface are almost independent of PFC and PS. In contrast, the particle front position strongly depends on PFC and PS. The particle front velocity showed two pattern according to time. First pattern is the increasing and then saturated behavior. Second pattern has increasing and decreasing behavior. These two patterns are mainly governed by the PFC. Generally, a smaller PFC and a larger PFC lead to first and second patterns, respectively. We could classify all cases of PFC and PS by the arrival time of particles into earlier and later arrivals than liquid front.

## I. INTRODUCTION

The presence of multiple phases such as liquid, gas, and particles can influence the behavior of a flow due to their interaction. The interaction is naturally complex and nonlinear, which makes the investigation of liquid-gas-particle flows extremely challenging. Nevertheless, these types of flows have attracted numerous researchers because they are widely observed and play an important role in various industries and nature.^{1–17}

It is difficult to study such flows through field observations and simple physical experiments. Therefore, most studies on these flows make significant simplifications, but they include the essential material properties and boundary conditions.^{18} In this view, the liquid-gas dam break flows with particles can be regarded as a simplified problem that involves essential elements of complex multiphase flows.

The numerical studies related to the liquid-gas-particle flows in dam breaks have mainly developed and validated numerical methods of such flows. Few researchers have investigated the behavior and characteristics of liquid-gas-particle flows. Guo and Morita^{17} compared the behavior of liquid-gas-particle flows according to height of the particle bed. They showed that the front position of the liquid became faster as the particle bed became higher. Xia et al.^{19} and Ran et al.^{20} emphasized that the sediment bed scouring and grain movement in a dam break are interesting subjects for researchers. Therefore, they investigated the propagation and physics of scouring near the interface between particles and the liquid.

Lin and Chen^{21} investigated the effect of the number of particles on the behavior of liquid-gas-particle flows and the impulse on a structure. The flow speed increased with decreasing concentration, while the impulse decreased. In addition, they examined the effects of parameters such as the stiffness, damping ratio, and friction coefficient of particles. The results showed that the friction coefficient has an obvious influence on the flow, while the stiffness and damping ratio did not. However, the parameters considered in that study were limited, and only the front position of the mixture flow was compared to investigate their effects.

Many numerical models for the computational fluid dynamics (CFD) have been developed for fluid flow with a free surface based on Eulerian or Lagrangian frameworks. For particle motion, the discrete element method (DEM) calculates the collision forces between particles and has been used to capture their motion.^{22} This method has been utilized by some researchers to model particle motion.^{23,24} Sun and Sakai^{25} and Guo et al.^{26} successively computed the behavior of a liquid-gas-particle mixture flow by adopting CFD-DEM which was validated by comparing with results of the experiments.

Also, the fluid flow has been simulated as particle motion in a Lagrangian framework using smoothed-particle hydrodynamics (SPH)^{27} and the moving particle semi-implicit (MPS) method.^{28} For instance, the sloshing was studied by using MPS-DEM which realized the liquid and the particles, respectively (Guo and Morita^{17}). The SPH-DEM was used to numerically simulate two phase flows with liquid-particle for a dam break and a rotating tank.^{29}

In this study, we used CFD-DEM to simulate the liquid-gas-particle flow because recent studies showed reasonable results using this method.^{25,26} The numerical method used in this study consisted of OpenFOAM for CFD^{30} and CFDEMproject.^{31,32} The present numerical methods were successively used in our previous study for the multiphase flow of a liquid-gas-particle mixture in dam break.^{33}

Previous studies related to liquid-gas-particle flows have not been thoroughly discussed compared to liquid-gas flows in dam breaks.^{25} Since the general particles in a granular flow have different characteristics and shapes, they can affect the flow. Their coefficients of friction and elastic restitution may affect the flow behavior due to interactions with other particles or the walls.^{34} In addition, the effect of grain size on the hydraulic permeability is also an important characteristic of particles.^{18} In this sense that particle’s properties have strong influence on the multiphase flow, we had investigated dependence of multiphase flow on particle density.^{33}

Park et al.^{33} studied the effect of particle density on three phases of liquid-gas-particle flow which occurs in a dam break. The adopted numerical method of CFD-DEM successively realized the interaction of two-phase fluids and particle. They systematically analyzed the dependence of the front head speeds of the free-surface and particles on the density of the particle. They reported that the arrival time of the front head of the free surface to the wall is almost the same for all particle densities. Otherwise, the arrival time of the front head of particles becomes shorter with increasing the particle density. The decrease rate of the arrival time of the front head of particles is about the linear relation to the particle density. Finally, they provided the motion map.

Authors’ previous study^{33} considered only one particle’s property, finding the effect of the particle density on the complex multiphase flow of liquid-gas-particle mixture in a dam-break. The present study extents our previous study^{33} by considering more particle’s properties. Therefore, the purpose of the present study investigates the effects of the particles’ friction coefficient and size on their behaviors and the free surface when the mixture moves in a dam break problem. We use numerical methods to simulate the liquid-gas-particle flow. This may obtain more useful information than experimental approaches because they are still limited in observing complicated phenomena.^{25} The governing equations used in this study are presented in Section II. The behavior of particles and the free surface with different particle friction coefficients and diameters are compared in Section III to investigate their effects.

## II. NUMERICAL MODEL

The present study uses the same numerical methods as our previous paper of Park et al.^{33} which provides a very detailed explanation for numerical methods and used OpenFOAM. Therefore, we very concisely describe the present numerical methods.

### A. Governing equations and numerical methods

The motion of particles is realized by DEM. The translational motion is governed by the following equation:

where the subscripts of *i* and *j* identify the particles. Thus, *m*_{i} and *v*_{i} are a mass and velocity of particle *i*, respectively. *F*_{n,ij}, *g*, *F*_{fluid} (*F*_{fluid} = *F*_{d} + *F*_{p}) are the normal contact force, gravity, and the fluid force, fluid drag force, and fluid pressure gradient force, respectively.

The governing equation of the rotational motion of particles is as follows;

where *I*_{i}, *ω*_{i}, *R*_{i}, and *F*_{t,ij} are the moment of inertia, angular velocity, distance from the center of mass to the contact point, and contact force in the tangential direction, respectively. The definitions of each forces in equations (1) and (2) can be refereed in Park et al.^{33} Note, the PFC (particle friction coefficient) is a key parameter in this study. The present PFC is the static friction coefficient which is defined as the ratio of the tangential force (*F*_{t,ij}) to the normal force (*F*_{n,ij}).^{33} In this study, the static friction coefficient is symbolized as *μ*_{S} which corresponds to *C*_{fs} in Park et al.^{33}

The continuity and the momentum equations govern the present unsteady viscous three-dimensional fluid flow as follows;

where *x*_{i}, *u*_{i}, and *ε* are the is a Cartesian coordinate corresponding velocity component, and the porosity, respectively. In equation (4), *p*, *ρ*, and *μ* are the pressure, the density, and the viscosity, respectively. The two-way coupling in fluids and particles is performed by the source term *S*_{m} in equation (4).

The VOF method is utilized to track the free surfaces as the interface between the liquid and gas. The free surface in each cell is defined when a volume fraction (*Q*) corresponding to each fluid is between 0 and 1. The transport of *Q* is governed by following equation;

The porosity (*ε*) in equation (5) is considered in the conventional VOF equation, which represents the fluid displaced by particles.^{35} The porosity is a fraction of void volume over the total volume.

The spatial discretization for the convection and diffusion terms are carried out by using a second-order upwind scheme and second-order central differencing scheme, respectively. The temporal discretization of the time derivative terms are performed by a second-order accurate backward implicit scheme. The PISO (Pressure Implicit with Splitting of Operators) algorithm^{36} is used for the velocity–pressure coupling and overall solution procedure. In addition, the schematics of a flow chart for the algorithm of the CFD-DEM can be refereed in Park et al.^{33}

### B. Computational conditions

Figure 1 shows that the three-dimensional computational domain is composed of a length, a breadth and a height (*L*, *B*, *H*) with dimension of (0.2 m, 0.1 m, 0.3 m). The initial liquid column is formed with a length, a breadth and a height (*L*_{l}, *B*_{l}, *H*_{l}) with dimension of (0.05 m, 0.1 m, 0.1 m). The number of particles is 3,883 which is stacked on the left side bottom. The air occupies the remained space. The properties of the particles and fluids are recorded in Table I and Table II, respectively.

Type . | Value . |
---|---|

Baseline size, d_{o} (m) | 0.0027 |

Diameter, d_{p} (m) | 0.00135, 0.00162, 0.00189, |

0.00216, 0.00243, 0.0027 | |

Density, ρ_{p} (kg/m^{3}) | 2,500 |

The number of particles | 3,883 |

Spring constant (N/m) | 1,000 |

Restitution coefficient | 0.9 |

Friction coefficient, μ_{s} | 0.1, 0.2, 0.3, 0.4, 0.5, |

0.6, 0.7, 0.8, 0.9 |

Type . | Value . |
---|---|

Baseline size, d_{o} (m) | 0.0027 |

Diameter, d_{p} (m) | 0.00135, 0.00162, 0.00189, |

0.00216, 0.00243, 0.0027 | |

Density, ρ_{p} (kg/m^{3}) | 2,500 |

The number of particles | 3,883 |

Spring constant (N/m) | 1,000 |

Restitution coefficient | 0.9 |

Friction coefficient, μ_{s} | 0.1, 0.2, 0.3, 0.4, 0.5, |

0.6, 0.7, 0.8, 0.9 |

Type . | Value . |
---|---|

Liquid density, ρ_{l} (kg/m^{3}) | 1,000 |

Liquid viscosity (Pa·s) | 10^{-3} |

Gas density, ρ_{g} (kg/m^{3}) | 1 |

Gas viscosity (Pa·s) | 10^{-5} |

Type . | Value . |
---|---|

Liquid density, ρ_{l} (kg/m^{3}) | 1,000 |

Liquid viscosity (Pa·s) | 10^{-3} |

Gas density, ρ_{g} (kg/m^{3}) | 1 |

Gas viscosity (Pa·s) | 10^{-5} |

For the case of a dam break flow, the Reynolds number is generally defined as $Re=HgH\rho f/\mu f$ where *ρ*_{f}, *μ*_{f} and *g* are the density and viscosity of fluid, and the gravity, respectively. This definition of Re leads to about 1.0×10^{5} as the value of the present Reynolds number. Like our previous study,^{33} we concentrate on the fluid flow and behavior of particles on a dry bed. From the collapse to the collision to the opposite wall, the dam break flow on a dry bed does not exhibit strong turbulence.^{37} Thus, this period on a dry bed is mainly exposed to the laminar flow rather than the turbulence. Therefore, the turbulence phenomena and turbulence model are not considered in this study.

In order to evaluate the drag force action on each particles in DEM, the particle Reynolds number (Re_{p}) is needed.^{33} The definition of Re_{p} is $Rep=\epsilon dpvf\u2212vp\rho f/\mu f$, where *d*_{p}, *v*_{f} and *v*_{p} are the particle diameter, fluid velocity and the particle velocity, respectively. The values of the present Re_{p} according to the particle diameter are in the range of about 1-4×10^{2}.

The structured grid is used and the grid size for the fluid flow is bigger than the particle size to avoid numerical errors in CFD-DEM. The vertical and horizontal walls of the computational domain were considered to have no-slip conditions. The dimensionless parameters to identify the effect of the particle friction coefficient and diameter are achieved by the characteristics length and time of *L*_{l} and $2g/Ll\u221212$, respectively.

### C. Validation of the numerical method

The validation of the present numerical method has been perform by comparing the present front head positions of the free surface and the particles with the previous computational and experimental results of Sun and Sakai.^{25} Figure 2 shows that the present results give a great agreement with both numerical and experimental results, which supports the reliability of the present numerical methods. In addition, the qualitative comparisons are shown in Fig. 3 where the time evolution of the formation of the free-surface and particles is presented with that of Sun and Sakai.^{25} The thick blue line indicates the free surface. The particles are colored individually with their velocity.

Sun and Sakai^{25} considered the effects of a rising gate by adding an obstacle in front of the liquid column. However, the present study ignores the gate effect because the focus is on the mixture flow. The overall process of the free surface and particle motion was similar to that of Sun and Sakai^{25} except for some differences in the free surface in the initial stages. After the collapse of the liquid column, the free surface advances against the wall, the particle bed expands over the bottom of the domain, and the particle front moves more slowly than the liquid front at *t*^{*}< 4 in Figs. 3(a) and (b) Even after hitting the opposite wall, the splash of the particles and the upward motion of the free surface are reasonably represented by the present results at *t*^{*}= 6 in Fig. 3(c).

## III. RESULTS AND DISCUSSION

### A. Motion of the front heads of the liquid and particles

We identified the effects of the particle friction coefficients (PFC) and particle sizes (PS) on the three-phase liquid-gas-particle flow. Figures 4 and 5 show typical three-dimensional perspective views and side views of the particles and the free surface between the liquid and the gas for different PFCs and PSs at *t*^{*}=1.0, 2.0, 3.0, and 3.5. In these cases, PFC=0.1 and 0.9, and PS=0.5 and 1.0.

When the liquid column starts to collapse due to gravity, the front of the liquid and particles advances against the wall. At a fixed time, the liquid front position and the formation of the free-surface are almost independent of the variation of PFC and PS, which is clarified by the side views in Fig. 5. This independent behavior of the liquid front position and the formation of the free-surface can be also founded in Park et al.^{33} where the variation of the particle density is considered to investigate a multiphase liquid-gas-particle flow from a dam break using CFD-DEM.

In contrast, the particle front positions differ according to PFC and PS at a fixed time, which can be clarified by comparing the liquid front independently of PFC and PS. Briefly, a smaller PFC contributes to the fast advance of particles when PS is the same, as shown in Fig. 4. A smaller particle advances quickly when the PFC is the same, as shown in the figure. Eventually, as both PFC and PS become smaller, the particles become faster and finally reach the wall earlier, as shown in Fig. 4 and Fig. 5. Particularly, the largest PS and PFC lead to slower travel of the particle front than the liquid front up to the wall, as shown in Fig. 4(d) and Fig. 5(d).

To clarify the time evolution of the front position of the liquid and particles, the spanwise-averaged dimensionless front position of the liquid and particles ($<xl*>$ and $<xp*>$) along the time are presented in Fig. 6. Additionally, only two phases of liquid-gas are considered and then a dam break flow without particles is numerically simulated to identify the effect of the particles on the dam break flow by comparing with those of the liquid-gas-particle flow.

In general, the liquid front in the liquid-gas-particle flows advances slower than that in the liquid-gas flow against the wall. Irrespective of the existence of particles, the liquid front position is almost the same at the beginning of the dam break flow, but the particles give an effect on the front head position of the liquid after *t*^{*}=0.5. As time goes by, the difference in the liquid front position increases until the liquid front hits the opposite wall, as observed in Fig. 6(a).

Irrespective of PFC and PS, the profile of the liquid front position in the liquid-gas-particle flows are almost the same, as shown in Fig. 6(a), which confirms the front formations of the liquid in Figs. 4 and 5. As a result, it can be concluded that the PFC and PS give a negligible effect on the liquid front.

The effect of PFC and PS on the time evolution of the particle front position is identified in Fig. 6(b). In order to direct compare the liquid front head, the positions of liquid front head without and with particles are put together with the particle front positions in Fig. 6(b). The particles move faster when simultaneously decreasing the size and friction coefficient. Thus, the smallest values of PS=0.5 and PFC= 0.1 lead to the earliest arrival of the particle front at the opposite wall. Particularly, this particle front runs ahead of the liquid front at all times until the impingement. The largest values of PS=1.0 and PFC=0.9 lead to the slowest advancing speed of the particle front over the time. Consequently, the arrival time of the particle front at the opposite wall is strongly dependent on PS and PFC. This will be discussed later for Fig. 11 where we cover all ranges of PS and PFC and classify the cases of PS and PFC into earlier and later arrivals of the particle front than the liquid front head.

Recently, authors’ previous study^{33} considered the effect of the particle density on the multiphase liquid-gas-particle flow from a dam break using CFD-DEM. We showed almost the invariance of the liquid front position and the formation of the free-surface according to the particle density. However, the particle front positions strongly depended on the particle density. The PFC and PS also gives the same effect on the liquid front position, the formation of the free-surface, and the particle front position as the particles density.

Figure 7 presents the dimensionless velocity and acceleration of the particle front for different PFC and PS. In addition, the liquid fronts in the liquid-gas and the liquid-gas-particle flows are considered in Fig. 7. When particles are present, the liquid front velocity reaches the maximum value at *t*^{*}=2.58. The time obtained without particles is about 2.5, as shown in Fig. 7(a).

The velocity of the liquid front without particles is always faster than that of particles until the liquid front velocity reached the maximum velocity. The liquid front velocity reduces irrespective of the existence of particles after the maximum velocity of the liquid front. The velocity of the liquid in case of the non-existence of particles diminishes more quickly than that in case of the existence of particles.

In contrast to the liquid front velocity in the liquid-gas-particle flows, those of the particle front are affected by PFC and PS. The profile of the particle front velocity can be divided into two groups. The velocity profile of the first group is composed of two sections. The first section has increasing velocity, while the second section has almost constant velocity. The velocity profile of the second group has increasing behavior, and reaches to the maximum. The second section has decreasing velocity.

Figure 7(b) shows the acceleration which correspond to the velocities in Fig. 7(a). The largest PS with the smallest PFC results in the largest inertia. The deceleration for the decreasing rate of the velocity is the smallest among the cases of PFC and PS considered in this study. As a result, for a smaller PFC, the inertia of the particles is preserved further downstream until the vertical wall, which leads to the increasing and constant behavior of the velocity. Otherwise, in case of the largest PFC, the magnitude of the deceleration becomes quickly increasing, which could be contributed to a larger friction force.

### B. Isometric views of the front heads of the liquid and particles

To cover all ranges of PFC and PS considered in this study, isometric views of the front positions of the liquid and particles are presented in Fig. 8. The coordinates for the liquid and the particles consist of (*d*_{p}, *μ*_{s}, $xl*$) and (*d*_{p}, *μ*_{s}, $xp*$), respectively. The liquid front position shows a plane surface at each time, as shown in Fig. 8(a), which means that the liquid front is independent of PFC and PS. However, the effect of PFC and PS on the particle front position is apparent, as shown in Fig. 8(b). The difference in the particle front position for different parameters is minor at *t*^{*}=1, and the surface is nearly planar. However, as time goes by, the effect of PFC and PS on the particle front position becomes bigger, resulting in the augmentation of the difference of each particle front position according to PFC and PS. The maximum value appears at a PS of *d*_{p}=0.5*d*_{o} and PFC of *μ*_{s}=0.1 at every time, which means that the particle front under these conditions has the fastest velocity among properties considered in this study. The early plane shape become a 3-D rhombus, which means that the discrepancy of the upper vertex (the maximum value) and the lower vertex (the minimum value) becomes bigger as time goes by.

The difference between the particle front position and the liquid front position is plotted as 3-D as function of PFC and PS in Fig. 9. The dashed line indicates that the front positions of the liquid and particles are the same. The figures correspond to *t*^{*}= 1, 2, 3, and 3.5. In the early stage, except for PS of *d*_{p}=0.5*d*_{o} and regardless of PFC, the difference of the front position of the liquid and particles is close to 0% for different PFCs and PSs at *t*^{*}=1.

The particle front was ahead of the liquid front for PS of *d*_{p}=0.5*d*_{o}, regardless of PFC (Fig. 9(a)). As time goes by, for more cases considered in this study, the particle front lags behind the liquid front at *t*^{*}=2 except at PS=0.5, regardless of PFC. However, at *t*^{*}=3, the difference of the liquid and particle front positions is close to 0%, which means that the particle front position overlapped with the liquid front position in most cases. In addition, the difference between the maximum value at *d*_{p}=0.5*d*_{o} and *μ*_{s}=0.1 and the minimum value at *d*_{p}=1.0*d*_{o} and *μ*_{s}=1.0 is the biggest at *t*^{*}=3.

Before the particle front collides against the other side walls, the trends lasted for *t*^{*}=3.5. However, the degree of difference between the maximum and minimum values is somewhat alleviated. In addition, when PFC is 0.5 and the range of PS is 0.1–0.3, the particle front is always ahead of the liquid front, while the particle front always lags behind the liquid front when the ranges of PFC and PS are 0.6–1.0 and 0.3–0.9, respectively.

The liquid has an effect of the reflection force against the approaching wall, leading to the decrease of the velocity. However, the particles keep the inertia close to the wall. Thus, the smaller particles with smaller friction coefficient keep maintains the velocity, which support the tendency of increasing and constant velocity, as discussed in Fig. 6.

### C. Map of the arrival time

In order to identify the material that arrives first on the opposite wall, the formation of the free-surface and particle distribution for representative cases are presented in Fig. 10 when particle or liquid arrive at opposite wall. To cover all ranges of PFC and PS, the arrival times of the liquid and particle fronts are quantitatively presented in Fig. 11.

At the lowest particle size of *d*_{p}=0.5*d*_{o} among the particle sizes considered in this study, the particle front arrives earlier than the liquid front head, which is independent of the particle friction coefficient, as shown in Fig. 10(a). This early arrival of the particle fronts at the lowest particle size of *d*_{p}=0.5*d*_{o} can be clarified by Fig. 11 where the arrival time of the liquid front is about $ta*$=4.0, regardless of the parameters of the particles.

Like the lowest particle size of *d*_{p}=0.5*d*_{o}, the smallest friction coefficient of for *μ*_{s}=0.1 leads to earlier arrival of the particle front than the liquid front head, regardless of the particles size, as shown in the first column of Fig. 10(a-f). Thus, Fig. 11 shows smaller arrival time of the particles front for *μ*_{s}=0.1 than the liquid front head.

In general, at a fixed PS, The arrival time of the particle front increases with increasing PFC. However, as PS increases, the increasing rate of arrival time of the particle front with increasing PFC becomes larger. When PS was 0.5, the difference in the maximum and minimum arrival times depending on PFC WAS about 3%. However, when PS was 1.0, the difference was about 15%. This dependence of the arrival time of the particle front on PFC at a fixed PS is the same as the variation of PS at a fixed PFC. Consequently, a bigger particle size corresponded with a bigger effect of PFC on the arrival time.

Finally, we classify the particle size and the particle friction coefficient into two groups where the particle front arrives faster or not. One is liquid front first arrival and another is liquid front first arrival. These two groups are mapped in Fig. 12 where the x- and y-axes represent PFC and PS, respectively.

## IV. CONCLUSION

The present study numerically investigated the liquid-gas-particles mixture flow during dam break using the CFD-DEM numerical method. The adopted numerical methods were validated through comparison with previous researches of computation and experiments. The main purpose is to investigate the effects of particle size (PS) and the particle friction coefficient (PFC) on the behavior of the particles and the free surface. To identify the effects of PFC and PS, the motion of the particles and free surface, the front position of the liquid and particles, and the arrival times were analyzed and compared for different values of PFC and PS.

The liquid front was almost independent of PFC and PS, but the particle front was strongly affected by them. The particles moved faster when simultaneously decreasing the size and friction coefficient. Thus, the smallest values of PS and PFC lead to the earliest arrival of the particle front at the opposite wall. Recently, authors’ previous study^{33} considered the effect of the particle density on the multiphase liquid-gas-particle flow from a dam break using CFD-DEM. We showed almost the invariance of the liquid front position and the formation of the free-surface according to the particle density. However, the particle front positions strongly depended on the particle density. The PFC and PS also gives the same effect on the liquid front position, the formation of the free-surface, and the particle front position as the particles density.

The liquid front velocity reduces irrespective of the existence of particles after the maximum velocity of the liquid front. The velocity of the liquid in case of the non-existence of particles diminishes more quickly than that in case of the existence of particles. Especially, the liquid front velocity was independent of the PS and PFC, which is consistent with the time variation of the liquid front head. The particle front velocity showed two pattern according to time. First pattern is the increasing and then saturated. Second pattern has increasing and decreasing behavior. These two patterns are mainly governed by the PFC. The smaller PFC leads to first pattern, since the inertia overcomes the weak friction, regardless of the PS. The larger PFC leads to second pattern, since the inertia becomes weak due to relatively the strong friction, regardless of the PS.

## ACKNOWLEDGMENTS

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) through GCRC-SOP (No. 2011- 0030662).