- Email: [email protected]
CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break
Accepted Manuscript
CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break Kyung Min Park , Hyun Sik Yoon , Min Il Kim PII: DOI: Reference:
S1007-5704(17)30383-0 10.1016/j.cnsns.2017.11.010 CNSNS 4367
To appear in:
Communications in Nonlinear Science and Numerical Simulation
Received date: Revised date: Accepted date:
25 April 2017 13 September 2017 9 November 2017
Please cite this article as: Kyung Min Park , Hyun Sik Yoon , Min Il Kim , CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break, Communications in Nonlinear Science and Numerical Simulation (2017), doi: 10.1016/j.cnsns.2017.11.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
Complex multiphase flow of liquid-gas-particle mixture in a dam-break is numerically investigated according to particle density.
CFD-DEM coupled method is successively adopted to find the natural phenomena involved in two-phase fluids and particle interaction.
The head speeds of the free-surface and particles show different patterns to particle density. The response of the free-surface and particles to particle density is identified
CR IP T
AC
CE
PT
ED
M
AN US
by three motion regimes of the advancing, overlapping and delaying motions.
1
ACCEPTED MANUSCRIPT
CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break Kyung Min Parka, Hyun Sik Yoonb*, Min Il Kimb, a
Global Core Research Center for Ships and Offshore Plants, Pusan National University, 2, Busandaehak-ro 63beon-gil, Gumjeong-Gu, Busan, 46241, Republic of Korea
b
CR IP T
Department of Naval Architecture and Ocean Engineering, Pusan National University, 2, Busandaehak-ro 63beon-gil, Gumjeong-Gu, Busan, 46241, Republic of Korea
ABSTRACT
AN US
This study investigates the multiphase flow of a liquid-gas-particle mixture in dam break. The open source codes, OpenFOAM and CFDEMproject, were used to reproduce the multiphase flow. The results of the present study are compared with those of previous results obtained by numerical and experimental methods, which guarantees validity of present numerical method
M
to handle the multiphase flow. The particle density ranging from 1100 to 2500 kg/m3 is considered to investigate the effect of the particle density on the behavior of the free-surface
ED
and the particles. The particle density has no effect on the liquid front, but it makes the particle
PT
front move with different velocity. The time when the liquid front reach at the opposite wall is independent of particle density. However, such time for particle front decrease as particle density
CE
increases, which turned out to be proportional to particle density. Based on these results, we classified characteristics of the movement by the front positions of the liquid and the particles.
AC
Eventually, the response of the free-surface and particles to particle density is identified by three motion regimes of the advancing, overlapping and delaying motions.
Keywords: DEM, CFD, Liquid-gas-particle mixture flow, Dam break, Three motion regimes
*
Corresponding author, Tel.: +82-51-510-3685, E-mail: [email protected] 2
ACCEPTED MANUSCRIPT
1. Introduction
The study of multiphase flows such as liquid-gas-particle mixture flow is challenging issue because the interactions between different phases result in complex behavior. Nevertheless, it is important to study multiphase flows because they are a common phenomenon in various areas in
CR IP T
academics, industry, and real life.
In geology, the multiphase flows on mountains are considered an important subject, where such flows are known as debris flows, which include mudslides, mudflow, and debris avalanches.
AN US
The compositions of a debris flow are mainly water, rock, and air, and the forces of particles and fluids greatly influence the flow [1]. Hilly areas and intense rainfall are particularly susceptible to debris flows. These flows are dangerous to life and property because of the high sediment concentrations and mobility. Many properties were inundated or crushed by the debris flows,
M
pastures and cornfields were buried, and livestock perished [2].
During volcanic eruptions, mixtures of lava, gas, and ash generate hot multiphase flows. There
ED
are two types of such flows, and the first type is lava flow. Lava flow occurs when molten or
PT
partially molten rock erupt from the surface of the earth. The melt spreads on the surface as a gravity current and forms a lava flow. The lava flow is not a simple liquid but generally a mixture
CE
that contains magma, rock, ash, crystals, and gas bubbles. Lava flows play important roles in the development of planetary crusts and ore deposits, and they pose dangers to people and property,
AC
which means the importance of research on such flows [3]. Another type of multiphase flow in volcanic eruptions is pyroclastic flow. When the density of
the ash particles is higher than the gas density, the eruption column may collapse and form a pyroclastic multiphase flow. A pyroclastic flow can propagate distances of several kilometers from its starting point [4]. Such flows therefore also pose dangers to populations and infrastructure [5-7], and prediction of their behavior is crucial for hazard assessment.
3
ACCEPTED MANUSCRIPT Also, in ocean the complex multiphase flow occurs. Littoral sediment transport is caused by the interactions of waves, currents, and tides. These phenomena determine the rate of beach erosion or accretion and act as boundary conditions for nearshore hydrodynamic and morph dynamic models. The study of nearshore multiphase flows is important for coastal engineering, environmental protection, recreation, and military operations. Some models have therefore been
CR IP T
developed to simulate and predict their behaviors [8-10]. Offshore multiphase flows play an important role in turbidity currents. Turbidity current occurs when the edge of the continental shelf fails and causes a submarine landslide. As the turbidity current propagates down the slope, it erodes the ocean floor and increases the amount of
AN US
suspended sediment. When the current reaches the abyssal plain, deposition becomes dominant and another sediment layer is added to the geological record of the ocean floor. The sediment layers deposited by sequential turbidity currents can form traps for oil. Sound understanding of
M
the dynamics of turbidity currents is instrumental for the reconstruction of the sediment layer thickness and potential location of oil reservoirs from drilling cores and seismological data [11].
ED
Multiphase flow also occurs in liquid-metal-cooled reactors such as mixture flow of liquid metal and resolidified fuel particles. The dynamics of the flow strongly affects the risk of
CE
phase [12].
PT
accidents, which makes it important to analyze the movements and interactions between each
As illustrated these relevant fields governed by the liquid-gas-particles mixture flow, the
AC
multiphase flows appear in many areas. However, these problems are very complex to be numerical simulated, since these problems are also considered as one of the challenging problems in the CFD and also simulation community. Therefore, as the simplified problem but involving essential elements of these problems, a dam break flow containing particles can be considered as reasonable problem to investigate the multiphase flows such as liquid-gas-particle mixture flows. In general, the classical dam breaks have been well studied in recent decades [13]. For
4
ACCEPTED MANUSCRIPT example, the classical dam break flow problem has been widely studied as a typical liquid-gas mixture flow problem to identify the interactions of each phase. Studies have examined the development of numerical methods [14-20], the effects of the density ratio of the fluids [19, 20], the impact on structures [21-23], the presence of obstacles [24], the effects of the bed slope [25] and turbulence intensity [26], and the physics of a two-phase flow [27, 28].
CR IP T
However, dam break flows containing particles have not been thoroughly discussed, there is a lack of numerical simulations and experiments dealing with flows of liquid-gas-particle mixtures in the dam break problem. Furthermore, the main focus of previous studies was to validate numerical methods by comparing simulation and experimental results rather than interpreting the
AN US
physics of the flow.
Such studies developed various numerical models in the Eulerian or Lagrangian framework to describe the motion of two-phase flows separated by an immiscible fluid interface. In terms of,
M
the particle motion was captured by the discrete element method (DEM), which has been widely applied in particle collision studies. The DEM model was essentially used to calculate the
ED
collision forces between particles or between particles and walls [29], and some researchers used the method to model particle motion [30, 31]
PT
In the Eulerian framework, the fluid flow is solved with a computational mesh by
CE
computational fluid dynamics (CFD), and the free surface of fluids was captured mainly by the volume of fraction (VOF) method [15]. For example, Lin and Chen [32] investigate the
AC
discrepancy of the particulate flow behavior and the impact of the particulate flow on the structure according to the number of particles using CFD-DEM. Guo et al. [33] reasonably simulated the transient behavior of the multiphase flow of a liquid-gas-particle mixture flow by CFD-DEM. The method was validated based on experiments of a water-particle dam break in tems of gas-liquidparticle flows. Sun and Sakai [13] developed CFD-DEM to simulate various three-dimensional flows of liquid-gas-particle mixtures. The numerical results were validated through comparison
5
ACCEPTED MANUSCRIPT with an experiment. In contrast, fluid flow in the Lagrangian framework considers particle motion based on smoothed-particle hydrodynamics (SPH) [34] and the moving particle semi-implicit (MPS) method [35]. For example, Sun et al. [36] used SPH-DEM to simulate liquid-particle flows such as those in a dam break and in a rotating tank. Guo and Morita [12] investigated the sloshing in a
and they ignored the effects of gas on the liquid and particles.
CR IP T
liquid-particle mixture flow using the MPS-DEM. However, they did not consider the gas phase
None of these studies provide a physical interpretation for various conditions, in terms of the fluid and particle properties and the liquid column. Especially, the multiphase flows such as debris
AN US
flow, turbidity flow, snow avalances, and pyroclastic flow have different particle densities ranging to 1000-3000 kg/m3 [37]. Therefore, the present study investigates the effects of the density of particles on their behaviors and the free surface in the dam break problem.
M
The CFD-DEM method was used for the simulation with the open source code in OpenFOAM for CFD [38] and CFDEMproject, which includes LIGGGHTS for DEM simulation
ED
and calculates the interaction data between CFD and DEM [39, 40]. The numerical method is effective for investigating the interaction of the different phases. A dam break flow with particles
PT
has intrinsic complexity and instability because of the interfaces between the phases. However, a
CE
few recent studies have reasonably simulated such flows using CFD-DEM [13, 33]. This paper is structured as follows. Section 2 presents the governing equations for the motion
AC
of each phase, the numerical scheme, the problem definition, and a validation case for the numerical method. The behavior of particles and the free surface between the liquid and gas at different particle densities was compared to identify the effects. Section 3 presents the different regimes of the particle and free-surface motion as a function of the particle/liquid density ratio and position.
6
ACCEPTED MANUSCRIPT
2. Numerical methodology and problem definition
Models to describe a dam break flow with particles can be divided into Eulerian models (in OpenFOAM) and Lagrangian models (in LIGGGHTS). Eulerian models calculate the fluid flow
CR IP T
based on CFD while capturing the liquid-gas interface provided by VOF. In the Lagrangian models, DEM calculates the path and motion of each particle. Particle interactions are modeled through specified collision rules. The detailed governing equations for each phase are described in
2.1 Governing equations of particle motion
AN US
the following.
M
In DEM, a particle has translational and rotational motion that depend on the forces and torques acting on it [29]. The forces and torques result from particle interactions with other
ED
particles, the walls, and the surrounding fluids. As shown in Fig. 1, there are three type forces: contact force, gravity force, and fluid force.
PT
The contact forces are considered based on a spring, dashpot, and friction-slider model. For a
CE
particle i with mass mi and velocity vi is in contact with another particle or geometric
AC
element j . The translation motion equation is expressed by the following: mi
dvi Fn, ij mi g Ffluid ( Fd Fp ) dt j
(1)
where Fn ,ij , g, and Ffluid are the contact force in the normal direction, gravity, and the fluid force due to the presence of fluid. The fluid force is a combination of the fluid drag force ( Fd
and
fluid pressure gradient force ( Fp ). The fluid force is also considered as fluid-particle interaction force based on Newton’s third law.
7
ACCEPTED MANUSCRIPT In Eq. (1), the normal-direction contact forces are calculated using a linear spring contact model:
Fn, ij K n dn Nn vn, ij
(2)
where K n , d n , N n , and v n ,ij are the normal spring stiffness, the particle overlap in the normal
CR IP T
direction at the contact point, the normal damping, and the normal velocity components of the relative sphere-surface velocity at the contact point. N n is defined by the following: N n 2 N n damp K n meq
(3)
where meq and N n damp are the equivalent particle mass and the normal damping coefficient: 1 1 1 mi m j
AN US
meq
ln(Cn rest )
2 ln(Cn rest )2
(5)
M
N n damp
(4)
coefficient of restitution.
ED
where mi , m j , and Cn rest are the mass of particle i , that of particle j , and the normal
PT
The drag force on the individual particles was calculated by the Gidaspow empirical model [41], which takes into account the effects of the surrounding particles. The drag force is defined
Fd 0.5Cd f A v f vi v f vi
(6)
AC
CE
by the following:
where Cd , f , A , and v f are the drag coefficient, fluid density, projected area of particle i , and fluid velocity. The drag coefficient is defined as: 4 1 CD 150 1.75 if min 3 Re p
Otherwise:
8
(7)
ACCEPTED MANUSCRIPT
CD
24 3.6 Re
0.687 p
Re p
3.65
(8)
where is the porosity, vmin is a cutoff porosity of 0.8, and Re p is the particle Reynolds number.
FP Vi pstatic
(9)
CR IP T
Eq. (9) represents the pressure gradient force, where Vi and pstatic are the volume of the particle i and the gradient of the static pressure in the fluid, respectively.
The particles also undergo rotational motion, which is governed by Eq. (10):
di R i Ft ,ij dt
AN US
Ii
(10)
where I i , i , Ri , and Ft ,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
M
tangential force is defined by the following:
ED
Ft ,ij Kt dt Nt vt ,ij
K n d n c fs dt dt
(12)
PT
Kt dt K n d n C fs ; Ft
(11)
where K t , d t , N t , vt , and C fs are the tangential spring stiffness, the overlap in the tangential
CE
direction at the contact point, the tangential damping, the tangential velocity components of the
AC
relative sphere-surface velocity at the contact point, and the static friction coefficient. N t is defined by the following: Nt 2 Nt damp Kt meq
(13)
where meq and Nt damp are the equivalent particle mass and the tangential damping coefficient defined by the following:
9
ACCEPTED MANUSCRIPT In Ct rest
Nt damp
2 In Ct rest
(14)
2
where Ct rest is the tangential coefficient of restitution.
CR IP T
2.2 Governing equations of fluid flow
The three-dimensional two-phase flow is governed by the continuity equation in Eq. (15) and the Navier–Stokes equation in Eq. (16):
AN US
ui 0 xi
p ui uiu j t x j xi x j
u u j i gi Sm x j xi
(15)
(16)
M
where x i is a Cartesian coordinate, ui is the corresponding velocity component, p is the
ED
pressure, is the density, is the viscosity, is the porosity, and g is gravity. The source term S m takes into account the interaction force between the fluid and particles by two-way
PT
coupling. Based on Newton’s third law, S m in a fluid cell is determined by adding up the fluid
Sm
DEM iteration Particle F Vcell
fluid
(17)
AC
CE
force of the particles located in the fluid cell:
where Vcell is the volume of the fluid cell. The VOF method was used to capture the free surface between the liquid and gas. This
method was designed for immiscible fluids, and each cell has a volume fraction ( Q ) corresponding to each fluid. A cell with a scalar Q value 0 indicates air, and a value of 1 represents liquid. The free surface in each cell is defined when the Q value is between 0 and 1.
10
ACCEPTED MANUSCRIPT At each time step, a transport equation is solved to find the distribution of the fluid [15]. In conventional VOF method, the porosity was additionally considered to represent the fluid displaced by particles [42]. The porosity is a fraction of void volume over the total volume such as
1
V
particle
Vcell
, where Vparticle is particle volume corresponding fluid cell. The distribution is
Q Q ui 0 t xi
CR IP T
found by solving the following transport equation of Q : (18)
The second-order upwind scheme and the second-order central differencing scheme were used
AN US
for the convection and diffusion terms, respectively. To the calculate the unsteady flow, the time derivative terms are discretized using a second-order accurate backward implicit scheme. The velocity–pressure coupling and overall solution procedure are based on a PISO (Pressure Implicit
ED
2.3 Algorithm of CFD-DEM scheme
M
with Splitting of Operators) algorithm [43].
PT
The overall algorithm of the CFD-DEM method is illustrated in Fig. 2. The calculation of
CE
CFD-DEM starts with DEM. The position and velocity of each particle are calculated by solving the governing equations of particle motion (1) and (10) based on the contact force and fluid force
AC
in DEM. The porosity in each cell is then calculated. Simultaneously, the fluid-particle interaction force is accumulated from the particles and transferred to the fluid phase. To calculate the governing equation of fluid flow in CFD, the VOF equation (18) is first calculated. Through this process, the volume fraction Q is obtained, which updates density and viscosity for each fluid in each fluid cell. The governing equations of the fluid flow (15) and (16) are then solved iteratively based on the PISO algorithm [43]. When the iteration converges, a new fluid velocity and pressure are obtained. After that, the fluid forces acting on each particle are calculated from 11
ACCEPTED MANUSCRIPT the results of CFD and are passed to the DEM solver using the CFD-DEM scheme. The time step in the DEM solver tends to be more limited than that in the CFD solver (it is 50 times smaller than that of the CFD solver in this study). This occurs because the total stability conditions are dominated by DEM [13]. This condition satisfies the qualification for particle
CR IP T
collision modeling in an extreme situation and ensures overall computational efficiency.
2.4 Computational conditions
The computational domain is shown in Fig. 3, where the length, breadth, and height are 0.2,
AN US
0.1, and 0.3 m respectively. The domain contains a liquid column and 3,883 particles. The liquid column’s initial length ( Ll ), breadth ( Bl ), and height ( H l ) are 0.05, 0.1, and 0.1 m, respectively. The rest of the domain is filled with air. Tables 1 and 2 show the properties of the particles and
DR, which ranges from 1.1 to 2.5.
M
fluids, respectively. In this study, the density ratio of particle density to liquid density is defined as
ED
In general, the Reynolds number in dam break flow is calculated as Re H gH f / f
PT
where H , g , f , and f are the height of initial liquid column, the gravity, the fluid density,
CE
and the fluid viscosity. According to this definition, the Reynolds number of the subject considered in this study is about 1.0 105. In terms of flow state, because the present study
AC
focused on dam break flow on the dry bed from standstill before hitting the opposite wall, fluids experiences developing and transient flow from the stationary. However in dam break flow on dry bed, the signs of strong turbulence did not exist generally [44]. It can be understood that laminar flow is predominant to this period before hitting the opposite wall. Therefore, the present study does not account for the turbulence phenomenon. However, even not for the dam break flow, there are a few studies for the effect of the particles on the turbulence flow. They reported that particles can suppress the turbulence intensity, which 12
ACCEPTED MANUSCRIPT can lead to tremendous acceleration of the flow due to effect of drag reduction [45, 46]. In addition, the drag reduction caused by particles in turbulence flow is substantial when the size of the suspended particles is subcritical, which shows the effect of particle type on the turbulence of the flow [47].
The particle Reynolds number is defined as Re p d p v f v p f / f where , d p , v f ,
CR IP T
v p , f , and f are the porosity, the particle diameter, the fluid velocity, the particle velocity,
the fluid density, and the fluid viscosity. The particle Reynolds number in this study is about 24 102, which is utilized to calculate drag coefficient of the particle in equations (7) and (8).
AN US
The grid size for computing the fluid motion is 0.005 m, which should be larger than the particle size to avoid numerical errors. Discrete particles would otherwise cause discontinuities in the cell of the fluid in the porosity calculation [48, 49]. The time steps for CFD and DEM are 5 10-4 and 1 10-5 s, respectively. A no-slip boundary condition is imposed on the vertical and
M
horizontal walls. The dimensionless variables used in this study are the dimensionless time ( t * ),
1 * 2 t t 2 g / L l x * x / L l l l * yl yl / Ll * x p x p / Ll
(19)
AC
CE
PT
ED
liquid front position ( xl* ), liquid height ( yl* ), and particle front position ( x*p ):
2.5 Validation of the numerical method
To validate the numerical methods, the time histories of dimensionless position of liquid front and liquid height were quantitatively compared with those of Sun and Sakai [13] who utilized numerical and experimental methods to study dam break flow containing particles under similar conditions we considered in this study, as shown in Fig. 4. Consequently, the comparison of the
13
ACCEPTED MANUSCRIPT present results with the experimental and numerical ones of Sun and Sakai [13] gives a good validation, which guarantees the reliability of the present numerical methods. Additionally, comparison in terms of the particle distribution and free-surface formation at each time was carried out with Sun and Sakai’s results [13], as shown in Fig. 5. The free-surface was indicated by blue line and the particle motions were individually captured with their velocity.
CR IP T
Sun and Sakai [13] considered the effects of a rising gate by adding an obstacle in front of the liquid column, which was not done in the present study. Although there are some discrepancies on the free surface in the early stage, the overall process of fluid flow and particle motion showed similar results. After collapse of liquid column, free-surface advances against the wall and the
AN US
particle bed expands over the bottom of the domain. In addition, the particle front moves more slowly than the liquid front at t * < 4.
M
3. Results and discussion
PT
ED
3.1 Effect of particle density on liquid front and particle front
To verify the effects of the particles on the free surface and particles, an additional simulation
CE
of a classical liquid-gas flow in a dam break was performed independently with the same conditions, and the results were compared. Fig. 6 presents typical three-dimensional perspective
AC
views and side views, which show the time evolution of the particles and the free surface between the liquid and the gas for the liquid-gas flow and liquid-gas-particle flows. In the case of the liquid-gas-particle flows, DR values of 1.1, 1.5, and 2.5 were representatively considered. To clarify the dependence of the front positions of the free surface, the temporal evolution of the freesurface profiles for the different flows on a two-dimensional (x-y) plane is plotted in Fig. 7. When the liquid column starts to collapse due to gravity, the front of the liquid and particles advances
14
ACCEPTED MANUSCRIPT against the wall. In this process, the liquid front behaves differently when the particles are present. In the initial stage of the dam breaking, the front positions of the liquid are almost identical, regardless of the presence of particles. However, over time, the front position of the liquid in the liquid-gas-particle flow lags behind that in the liquid-gas flow, as shown in Fig. 6 (a). At t * =2, the difference between the front position of the liquid according to presence of particles was greater
CR IP T
than that at t * =1. This tendency remained until t * =3, as shown Figs. 6 (b) and (c). Above results indicating the effect by particle can be clearly identified according to time in Fig. 7.
These observations similar to findings by Sun and Sakai [13] and Li and Zhao [50], who observed that the front position of the liquid in a liquid-gas-particle dam break flow seems to
AN US
move slower than that in a liquid-gas dam break flow.
We also identified the effects of DR on the front position of the liquid in liquid-gas-particle flows. Regardless of DR, the front position of the liquid is almost the same at each time, as shown
M
in Fig. 7. However, the front position of the particles was greatly influenced by DR. In the initial stage, the front position of the particles was almost the same, like that of the liquid. As shown in
ED
Fig. 6 (a), the font position of the particles at larger DR lagged behind that at smaller DR at t*=1.
PT
Over time, the difference of the front position of particles at different DR was clear, as shown in Figs. 6 (b) and (c).
CE
Consequently, the DR has no effect on the liquid front, but it makes the particle front move with different velocity. Based on these results, we can classify characteristics of the movement in
AC
terms of the front positions of the liquid and the particles. For DR=1.1, the particle front is always ahead of the particle front after t * =1. In the case of DR=1.5, the front positions of the liquid and particles are almost the same every time. However, the particle front for DR=2.5 lagged behind the liquid front after t * =1. Fig. 8 shows the dimensionless positions and velocity of the liquid front for the liquid-gas flow and liquid-gas-particle flows, for which DR values ranging from 1.1 to 2.5 were considered.
15
ACCEPTED MANUSCRIPT The dimensionless liquid front position was averaged in the spanwise direction, xl* , to observe the variation in this direction. As shown in Figs. 6 and 7, the liquid front in the liquid-gas flow generally moves faster than that in the liquid-gas-particle flows against the wall. Early on, the liquid front position is almost the same, regardless of the particles. However, the particles
liquid front hits the wall on the other side, as shown in Fig. 8 (a).
CR IP T
affect the liquid front position after t * =0.5. This difference can be observed clearly until the
The results were also compared with those for a liquid-gas flow from Sun and Sakai [13], Martin and Moyce [51], and Koshizuka and Oka [17], as shown in Fig. 8 (a). This was done to
AN US
emphasize the effects of the particles. The maximum velocity of the liquid front with the particles was identified at t * =2.58, and that without particles occurred at t * =2.5. Before the maximum velocity, the liquid front velocity without particles was always greater than that with particles. After the maximum velocity of the liquid front, the liquid front velocity decreases, regardless of
rapidly rather than that with particles
M
the presence of particles. Particularly, the liquid front velocity without particles decreases more
ED
The liquid front in the liquid-gas-particles flows was also compared at different DR. As shown
PT
in Fig. 8 (a), the liquid front positions for different DR generally have the same profile. This means that the effect of DR on the liquid front is almost negligible. As with the liquid front
CE
position, the time history of the liquid front velocity for different DR is almost the same. However, the particle front is indeed affected by DR. Fig. 9 shows top views of the free
AC
surface and particle formation for different DR. As expected, greater DR impedes the particle motion, and the particle velocity decreases as DR increases. Furthermore, instability of the particles was observed at DR=1.1. However, the particle formation along the bottom for DR=2.5 does not display significant instabilities and stays perpendicular to the direction of the flow. Fig. 10 shows the spanwise-averaged dimensionless front position of particles, x*p , and the streamwise velocity of the particle front, u p , for different particle densities. Regardless of DR,
16
ACCEPTED MANUSCRIPT the profile of the particle front position is composed of three sections. The first section has the smallest slope in particle front position in Fig. 10 (a) and corresponds to the front velocity at t * < 0.2, as shown in Fig. 10 (b). In this section, the front position is independent of DR. The second section has a sharp slope in particle front position in Fig. 10 (a) induced by the higher velocity than the first section. In this section, the front position depends on DR. As a result, the particle
CR IP T
front position decreases as the DR increases, as shown in Fig. 10 (a). This suggests that the front velocity decreases with increasing DR at t * > 0.2, as shown in Fig. 10 (b). Finally, before the particles front collide against the wall, the front velocity in near the wall drops dramatically and reaches zero.
AN US
Fig. 11 shows the moments when the liquid and particle fronts arrive at the opposite wall ( ta* ) for different DR. ta* of the liquid front for different DR is consistent, while the particle front takes more time to arrive at the opposite wall as DR increases. Furthermore, ta* of the particle
M
front turned out to be proportional to DR.
ED
The maximum velocities of the liquid front and particle front are plotted according to DR in Fig. 12. The maximum velocity of the liquid front is independent of DR, while that of the particle
PT
front decreases as DR increases, as shown in Fig. 12 (a). The difference between the maximum velocity of the liquid front and particle front, u*max , is also compared according to DR in Fig.
CE
12 (b). In the case of DR below 1.7, u*max is positive. For DR=1.7, u*max is almost zero,
AC
and after that u*max is negative as DR increases.
3.2 Motion regimes
As discussed in section 3.1. The front speed of the liquid and particles according to DR changes differently over time. Therefore, differences occur in the front position of the liquid and particles. Fig. 13 presents the time histories of the spanwise-averaged front position of the liquid
17
ACCEPTED MANUSCRIPT and particles for different DR. In the initial stage of the dam breaking, the front positions of the liquid and the particles are almost identical for different DR at t * 0.4. However, at t* > 0.4 and DR = 1.1, the particle front position is larger than the liquid front position at t* > 0.4. On the other hand, the particle front position at DR = 2.5 is smaller than the liquid front position at t* > 0.4. These results were attributed to the varying speed of the particle front while liquid front
overlapped position before impingement of the mixture flow.
CR IP T
speed hardly changed. At DR = 1.5, the front positions of the liquid and particle keep up the
Based on these results, three different regimes in the motion of the dam break flow with particles were identified. The regimes were simply classified based on the distance between
AN US
the particle front and the liquid front, which was strongly related to DR. The motion regimes identified are the advancing motion regime, overlapping motion regime, and delaying motion regime. The distance between the particle front and liquid front, Dp l x p xl , was
M
used to classify the regimes. The advancing motion regime occurs when the particle front is in
ED
front of that of the liquid ( 0.5d Dp l , where d is the particle diameter). The overlapping motion regime is in the region of 0.5d Dp l 0.5d , where the front positions of the liquid and
PT
particles are almost identical. The delaying motion regime occurs when the particle front lags
CE
behind the liquid front ( Dp l 0.5d ). Illustrations of the regimes at the same instance are shown in Fig. 14. The three motion
AC
regimes are greatly influenced by DR and the travel distance of dam break flow. The three intervals of DR can be demarcated in terms of the three motion regimes. For each computation, the travel distance was considered before the impingement of the mixture flow to confirm that the results present distinct flow regimes. The map of the motion regimes is shown in Fig. 15, where the x- and y-axis are the dimensionless travel distance and DR, respectively. In the region of 1.1 DR 1.4, the overlapping motion regime occurred in an early stage. After that, the advancing motion regime dominated the remaining domain. In addition, the area of the advancing regime 18
ACCEPTED MANUSCRIPT decreased with increasing DR. The series of changing regimes in the ranger of 1.1 DR 1.4 was named the overlapping-advancing (O-A) regime. At 1.5 DR 1.6, a transition regime from the advancing to the delaying motion regime was observed. In this stage, the delaying motion regime first appears at DR = 1.6. In addition, the overlapping motion regime was more dominant than the other regimes. At 1.6 DR 2.5, the
CR IP T
liquid and particle fronts moved in an overlapping-delaying (O-D) regime. Early on, the overlapping motion regime was observed, similar to the O-A regime. However, the delaying motion regime dominated the remaining domain. In addition, the area of the delaying motion regime expanded as DR increased. The map thus showed an O-A regime, transition regime, and
AN US
O-D regime. The map could be used to predict which material (liquid or particles) will arrive more quickly at a target point.
M
4. Conclusion
ED
This study investigated the liquid-gas-particles mixture flow during dam break. The CFDDEM numerical method was used for the simulation. Such numerical method consist of the open
PT
source codes in OpenFOAM for CFD and CFDEMproject, which includes LIGGGHTS for DEM
CE
simulation and calculates the interaction data between CFD and DEM. The computed results are in good agreement with those from previous studies. The main purpose is to investigate the effects
AC
of particle density on the behavior of the particles and the free surface. Basically, if the particles were added in dam break flow, they interfered with the movement of
the liquid front. Therefore, the liquid front in a liquid-gas-particle dam break flow seems to move slower than that in a liquid-gas dam break flow. In addition, the liquid front in the liquid-gasparticles flows was also compared at different particle density. The time history of liquid front positions for different DR generally have the same profile. Which means that the effect of particle
19
ACCEPTED MANUSCRIPT density on the liquid front is almost negligible. However, as the particle density grows, particle motion was impeded, and the particle velocity decreased. Consequently, the DR has no effect on the liquid front, but it makes the particle front move with different velocity. Furthermore, instability of the particles formation was observed as particle density was the smallest. However, the particle formation along the bottom for the biggest particle density does not display
CR IP T
significant instabilities and stays perpendicular to the direction of the flow. The time when the liquid front reach at the opposite wall is consistent according to particle density. However, such time for particle front decrease as particle density increases. Furthermore, the time of the particle front turned out to be proportional to particle density.
AN US
Based on these results, we can classify characteristics of the movement in terms of the front positions of the liquid and the particles. The characteristic motion was classified into an advancing, overlapping, and delaying motion regime. In the initial stage of the liquid collapse, the
M
overlapping motion regime governs the front positions of the liquid and the particles, regardless of the particle density. This results in almost identical front positions of the liquid and the particles.
ED
In the advancing motion regime, the particle front moves more quickly than the liquid front does. For this reason, the particles front is ahead of the liquid front after penetration of the particless
PT
across the free surface. In the delaying motion regime, the particle front lags behind the liquid
CE
front. The three regimes are greatly influenced by particle density and the travel distance of the dam break flow. The map of the motion regimes was created as afunction of DR and streamwise
AC
position. The map can be used to predict which phase will arrive more quickly at a target point.
Acknowledgment
This work was supported by National Research Foundation of Korea (NRF) grant founded by the Korea government (MSIP) through GCRC-SOP (No. 2011–0030013) and (NRF2015R1D1A3A01020867).
20
ACCEPTED MANUSCRIPT References [1] Iverson RM. The physics of debris flows. Rev Geophs 1997; 35(3):245-96. [2] Highland LM, Ellen SD, Christian SB, Brown III WM. Debris-flow hazards in the United States. https://pubs.usgs.gov/fs/fs-176-97/fs-176-97.html; 1997
CR IP T
[3] Griffiths RW. The dynamics of lava flows. Annu Rev Fluid Mech 2000;32(1):477-518. [4] Roche O, Montserrat S, Ni o Y, Tamburrino A. Experimental observations of water-like behavior of initially fluidized, dam break granular flows and their relevance for the propagation of ash-rich pyroclastic flows. J Geophys Res-Solid Earth 2008;113(B12).
AN US
[5] Ancey C. Plasticity and geophysical flow: A review. J Non-Newton Fluid 2007;142(1):435.
[6] Balmforth NJ, Kerswell RR. Granular collapse in two dimensions. J Fluid Mech 2005;538:399-428.
M
[7] Bagnold RA. Experiments on a gravity-free dispersion of large solid spheres in a
ED
Newtonian fluid under shear. P Roy Soc A-Math Phy 1954;225(1160):49-63. [8] Bakhtyar R, Ghaheri A, Yeganeh-Bakhtiary A, Barry DA. Process-based model for
PT
nearshore hydrodynamics, sediment transport and morphological evolution in the surf and swash zones. Appl Ocean Res 2009;31(1):44-56.
CE
[9] Tang HS, Keen TR, Khanbilvardi R. A model-coupling framework for nearshore waves,
AC
currents, sediment transport, and seabed morphology. Commun Nonlinear Sci 2009;14(7):2935-47.
[10] Wang JY, Tang HS, Fang HW. A fully coupled method for simulation of wave-currentseabed systems. Commun Nonlinear Sci 2013;18(7):1694-709. [11] De Rooij F. Sedimenting particle-laden flows in confined geometries. Diss. University of Cambridge; 2000.
21
ACCEPTED MANUSCRIPT [12] Guo L, Morita K. Numerical simulation of 3D sloshing in a liquid–solid mixture using particle methods. Int J Numer Meth Eng 2013;95(9):771-90. [13] Sun X, Sakai M. Three-dimensional simulation of gas–solid–liquid flows using the DEM–VOF method. Chem Eng Sci 2015;134:531-48. [14] Shigematsu T, Liu PLF, Oda K. Numerical modeling of the initial stages of dam-break
CR IP T
waves. J Hydraul Res 2004;42(2):183-95. [15] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39(1):201-25.
[16] Harlow FH, Welch JE. Numerical calculation of time-dependent viscous incompressible
AN US
flow of fluid with free surface. Phys Fluids 1965;8(12):2182-9.
[17] Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng 1996;123(3):421-34.
M
[18] Shao S, Lo EY. Incompressible SPH method for simulating Newtonian and nonNewtonian flows with a free surface. Adv Water Resour 2003;26(7):787-800.
ED
[19] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J Comput Phys 2003;191(2):448-75.
PT
[20] Monaghan JJ, Rafiee A. A simple SPH algorithm for multi-fluid flow with high density
CE
ratio. Int J Numer Meth Fl 2013;71(5):537-61. [21] Abdolmaleki K, Thiagarajan KP, Morris-Thomas MT. Simulation of the dam break
AC
problem and impact flows using a Navier-Stokes solver. In: 15th Australasian Fluid Mechanics Conference The University of Sydney; 2004.
[22] Lobovský L, Botia-Vera E, Castellana F, Mas-Soler J, Souto-Iglesias A. Experimental investigation of dynamic pressure loads during dam break. J Fluid Struct 2014;48:407-34. [23] Kleefsman KMT, Fekken G, Veldman AEP, Iwanowski B, Buchner B. A volume-of-fluid based simulation method for wave impact problems. J Comput Phy 2005;206(1):363-93.
22
ACCEPTED MANUSCRIPT [24] Ozmen-Cagatay H, Kocaman S. Dam-break flow in the presence of obstacle: experiment and CFD simulation. Eng Appl Comp Fluid 2011;5(4):541-52. [25] Nsom B, Debiane K, Piau JM. Bed slope effect on the dam break problem. J Hydraul Res 2000;38(6):459-64. [26] Park IR, Kim KS, Kim J, Van SH. Numerical investigation of the effects of turbulence
CR IP T
intensity on dam-break flows. Ocean Eng 2012;42:176-87. [27] Jiang CB, Chen J, Tang HS, Cheng YZ. Hydrodynamic processes on beach: Wave breaking, up-rush, and backwash. Commun Nonlinear Sci 2011;16(8):3126-39.
Sci 2012;17(12):5273-85
AN US
[28] Zamankhan P. Rollers in low-head dams – challenges and solutions. Commun Nonlinear
[29] Cundall PA, Strack OD. A discrete numerical model for granular assemblies. Geotechnique 1979;29(1):47-65.
M
[30] Stevens AB, Hrenya CM. Comparison of soft-sphere models to measurements of collision properties during normal impacts. Powder Technol 2005;154(2):99-109.
ED
[31] Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol, 1993;77(1):79-87.
PT
[32] Lin SY, Chen YC. A pressure correction-volume of fluid method for simulations of
CE
fluid–particle interaction and impact problems. Int J Multiphas Flow 2013;49:31-48. [33] Guo L, Morita K, Tobita Y. Numerical simulations of gas-liquid-particle three-phase
AC
flows using a hybrid method. J Nucl Sci Technol 2016;53(2):271-80.
[34] Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J 1977;82(12):1013-24.
[35] Yamada Y, Sakai M, Mizutani S, Koshizuka S, Oochi M, Murozono K. Numerical simulation of three-dimensional free-surface flows with explicit moving particle simulation method. J Atom Energ Soc Jpn 2011;10(3):185–93.
23
ACCEPTED MANUSCRIPT [36] Sun X, Sakai M, Yamada Y. Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method. J Comput Phys 2013;248:147-76. [37] Delannay R, Valance A, Mangeney A, Roche O, Richard P. Granular and particle-laden flows: from laboratory experiments to field observations. J Phys D Appl Phys 2017;50(5):053001. Ltd.
OpenFOAM
-
The
Open
Source
CFD
Toolbox.
CR IP T
[38] OpenCFD
http://www.openfoam.com; 2009.
[39] Kloss C, Goniva C, Hager A, Amberger S, Pirker S. Models, algorithms and validation for opensource DEM and CFD–DEM. Prog Comput Fluid Dy, an International Journal
AN US
2012; 12(2-3):140–52.
[40] Goniva C, Kloss C. Hager A, Pirker S. An open source CFD-DEM perspective. In: Proceedings of OpenFOAM Workshop; 2010. p. 22-4.
M
[41] Gidaspow D. Multiphase flow and fluidization-continuum and kinetic theory descriptions. Academic Press; 1994. ISBN 0122824709.
ED
[42] Goniva C, Gruber K, Kloss C. Sediment erosion a numerical and experimental study. River Flow 2012; 2012. p. 415-22. ISBN 9780415621298.
PT
[43] Issa, R. Solution of the implicitly discretized fluid flow equations by operator- splitting. J
CE
Comput Phys 1986;62(1):40-65. [44] Hanosi IM, Jan D, Szabo KG, Tel T. Turbulent drag reduction in dam-break flows. Exp
AC
Fluids 2004;37(2):219-29.
[45] Barenblatt GI. Scaling, self-similarity, and intermediate asymptotics: dimensional analysis and intermediate asymptotics. Cambridge Press; 1996. ISBN 0521435226
[46] Barenblatt GI, Chorin AJ, Prostokishin VM. A note concerning the Lighthill “sandwich model” of tropical cyclones. P Natl Acad Sci USA 2005;102(32):11148-50.
24
ACCEPTED MANUSCRIPT [47] Bersch M, Hulshof J, Prostokishin VM. Flow laminarization and acceleration by suspended particles. SIAM J Appl Math 2015;75(4):1852-83. [48] Kawaguchi T, Sakamoto M, Tanaka T, Tsuji Y. Quasi-three-dimensional numerical simulation of spouted beds in cylinder. Powder Technol 2000;109(1):3–12. [49] Link JM, Cuypers LA, Deen NG, Kuipers JAM. Flow regimes in a spout–fluid bed: A
CR IP T
combined experimental and simulation study. Chem Eng Sci 2005;60(13):3425–42 [50] Li X, Zhao J. Numerical simulation of dam break by a coupled CFD-DEM approach. In: The 15th Asian Regional Conference on Soil Mechanics and Geotechnical Engineering; 2016;2(18):691-6.
AN US
[51] Martin JC, Moyce WJ. An experimental study of the collapse of liquid columns on a rigid
AC
CE
PT
ED
M
horizontal plane. Philos Trans R Soc Lond A 1952;244(882):312–24.
25
ACCEPTED MANUSCRIPT
LIST OF TABLES TABLE 1. Properties of particles TABLE 2. Properties of fluids
CR IP T
Table 1 Properties of particles Type
Value
Diameter, d (m)
0.0027 1100,1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2,500 3,883
The number of particles Spring constant (N/m) Restitution coefficient
1,000 0.9 0.3
ED
M
Friction coefficient
AN US
Density, ρp (kg/m3)
PT
Table 2 Properties of fluids
Value
Liquid density, ρl (kg/m3)
1000
AC
CE
Type
Liquid viscosity (Pa·s)
10-3
Gas density, ρg (kg/m3)
1
Gas viscosity (Pa·s)
10-5
26
ACCEPTED MANUSCRIPT
LIST OF FIGURES FIGURE 1. Illustration of force acting on particles in DEM. FIGURE 2. The algorithm of CFD-DEM.
CR IP T
FIGURE 3. Schematic of the three-dimensional computational domain. FIGURE 4. Comparison of dimensionless positions of (a) liquid front and (b) liquid height with the present results and Sun and Sakai’s results [13] of computation (Com) and experiment (Exp).
AN US
FIGURE 5. Comparison between Sun and Sakai’s results [13] (left column) and the present results (right column) for the particle distribution and free-surface formation at (a) t*=2 and (b) t*=4.
FIGURE 6. Illustrations of dam break flow without particles (1st column) and with particles
M
for DR= 1.1 (2nd column), 1.5 (3rd column), and 2.5 (4th column) at (a) t*=1.0, (b)
ED
t*=2.0, and (c) t*=3.0. The white surface and line indicate the free surface. FIGURE 7. Free-surface profile for dam break according to presence of particle at t*= 0.5, 1.0,
PT
2.0, and 3.0.
CE
FIGURE 8. Time histories of (a) dimensionless position of liquid front with data of liquid-gas flow from Martin and Moyce [40], Koshizuka and Oka [17], and Sun and Sakai
AC
[13]; (b) dimensionless streamwise velocity of liquid front for different DR.
FIGURE 9. Free-surface and particle formation for different DR at t*=3. FIGURE 10. Time histories of (a) dimensionless position and (b) dimensionless streamwise velocity of particle front for different DR. FIGURE 11. Arrival time (ta*) at the opposite wall versus DR. FIGURE 12. (a) Maximum front velocity of liquid and particles and (b) difference of maximum velocity between liquid and particles according to DR. 27
ACCEPTED MANUSCRIPT FIGURE 13. Comparison of dimensionless front position of liquid and particles. FIGURE 14. The front position lines of typical motion regimes for (a) advancing motion regime, (b) overlapping motion regime, and (c) delaying motion regime (top view). The blue and orange lines represent the liquid front and particle front, respectively.
CR IP T
FIGURE 15. Classification of three motion regimes as a function of DR and dimensionless
AC
CE
PT
ED
M
AN US
streamwise position (x*).
28
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 1. Illustration of force acting on particles in DEM.
29
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 2. The algorithm of CFD-DEM.
30
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
Fig. 3. Schematic of the three-dimensional computational domain.
31
CE
PT
ED
M
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
(b)
AC
Fig. 4. Comparison of dimensionless positions of (a) liquid front and (b) liquid height with the present results and Sun and Sakai’s results [13] of computation (Com) and experiment (Exp).
32
ACCEPTED MANUSCRIPT
AN US
CR IP T
(a)
(b)
Fig. 5. Comparison between Sun and Sakai’s results [13] (left column) and the present results
AC
CE
PT
ED
M
(right column) for the particle distribution and free-surface formation at (a) t*=2 and (b) t*=4.
33
M
AN US
(a) t*= 1.0
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
(b) t*= 2.0
(c) t*= 3.0
Fig. 6. Illustrations of dam break flow without particles (1st column) and with particles for DR= 1.1 (2nd column), 1.5 (3rd column), and 2.5 (4th column) at (a) t*=1.0, (b) t*=2.0, and (c) t*=3.0. The yellow surface and line indicate the free surface.
34
CR IP T
ACCEPTED MANUSCRIPT
AN US
Fig. 7. Free-surface profile for dam break according to presence of particle at t*= 0.5, 1.0, 2.0,
AC
CE
PT
ED
M
and 3.0.
35
CR IP T
ACCEPTED MANUSCRIPT
CE
PT
ED
M
AN US
(a)
(b)
AC
Fig. 8. Time histories of (a) dimensionless position of liquid front with data of liquid-gas flow from Martin and Moyce [40], Koshizuka and Oka [17], and Sun and Sakai [13]; (b) dimensionless streamwise velocity of liquid front for different DR.
36
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
Fig. 9. Free-surface and particle formation for different DR at t*=3.
37
CR IP T
ACCEPTED MANUSCRIPT
CE
PT
ED
M
AN US
(a)
(b)
AC
Fig. 10. Time histories of (a) dimensionless position and (b) dimensionless streamwise velocity of particle front for different DR.
38
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 11. Arrival time (ta*) at the opposite wall versus DR.
39
CR IP T
ACCEPTED MANUSCRIPT
PT
(b)
ED
M
AN US
(a)
CE
Fig. 12. (a) Maximum front velocity of liquid and particles and (b) difference of maximum
AC
velocity between liquid and particles according to DR.
40
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 13. Comparison of dimensionless front position of liquid and particles.
41
ACCEPTED MANUSCRIPT
AN US
CR IP T
(a)
PT
(c)
ED
M
(b)
Fig. 14. The front position lines of typical motion regimes for (a) advancing motion regime,
CE
(b) overlapping motion regime, and (c) delaying motion regime (top view). The blue and
AC
orange lines represent the liquid front and particle front, respectively.
42
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 15. Classification of three regimes as a function of DR and dimensionless streamwise
AC
CE
PT
ED
position (x*).
43
CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break Kyung Min Park , Hyun Sik Yoon , Min Il Kim PII: DOI: Reference:
S1007-5704(17)30383-0 10.1016/j.cnsns.2017.11.010 CNSNS 4367
To appear in:
Communications in Nonlinear Science and Numerical Simulation
Received date: Revised date: Accepted date:
25 April 2017 13 September 2017 9 November 2017
Please cite this article as: Kyung Min Park , Hyun Sik Yoon , Min Il Kim , CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break, Communications in Nonlinear Science and Numerical Simulation (2017), doi: 10.1016/j.cnsns.2017.11.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
Complex multiphase flow of liquid-gas-particle mixture in a dam-break is numerically investigated according to particle density.
CFD-DEM coupled method is successively adopted to find the natural phenomena involved in two-phase fluids and particle interaction.
The head speeds of the free-surface and particles show different patterns to particle density. The response of the free-surface and particles to particle density is identified
CR IP T
AC
CE
PT
ED
M
AN US
by three motion regimes of the advancing, overlapping and delaying motions.
1
ACCEPTED MANUSCRIPT
CFD-DEM based numerical simulation of liquid-gas-particle mixture flow in dam break Kyung Min Parka, Hyun Sik Yoonb*, Min Il Kimb, a
Global Core Research Center for Ships and Offshore Plants, Pusan National University, 2, Busandaehak-ro 63beon-gil, Gumjeong-Gu, Busan, 46241, Republic of Korea
b
CR IP T
Department of Naval Architecture and Ocean Engineering, Pusan National University, 2, Busandaehak-ro 63beon-gil, Gumjeong-Gu, Busan, 46241, Republic of Korea
ABSTRACT
AN US
This study investigates the multiphase flow of a liquid-gas-particle mixture in dam break. The open source codes, OpenFOAM and CFDEMproject, were used to reproduce the multiphase flow. The results of the present study are compared with those of previous results obtained by numerical and experimental methods, which guarantees validity of present numerical method
M
to handle the multiphase flow. The particle density ranging from 1100 to 2500 kg/m3 is considered to investigate the effect of the particle density on the behavior of the free-surface
ED
and the particles. The particle density has no effect on the liquid front, but it makes the particle
PT
front move with different velocity. The time when the liquid front reach at the opposite wall is independent of particle density. However, such time for particle front decrease as particle density
CE
increases, which turned out to be proportional to particle density. Based on these results, we classified characteristics of the movement by the front positions of the liquid and the particles.
AC
Eventually, the response of the free-surface and particles to particle density is identified by three motion regimes of the advancing, overlapping and delaying motions.
Keywords: DEM, CFD, Liquid-gas-particle mixture flow, Dam break, Three motion regimes
*
Corresponding author, Tel.: +82-51-510-3685, E-mail: [email protected] 2
ACCEPTED MANUSCRIPT
1. Introduction
The study of multiphase flows such as liquid-gas-particle mixture flow is challenging issue because the interactions between different phases result in complex behavior. Nevertheless, it is important to study multiphase flows because they are a common phenomenon in various areas in
CR IP T
academics, industry, and real life.
In geology, the multiphase flows on mountains are considered an important subject, where such flows are known as debris flows, which include mudslides, mudflow, and debris avalanches.
AN US
The compositions of a debris flow are mainly water, rock, and air, and the forces of particles and fluids greatly influence the flow [1]. Hilly areas and intense rainfall are particularly susceptible to debris flows. These flows are dangerous to life and property because of the high sediment concentrations and mobility. Many properties were inundated or crushed by the debris flows,
M
pastures and cornfields were buried, and livestock perished [2].
During volcanic eruptions, mixtures of lava, gas, and ash generate hot multiphase flows. There
ED
are two types of such flows, and the first type is lava flow. Lava flow occurs when molten or
PT
partially molten rock erupt from the surface of the earth. The melt spreads on the surface as a gravity current and forms a lava flow. The lava flow is not a simple liquid but generally a mixture
CE
that contains magma, rock, ash, crystals, and gas bubbles. Lava flows play important roles in the development of planetary crusts and ore deposits, and they pose dangers to people and property,
AC
which means the importance of research on such flows [3]. Another type of multiphase flow in volcanic eruptions is pyroclastic flow. When the density of
the ash particles is higher than the gas density, the eruption column may collapse and form a pyroclastic multiphase flow. A pyroclastic flow can propagate distances of several kilometers from its starting point [4]. Such flows therefore also pose dangers to populations and infrastructure [5-7], and prediction of their behavior is crucial for hazard assessment.
3
ACCEPTED MANUSCRIPT Also, in ocean the complex multiphase flow occurs. Littoral sediment transport is caused by the interactions of waves, currents, and tides. These phenomena determine the rate of beach erosion or accretion and act as boundary conditions for nearshore hydrodynamic and morph dynamic models. The study of nearshore multiphase flows is important for coastal engineering, environmental protection, recreation, and military operations. Some models have therefore been
CR IP T
developed to simulate and predict their behaviors [8-10]. Offshore multiphase flows play an important role in turbidity currents. Turbidity current occurs when the edge of the continental shelf fails and causes a submarine landslide. As the turbidity current propagates down the slope, it erodes the ocean floor and increases the amount of
AN US
suspended sediment. When the current reaches the abyssal plain, deposition becomes dominant and another sediment layer is added to the geological record of the ocean floor. The sediment layers deposited by sequential turbidity currents can form traps for oil. Sound understanding of
M
the dynamics of turbidity currents is instrumental for the reconstruction of the sediment layer thickness and potential location of oil reservoirs from drilling cores and seismological data [11].
ED
Multiphase flow also occurs in liquid-metal-cooled reactors such as mixture flow of liquid metal and resolidified fuel particles. The dynamics of the flow strongly affects the risk of
CE
phase [12].
PT
accidents, which makes it important to analyze the movements and interactions between each
As illustrated these relevant fields governed by the liquid-gas-particles mixture flow, the
AC
multiphase flows appear in many areas. However, these problems are very complex to be numerical simulated, since these problems are also considered as one of the challenging problems in the CFD and also simulation community. Therefore, as the simplified problem but involving essential elements of these problems, a dam break flow containing particles can be considered as reasonable problem to investigate the multiphase flows such as liquid-gas-particle mixture flows. In general, the classical dam breaks have been well studied in recent decades [13]. For
4
ACCEPTED MANUSCRIPT example, the classical dam break flow problem has been widely studied as a typical liquid-gas mixture flow problem to identify the interactions of each phase. Studies have examined the development of numerical methods [14-20], the effects of the density ratio of the fluids [19, 20], the impact on structures [21-23], the presence of obstacles [24], the effects of the bed slope [25] and turbulence intensity [26], and the physics of a two-phase flow [27, 28].
CR IP T
However, dam break flows containing particles have not been thoroughly discussed, there is a lack of numerical simulations and experiments dealing with flows of liquid-gas-particle mixtures in the dam break problem. Furthermore, the main focus of previous studies was to validate numerical methods by comparing simulation and experimental results rather than interpreting the
AN US
physics of the flow.
Such studies developed various numerical models in the Eulerian or Lagrangian framework to describe the motion of two-phase flows separated by an immiscible fluid interface. In terms of,
M
the particle motion was captured by the discrete element method (DEM), which has been widely applied in particle collision studies. The DEM model was essentially used to calculate the
ED
collision forces between particles or between particles and walls [29], and some researchers used the method to model particle motion [30, 31]
PT
In the Eulerian framework, the fluid flow is solved with a computational mesh by
CE
computational fluid dynamics (CFD), and the free surface of fluids was captured mainly by the volume of fraction (VOF) method [15]. For example, Lin and Chen [32] investigate the
AC
discrepancy of the particulate flow behavior and the impact of the particulate flow on the structure according to the number of particles using CFD-DEM. Guo et al. [33] reasonably simulated the transient behavior of the multiphase flow of a liquid-gas-particle mixture flow by CFD-DEM. The method was validated based on experiments of a water-particle dam break in tems of gas-liquidparticle flows. Sun and Sakai [13] developed CFD-DEM to simulate various three-dimensional flows of liquid-gas-particle mixtures. The numerical results were validated through comparison
5
ACCEPTED MANUSCRIPT with an experiment. In contrast, fluid flow in the Lagrangian framework considers particle motion based on smoothed-particle hydrodynamics (SPH) [34] and the moving particle semi-implicit (MPS) method [35]. For example, Sun et al. [36] used SPH-DEM to simulate liquid-particle flows such as those in a dam break and in a rotating tank. Guo and Morita [12] investigated the sloshing in a
and they ignored the effects of gas on the liquid and particles.
CR IP T
liquid-particle mixture flow using the MPS-DEM. However, they did not consider the gas phase
None of these studies provide a physical interpretation for various conditions, in terms of the fluid and particle properties and the liquid column. Especially, the multiphase flows such as debris
AN US
flow, turbidity flow, snow avalances, and pyroclastic flow have different particle densities ranging to 1000-3000 kg/m3 [37]. Therefore, the present study investigates the effects of the density of particles on their behaviors and the free surface in the dam break problem.
M
The CFD-DEM method was used for the simulation with the open source code in OpenFOAM for CFD [38] and CFDEMproject, which includes LIGGGHTS for DEM simulation
ED
and calculates the interaction data between CFD and DEM [39, 40]. The numerical method is effective for investigating the interaction of the different phases. A dam break flow with particles
PT
has intrinsic complexity and instability because of the interfaces between the phases. However, a
CE
few recent studies have reasonably simulated such flows using CFD-DEM [13, 33]. This paper is structured as follows. Section 2 presents the governing equations for the motion
AC
of each phase, the numerical scheme, the problem definition, and a validation case for the numerical method. The behavior of particles and the free surface between the liquid and gas at different particle densities was compared to identify the effects. Section 3 presents the different regimes of the particle and free-surface motion as a function of the particle/liquid density ratio and position.
6
ACCEPTED MANUSCRIPT
2. Numerical methodology and problem definition
Models to describe a dam break flow with particles can be divided into Eulerian models (in OpenFOAM) and Lagrangian models (in LIGGGHTS). Eulerian models calculate the fluid flow
CR IP T
based on CFD while capturing the liquid-gas interface provided by VOF. In the Lagrangian models, DEM calculates the path and motion of each particle. Particle interactions are modeled through specified collision rules. The detailed governing equations for each phase are described in
2.1 Governing equations of particle motion
AN US
the following.
M
In DEM, a particle has translational and rotational motion that depend on the forces and torques acting on it [29]. The forces and torques result from particle interactions with other
ED
particles, the walls, and the surrounding fluids. As shown in Fig. 1, there are three type forces: contact force, gravity force, and fluid force.
PT
The contact forces are considered based on a spring, dashpot, and friction-slider model. For a
CE
particle i with mass mi and velocity vi is in contact with another particle or geometric
AC
element j . The translation motion equation is expressed by the following: mi
dvi Fn, ij mi g Ffluid ( Fd Fp ) dt j
(1)
where Fn ,ij , g, and Ffluid are the contact force in the normal direction, gravity, and the fluid force due to the presence of fluid. The fluid force is a combination of the fluid drag force ( Fd
and
fluid pressure gradient force ( Fp ). The fluid force is also considered as fluid-particle interaction force based on Newton’s third law.
7
ACCEPTED MANUSCRIPT In Eq. (1), the normal-direction contact forces are calculated using a linear spring contact model:
Fn, ij K n dn Nn vn, ij
(2)
where K n , d n , N n , and v n ,ij are the normal spring stiffness, the particle overlap in the normal
CR IP T
direction at the contact point, the normal damping, and the normal velocity components of the relative sphere-surface velocity at the contact point. N n is defined by the following: N n 2 N n damp K n meq
(3)
where meq and N n damp are the equivalent particle mass and the normal damping coefficient: 1 1 1 mi m j
AN US
meq
ln(Cn rest )
2 ln(Cn rest )2
(5)
M
N n damp
(4)
coefficient of restitution.
ED
where mi , m j , and Cn rest are the mass of particle i , that of particle j , and the normal
PT
The drag force on the individual particles was calculated by the Gidaspow empirical model [41], which takes into account the effects of the surrounding particles. The drag force is defined
Fd 0.5Cd f A v f vi v f vi
(6)
AC
CE
by the following:
where Cd , f , A , and v f are the drag coefficient, fluid density, projected area of particle i , and fluid velocity. The drag coefficient is defined as: 4 1 CD 150 1.75 if min 3 Re p
Otherwise:
8
(7)
ACCEPTED MANUSCRIPT
CD
24 3.6 Re
0.687 p
Re p
3.65
(8)
where is the porosity, vmin is a cutoff porosity of 0.8, and Re p is the particle Reynolds number.
FP Vi pstatic
(9)
CR IP T
Eq. (9) represents the pressure gradient force, where Vi and pstatic are the volume of the particle i and the gradient of the static pressure in the fluid, respectively.
The particles also undergo rotational motion, which is governed by Eq. (10):
di R i Ft ,ij dt
AN US
Ii
(10)
where I i , i , Ri , and Ft ,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
M
tangential force is defined by the following:
ED
Ft ,ij Kt dt Nt vt ,ij
K n d n c fs dt dt
(12)
PT
Kt dt K n d n C fs ; Ft
(11)
where K t , d t , N t , vt , and C fs are the tangential spring stiffness, the overlap in the tangential
CE
direction at the contact point, the tangential damping, the tangential velocity components of the
AC
relative sphere-surface velocity at the contact point, and the static friction coefficient. N t is defined by the following: Nt 2 Nt damp Kt meq
(13)
where meq and Nt damp are the equivalent particle mass and the tangential damping coefficient defined by the following:
9
ACCEPTED MANUSCRIPT In Ct rest
Nt damp
2 In Ct rest
(14)
2
where Ct rest is the tangential coefficient of restitution.
CR IP T
2.2 Governing equations of fluid flow
The three-dimensional two-phase flow is governed by the continuity equation in Eq. (15) and the Navier–Stokes equation in Eq. (16):
AN US
ui 0 xi
p ui uiu j t x j xi x j
u u j i gi Sm x j xi
(15)
(16)
M
where x i is a Cartesian coordinate, ui is the corresponding velocity component, p is the
ED
pressure, is the density, is the viscosity, is the porosity, and g is gravity. The source term S m takes into account the interaction force between the fluid and particles by two-way
PT
coupling. Based on Newton’s third law, S m in a fluid cell is determined by adding up the fluid
Sm
DEM iteration Particle F Vcell
fluid
(17)
AC
CE
force of the particles located in the fluid cell:
where Vcell is the volume of the fluid cell. The VOF method was used to capture the free surface between the liquid and gas. This
method was designed for immiscible fluids, and each cell has a volume fraction ( Q ) corresponding to each fluid. A cell with a scalar Q value 0 indicates air, and a value of 1 represents liquid. The free surface in each cell is defined when the Q value is between 0 and 1.
10
ACCEPTED MANUSCRIPT At each time step, a transport equation is solved to find the distribution of the fluid [15]. In conventional VOF method, the porosity was additionally considered to represent the fluid displaced by particles [42]. The porosity is a fraction of void volume over the total volume such as
1
V
particle
Vcell
, where Vparticle is particle volume corresponding fluid cell. The distribution is
Q Q ui 0 t xi
CR IP T
found by solving the following transport equation of Q : (18)
The second-order upwind scheme and the second-order central differencing scheme were used
AN US
for the convection and diffusion terms, respectively. To the calculate the unsteady flow, the time derivative terms are discretized using a second-order accurate backward implicit scheme. The velocity–pressure coupling and overall solution procedure are based on a PISO (Pressure Implicit
ED
2.3 Algorithm of CFD-DEM scheme
M
with Splitting of Operators) algorithm [43].
PT
The overall algorithm of the CFD-DEM method is illustrated in Fig. 2. The calculation of
CE
CFD-DEM starts with DEM. The position and velocity of each particle are calculated by solving the governing equations of particle motion (1) and (10) based on the contact force and fluid force
AC
in DEM. The porosity in each cell is then calculated. Simultaneously, the fluid-particle interaction force is accumulated from the particles and transferred to the fluid phase. To calculate the governing equation of fluid flow in CFD, the VOF equation (18) is first calculated. Through this process, the volume fraction Q is obtained, which updates density and viscosity for each fluid in each fluid cell. The governing equations of the fluid flow (15) and (16) are then solved iteratively based on the PISO algorithm [43]. When the iteration converges, a new fluid velocity and pressure are obtained. After that, the fluid forces acting on each particle are calculated from 11
ACCEPTED MANUSCRIPT the results of CFD and are passed to the DEM solver using the CFD-DEM scheme. The time step in the DEM solver tends to be more limited than that in the CFD solver (it is 50 times smaller than that of the CFD solver in this study). This occurs because the total stability conditions are dominated by DEM [13]. This condition satisfies the qualification for particle
CR IP T
collision modeling in an extreme situation and ensures overall computational efficiency.
2.4 Computational conditions
The computational domain is shown in Fig. 3, where the length, breadth, and height are 0.2,
AN US
0.1, and 0.3 m respectively. The domain contains a liquid column and 3,883 particles. The liquid column’s initial length ( Ll ), breadth ( Bl ), and height ( H l ) are 0.05, 0.1, and 0.1 m, respectively. The rest of the domain is filled with air. Tables 1 and 2 show the properties of the particles and
DR, which ranges from 1.1 to 2.5.
M
fluids, respectively. In this study, the density ratio of particle density to liquid density is defined as
ED
In general, the Reynolds number in dam break flow is calculated as Re H gH f / f
PT
where H , g , f , and f are the height of initial liquid column, the gravity, the fluid density,
CE
and the fluid viscosity. According to this definition, the Reynolds number of the subject considered in this study is about 1.0 105. In terms of flow state, because the present study
AC
focused on dam break flow on the dry bed from standstill before hitting the opposite wall, fluids experiences developing and transient flow from the stationary. However in dam break flow on dry bed, the signs of strong turbulence did not exist generally [44]. It can be understood that laminar flow is predominant to this period before hitting the opposite wall. Therefore, the present study does not account for the turbulence phenomenon. However, even not for the dam break flow, there are a few studies for the effect of the particles on the turbulence flow. They reported that particles can suppress the turbulence intensity, which 12
ACCEPTED MANUSCRIPT can lead to tremendous acceleration of the flow due to effect of drag reduction [45, 46]. In addition, the drag reduction caused by particles in turbulence flow is substantial when the size of the suspended particles is subcritical, which shows the effect of particle type on the turbulence of the flow [47].
The particle Reynolds number is defined as Re p d p v f v p f / f where , d p , v f ,
CR IP T
v p , f , and f are the porosity, the particle diameter, the fluid velocity, the particle velocity,
the fluid density, and the fluid viscosity. The particle Reynolds number in this study is about 24 102, which is utilized to calculate drag coefficient of the particle in equations (7) and (8).
AN US
The grid size for computing the fluid motion is 0.005 m, which should be larger than the particle size to avoid numerical errors. Discrete particles would otherwise cause discontinuities in the cell of the fluid in the porosity calculation [48, 49]. The time steps for CFD and DEM are 5 10-4 and 1 10-5 s, respectively. A no-slip boundary condition is imposed on the vertical and
M
horizontal walls. The dimensionless variables used in this study are the dimensionless time ( t * ),
1 * 2 t t 2 g / L l x * x / L l l l * yl yl / Ll * x p x p / Ll
(19)
AC
CE
PT
ED
liquid front position ( xl* ), liquid height ( yl* ), and particle front position ( x*p ):
2.5 Validation of the numerical method
To validate the numerical methods, the time histories of dimensionless position of liquid front and liquid height were quantitatively compared with those of Sun and Sakai [13] who utilized numerical and experimental methods to study dam break flow containing particles under similar conditions we considered in this study, as shown in Fig. 4. Consequently, the comparison of the
13
ACCEPTED MANUSCRIPT present results with the experimental and numerical ones of Sun and Sakai [13] gives a good validation, which guarantees the reliability of the present numerical methods. Additionally, comparison in terms of the particle distribution and free-surface formation at each time was carried out with Sun and Sakai’s results [13], as shown in Fig. 5. The free-surface was indicated by blue line and the particle motions were individually captured with their velocity.
CR IP T
Sun and Sakai [13] considered the effects of a rising gate by adding an obstacle in front of the liquid column, which was not done in the present study. Although there are some discrepancies on the free surface in the early stage, the overall process of fluid flow and particle motion showed similar results. After collapse of liquid column, free-surface advances against the wall and the
AN US
particle bed expands over the bottom of the domain. In addition, the particle front moves more slowly than the liquid front at t * < 4.
M
3. Results and discussion
PT
ED
3.1 Effect of particle density on liquid front and particle front
To verify the effects of the particles on the free surface and particles, an additional simulation
CE
of a classical liquid-gas flow in a dam break was performed independently with the same conditions, and the results were compared. Fig. 6 presents typical three-dimensional perspective
AC
views and side views, which show the time evolution of the particles and the free surface between the liquid and the gas for the liquid-gas flow and liquid-gas-particle flows. In the case of the liquid-gas-particle flows, DR values of 1.1, 1.5, and 2.5 were representatively considered. To clarify the dependence of the front positions of the free surface, the temporal evolution of the freesurface profiles for the different flows on a two-dimensional (x-y) plane is plotted in Fig. 7. When the liquid column starts to collapse due to gravity, the front of the liquid and particles advances
14
ACCEPTED MANUSCRIPT against the wall. In this process, the liquid front behaves differently when the particles are present. In the initial stage of the dam breaking, the front positions of the liquid are almost identical, regardless of the presence of particles. However, over time, the front position of the liquid in the liquid-gas-particle flow lags behind that in the liquid-gas flow, as shown in Fig. 6 (a). At t * =2, the difference between the front position of the liquid according to presence of particles was greater
CR IP T
than that at t * =1. This tendency remained until t * =3, as shown Figs. 6 (b) and (c). Above results indicating the effect by particle can be clearly identified according to time in Fig. 7.
These observations similar to findings by Sun and Sakai [13] and Li and Zhao [50], who observed that the front position of the liquid in a liquid-gas-particle dam break flow seems to
AN US
move slower than that in a liquid-gas dam break flow.
We also identified the effects of DR on the front position of the liquid in liquid-gas-particle flows. Regardless of DR, the front position of the liquid is almost the same at each time, as shown
M
in Fig. 7. However, the front position of the particles was greatly influenced by DR. In the initial stage, the front position of the particles was almost the same, like that of the liquid. As shown in
ED
Fig. 6 (a), the font position of the particles at larger DR lagged behind that at smaller DR at t*=1.
PT
Over time, the difference of the front position of particles at different DR was clear, as shown in Figs. 6 (b) and (c).
CE
Consequently, the DR has no effect on the liquid front, but it makes the particle front move with different velocity. Based on these results, we can classify characteristics of the movement in
AC
terms of the front positions of the liquid and the particles. For DR=1.1, the particle front is always ahead of the particle front after t * =1. In the case of DR=1.5, the front positions of the liquid and particles are almost the same every time. However, the particle front for DR=2.5 lagged behind the liquid front after t * =1. Fig. 8 shows the dimensionless positions and velocity of the liquid front for the liquid-gas flow and liquid-gas-particle flows, for which DR values ranging from 1.1 to 2.5 were considered.
15
ACCEPTED MANUSCRIPT The dimensionless liquid front position was averaged in the spanwise direction, xl* , to observe the variation in this direction. As shown in Figs. 6 and 7, the liquid front in the liquid-gas flow generally moves faster than that in the liquid-gas-particle flows against the wall. Early on, the liquid front position is almost the same, regardless of the particles. However, the particles
liquid front hits the wall on the other side, as shown in Fig. 8 (a).
CR IP T
affect the liquid front position after t * =0.5. This difference can be observed clearly until the
The results were also compared with those for a liquid-gas flow from Sun and Sakai [13], Martin and Moyce [51], and Koshizuka and Oka [17], as shown in Fig. 8 (a). This was done to
AN US
emphasize the effects of the particles. The maximum velocity of the liquid front with the particles was identified at t * =2.58, and that without particles occurred at t * =2.5. Before the maximum velocity, the liquid front velocity without particles was always greater than that with particles. After the maximum velocity of the liquid front, the liquid front velocity decreases, regardless of
rapidly rather than that with particles
M
the presence of particles. Particularly, the liquid front velocity without particles decreases more
ED
The liquid front in the liquid-gas-particles flows was also compared at different DR. As shown
PT
in Fig. 8 (a), the liquid front positions for different DR generally have the same profile. This means that the effect of DR on the liquid front is almost negligible. As with the liquid front
CE
position, the time history of the liquid front velocity for different DR is almost the same. However, the particle front is indeed affected by DR. Fig. 9 shows top views of the free
AC
surface and particle formation for different DR. As expected, greater DR impedes the particle motion, and the particle velocity decreases as DR increases. Furthermore, instability of the particles was observed at DR=1.1. However, the particle formation along the bottom for DR=2.5 does not display significant instabilities and stays perpendicular to the direction of the flow. Fig. 10 shows the spanwise-averaged dimensionless front position of particles, x*p , and the streamwise velocity of the particle front, u p , for different particle densities. Regardless of DR,
16
ACCEPTED MANUSCRIPT the profile of the particle front position is composed of three sections. The first section has the smallest slope in particle front position in Fig. 10 (a) and corresponds to the front velocity at t * < 0.2, as shown in Fig. 10 (b). In this section, the front position is independent of DR. The second section has a sharp slope in particle front position in Fig. 10 (a) induced by the higher velocity than the first section. In this section, the front position depends on DR. As a result, the particle
CR IP T
front position decreases as the DR increases, as shown in Fig. 10 (a). This suggests that the front velocity decreases with increasing DR at t * > 0.2, as shown in Fig. 10 (b). Finally, before the particles front collide against the wall, the front velocity in near the wall drops dramatically and reaches zero.
AN US
Fig. 11 shows the moments when the liquid and particle fronts arrive at the opposite wall ( ta* ) for different DR. ta* of the liquid front for different DR is consistent, while the particle front takes more time to arrive at the opposite wall as DR increases. Furthermore, ta* of the particle
M
front turned out to be proportional to DR.
ED
The maximum velocities of the liquid front and particle front are plotted according to DR in Fig. 12. The maximum velocity of the liquid front is independent of DR, while that of the particle
PT
front decreases as DR increases, as shown in Fig. 12 (a). The difference between the maximum velocity of the liquid front and particle front, u*max , is also compared according to DR in Fig.
CE
12 (b). In the case of DR below 1.7, u*max is positive. For DR=1.7, u*max is almost zero,
AC
and after that u*max is negative as DR increases.
3.2 Motion regimes
As discussed in section 3.1. The front speed of the liquid and particles according to DR changes differently over time. Therefore, differences occur in the front position of the liquid and particles. Fig. 13 presents the time histories of the spanwise-averaged front position of the liquid
17
ACCEPTED MANUSCRIPT and particles for different DR. In the initial stage of the dam breaking, the front positions of the liquid and the particles are almost identical for different DR at t * 0.4. However, at t* > 0.4 and DR = 1.1, the particle front position is larger than the liquid front position at t* > 0.4. On the other hand, the particle front position at DR = 2.5 is smaller than the liquid front position at t* > 0.4. These results were attributed to the varying speed of the particle front while liquid front
overlapped position before impingement of the mixture flow.
CR IP T
speed hardly changed. At DR = 1.5, the front positions of the liquid and particle keep up the
Based on these results, three different regimes in the motion of the dam break flow with particles were identified. The regimes were simply classified based on the distance between
AN US
the particle front and the liquid front, which was strongly related to DR. The motion regimes identified are the advancing motion regime, overlapping motion regime, and delaying motion regime. The distance between the particle front and liquid front, Dp l x p xl , was
M
used to classify the regimes. The advancing motion regime occurs when the particle front is in
ED
front of that of the liquid ( 0.5d Dp l , where d is the particle diameter). The overlapping motion regime is in the region of 0.5d Dp l 0.5d , where the front positions of the liquid and
PT
particles are almost identical. The delaying motion regime occurs when the particle front lags
CE
behind the liquid front ( Dp l 0.5d ). Illustrations of the regimes at the same instance are shown in Fig. 14. The three motion
AC
regimes are greatly influenced by DR and the travel distance of dam break flow. The three intervals of DR can be demarcated in terms of the three motion regimes. For each computation, the travel distance was considered before the impingement of the mixture flow to confirm that the results present distinct flow regimes. The map of the motion regimes is shown in Fig. 15, where the x- and y-axis are the dimensionless travel distance and DR, respectively. In the region of 1.1 DR 1.4, the overlapping motion regime occurred in an early stage. After that, the advancing motion regime dominated the remaining domain. In addition, the area of the advancing regime 18
ACCEPTED MANUSCRIPT decreased with increasing DR. The series of changing regimes in the ranger of 1.1 DR 1.4 was named the overlapping-advancing (O-A) regime. At 1.5 DR 1.6, a transition regime from the advancing to the delaying motion regime was observed. In this stage, the delaying motion regime first appears at DR = 1.6. In addition, the overlapping motion regime was more dominant than the other regimes. At 1.6 DR 2.5, the
CR IP T
liquid and particle fronts moved in an overlapping-delaying (O-D) regime. Early on, the overlapping motion regime was observed, similar to the O-A regime. However, the delaying motion regime dominated the remaining domain. In addition, the area of the delaying motion regime expanded as DR increased. The map thus showed an O-A regime, transition regime, and
AN US
O-D regime. The map could be used to predict which material (liquid or particles) will arrive more quickly at a target point.
M
4. Conclusion
ED
This study investigated the liquid-gas-particles mixture flow during dam break. The CFDDEM numerical method was used for the simulation. Such numerical method consist of the open
PT
source codes in OpenFOAM for CFD and CFDEMproject, which includes LIGGGHTS for DEM
CE
simulation and calculates the interaction data between CFD and DEM. The computed results are in good agreement with those from previous studies. The main purpose is to investigate the effects
AC
of particle density on the behavior of the particles and the free surface. Basically, if the particles were added in dam break flow, they interfered with the movement of
the liquid front. Therefore, the liquid front in a liquid-gas-particle dam break flow seems to move slower than that in a liquid-gas dam break flow. In addition, the liquid front in the liquid-gasparticles flows was also compared at different particle density. The time history of liquid front positions for different DR generally have the same profile. Which means that the effect of particle
19
ACCEPTED MANUSCRIPT density on the liquid front is almost negligible. However, as the particle density grows, particle motion was impeded, and the particle velocity decreased. Consequently, the DR has no effect on the liquid front, but it makes the particle front move with different velocity. Furthermore, instability of the particles formation was observed as particle density was the smallest. However, the particle formation along the bottom for the biggest particle density does not display
CR IP T
significant instabilities and stays perpendicular to the direction of the flow. The time when the liquid front reach at the opposite wall is consistent according to particle density. However, such time for particle front decrease as particle density increases. Furthermore, the time of the particle front turned out to be proportional to particle density.
AN US
Based on these results, we can classify characteristics of the movement in terms of the front positions of the liquid and the particles. The characteristic motion was classified into an advancing, overlapping, and delaying motion regime. In the initial stage of the liquid collapse, the
M
overlapping motion regime governs the front positions of the liquid and the particles, regardless of the particle density. This results in almost identical front positions of the liquid and the particles.
ED
In the advancing motion regime, the particle front moves more quickly than the liquid front does. For this reason, the particles front is ahead of the liquid front after penetration of the particless
PT
across the free surface. In the delaying motion regime, the particle front lags behind the liquid
CE
front. The three regimes are greatly influenced by particle density and the travel distance of the dam break flow. The map of the motion regimes was created as afunction of DR and streamwise
AC
position. The map can be used to predict which phase will arrive more quickly at a target point.
Acknowledgment
This work was supported by National Research Foundation of Korea (NRF) grant founded by the Korea government (MSIP) through GCRC-SOP (No. 2011–0030013) and (NRF2015R1D1A3A01020867).
20
ACCEPTED MANUSCRIPT References [1] Iverson RM. The physics of debris flows. Rev Geophs 1997; 35(3):245-96. [2] Highland LM, Ellen SD, Christian SB, Brown III WM. Debris-flow hazards in the United States. https://pubs.usgs.gov/fs/fs-176-97/fs-176-97.html; 1997
CR IP T
[3] Griffiths RW. The dynamics of lava flows. Annu Rev Fluid Mech 2000;32(1):477-518. [4] Roche O, Montserrat S, Ni o Y, Tamburrino A. Experimental observations of water-like behavior of initially fluidized, dam break granular flows and their relevance for the propagation of ash-rich pyroclastic flows. J Geophys Res-Solid Earth 2008;113(B12).
AN US
[5] Ancey C. Plasticity and geophysical flow: A review. J Non-Newton Fluid 2007;142(1):435.
[6] Balmforth NJ, Kerswell RR. Granular collapse in two dimensions. J Fluid Mech 2005;538:399-428.
M
[7] Bagnold RA. Experiments on a gravity-free dispersion of large solid spheres in a
ED
Newtonian fluid under shear. P Roy Soc A-Math Phy 1954;225(1160):49-63. [8] Bakhtyar R, Ghaheri A, Yeganeh-Bakhtiary A, Barry DA. Process-based model for
PT
nearshore hydrodynamics, sediment transport and morphological evolution in the surf and swash zones. Appl Ocean Res 2009;31(1):44-56.
CE
[9] Tang HS, Keen TR, Khanbilvardi R. A model-coupling framework for nearshore waves,
AC
currents, sediment transport, and seabed morphology. Commun Nonlinear Sci 2009;14(7):2935-47.
[10] Wang JY, Tang HS, Fang HW. A fully coupled method for simulation of wave-currentseabed systems. Commun Nonlinear Sci 2013;18(7):1694-709. [11] De Rooij F. Sedimenting particle-laden flows in confined geometries. Diss. University of Cambridge; 2000.
21
ACCEPTED MANUSCRIPT [12] Guo L, Morita K. Numerical simulation of 3D sloshing in a liquid–solid mixture using particle methods. Int J Numer Meth Eng 2013;95(9):771-90. [13] Sun X, Sakai M. Three-dimensional simulation of gas–solid–liquid flows using the DEM–VOF method. Chem Eng Sci 2015;134:531-48. [14] Shigematsu T, Liu PLF, Oda K. Numerical modeling of the initial stages of dam-break
CR IP T
waves. J Hydraul Res 2004;42(2):183-95. [15] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39(1):201-25.
[16] Harlow FH, Welch JE. Numerical calculation of time-dependent viscous incompressible
AN US
flow of fluid with free surface. Phys Fluids 1965;8(12):2182-9.
[17] Koshizuka S, Oka Y. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng 1996;123(3):421-34.
M
[18] Shao S, Lo EY. Incompressible SPH method for simulating Newtonian and nonNewtonian flows with a free surface. Adv Water Resour 2003;26(7):787-800.
ED
[19] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J Comput Phys 2003;191(2):448-75.
PT
[20] Monaghan JJ, Rafiee A. A simple SPH algorithm for multi-fluid flow with high density
CE
ratio. Int J Numer Meth Fl 2013;71(5):537-61. [21] Abdolmaleki K, Thiagarajan KP, Morris-Thomas MT. Simulation of the dam break
AC
problem and impact flows using a Navier-Stokes solver. In: 15th Australasian Fluid Mechanics Conference The University of Sydney; 2004.
[22] Lobovský L, Botia-Vera E, Castellana F, Mas-Soler J, Souto-Iglesias A. Experimental investigation of dynamic pressure loads during dam break. J Fluid Struct 2014;48:407-34. [23] Kleefsman KMT, Fekken G, Veldman AEP, Iwanowski B, Buchner B. A volume-of-fluid based simulation method for wave impact problems. J Comput Phy 2005;206(1):363-93.
22
ACCEPTED MANUSCRIPT [24] Ozmen-Cagatay H, Kocaman S. Dam-break flow in the presence of obstacle: experiment and CFD simulation. Eng Appl Comp Fluid 2011;5(4):541-52. [25] Nsom B, Debiane K, Piau JM. Bed slope effect on the dam break problem. J Hydraul Res 2000;38(6):459-64. [26] Park IR, Kim KS, Kim J, Van SH. Numerical investigation of the effects of turbulence
CR IP T
intensity on dam-break flows. Ocean Eng 2012;42:176-87. [27] Jiang CB, Chen J, Tang HS, Cheng YZ. Hydrodynamic processes on beach: Wave breaking, up-rush, and backwash. Commun Nonlinear Sci 2011;16(8):3126-39.
Sci 2012;17(12):5273-85
AN US
[28] Zamankhan P. Rollers in low-head dams – challenges and solutions. Commun Nonlinear
[29] Cundall PA, Strack OD. A discrete numerical model for granular assemblies. Geotechnique 1979;29(1):47-65.
M
[30] Stevens AB, Hrenya CM. Comparison of soft-sphere models to measurements of collision properties during normal impacts. Powder Technol 2005;154(2):99-109.
ED
[31] Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol, 1993;77(1):79-87.
PT
[32] Lin SY, Chen YC. A pressure correction-volume of fluid method for simulations of
CE
fluid–particle interaction and impact problems. Int J Multiphas Flow 2013;49:31-48. [33] Guo L, Morita K, Tobita Y. Numerical simulations of gas-liquid-particle three-phase
AC
flows using a hybrid method. J Nucl Sci Technol 2016;53(2):271-80.
[34] Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J 1977;82(12):1013-24.
[35] Yamada Y, Sakai M, Mizutani S, Koshizuka S, Oochi M, Murozono K. Numerical simulation of three-dimensional free-surface flows with explicit moving particle simulation method. J Atom Energ Soc Jpn 2011;10(3):185–93.
23
ACCEPTED MANUSCRIPT [36] Sun X, Sakai M, Yamada Y. Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method. J Comput Phys 2013;248:147-76. [37] Delannay R, Valance A, Mangeney A, Roche O, Richard P. Granular and particle-laden flows: from laboratory experiments to field observations. J Phys D Appl Phys 2017;50(5):053001. Ltd.
OpenFOAM
-
The
Open
Source
CFD
Toolbox.
CR IP T
[38] OpenCFD
http://www.openfoam.com; 2009.
[39] Kloss C, Goniva C, Hager A, Amberger S, Pirker S. Models, algorithms and validation for opensource DEM and CFD–DEM. Prog Comput Fluid Dy, an International Journal
AN US
2012; 12(2-3):140–52.
[40] Goniva C, Kloss C. Hager A, Pirker S. An open source CFD-DEM perspective. In: Proceedings of OpenFOAM Workshop; 2010. p. 22-4.
M
[41] Gidaspow D. Multiphase flow and fluidization-continuum and kinetic theory descriptions. Academic Press; 1994. ISBN 0122824709.
ED
[42] Goniva C, Gruber K, Kloss C. Sediment erosion a numerical and experimental study. River Flow 2012; 2012. p. 415-22. ISBN 9780415621298.
PT
[43] Issa, R. Solution of the implicitly discretized fluid flow equations by operator- splitting. J
CE
Comput Phys 1986;62(1):40-65. [44] Hanosi IM, Jan D, Szabo KG, Tel T. Turbulent drag reduction in dam-break flows. Exp
AC
Fluids 2004;37(2):219-29.
[45] Barenblatt GI. Scaling, self-similarity, and intermediate asymptotics: dimensional analysis and intermediate asymptotics. Cambridge Press; 1996. ISBN 0521435226
[46] Barenblatt GI, Chorin AJ, Prostokishin VM. A note concerning the Lighthill “sandwich model” of tropical cyclones. P Natl Acad Sci USA 2005;102(32):11148-50.
24
ACCEPTED MANUSCRIPT [47] Bersch M, Hulshof J, Prostokishin VM. Flow laminarization and acceleration by suspended particles. SIAM J Appl Math 2015;75(4):1852-83. [48] Kawaguchi T, Sakamoto M, Tanaka T, Tsuji Y. Quasi-three-dimensional numerical simulation of spouted beds in cylinder. Powder Technol 2000;109(1):3–12. [49] Link JM, Cuypers LA, Deen NG, Kuipers JAM. Flow regimes in a spout–fluid bed: A
CR IP T
combined experimental and simulation study. Chem Eng Sci 2005;60(13):3425–42 [50] Li X, Zhao J. Numerical simulation of dam break by a coupled CFD-DEM approach. In: The 15th Asian Regional Conference on Soil Mechanics and Geotechnical Engineering; 2016;2(18):691-6.
AN US
[51] Martin JC, Moyce WJ. An experimental study of the collapse of liquid columns on a rigid
AC
CE
PT
ED
M
horizontal plane. Philos Trans R Soc Lond A 1952;244(882):312–24.
25
ACCEPTED MANUSCRIPT
LIST OF TABLES TABLE 1. Properties of particles TABLE 2. Properties of fluids
CR IP T
Table 1 Properties of particles Type
Value
Diameter, d (m)
0.0027 1100,1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2,500 3,883
The number of particles Spring constant (N/m) Restitution coefficient
1,000 0.9 0.3
ED
M
Friction coefficient
AN US
Density, ρp (kg/m3)
PT
Table 2 Properties of fluids
Value
Liquid density, ρl (kg/m3)
1000
AC
CE
Type
Liquid viscosity (Pa·s)
10-3
Gas density, ρg (kg/m3)
1
Gas viscosity (Pa·s)
10-5
26
ACCEPTED MANUSCRIPT
LIST OF FIGURES FIGURE 1. Illustration of force acting on particles in DEM. FIGURE 2. The algorithm of CFD-DEM.
CR IP T
FIGURE 3. Schematic of the three-dimensional computational domain. FIGURE 4. Comparison of dimensionless positions of (a) liquid front and (b) liquid height with the present results and Sun and Sakai’s results [13] of computation (Com) and experiment (Exp).
AN US
FIGURE 5. Comparison between Sun and Sakai’s results [13] (left column) and the present results (right column) for the particle distribution and free-surface formation at (a) t*=2 and (b) t*=4.
FIGURE 6. Illustrations of dam break flow without particles (1st column) and with particles
M
for DR= 1.1 (2nd column), 1.5 (3rd column), and 2.5 (4th column) at (a) t*=1.0, (b)
ED
t*=2.0, and (c) t*=3.0. The white surface and line indicate the free surface. FIGURE 7. Free-surface profile for dam break according to presence of particle at t*= 0.5, 1.0,
PT
2.0, and 3.0.
CE
FIGURE 8. Time histories of (a) dimensionless position of liquid front with data of liquid-gas flow from Martin and Moyce [40], Koshizuka and Oka [17], and Sun and Sakai
AC
[13]; (b) dimensionless streamwise velocity of liquid front for different DR.
FIGURE 9. Free-surface and particle formation for different DR at t*=3. FIGURE 10. Time histories of (a) dimensionless position and (b) dimensionless streamwise velocity of particle front for different DR. FIGURE 11. Arrival time (ta*) at the opposite wall versus DR. FIGURE 12. (a) Maximum front velocity of liquid and particles and (b) difference of maximum velocity between liquid and particles according to DR. 27
ACCEPTED MANUSCRIPT FIGURE 13. Comparison of dimensionless front position of liquid and particles. FIGURE 14. The front position lines of typical motion regimes for (a) advancing motion regime, (b) overlapping motion regime, and (c) delaying motion regime (top view). The blue and orange lines represent the liquid front and particle front, respectively.
CR IP T
FIGURE 15. Classification of three motion regimes as a function of DR and dimensionless
AC
CE
PT
ED
M
AN US
streamwise position (x*).
28
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 1. Illustration of force acting on particles in DEM.
29
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 2. The algorithm of CFD-DEM.
30
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
Fig. 3. Schematic of the three-dimensional computational domain.
31
CE
PT
ED
M
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
(b)
AC
Fig. 4. Comparison of dimensionless positions of (a) liquid front and (b) liquid height with the present results and Sun and Sakai’s results [13] of computation (Com) and experiment (Exp).
32
ACCEPTED MANUSCRIPT
AN US
CR IP T
(a)
(b)
Fig. 5. Comparison between Sun and Sakai’s results [13] (left column) and the present results
AC
CE
PT
ED
M
(right column) for the particle distribution and free-surface formation at (a) t*=2 and (b) t*=4.
33
M
AN US
(a) t*= 1.0
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
(b) t*= 2.0
(c) t*= 3.0
Fig. 6. Illustrations of dam break flow without particles (1st column) and with particles for DR= 1.1 (2nd column), 1.5 (3rd column), and 2.5 (4th column) at (a) t*=1.0, (b) t*=2.0, and (c) t*=3.0. The yellow surface and line indicate the free surface.
34
CR IP T
ACCEPTED MANUSCRIPT
AN US
Fig. 7. Free-surface profile for dam break according to presence of particle at t*= 0.5, 1.0, 2.0,
AC
CE
PT
ED
M
and 3.0.
35
CR IP T
ACCEPTED MANUSCRIPT
CE
PT
ED
M
AN US
(a)
(b)
AC
Fig. 8. Time histories of (a) dimensionless position of liquid front with data of liquid-gas flow from Martin and Moyce [40], Koshizuka and Oka [17], and Sun and Sakai [13]; (b) dimensionless streamwise velocity of liquid front for different DR.
36
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
Fig. 9. Free-surface and particle formation for different DR at t*=3.
37
CR IP T
ACCEPTED MANUSCRIPT
CE
PT
ED
M
AN US
(a)
(b)
AC
Fig. 10. Time histories of (a) dimensionless position and (b) dimensionless streamwise velocity of particle front for different DR.
38
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 11. Arrival time (ta*) at the opposite wall versus DR.
39
CR IP T
ACCEPTED MANUSCRIPT
PT
(b)
ED
M
AN US
(a)
CE
Fig. 12. (a) Maximum front velocity of liquid and particles and (b) difference of maximum
AC
velocity between liquid and particles according to DR.
40
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 13. Comparison of dimensionless front position of liquid and particles.
41
ACCEPTED MANUSCRIPT
AN US
CR IP T
(a)
PT
(c)
ED
M
(b)
Fig. 14. The front position lines of typical motion regimes for (a) advancing motion regime,
CE
(b) overlapping motion regime, and (c) delaying motion regime (top view). The blue and
AC
orange lines represent the liquid front and particle front, respectively.
42
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Fig. 15. Classification of three regimes as a function of DR and dimensionless streamwise
AC
CE
PT
ED
position (x*).
43