- Email: [email protected]

Contents lists available at ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

Large-Eddy Simulation (LES) of settling particle cloud dynamics Ruo-Qian Wang a, Adrian Wing-Keung Law b, E. Eric Adams a,⇑ a b

Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, MA 02139, USA School of Civil and Environmental Engineering, Nanyang Technological University, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 17 October 2013 Received in revised form 1 August 2014 Accepted 8 August 2014 Available online 19 August 2014 Keywords: Euler–Lagrangian Four-way coupling CFD–DEM Sediment disposal Entrainment coefﬁcient

a b s t r a c t A series of Euler–Lagrangian four-way coupling Large-Eddy Simulations (LES) are performed to study the dynamics of settling particle clouds. The numerical method is ﬁrst validated by comparing with existing experimental results. Then, a parametric study is performed with various initial release shapes, particle sizes, and cloud buoyancies. Three issues are examined in detail, including the initial release aspect ratio, two-phase interactions and polydispersion. Incorporating these factors, empirical relationships are developed to predict the phase separation time and height, the position of the cloud front and the edge radius based on numerical results. The entrainment rate and deposition patterns are also analyzed. These relationships should be useful for engineering applications. Ó 2014 Elsevier Ltd. All rights reserved.

Introduction A particle cloud settling through a water column is inherently a two-phase interactive process that is closely coupled. Once released, the particle cloud ﬁrst accelerates due to the density difference (so called ‘‘acceleration stage’’), and then forms a circulating system with the solid particles and a rotating ﬂuid vortex ring (so called ‘‘thermal stage’’). Finally, individual particles separate from the vortex ring and descend with their terminal settling velocity, while the ﬂuid vortex ring is left behind due to the removal of driving force by the heavier particles (so called ‘‘dispersive stage’’) (Rahimipour and Wilkinson, 1992; Lai et al., 2013). In this study, we focus on the thermal and dispersive stages of the particle cloud with the initial release of particles in the middle of the water column. The ‘‘acceleration stage’’ is very short with this release conﬁguration, and can therefore be ignored. The key in understanding the coupled system is the two-phase interactions that govern the settling dynamics. The route of momentum transfer differs between the thermal and dispersive stages. In the thermal stage, momentum is transferred from the particles to the ﬂuid generating the vortex ﬂow, which then redistributes the solid phase and feeds the momentum back to the particles. Contrary to this interactive loop, the dispersive stage has an open route – the momentum is primarily transferred from the solid to the ﬂuid phase through the settling process, which is ⇑ Corresponding author. Tel.: +1 617 253 6595. E-mail addresses: [email protected] (R.-Q. (A.W.-K. Law), [email protected] (E.E. Adams).

Wang),

http://dx.doi.org/10.1016/j.ijmultiphaseﬂow.2014.08.004 0301-9322/Ó 2014 Elsevier Ltd. All rights reserved.

[email protected]

then dissipated by viscous diffusion. In other words, momentum is lost after the transfer. Beyond the importance of basic understanding, including the ﬂuid phase in the engineering analysis has practical implications. Sediment particles released in the water body may possess dissolvable or volatile chemical substance, e.g. organic pollutants in disposed dredged sediment. The ﬂuid-phase pollution can pose a secondary impact on the surrounding environment alongside the solid phase. Detailed knowledge of the solid–ﬂuid phase interactions can help control and minimize the consequences. Until recently, the tracking of both phases has received little attention. Lai et al. (2013) presented the ﬁrst study that addresses the dynamics of both phases using the simpliﬁed approach of the Hill’s vortex for the ﬂuid phase. We aim to address this issue in fuller detail using the Large-Eddy Simulation (LES) numerical approach as the ﬁrst objective of the present study. From a physical perspective, the process of particle settlement is an initial value problem where the boundaries are relatively far away compared to the cloud sizes. Thus, the computational domain is relatively simple in geometry, and the initial release conditions set the tone for the settlement process. In previous investigations, the initial conditions often covered a selected range of several basic parameters, including particle sizes, materials, densities, and total mass (Rahimipour and Wilkinson, 1992; Bühler and Papantoniou, 1991; and Nakatrsuji et al., 1990). The typical particle cloud characteristics that were measured included the penetration depth of the particle cloud front (or centroid) and the cloud radius. Ruggaber (2000) was among the ﬁrst to systematically vary all the basic parameters mentioned and generalize his results into an

66

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

integral model. In addition, he also explored the initial condition of moisture content and quantiﬁed the mass in the trailing stems. Subsequently, additional ambient and initial conditions were studied, including the ambient stratiﬁcation (Bush et al., 2003), currents (Gensheimer et al., 2013), waves (Zhao et al., 2013), and release height (Zhao et al., 2012). The previous studies have largely ignored the effect of the release shape, which is highly relevant, e.g., in the laboratory design of particle cloud generators, buckets used for open water sediment disposal, or injection chambers used in plant seed dispersions. In addition, due to facility constraints, past researchers were mostly preoccupied by a particular release mechanism, e.g. a bowl shape container in Rahimipour and Wilkinson (1992), a cylindrical pipe in Ruggaber (2000), and two sizes of funnels in Bush et al. (2003). This restricted the research to certain release aspect ratios, precluding a systematic parametric study. As the second objective of the present study, we aim to complete a parametric study of the initial release conditions numerically as described in the following sections. Finally, the third objective of the present study concerns polydispersion: the particles in the release may be non-uniform and possess different sizes. Again, there have been very few systematic studies involving multi-size particles. Bühler and Papantoniou (1991, 2001) observed qualitatively that the particle distribution has a signiﬁcant effect on the growth and settlement of particle clouds. Here, we aim to clarify the quantitative details of the polydispersion dynamics also through the parametric study. As discussed above, we adopt LES for the present investigation. LES has great ﬂexibility in terms of simulating the particle cloud behavior with respect to the initial release shape and particle size distributions. In contrast, most classical integral models originally developed by Koh and Chang (1973) and extended by others (Ruggaber, 2000; Bush et al., 2003) assumed a speciﬁc cloud shape (e.g. a sphere or ellipse). The novel semi-integral method recently developed by Lai et al. (2013) approximated the ﬂuid phase motion as an expanding Hill’s spherical vortex, and tracked the motion of individual particle groups. The simpliﬁed approach is ﬂexible for various applications, but is also limited by the approximate nature of the assumed ﬂuid ﬂow ﬁeld. More advanced CFD methods can potentially provide more accurate predictions with various initial release conditions, if they can be validated. Previously, Li (1997) performed a two-ﬂuid simulation with a mixing length model as turbulence closure. Gu and Li (2004) simulated the two-phase ﬂows directly. The ﬂuid phase motion was solved by the Reynolds Averaged Navier–Stokes method with k— turbulence model, and the solid phase motion was computed with Lagrangian tracking. Harada et al. (2013) replaced the ﬂuid phase computational approach with LES, and obtained more detailed turbulence structures and better agreement with experimental data. Generally speaking, a relatively small number of particles were used in these previous studies to reduce computational resources, and particle– particle collision was typically ignored which might affect the computational results in the high particle concentration region near the point of release. In this study, we develop a new higher order numerical scheme that can circumvent the limitations and accommodate higher particle concentrations. In the following, Section ‘Numerical methods’ introduces the numerical method and computational domain; Section ‘Validation of numerical method’ validates the code by comparing with existing experiments; Sections ‘Phase separation and Penetration and growth’ analyze the monodispersion results to address the issues of phase separation and the effect of initial release condition; Section ‘Entrainment and deposition’ focuses on the entrainment rate and deposition patterns; the particle size distribution effect of polydispersion is shown in Section ‘Polydispersion’; and conclusions are drawn in Section ‘Summary’.

Numerical methods A series of 3-D numerical simulations was performed with Euler–Lagrangian LES using the software CFDEM, which couples the open source CFD toolbox (OpenFOAM) and the Discrete Element Method (DEM) package (LIGGGHTS) (Kloss et al., 2012). For the ﬂuid phase, the void fractioned Navier–Stokes equations,

@ af þ r af uf ¼ 0 @t ! @ af uf af p þ r af uf uf ¼ af r Rsl þ af mr ruf þ r s ; @t qf qf ð1Þ are solved numerically, where af is the volume fraction of the ﬂuid, uf is the velocity of the ﬂuid phase, t is time, p is the dynamic pressure, qf is the density of the ﬂuid, m is the kinematic viscosity, and Rsl is the inter-phase momentum transfer term to be discussed later. A ﬁnite volume method is adopted and all variables are ﬁltered except the sub-grid stress tensor s, which is resolved by the local dynamic one equation eddy viscosity sub-grid model (locDynOneEqEddy) (Fureby et al., 1997). Note that the ﬂuid phase solver is the same as in the CFDEM ofﬁcial package. The solid phase motion is resolved using the ODE by LIGGGHTS, i.e.

mp

dup ¼ RF; dt

ð2Þ

where mp is the particle mass, up is the particle velocity, and RF is the total force acting on the particle, which can be speciﬁed as

RF ¼ FG þ FS þ FD þ FL þ FA

ð3Þ

where FG is the gravity force, FS is the ﬂuid stress, FD is the drag force, FL is the lift force, and FA is the added mass force. These forces can be modeled as follow:

FG ¼ mp ð1 cÞg;

ð4Þ

which includes the gravitational and buoyancy forces;

FS ¼ mp rp=qp ;

ð5Þ

which is in the same form of Zhou et al. (2010); and

3 FD ¼ C D mp cjup ujðup uÞ=dp ; 4 Du dup ; FA ¼ C A cmp Dt dt 3 FL ¼ C L mp cjup ujðup uÞ x=ðdp jxjÞ; 4

ð6Þ ð7Þ ð8Þ

where the density ratio c ¼ qf =qp and u is the ﬂuid velocity at the position of the particle center. The last three force models are equivalent to Loth and Dorgan (2009). Speciﬁcally, the drag coefﬁcient is

( CD ¼

24 Rep

ð1 þ 0:15Re0:687 Þ; Rep 6 800 p Rep > 800:

0:44;

ð9Þ

The added mass coefﬁcient is C A ¼ 0:5, and the lift coefﬁcient is C L ¼ J

12:92

p

n

sﬃﬃﬃﬃﬃﬃﬃﬃ

x

Rep

þ Xp;eq C L;X h

ð10Þ

io C L;X ¼ 1 0:825 þ 0:15tanh 0:28 Xp;eq 2 tanh 0:18Re1=2 ; ð11Þ p ( " !#)( " sﬃﬃﬃﬃﬃﬃﬃﬃ #) sﬃﬃﬃﬃﬃﬃﬃﬃ x 2 x þ 0:191 1:92 ;ð12Þ J 0:3 1 þ tanh 2:5 log 10 þ tanh 6 3 Rep Rep x Xp;eq ¼ ð1 0:0075Rex Þ 1 0:062Re1=2 ð13Þ p 0:001Rep ; 2 x j jd p ; x ¼ ð14Þ up u

67

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

dp is the diameter of the particle, Rep ¼ dp jup uj=m, and 2 Rex ¼ dp jxj=m. Coding for the gravitational force and the ﬂuid stress is unchanged from the original ofﬁcial version, but the drag, added mass, and lift force models are newly coded in this study. The new codes are included in the most updated version of CFDEM (version 2.7.1) and the source code is open to the public. As implemented in the ofﬁcial release of CFDEM, the gravitational force and the ﬂuid stress are not fed back to the ﬂuid phase, whereas the newly developed drag, added mass, and lift forces are fed back as discussed later. In each time step, the solver starts a loop over all the particles. For each individual particle, the ﬂuid velocity is interpolated from the ﬂuid ﬂow ﬁeld based on the particle position. Then, the forces are calculated using the particle’s information according to the above equations, and the added mass, drag, and lift forces are fed back to the ﬂuid phase in order to conserve the total momentum. The calculated forces are also spread to LIGGGHTS and used to track the particle’s next position. The tracking process is divided into 100 small time steps and the particle pathlines would change if collision occurs. The inter-phase momentum term Rsl needs to be adjusted every time step to conserve the total momentum by the interface of

Ps P CFDEM, i.e. Rsl ¼ si¼1 ðFD;i þ FA;i þ FL;i Þ i¼1 mp;i , where s is the number of particles within a single grid cell. The original developers of CFDEM established the function Rsl to feed back the solid phase forces to the ﬂuid, and we have retained this function with our own models for FD ; FA , and FL . In the present study, the four-way inter-phase coupling is realized by LIGGGHTS, i.e. particle–particle interaction is incorporated. We assume that the particles can overlap by a small distance d when they collide. Then, the normal force and tangential force can be computed by

Fn ¼ kn d þ C n DVn ;

ð15Þ

and

Z t Ft ¼ kt DVt ðsÞds t þ C t DVt ; tc;0

ð16Þ

where DVn is the normal relative velocity at the contact point, DVt is the relative tangential velocity, t is the tangential vector at the contact point, t c;0 is the time when the contact begins, and C t and C n are model constants. As a result, the particles’ motion after collision can be modeled by

€p ¼ Fn þ Ft þ mp g mp x

ð17Þ

After every simulation, a passive tracer is applied by solving the scalar transport equation,

@ af T þ r af uf T ¼ 0; @t

ð18Þ

where T is the scalar concentration. The scalar transport reveals the ﬂuid phase motion. The computational domain was discretized with an unstructured mesh in a cylindrical shape as shown in Fig. 1. Finer grid cells were concentrated along the axis of symmetry to better capture the shear layer. The domain extended to Hdepth ¼ 1:2 m in depth and 0.567 m in diameter (two exceptions were BL4 and BL16, whose length scales were increased by factors of 4 and 16). Zero-gradient open boundary conditions were applied to all the surfaces. The particles were released in the middle of the quiescent domain with the same upper bound at 4.33% of Hdepth from the top. The initial release was contained in a predeﬁned cylinder with radius R0 and height H, the geometry of which was varied to obtain different aspect ratios of A ¼ H=R0 . The particles in the initial release were compacted with a volume fraction as ¼ 1 af ¼ p=6 ¼ 0:524. The non-buoyant

Fig. 1. Mesh of the computational domain.

passive tracer was released at a concentration of 1 g/cm3 in the ﬂuid inside the particles volume. The boundary of the ﬂuid phase was outlined by 1% of the maximum scalar concentration. All the simulations were executed with 32 processors in the High Performance Cluster in the Nanyang Technological University, Singapore. Validation of numerical method A non-trivial issue in multiphase simulation meshing is the relationship between the particles and grid cells. In principle, LES should pass the grid convergence test, which guarantees that the grid length scale is within the turbulence cascade range and thus that the ﬁnal results are independent of the grid resolution. This requires a relatively small grid length. At the same time, in the present numerical strategy, the grid cells have to be bigger than the particles to properly apply the point model of particles in the governing equations. If the particles are relatively small or the Reynolds number is fairly low, there is no conﬂict. However, this is not the typical situation and thus the conﬂict needs to be reconciled. Here, we propose a new numerical scheme: instead of using small grid cells, large grids are adopted but with the spatial discretization enhanced to high order accuracy, which also leads to a quicker grid convergence with the reduced number of grid cells. Speciﬁcally, the convection terms are discretized using the Total Variation Diminishing (TVD) scheme with a cubic limiter, and the cubic discretization schemes are applied to second-order differential terms. The success of this new numerical scheme for two-phase simulations is demonstrated by the following grid convergence test. The grid convergence test was performed by systematically increasing the resolution of the grid. The conditions are summarized in Table 1. Five different grids with 6098, 13,186, 17,695 and 26,931 cells were used. From Fig. 2, we can observe that the different grids yielded almost the same penetration rate of the particle cloud, and only the grid with 6098 cells was different from the others in cloud growth. Therefore, for all the following studies, we used the grid with 26,931 cells, for which results can be considered to be grid-independent. It should be noted that the smallest grid cell still had a volume that was 1.5 times that of the particles. Note that the penetration rate is deﬁned as the time derivative of the cloud front position. The solid phase front is deﬁned by the lowest

Table 1 Flow conﬁguration for grid convergence. Case

dp (mm)

R0 (mm)

Total mass (g)

H (mm)

GC

0.513

8.4

8.4

50.4

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

45

10

40

9

35

8

cloud radius (cm)

cloud front (cm)

68

30 25

grid size of 6098 grid size of 13186 grid size of 17695 grid size of 26931

20 15 10

6 5

grid size of 6098 grid size of 13186 grid size of 17695 grid size of 26931

4 3 2

5 0

7

1 0

0.5

1

1.5

2

2.5

3

3.5

4

0

0

0.5

1

1.5

2

2.5

3

3.5

4

time (s)

time (s)

Fig. 2. Grid convergence test (the number in the legend indicates the total number of grid cells used).

particle within a cloud, and the ﬂuid phase is the lowest point of the 1% tracer concentration. Three experiments from Lai et al. (2013) were used to validate the present numerical method. The ﬂow cases are speciﬁed in the ﬁrst three rows of Table 2, the only difference among them being the particle size. The comparison for Case A is shown in Fig. 3. The green dots are the particle group (solid phase) and the bright yellow region indicates the tracer distribution (ﬂuid phase). At t ¼ 0:2 s, the ﬂuid phase still remained intertwined with the solid phase, with the particle cloud starting to accelerate and a vortex ring being formed in the ﬂuid phase. The numerical simulation results generally matched the experimental measurements – overlap can be observed in both the left and right columns in Fig. 3. At t ¼ 0:4 s, the particles were half inside the tracer cloud and half outside, which indicated that the particle cloud and the ﬂuid were separating. Again, the numerical results well captured the transient separation process shown in the experiments. At t ¼ 0:8 s, the phase separation was nearly complete, and by t ¼ 1:2 s, the solid particles began to settle individually. The numerical simulations were able to predict the position and bowl shape of the particle cloud. Similar observations can also be made in cases B and C (omitted here for brevity). Note that there was a small difference in size comparison, because in the experiments, the particles had a distribution in sizes, but the particles in the numerical simulations were strictly uniform. In the dispersive stage, particles with nonuniform size would settle with different velocities. As a result, the experimental particle cloud was stretched vertically with smaller particles settling slowly and bigger particles settle faster.

This is consistent with the observation in the following section of poly-dispersion. Quantitative comparison further validated the present numerical approach. We focused on four variables – the frontal depth and cloud radius of both the solid and ﬂuid phases which are shown in Fig. 4. The solid lines denote the numerical results, and the error bars bound the lower and upper limits of ﬁve realizations of the same experiment. As shown in Fig. 4a, good agreement was found in the frontal position of the solid phase (zs ), despite a small overestimation in cases A and C and a slight underestimation in case B. A better comparison was noted in the ﬂuid phase front (zf ), where the numerical predictions were generally within the ranges of experimental uncertainties (Fig. 4b). Furthermore, the cloud radius of the solid phase (r s ) shown in Fig. 4c was well predicted, although a slight discrepancy can be identiﬁed at earlier times in case C. The most challenging comparison was in the ﬂuid phase radius (rf ) shown in Fig. 4d. The numerical predictions for cases A and B generally followed the experimental measurements, while case C had a signiﬁcant discrepancy at the beginning, which was removed at later times. In summary, the numerical results compared well to the experimental data which veriﬁed the applicability of the numerical approach. Phase separation After the code validation, a parametric study with 14 simulations was implemented covering three sizes of particles, with a range of total mass and aspect ratios. The ﬂow cases are summarized in Table 2, which also lists the corresponding

Table 2 Flow conﬁguration for monodispersion studies. Case

dp (mm)

R0 (mm)

Total mass (g)

H (mm)

Aspect ratio

Ra

Re

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16

0.73 0.51 0.26 0.73 0.51 0.26 0.51 0.26 0.26 0.73 0.51 0.26 0.99 2.20

9 9 9 6.7 6.7 6.7 5.6 5.6 3.9 3.9 3.9 3.9 36 144

3 3 3 3 3 3 3 3 2 3 3 3 192 12,288

11 11 11 20 20 20 28 28 39 58 58 58 43 171

1.2 1.2 1.2 3 3 3 5 5 10 15 15 15 1.2 1.2

19 41 214 36 76 398 106 557 772 104 220 1158 41 41

2667 2667 2667 1956 1956 1956 1654 1654 1072 1147 1147 1147 21,339 170,709

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

69

of an individual particle. Ra can also be used to scale up the results from lab-scale experiments to ﬁeld scale (Ruggaber, 2000). The most pronounced phenomenon in such two-phase ﬂows is phase separation. While it is inherently a continuous process, for analytical convenience we deﬁne a phase separation time T sp as the time that the centroid of the particle cloud reaches the front of the ﬂuid phase. Because the front of the ﬂuid phase is usually easy to observe in the experiments and the particle centroid is normally available after image processing, this criterion can be easily applied to the experimental measurements. Fig. 5a sketches the separation time as discussed above. Fig. 5b further illustrates this time for case B through plotting the positions of the particle cloud centroid and ﬂuid front. In the ﬁgure, it can be observed that the particle centroid started slightly behind the ﬂuid front, but eventually accelerated ahead of the ﬂuid front (a red dash line marks the time of crossing over). Also shown in the ﬁgure is the portion of particles still held inside the ﬂuid phase cloud. As the phase separation developed, the portion reduced from 1 (full overlap) to around zero. Using this deﬁnition, the phase separation time for different cases was determined and shown in Fig. 6 normalized by a time scale of

R20 ﬃ: T s ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ B=q0

ð20Þ

Including the aspect ratio A in the normalization, all the separation times collapse to a straight line, as follow:

T sp

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ B=q0 1=4 ¼ 0:85 A : w2s

ð21Þ

For comparison, the phase separation time determined alternatively by the instant that the particle front velocity approaches the settling velocity is also shown in Fig. 6. Here, 1.4ws is deﬁned to be the threshold, below which the particle cloud is considered in the dispersive phase (Bühler and Papantoniou, 2001). The ﬁgure shows that these two methods are roughly the same and can be ﬁtted by the same Eq. (21). In a similar fashion, the ratio of the phase separation height over the separation time collapses to a straight line in Fig. 7. Using Eqs. (21) and (22), the phase separation time and heights can be predicted according to the ﬂow condition. Furthermore, the fact that the aspect ratio A is raised to a small power (1/4) in the equation of T sp indicates that the separation has only a weak dependence on the initial aspect ratio. The ﬁt in Fig. 7 is given by

Hsp ¼ 2:4ws T sp :

Fig. 3. Visual comparison between numerical simulations (left) and experiments (right).

Rayleigh number, Ra, and the characteristic Reynolds number, 1=6 R0 =m, for each case. The non-dimensional Rayleigh Re ¼ Bg 2 =qs

ð22Þ

Combining Eqs. (21) and (22), we can draw a comparison to Bush et al. (2003) as shown in Fig. 8. Bush et al. (2003) gave a lower and an upper limit for their empirical equation for the separation height (solid lines in Fig. 8). The initial release aspect ratio was not reported in their experiments, hence we assumed a range of 1–3 based on their experimental description. Furthermore, they deﬁned the phase separation time differently – as the instant that the cloud centroid reached the individual particles’ settling velocity – which should give a slightly larger phase separation height compared to ours. However, generally speaking, Fig. 8 shows that Eq. (22) is able to predict the trend and range of phase separation heights observed in their experiments.

number is deﬁned as

Ra ¼

B

q0 w2s R20

;

ð19Þ

and characterizes the ratio of buoyancy force to viscous drag force, P where B ¼ g ð1 cÞ Ni¼1 mpi is the total buoyancy, N is the total number of particles within a particle cloud, and ws is the settling velocity

Penetration and growth The phase separation results derived above can help understand the general settling characteristics including the penetration and growth, which are represented by the time evolution of the front position and the radius of the solid and ﬂuid phases of the settling

70

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

0.5 0.9

0.45

0.8

0.4

0.7

0.35

f

0.5 A B C

0.4 0.3

0.3 0.25 0.2 0.15

0.2

0.1

0.1

0.05

0

0

2

4

6

8

0

10

0

2

4

6

8

10

time (s)

(a) Solid phase front

(b) Fluid phase front

0.12

0.12

0.1

0.1

0.08

0.08

0.06

A B C

0.02

0

2

4

6

8

12

0.06 A B C

0.04

0.04

0

A B C

time (s)

rf (m)

rs (m)

z (m)

zs (m)

0.6

0.02

0

10

0

2

4

6

8

10

time (s)

time (s)

(c) Solid phase radius

(d) Fluid phase radius

12

14

Fig. 4. Penetration and growth of particle clouds compared to experiments (lines: numerical results).

Particle portion in fluid phase Particle cloud center depth (m) Fluid cloud front depth (m)

1

0.8

0.6

0.4

0.2

0

0

2

4

6

8

10

time (s)

(a) A sketch of phase separation

(b) Definition of phase separation

Fig. 5. Illustration of phase separation.

particle cloud. They were extracted from the numerical predictions in Figs. 9, 10, and 13. As shown below, the proper scaling is able to simplify the present numerical results andﬃ collapse all the data. pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Deﬁning the length scale of Lm ¼ B=q0 =ws , the penetration distance can be normalized and generalized into simple relationships.

For dimensional scaling, the penetration distance before and after phase separation has to be treated separately. The difference between the two normalizations is the presence of the aspect ratio A. Eqs. (23) and (24) are best ﬁt equations of the penetration distance shown in Figs. 9 and 10 (see Fig. 11).

71

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

1200 1000

Tsp /(TsA1/4 )

800 600 400

by 1.4ws

200

by center−front cross 0

0

500

1000

1500

Ra Fig. 6. Phase separation time determined by two methods (symbols: determined separation time; dash line: Eq. (21)).

160 140 120

sp

H /R0

100 80 60 40

The ﬁtting results reveal that the settling dynamics before and after the phase separation (i.e. in the thermal and dispersive stages) are different, as anticipated. Before phase separation, particle clouds develop as ‘‘thermals’’, whereas after phase separation, particles settle individually. This was conﬁrmed in the ﬁtting exponents of the non-dimensional time in the solid phase, i.e. in Eq. (23). For t < T sp , the exponent was 2/3, which was consistent with the penetration rate of starting plume (Diez et al., 2003); for t > T sp the penetration was linear in time, and close to the settling velocity of the particles. Although theoretically the coefﬁcient of the second term should be equal to one, a value of 1.1 was obtained in the present observation. This may be attributed to our deﬁnition of the dispersive stage, which covers a section where some particles still remain in the ﬂuid phase and contribute buoyancy to the penetration. The difference in the settling dynamics in the two stages can also be observed in the ﬂuid phase as demonstrated in Eq. (24), but no direct analogy can be drawn to known phenomenon (e.g. thermals/puffs or starting plumes/jets (Wang et al., 2010)). Finally, we address the effect of the aspect ratio A, whose inﬂuence is only seen before phase separation. In other words, particle clouds only ‘‘remember’’ the initial shape in the thermal stage, and the dispersive stage is not directly affected. Eqs. (23) and (24) also clarify the difference between the solid and ﬂuid phases. First, as expected, the ﬂuid phase penetrates slower than the solid phase, which results in the phase separation. Second, the behavior of the ﬂuid phase is more sensitive to the initial release shape than the solid phase, as indicated by the exponents in these equations. Similarly, the growth of the particle cloud can also be extracted in Fig. 13 and ﬁtted by

20

rf rs t ¼ ¼ 0:7 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Lm Lm B=q0 =w2s

0 0

10

20

30

40

50

60

T /(R /w ) sp

0

s

Fig. 7. Phase separation height (symbols: determined separation height; dash line: Eq. (22)).

3500 Bush et al. upper limit Bush et al. lower limit Present formula A=1 Present formula A=3

3000

Hsp/R0

2500

!0:3 :

ð25Þ

Different from the penetration results, the growth cannot clearly be separated into two stages. The growth of the solid and ﬂuid phases followed the same temporal rate. The effect of the initial aspect ratio A was also not observed, which suggests that the growth was independent of the initial release shape. To summarize, Eqs. (23)–(25) can be used to predict the position and size of the particle cloud released in the water column, which could be useful in engineering applications. Entrainment and deposition

2000 1500 1000 500 0 0

200

400

600

800

1000

1200

Ra Fig. 8. Comparison of the phase separation height between Bush et al. (2003) (solid lines) and the present Eq. (22) (dash lines).

8 > > tﬃ < 2:1A0:17 pﬃﬃﬃﬃﬃﬃﬃ

2=3

zs B=q0 =w2s ¼ Lm > tT 0:34 > : 1:89A þ 1:1 pﬃﬃﬃﬃﬃﬃﬃspﬃ

B=q0 =w2s

t < T sp

t > T sp

ð23Þ

;

and

8 > 0:3 tﬃ > pﬃﬃﬃﬃﬃﬃﬃ > 1:5A <

3=5

B=q0 =w2s zf ¼ 4=5 Lm > > tT sp > : 1:36A0:45 þ 0:75 pﬃﬃﬃﬃﬃﬃﬃ ﬃ 2 B=q0 =ws

t < T sp

t > T sp

:

ð24Þ

One of the key ingredients in integral modeling is the entrainment rate, which determines the general behavior of the model in terms of the particle cloud settlement and growth. Under the self-similar assumption, the entrainment rate is often represented by the entrainment coefﬁcient, which is the ratio of the velocities of entrainment over the cloud settling. In this section, the entrainment coefﬁcient is extracted by post-processing the present twophase Large-Eddy Simulation results aiming to further reveal the basic settling dynamics. The entrainment ﬂux can be derived by integrating the velocity normal to the sphere that contains the whole particle cloud. Specifically, we ﬁnd the particle cloud center and the distance from the center to the farthest particle (Fig. 12). Then, the containing sphere is built with the center at the determined particle cloud center and the radius equal to the determined distance. At last, the normal ﬂux to the containing sphere is numerically integrated to yield the entrainment ﬂux. The entrainment coefﬁcients with different particle sizes and initial aspect ratios are shown in Fig. 13a and b, respectively. The entrainment coefﬁcient has a sharp drop in the beginning; then, it quickly climbs up to the peak; and it slowly decays to zero. As observed from all the simulations, the entrainment

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

4 3 2 1

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16 formula

20 s

5

25

15

0

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16 formula

10

s

zs/A0.17 /((B/ρ0)1/2 /ws)

6

z /((B/ρ )1/2 /w )

72

5

0

0 0

1

2

3

4

0

5

10

15

20

t/(R2/(B/ρ )1/2 )

t/((B/ρ )1/2 /w2)

(a) Before phase separation (t < Tsp )

(b) After phase separation (t > Tsp )

0

0

0

s

Fig. 9. Normalized penetration of solid phase front.

zf /A

0.3

3 2.5 2 1.5 1 0.5

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16 formula

8

/ws)

/((B/ρ0)1/2 /ws)

3.5

10

1/2

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16 formula

6

zf /((B/ρ0)

4

4 2

0

0

0

2

4 2

1/2

t/(R0/(B/ρ0)

0

5

10

15

20

t/((B/ρ0)1/2 /w 2s )

)

(a) Before phase separation (t < Tsp )

(b) After phase separation (t > Tsp )

Fig. 10. Normalized penetration of ﬂuid phase front.

1

0.5

0

0

10

2 s

1.5 1 0.5 0

20

2 1/2 t/(R /(B/ρ ) ) 0 0

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16 formula

0

1.5

2.5

f

A B C A3 B3 C3 B5 C5 C10 A15 B15 C15 BL4 BL16 formula

r /((B/ρ )1/2 /w )

0

rs /((B/ρ )1/2 /ws)

2

0

10

20

2 1/2 t/(R /(B/ρ ) ) 0 0

(a) Solid phase

(b) Fluid phase

Fig. 11. Normalized growth of particle clouds.

coefﬁcient always converges to a curve, which can be ﬁtted by pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ae ¼ 0:45 exp t 16 B=q0 =w2s . As shown in Fig. 13a, the cloud with smaller particles has a lower maximum entrainment coefﬁcient and a smaller overshoot. In Fig. 13b, the cloud with a higher initial aspect ratio has a larger increase in entrainment coefﬁcient. The overshoot is more apparent in low aspect ratios and case B15 has no overshoot. The deposition pattern is relevant in certain engineering applications; however, it’s not the main topic here. As a brief discussion,

we focus on a particular case – case B. Fig. 14 shows the radial distribution of particle projections onto the bottom at different time. A clear ring shape pattern can be observed, which is enlarging with few particles in the center and a peak at the boundary. Beyond the peak, the particle number quickly decreases to zero. Note that this projection pattern can reﬂect the ﬁnal deposition pattern if the cloud reaches the bottom in the dispersive stage. However, the ﬁnal deposition pattern might be quite different from the projection pattern mentioned above, if the particle cloud impacts a soft

73

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75 7

number of particles

8

4

2

0

0.02

0.04

0.06

0.08

0.1

0.12

Radius (m)

bottom and causes the sediment to surge. More experimental studies on the deposition can be found in Gensheimer et al. (2013).

Fig. 14. Radial distributions of particle projections for case B.

Table 3 Flow conﬁguration for polydispersion studies.

Polydispersion In real applications, the released particles are not always uniform in size or ‘‘monodispersed’’. A spectrum of particle size is much more likely, yet there are limited studies on this issue. Questions remain such as how the ‘‘polydispersion’’ changes the dynamics; how signiﬁcant is the effect; and whether there is a simple way to predict the cloud size and position. To address these questions, several ‘‘polydispersion’’ releases were simulated numerically. Their conﬁgurations can be found in Table 3. The set of simulations covered two types of particle distribution by number, i.e. uniform and Gaussian. The median particle diameter by number was around 0.51 mm; the speciﬁc range was listed in the table and their distributions can be found in Fig. 15. For engineering applications, the particle diameter d50 , deﬁned by the size that 50% of the mass is ﬁner, is also listed. Corresponding to d50 and the maximum particle diameter respectively, the median (w50 ) and the maximum (wmax ) settling velocities were calculated using the formula offered by Brown and Lawler (2003) and are shown in the same table. All the other ﬂow parameters including initial particle cloud radius, total mass and initial shape are the same as case B in Table 2. As shown in Fig. 16, the cloud patterns were similar before phase separation but different afterwards. In the thermal stage (Fig. 16a), the particles were well mixed and distributed similarly to the monodispersion cases. In contrast, the particles settled according to their sizes or individual terminal settling velocities in the dispersive stage (Fig. 16b). As a result, the particles were

Distr. by num.

dp (mm)

d50 (mm)

w50 (cm/s)

wmax (cm/s)

Bu1 Bu2 Bu3 Bg1 Bg2

Uniform Uniform Uniform Gaussian Gaussian

0.26–0.73 0.26–0.76 0.43–0.60 0.45–0.57 0.22–0.81

0.61 0.64 0.53 0.52 0.57

8.8 9.3 7.6 7.4 8.2

10.6 11.2 8.7 8.3 12.0

0.15

0.05 0 −0.05 −0.1

Entrainment Coefficient

A B C

0.1

0.1 0.05 0

B B3 B5 B15

−0.05 −0.1 −0.15

−0.15 −0.2

Case

sorted by their sizes, with bigger particles settling faster. The vertical distribution of size was also consistent with the particle distribution by number at release. This observation is easy to understand. In the thermal stage, the two phases are closely coupled and the solid phase is continuously re-distributed by the ﬂuid phase. This offers a good stirring mechanism keeping the particles well-mixed. Conversely, in the dispersive stage, particles are more separated and able to settle individually. This analysis thus suggests that the thermal stage dynamics for polydispersed particles is similar to that for monodispersed particles, and that the different particle sizes can be modeled using the same formulae developed for monodispersion. In addition, the phase separation results of monodispersion can be applied in a similar manner to the polydispersion. To test the above hypothesis, we develop an empirical model to compare with the polydispersion results. First, we deﬁne pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Lp1 ¼ B=q0 =w50 and Lp2 ¼ B=q0 =wmax . Dictating that Lp ¼ Lp1 at t < T sp and Lp ¼ Lp2 at t > T sp , the polydispersion data are normalized

0.15

Entrainment Coefficient

0.02 s 0.06 s 0.1 s 0.3 s 0.5 s 1s 1.5 s 2s

6

0

Fig. 12. The entrainment ﬂow of the particle cloud in case B (the dash line marks the surface of the containing sphere).

x 10

0

5

10

15 1/2

t/((B/ρ0)

(a)

2 /w s )

20

25

−0.2

0

5

10 1/2

t/((B/ρ0)

2 /w s )

(b)

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Fig. 13. Entrainment coefﬁcients with (a) different particle sizes and (b) initial aspect ratios (Dash lines are given by ae ¼ 0:45 exp t 16 B=q0 =w2s ).

74

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

3000 Bu1 Bu2 Bu3 Bg1 Bg2

Particle number

2500 2000 1500 1000 500 0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Particle diameter (mm) Fig. 15. Particle size distribution.

(a) Thermal stage (t = 0.1s)

(b) Dispersive stage (t = 6.0s) Fig. 16. Comparison of polydispersion (from left to right: cases Bu1, Bu2, Bu3, Bg1, Bg2 as in Table 3).

35

14 Bu1 Bu2 Bu3 Bg1 Bg2 empirical model

30

20

10

z f /L p

zs /Lp

25

12

15

8 6

10

4

5

2

Bu1 Bu2 Bu3 Bg1 Bg2 empirical model

0

0 0

5

10

15

20

25

0

5

10

15

20

t/(R20/(B/ρ0)1/2 )

t/(R2/(B/ρ )1/2 ) 0 0

(a) Solid phase

(b) Fluid phase

Fig. 17. Normalized penetration of particle clouds in polydispersion.

25

30

75

R.-Q. Wang et al. / International Journal of Multiphase Flow 67 (2014) 65–75

1.8

1.5

1.6 1.4 1.2 p2

1

f

0.8

s

r /L

r /Lp2

1

Bu1 Bu2 Bu3 Bg1 Bg2 empirical model

0.5

Bu1 Bu2 Bu3 Bg1 Bg2 empirical model

0.6 0.4 0.2 0

0 0

2

4

6

8

10

12

0

5

10

t/(R20/(B/ρ0)1/2 )

t/(R20/(B/ρ0)1/2 )

(a) Solid phase

(b) Fluid phase

15

20

Fig. 18. Normalized growth of particle clouds in polydispersion.

and the results callapse as shown in Figs. 17 and 18 using Lp in penetration and Lp2 in growth. Correspondingly, replacing Lm by Lp in Eqs. (23) and (24) and Lp2 in (25), the collapsed data can be ﬁtted well by the current empirical model. The most signiﬁcant discrepancy is observed in the radius of the ﬂuid phase, where the present empirical model overestimates the numerical results by around 20%. Summary In the present paper, a new numerical scheme with LES is developed for two-phase simulations to address three speciﬁc issues related to the settlement of particle clouds in the water column: the two-phase separation, the effect of the initial cloud shape and polydispersion dynamics. Three series of simulations were performed. The ﬁrst series validated the current numerical approach. The second series provided quantitative results on the monodispersion dynamics including phase separation, the penetration and growth of the particle cloud, and the entrainment rate and deposition patterns. Finally, the last series clariﬁed the polydispersion dynamics. With proper scaling, we obtain best-ﬁt equations of the present numerical results. They show that the initial aspect ratio plays a weak role, being more important in penetration than in growth, and it only affects the behavior before phase separation. These equations can be easily applied in engineering analysis. Acknowledgements This research was supported by the National Research Foundation Singapore through the Singapore–MIT Alliance for Research and Technology’s CENSAM IRG research programme. References Brown, P., Lawler, D., 2003. Sphere drag and settling velocity revisited. J. Environ. Eng. 129 (3), 222–231.

Bühler, J., Papantoniou, D., 1991. Swarms of coarse particles falling through a ﬂuid. In: Lee, J.T., Cheung, T.K. (Eds.), Environmental Hydraulics. Balkema, Rotterdam, pp. 135–140. Bühler, J., Papantoniou, D.A., 2001. On the motion of suspension thermals and particle swarms. J. Hydraul. Res. 39 (6), 643–653. Bush, J.W.M., Thurber, B.A., Blanchette, F., 2003. Particle clouds in homogeneous and stratiﬁed environments. J. Fluid Mech. 489, 29–54. Diez, F.J., Bernal, L.P., Faeth, G.M., 2003. Round turbulent thermals, puffs, starting plumes and starting jets in uniform crossﬂow. J. Heat Transfer 125 (6), 1046. Fureby, C., Tabor, G., Weller, H.G., Gosman, A.D., 1997. A comparative study of subgrid scale models in homogeneous isotropic turbulence. Phys. Fluids 9 (5), 1416. Gensheimer, R., Adams, E., Law, A., 2013. Dynamics of particle clouds in ambient currents with application to open-water sediment disposal. J. Hydraul. Eng. 139 (2), 114–123. Gu, J., Li, C., 2004. Modeling instantaneous discharge of unsorted particle cloud in ambient water by an Eulerian Lagrangian method. J. Hydraul. Res. 42 (4), 399– 405. Harada, E., Tsuruta, N., Gotoh, H., 2013. Two-phase ﬂow LES of the sedimentation process of a particle cloud. J. Hydraul. Res. 51 (2), 186–194. Kloss, C., Goniva, C., Hager, A., Amberger, S., Pirker, S., 2012. Models, algorithms and validation for opensource DEM and CFD–DEM. Prog. Comput. Fluid Dynam. Int. J. 12 (2), 140–152. Koh, R.C.Y., Chang, Y.C., 1973. Mathematical Model for Barged Ocean Disposal of Wastes. Tech. rep., Washington, DC. Lai, A., Zhao, B., Law, A.-K., Adams, E., 2013. Two-phase modeling of sediment clouds. Environ. Fluid Mech. 13 (5), 435–463. Li, C., 1997. Convection of particle thermals. J. Hydraul. Res. 35 (3), 363–376. Loth, E., Dorgan, A.J., 2009. An equation of motion for particles of ﬁnite Reynolds number and size. Environ. Fluid Mech. 9 (2), 187–206. Nakatrsuji, K., Tamian, M., Murta, A., 1990. Dynamics behavior of sand cloud in water. In: Proceeding of Conference on Physical Modeling of Transport and Dispersion. MIT, Boston, 8C.1–8C.6. Rahimipour, H., Wilkinson, D., 1992. Dynamic behavior of particle clouds. In: 11th Australasian Fluid Mechanics Conference. Hobart, Australia, pp. 743–746. Ruggaber, G.J., 2000. Dynamics of Particle Clouds Related to Open-Water Sediment Disposal. Ph.D. Thesis, MIT. Wang, R.-Q., Law, A.W.-K., Adams, E.E., Fringer, O.B., 2010. Large-eddy simulation of starting buoyant jets. Environ. Fluid Mech. 11 (6), 591–609. Zhao, B., Law, A.W.K., Adams, E.E., Shao, D., Huang, Z., 2012. Effect of air release height on the formation of sediment thermals in water. J. Hydraul. Res. 50 (5), 532–540. Zhao, B., Law, A.W.K., Huang, Z., Adams, E.E., Lai, A.C.H., 2013. Behavior of sediment clouds in waves. J. Waterway Port Coastal Ocean Eng. 139 (1), 24–33. Zhou, Z.Y., Kuang, S.B., Chu, K.W., Yu, A.B., 2010. Discrete particle simulation of particle?ﬂuid ﬂow: model formulations and their applicability. J. Fluid Mech. 661, 482–510.