- Email: [email protected]
Numerical simulation of particle-laden gas flow by Vortex in Cell method
Powder Technology 235 (2013) 376–385
Contents lists available at SciVerse ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
Numerical simulation of particle-laden gas flow by Vortex in Cell method Tomomi Uchiyama ⁎ EcoTopia Science Institute, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
a r t i c l e
i n f o
Article history: Received 18 July 2012 Received in revised form 1 October 2012 Accepted 24 October 2012 Available online 31 October 2012 Keywords: Gas–particle two-phase flow CFD Vortex method Vorticity
a b s t r a c t This study proposes a simulation method for incompressible gas flow laden with small solid particles. It is based on a Vortex in Cell (VIC) method, which was originally developed to simulate incompressible single-phase flows. The proposed VIC method discretizes the gas vorticity field into vortex elements and computes the time evolution of the two-phase flow by calculating the behavior of the vortex element as well as the particle motion with the Lagrangian approach. This study also applies the VIC method to simulate a free fall of small solid particles in an unbounded air. The particles, initially arranged within a spherical region in a quiescent air, are made to fall, and their fall induces the air flow around them. The interactions between the particle motion and the air flow are favorably compared with the existing measured and simulated results, demonstrating the validity of the proposed VIC method. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Gas flows containing small solid particles are found in many industrial applications, such as energy conversion systems and chemical reactors. They are also closely related with atmospheric phenomena, such as snowstorm and microburst. Various simulation methods have thus far been proposed for such particle-laden gas flow [1]. Most of the methods are of a Lagrangian–Eulerian type: Lagrangian approach is applied to simulate each particle motion, while Eulerian methods such as a finite difference method are used for the simulation of gas flow. Vortex in Cell (VIC) method with redistribution of vortex element is one of the vortex methods to simulate incompressible flows [2]. It discretizes the vorticity field with vortex elements and computes the time evolution of the flow by tracing the convection of each vortex element through the Lagrangian approach. The Lagrangian calculation markedly reduces the numerical diffusion as well as ensures the higher numerical stability. Therefore, the VIC method is expected to be usefully employed for the direct numerical simulation (DNS) of turbulent flows and various results have been reported [3–5]. Cottet and Poncet [3] applied the VIC method for the wake simulation of a circular cylinder, and captured the streamwise vortices occurring behind the cylinder. Cocle et al. [4] analyzed the behavior of two vortex system near a solid wall, and made clear the interaction between two counter-rotating vortices and the eddies induced in the vicinity of the wall. Chatelain et al. [5] simulated trailing edge vortices, and visualized the unsteady phenomena caused by disturbances. These studies are concerned with the timedeveloping free shear flows. But the VIC method has not been applied to the turbulent flows bounded by solid walls, which are closely related with the turbulent friction and the heat transfer. Thus, the authors [6] ⁎ Tel./fax: +81 52 789 5187. E-mail address: [email protected]. 0032-5910/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.powtec.2012.10.051
performed the DNS for a turbulent channel flow, which is a representative example for the wall turbulent flows. When applying the classic VIC method, the oscillation of the flow increased with the progress of the computation, and eventually the computation collapsed. This was caused by the fact that the consistency among the discretized equations is not ensured. This was also because the solenoidal condition for the vorticity is not fully satisfied. To overcome such problems of the existing VIC method, the authors [6] proposed two improvements for VIC method to heighten the numerical accuracy and efficiency: A discretization method employing a staggered grid was presented to ensure the consistency among the discretized equations as well as to prevent the numerical oscillation of the solution, and a correction method for vorticity was also proposed to compute the vorticity field satisfying the solenoidal condition. The improved VIC simulation for the turbulent channel flow highlighted that the time evolution of the flow is fully performed and that the statistically steady turbulent flow is favorably obtained. It also demonstrated that the organized flow structures, such as the streaks and the streamwise vortices appearing in the near wall region, are successfully captured and that the turbulence statistics, such as the mean velocity and the Reynolds stress, agree well with the existing DNS results. The author [7] applied the method to simulate the motion of small air bubbles in water and the behavior of a vortex ring launched toward the bubbles. The simulation highlighted that the bubbles are entrained into the vortex ring and that the entrained bubbles are transported by the convection of the vortex ring. The objective of this study is to propose a simulation method, which is based on the improved VIC method, for incompressible gas flow laden with small solid particles. It computes the time evolution of the flow by calculating the behavior of the vortex element as well as the particle motion with the Lagrangian approach. Though Walther and Koumoutsakos [8] proposed a VIC method for particle-laden gas flow, the present method based on the improved VIC method is
T. Uchiyama / Powder Technology 235 (2013) 376–385
expected to have higher numerical accuracy. It should also be noted that the present method is considered to be usefully applied to the DNS for particle-laden turbulent flow. Because the applicability and validity of the method were demonstrated through the DNS for a single-phase turbulent channel flow [6]. This study also applies the proposed method to simulate a free fall of small spherical solid particles in air so as to demonstrate the validity of the method.
2.1. Governing equations for gas and particle When the gas-phase is incompressible, the mass and momentum conservation equations are expressed as: ð1Þ
1 1 2 þ ug ⋅∇ ug ¼ − ∇p þ ν∇ ug − F D ρg ρg ∂t
∂ug
ð2Þ
ð3Þ
where ω is the vorticity. It is postulated that the dominant forces on the particles are the drag and gravitational forces. It is also assumed that the virtual mass force, the Basset force, and the pressure gradient force are negligible. The equation of motion for a particle is written as: dup f ¼ u −up þ τp g dt
1−
! ρg g ρp
ð4Þ
where τp is the particle response time, ρpd 2/(18ρgν). f is the drag factor defined by the following equation using the drag coefficient CD and the particle Reynolds number Rep (=d|ug − up|/ν): f ¼ C D Rep =24:
ð5Þ
According to Schiller and Naumann [9], f at Rep ≤ 800 is given as: 0:687
f ¼ 1 þ 0:15Re
:
as: ð7Þ
The number of particles per unit volume, np, is written by the following equation: π 3 np ¼ α p = d 6
ug ¼ ∇ϕ þ ∇ ψ:
∇⋅ψ ¼ 0
ð10Þ
ð11Þ
Taking the curl of Eq. (10) and substituting Eq. (11) into the resultant equation, the vector Poisson equation for ψ is derived: ∇ ψ ¼ −ω:
ð12Þ
When substituting Eq. (10) into Eq. (1) and rewriting the resultant equation by using the relation ∇⋅ (∇× ψ) = 0, the Laplace equation for ϕ is obtained: 2
∇ ϕ ¼ 0:
ð13Þ
When ψ and ϕ are calculated by solving Eqs. (12) and (13) respectively, the velocity ug is computed from Eq. (10). The vorticity ω in Eq. (12) is estimated from Eq. (3). The Vortex in Cell (VIC) method discretizes the vorticity field with vortex elements to calculate the distribution of ω by tracing the convection of each vortex element. It is postulated that the position vector and vorticity for the vortex element v are xv = (xv, yv, zv) and ωv respectively. The Lagrangian form of the vorticity equation, Eq. (3), is written as: dxv ¼ uðxv Þ dt
ð14Þ
dωv 1 2 ¼ ∇⋅ðuðxv Þωðxv ÞÞ þ ν∇ ωðxv Þ− ∇ F D : ρg dt
ð15Þ
When the position and vorticity of vortex element are known at t = t, the values at t = t + Δt are computed from Eqs. (14) and (15). In the VIC method, the flow field is divided into computational grid cells to define ψ, ϕ and ω on the grids. If ω is defined at a position xk = (xk, yk, zk), the vorticity ω is given to xk, or a vortex element with ω is redistributed onto xk
ð6Þ
Considering Eq. (4), the drag force acting on a particle, Fp, is given π 3 f F p ¼ d ρp ug −up : 6 τp
According to the Helmholtz theorem, the gas velocity ug is represented as the summation of the gradient of a scalar potential ϕ and the curl of a vector potential ψ: the gas velocity ug is expressed as:
2
where FD is the force exerted by the particle acting on the gas-phase per unit volume. Taking the curl of Eq. (2) and substituting Eq. (1) into the resultant equation, the vorticity equation for the gas is derived: ∂ω 1 2 þ ∇⋅ ωug ¼ ∇⋅ ug ω þ ν∇ ω− ∇ F D ρg ∂t
2.2. Orthogonal decomposition of velocity and discretization of vorticity
The velocity calculated from Eq. (10) remains unaltered even when any gradient of scalar potential function is added to ψ. To remove this arbitrariness, a solenoidal condition is imposed on ψ:
2. Basic equations and numerical method
∇⋅ug ¼ 0
377
ωðxk Þ ¼
Nv X v
ωv W
x −x y −y z −z k v v W k v W k Δx Δy Δz
ð16Þ
where Nv is the number of vortex elements, and Δx, Δy and Δz are the grid widths. For the redistribution function W, the following equation is employed [10]: 8 2 3 jεjb1 < 1−2:5ε þ 1:5jεj W ðεÞ ¼ 0:5ð2−jεjÞ2 ð1−jεjÞ 1≤jε j≤2 : : 0 jε j > 2
ð17Þ
ð8Þ 2.3. Discretization by using staggered grid
where αp is the particle volume fraction. Since the force FD in Eqs. (2) and (3) is expressed as npFp, the following relation is derived from Eqs. (7) and (8). F D ¼ α p ρp
f ug −up τp
ð9Þ
For incompressible flow simulations, the MAC method and the SMAC method solve the Poisson equation for pressure, which is derived from the equation for pressure gradient and the continuity equation. They employ a staggered grid to ensure the consistency between the discretized equations as well as to prevent the numerical oscillation of the solution.
378
T. Uchiyama / Powder Technology 235 (2013) 376–385
In the VIC method, the Poisson equation for ψ and the Laplace equation for ϕ are derived. Therefore, the use of a staggered grid is considered to be indispensable for their discretization. But the existing VIC method has less explained the grids, and some simulations seem to employ a regular grid. The author proposed a VIC method, which uses a staggered grid for single-phase flow [6]. This simulation also employs the staggered grid illustrated in Fig. 1. The scalar potential ϕ and the particle volume fraction αp are defined at the center. The gas velocity ug is defined at the sides, while the vorticity ω and the vector potential ψ are defined on the edges.
D x y
2.4. Correction of vorticity In the VIC method, the vorticity field is discretized with vortex elements, and the field is expressed by superpositioning the vorticity distribution around each vortex element. The superposition is performed by Eq. (16). The vorticity field obtained by Eq. (16), denoted by ωr, does not always satisfy the solenoidal condition [10]. When the vorticity satisfying the condition is denoted by ωs (=∇× ug), ωr is expressed as [4]:
ωr ¼ ∇F þ ωs
z Fig. 2. Initial arrangement of particles.
ð18Þ
where F is a scalar function. Eq. (18) corresponds to the Helmholtz decomposition of ωr. Taking the divergence of Eq. (18), the Poisson equation for F is obtained: 2
∇ F ¼ ∇⋅ωr :
ð19Þ
If F is calculated from Eq. (19) and substituted into Eq. (18), the vorticity ωs is recalculated [4]. This correction for vorticity needs to solve the Poisson equation, causing the increase of the computational time. The author [6] has proposed a simplified correction method.
The vorticity before the correction, ωr, is connected with ψr through Eq. (12). When taking the curl of Eq. (12) and substituting Eq. (19) into the resultant equation, the following relations are obtained: 2
∇ ð∇⋅ψr Þ ¼ −∇⋅ωr ¼ −∇2 F:
ð20Þ
On the contrary to Eq. (11) assuming the solenoidal condition for ψ, the following equation is derived from Eq. (20) if the vorticity is not solenoidal. ∇⋅ψr ¼ −F:
ð21Þ
i-1/2
xi x i+1/2
x/D
y/D 0
3 -3 0
-3 -1
z k+1/2
0 Particles at t*=0
zk
Δz
z k-1/2
z x
Δx
y
5
Gravity
y j+1/2
yj
z/D
y j-1/2
Δy
Gas velocity Gas vorticity Vector potential Scalor potential Particle volume fraction
ug
vg
ω x , ψx
ωy , ψ y
Fig. 1. Staggered grid.
φ, αp
wg ω z , ψz 10 11 Fig. 3. Computational domain.
3
T. Uchiyama / Powder Technology 235 (2013) 376–385 Table 1 Simulation conditions.
379
t*=0
1.2 kg/m3 1.5 × 10−5 m2/s 1 mm 15 kg/m3 0.0463 s 0.216 m/s 100, 200, 400 0.0125, 0.025, 0.05 20d 0.5D/ut 1 6D × 6D × 12D 40 × 40 × 80 5.4 × 10−3
Gas density ρg Gas kinematic viscosity ν Particle diameter d Particle density ρp Particle velocity response time τp Particle terminal velocity ut Number of particles Np Initial particle volume fraction αp0 Diameter of particle-cluster D Time-scale of gas flow τg Stokes number St Computational domain Number of grids Time increment Δt/τp
2.7
x/D
1
0
-1 -1
0 y/D 1
4.59
6.88
9.18
11.6
When calculating ψr from Eq. (12) with the use of ωr and determining ϕr from Eq. (13), the curl of ur is transformed as follows: ∇ ur ¼ ∇ ð∇ϕr þ ∇ ψr Þ ¼ ∇ð∇⋅ψr Þ−∇2 ψr
ð22Þ
¼ −∇F þ ωr ¼ ωs :
-1
Fig. 5. Time evolution for particle distribution projected onto the horizontal plane (αp0 = 0.05).
t*=0 Eq. (22) demonstrates that the curl of the velocity ur calculated from ωr yields a vorticity ωs satisfying the solenoidal condition. If this recalculation of vorticity, Eq. (22), or the correction of vorticity is performed just after the calculation of Eq. (10), the discretization error of vorticity is completely removed and the time evolution of the flow is accurately simulated without solving the Poisson equation [Eq. (19)]. It should be noted that the staggered grid is indispensable to make the transformation in Eq. (22) applicable to the corresponding discretized equations.
z/D
0
2.7 2
4.59 2.5. Numerical procedure When the flow at time t is known, the flow at t + Δt is simulated by the following procedure:
6.88
4
(1) (2) (3) (4) (5) (6) (7) (8)
9.18 6
11.6
Calculate the particle motion from Eq. (4). Calculate the time variation of ω at every grid point from Eq. (15). Calculate the convection of vortex element from Eq. (14). Calculate ω from Eq. (16). Calculate ψ from Eq. (12). Calculate ϕ from Eq. (13). Calculate ug from Eq. (10). Correct the vorticity, or take the curl of ug to calculate the corrected vorticity.
3. Application to free falling of solid particles in air 3.1. Computational conditions
8 -2
-1
0
1
2
x/D Fig. 4. Time evolution for particle distribution projected onto the vertical plane (αp0 =0.05).
To confirm the validity of the proposed simulation method, the behavior of solid particles falling in an unbounded air is simulated by the method. At the initial time t ∗ (=t/τp) = 0, spherical solid particles are arranged within a spherical region of diameter D in a quiescent
380
T. Uchiyama / Powder Technology 235 (2013) 376–385
-1
air, as shown in Fig. 2. The particles are made to fall, and their fall induces the air flow around them, resulting in a particle-laden air flow. Fig. 3 shows the computational domain of 6D × 6D × 12D. It is resolved into 40× 40× 80 uniform grid cells. The particle position at the commencement of the fall (t∗ = 0) is given by using random numbers. The origin of coordinates is set at the center of the spherical region, and the z-axis is taken in the vertically downward direction. The particle diameter d and density ρp are 1 mm and 15 kg/m 3 respectively. The diameter of spherical region D is set at 20d. To investigate the effect of the number of particles Np, the simulations are performed for Np =100, 200 and 400. Accordingly, the particle volume fraction αp0 at t∗ = 0 ranges from 0.0125 to 0.05. The particle terminal velocity ut is 0.216 m/s. When the time scale of the induced air flow τg is assumed to be 0.5D/ut, the Stokes number St, defined as the ratio of the particle response time τp to τg, is 1. The time increment Δt is set at Δt/τp = 5.4 × 10−3. The computational conditions are listed in Table 1. Since the particles fall in an unbounded air, ϕ is zero and the analysis of Eq. (13) can be omitted. The second-order Adams–Bashforth method is applied to the Lagrangian calculation of the particle motion. The collision between the particles is ignored, because the particle volume fraction is low. For the calculation of the particle volume fraction, αp, the volume of the particles existing in each computational grid cell is computed, and the ratio of the volume against the cell volume is calculated. It should be mentioned that the trial simulations employing the existing vortex in cell method collapsed in the cases of Np = 200 and 400 due to the divergence of the gas velocity field.
t*=0
0
2.7 2
z/D
4.59
6.88
4
9.18 6
3.2. Results for particle volume fraction of αp0 = 0.05
11.6
Fig. 4 shows the time evolution for the particle distribution in the case of Np = 400 or αp0 = 0.05. The particle positions at six time points are projected onto a vertical plane (x–z plane). The particles arranged in the half upper part and half lower one of the spherical region at t∗ = 0 are colored by red and blue respectively. The particles fall with destroying their initial arrangement, namely they tend to disperse in the horizontal direction. Most of the particles initially arranged in the half upper part sediment relatively toward the bottom of the particlecluster, while the particles initially arranged in the half lower part rise relatively toward the top. It is found that the particle behavior depends on the particle initial position. The particle dispersion in the horizontal (x) direction becomes markedly when t∗ > 2.7. These particle motions are caused by the air flow induced by them, as explained later.
8 -2
-1
0
1
2
x/D Fig. 6. Time evolution for particle distribution near the central vertical section (αp0 =0.05).
0
z/D
4
8
t*=2.7 -3
t*=6.88 0
x/D
3
-3
0
x/D
t*=11.6 3
-3
0
x/D
Fig. 7. Time evolution for gas velocity on the central vertical section and particle distribution (αp0 = 0.05).
3
T. Uchiyama / Powder Technology 235 (2013) 376–385
381
z/D
-1
1.4
t*=2.7 3.8 6.2 8.6 11
-3
-3
y/D
0
0
x/D
3
Particle
Q
ωz
Q
ωz
Q
ωz
z/D
-1
1.4
t*=6.88 3.8 6.2 8.6 11
-3
-3
y/D
x/D
0
0 3
Particle z/D
-1
1.4
t*=11.6 3.8
6.2
8.6 11
-3
-3
y/D
0
0
x/D
3
Particle
Fig. 8. Particle distribution and gas vortical structures (αp0 = 0.05).
382
T. Uchiyama / Powder Technology 235 (2013) 376–385
-1 2
t*=0
up /ut
0
z/D
1
0
5
10
2.97
15
t*
2
Fig. 9. Time variation for falling velocity of particles (αp0 = 0.05).
5.13 The time evolution for the particle distribution projected onto a horizontal plane (x–y plane) is shown in Fig. 5, where the distributions at the same instance as Fig. 4 are plotted. Though the particles fall with dispersing in the horizontal direction, they distribute within a circular region. Most of the particles initially arranged in the half lower part, colored by blue, concentrate on the periphery of the circular region when t∗ ≥ 6.88. This is because the particles are whirled up by a vortex ring induced by them, as explained later. Fig. 6 shows the particle distribution near a vertical plane (x–z plane) passing through the center of the particle-cluster (x=y=0). The particles within a vertical slit region of |y/D|≤0.0625 are plotted. When t∗ = 2.7, the particles in the half upper part sediment relatively in the downward direction, and all particles distribute in a hemispherical region. The particle distributions of concave and thin layer types are observed at t∗ = 4.59. The distributions are attributable to the sedimentation of the particles in the half upper part. They are also caused by the fact that the particles near the periphery of the particle-cluster are relatively whirled up. Accordingly, the particles initially arranged in the half lower part of the spherical region, plotted by blue, concentrate on the periphery of the particle-cluster at t∗ ≥6.88. The particles disperse more in the horizontal direction with the passage of time. The abovementioned relative sedimentation, while up and dispersion are caused by the air flow induced by the falling particles. The time evolution for the air velocity on the vertical plane (x–z plane) passing through the center of the particle-cluster (x = y = 0) is
1.5
rp’/D
1
0.5
zp’/D
0 0.5
0.25
0
5
10
t* Fig. 10. Time variation for dispersion of particles (αp0 = 0.05).
15
4
7.69
10.1 6
12.7
8 -2
-1
0
1
2
x/D Fig. 11. Time evolution for particle distribution projected onto vertical plane (αp0 =0.025).
depicted in Fig. 7. The velocity at three time points and the corresponding particle distributions are plotted. When t∗ = 2.7, a downward air flow with higher velocity occurs just above the particle-cluster. The gravitational acceleration, which causes the particles to fall, induces the air flow through the drag force acting on the particles. The relative sedimentation of the particles initially arranged in the half upper part shown in Figs. 4 and 6 is due to this downward air flow. A pair of large-scale eddies appears around the downward air flow. Such organized vortical structure is the wake of the particle-cluster. It is a vortex ring, as explained later. At t ∗ = 6.88, the developed vortex ring whirls up the particles, and accordingly it produces the particle distribution of thin layer type. When t ∗ = 11.6, an air flow in the horizontal direction occurs just below the vortex ring, making the particles disperse in the direction. Experimental studies on a plane mixing layer [11] and a twodimensional wake [12] reported that the particles of the Stokes number St≃1 are focused on thin layers, outlining the boundaries of the largescale eddies. In this simulation, the particles distribute on the periphery of the vortex ring induced by them. Considering that St value of this simulation is 1, such particle distribution corresponds to the abovementioned preferential distribution. According to a two-dimensional simulation for the falling of particles of St=1 by Chen and Marshall [13], an air flow is induced behind the particles, and a concave distribution of the particles occurs due to the wake flow. These support the validity of the present simulation.
T. Uchiyama / Powder Technology 235 (2013) 376–385
383
0
z/D
4
8
t*=2.97 -3
0
t*=7.69 3
-3
x/D
t*=12.7
0
3
-3
x/D
0
3
x/D
Fig. 12. Time evolution for gas velocity on the central vertical section and particle distribution (αp0 = 0.025).
Fig. 8 shows the particle distribution and the air vortical structures. The second invariant of the velocity gradient tensor Q and the vertical (z) component 2of the vorticity ωz are presented with their iso-surfaces of Q = 1=τ p = 0.0016 and ωz/(1/τp) = ± 0.007. A vortex ring occurs just above the particles at t ∗ = 2.7. The scale of the vortex ring increases with the passage of time. Vortex tubes with positive and negative values of ωz entangle behind the particles. One can observe the three-dimensional vortical flow having vertical vorticities. The time evolution for the particle falling velocity ūp is shown in Fig. 9, where ūp is the velocity for the center of gravity of the particles. The velocity ūp increases rapidly just after the commencement of the particle fall. It reaches its maximum value at t ∗ = 3.1. This is attributable to the fact that the particles are accelerated by the downward air flow induced by them. The velocity ūp lessens monotonously at t ∗ > 3.1, and it tends towards the terminal velocity ut. Because the distance between the vortex ring and the particle-cluster increases with the time evolution as seen in Figs. 7 and 8, and accordingly the effect of the vortex ring on the particles reduces. A trial simulation using finer grids (48 × 48 × 96) was also performed, but the particle velocity shown in Fig. 9 was not affected by the grids. Therefore, the present grids are found to be appropriate.
To grasp the particle dispersion, the variances in the horizontal and vertical directions,r ′ p andz′ p respectively, are estimated by the following equations: 0 r
′
p
1 ¼@ Np
11=2 N p 2 2 X A X−xp þ Y−yp
0
11=2 Np 2 X 1 zp¼@ Z−zp A Np p ′
ð24Þ
where (X,Y,Z) is the position vector for the center of gravity of the particles. The time variations for r ′ p =D and z′ p =D are presented in Fig. 10. r ′ p =D increases with the passage of time. A marked increment is observed at t∗ ≥ 2.54. This is because the particles are transported in the horizontal direction by the developing vortex ring. z′ p =D takes its minimum value at t ∗ =3.24 because the particles sediment relatively within the particle-cluster, and therefore the vertical distance between the particles lessens, as found from the particle distribution at t∗ = 2.7 shown
2
up /u t
αp0 =0.05 1
0.025
0
5
ð23Þ
p
0.0125
t*
10
Fig. 13. Effect of αp0 on falling velocity of particles.
15
384
T. Uchiyama / Powder Technology 235 (2013) 376–385
1.5
1
, r p /D
αp0 =0.05 0.025
0.5
0.0125
0
, z p /D
0.5
αp0 =0.05 0.25
0.025 0.0125 0
5
10
15
t* Fig. 14. Effect of αp0 on particle dispersion.
in Fig. 4. z′ p =D increases at 3.24b t ∗ ≤9.85. This is because the particles are whirled up by the vortex ring behind them. 3.3. Effects of initial particle volume fraction αp0 In the case of Np = 200 or αp0 = 0.025, the particles distribute on the vertical plane (x–z plane) as shown in Fig. 11. When compared with the result for αp0 = 0.05 presented in Fig. 4, the particle sedimentation, whirl up and dispersion are smaller. Fig. 12 shows the distributions for particles and air velocity in the case of αp0 = 0.025. The air velocity is lower than that for αp0 = 0.05. The vortex ring is not observed. The effect of αp0 on the particle falling velocity ūp is demonstrated in Fig. 13. The velocity ūp becomes smaller with the decrement of αp0. This is because the drag force on the particles is smaller, and accordingly the downward air flow accelerating the particles in the downward direction has lower velocity. The velocity ūp for αp0 = 0.0125 is almost parallel to the particle terminal velocity ut. One can find that the air flow induced by the particles is markedly weak. The time variation for the particle dispersion is plotted in Fig. 14. The dispersion reduces with the decrement of αp0. In the case of αp0 = 0.0125, the dispersion is so small, indicating that the particles fall with maintaining their initial arrangement. 4. Concluding remarks A numerical simulation method for incompressible gas flow laden with small solid particles is proposed. It is based on a Vortex in Cell method, of which accuracy for single-phase flow simulation was successfully improved by the author in a prior study. The proposed simulation method is applied to compute a free fall of small solid particles in an unbounded air. The particles, initially arranged
within a spherical region in a quiescent air, are made to fall, and their fall induces the air flow around them, resulting in a particle-laden air flow. The particles are accelerated by the induced downward air flow just after the commencement of their fall. After the acceleration, they are whirled up by a vortex ring produced around the downward air flow. These interactions between the particle motion and the air flow are favorably compared with the existing measured and simulated results. In this study, the simulation for solid particles having a relatively low density was carried out. It should be noted that the proposed method is indeed applicable to the simulation for the higher density. This study applied the proposed method to the simulation for gas–particle two-phase flow in an unbounded space. It was confirmed in the author's previous study that the proposed method can perform the DNS for a singlephase turbulent channel flow. Therefore, the method promises to be usefully applied to wall turbulent flow laden with solid particles. Notations CD d D FD f g p Q t t⁎ u ut ūp W
drag coefficient of particle particle diameter diameter of particle cluster at t ∗ = 0 force exerted by particle acting on gas-phase per unit volume Stokes drag factor gravitational constant pressure second invariant of velocity gradient tensor time nondimensional time = t/τp velocity terminal velocity of particle falling velocity of particles redistribution function
T. Uchiyama / Powder Technology 235 (2013) 376–385
x, y, z orthogonal coordinates α volume fraction Δt time increment Δx, Δy, Δz grid width v kinematic viscosity of gas-phase ρ density τp particle response time ϕ scalar potential ψ vector potential ω vorticity = ∇ × ug Subscripts 0 g p x, y, z
initial condition gas particle components in directions of x, y, z
References [1] C.T. Crowe, M. Sommerfeld, Y. Tsuji, Multiphase Flow with Droplets and Particles, CRC Press, 1992. [2] I.P. Christiansen, Numerical simulation of hydrodynamics by the method of point vortices, Journal of Computational Physics 13 (1973) 363–379.
385
[3] G.-H. Cottet, P. Poncet, Advances in direct numerical simulations of 3D wall-bounded flows by Vortex-in-Cell methods, Journal of Computational Physics 193 (2003) 136–158. [4] R. Cocle, G. Winckelmans, G. Daeninck, Combining the vortex-in-cell and parallel fast multipole methods for efficient domain decomposition simulations, Journal of Computational Physics 227 (2008) 9091–9120. [5] P. Chatelain, A. Curioni, M. Bergdorf, D. Rossinelli, W. Andreoni, P. Koumoutsakos, Billion vortex particle direct numerical simulations of aircraft wakes, Computer Methods in Applied Mechanics and Engineering 197 (2008) 1296–1304. [6] T. Uchiyama, Y. Yoshii, H. Hamada, Direct numerical simulation of a turbulent channel flow by an improved vortex in cell method, International Journal of Numerical Methods for Heat and Fluid Flow (in press). [7] T. Uchiyama, Y. Yoshii, Numerical simulation of entrainment and transport of gas bubbles by vortex in cell method, in: Proc. ASME-JSME-KSME Joint Fluids Engineering Conference, Hamamatsu, 2011, (on USB flash drive). [8] J.H. Walther, P. Koumoutsakos, Three-dimensional vortex methods for particle-laden flows with two-way coupling, Journal of Computational Physics 167 (2001) 39–71. [9] L. Schiller, A.Z. Naumann, Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung, Zeitschrift des Vereins Deutscher Ingenieure 77 (1933) 318–321. [10] G.-H. Cottet, P. Koumoutsakos, Vortex Methods: Theory and Practice, Cambridge University Press, 2000. [11] F. Wen, N. Kamalu, J.N. Chung, C.T. Crowe, T.R. Troutt, Particle dispersion by vortex structures in plane mixing layers, Transactions of the American Society of Mechanical Engineers, Journal of Fluids Engineering 114 (1992) 657–666. [12] Y. Yang, C.T. Crowe, J.N. Chung, T.R. Troutt, Experiments on particle dispersion in a plane wake, International Journal of Multiphase Flow 26 (2000) 1583–1607. [13] H. Chen, J.S. Marshall, A Lagrangian vorticity method for two-phase particulate flows with two-way phase coupling, Journal of Computational Physics 148 (1999) 169–198.
Contents lists available at SciVerse ScienceDirect
Powder Technology journal homepage: www.elsevier.com/locate/powtec
Numerical simulation of particle-laden gas flow by Vortex in Cell method Tomomi Uchiyama ⁎ EcoTopia Science Institute, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
a r t i c l e
i n f o
Article history: Received 18 July 2012 Received in revised form 1 October 2012 Accepted 24 October 2012 Available online 31 October 2012 Keywords: Gas–particle two-phase flow CFD Vortex method Vorticity
a b s t r a c t This study proposes a simulation method for incompressible gas flow laden with small solid particles. It is based on a Vortex in Cell (VIC) method, which was originally developed to simulate incompressible single-phase flows. The proposed VIC method discretizes the gas vorticity field into vortex elements and computes the time evolution of the two-phase flow by calculating the behavior of the vortex element as well as the particle motion with the Lagrangian approach. This study also applies the VIC method to simulate a free fall of small solid particles in an unbounded air. The particles, initially arranged within a spherical region in a quiescent air, are made to fall, and their fall induces the air flow around them. The interactions between the particle motion and the air flow are favorably compared with the existing measured and simulated results, demonstrating the validity of the proposed VIC method. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Gas flows containing small solid particles are found in many industrial applications, such as energy conversion systems and chemical reactors. They are also closely related with atmospheric phenomena, such as snowstorm and microburst. Various simulation methods have thus far been proposed for such particle-laden gas flow [1]. Most of the methods are of a Lagrangian–Eulerian type: Lagrangian approach is applied to simulate each particle motion, while Eulerian methods such as a finite difference method are used for the simulation of gas flow. Vortex in Cell (VIC) method with redistribution of vortex element is one of the vortex methods to simulate incompressible flows [2]. It discretizes the vorticity field with vortex elements and computes the time evolution of the flow by tracing the convection of each vortex element through the Lagrangian approach. The Lagrangian calculation markedly reduces the numerical diffusion as well as ensures the higher numerical stability. Therefore, the VIC method is expected to be usefully employed for the direct numerical simulation (DNS) of turbulent flows and various results have been reported [3–5]. Cottet and Poncet [3] applied the VIC method for the wake simulation of a circular cylinder, and captured the streamwise vortices occurring behind the cylinder. Cocle et al. [4] analyzed the behavior of two vortex system near a solid wall, and made clear the interaction between two counter-rotating vortices and the eddies induced in the vicinity of the wall. Chatelain et al. [5] simulated trailing edge vortices, and visualized the unsteady phenomena caused by disturbances. These studies are concerned with the timedeveloping free shear flows. But the VIC method has not been applied to the turbulent flows bounded by solid walls, which are closely related with the turbulent friction and the heat transfer. Thus, the authors [6] ⁎ Tel./fax: +81 52 789 5187. E-mail address: [email protected]. 0032-5910/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.powtec.2012.10.051
performed the DNS for a turbulent channel flow, which is a representative example for the wall turbulent flows. When applying the classic VIC method, the oscillation of the flow increased with the progress of the computation, and eventually the computation collapsed. This was caused by the fact that the consistency among the discretized equations is not ensured. This was also because the solenoidal condition for the vorticity is not fully satisfied. To overcome such problems of the existing VIC method, the authors [6] proposed two improvements for VIC method to heighten the numerical accuracy and efficiency: A discretization method employing a staggered grid was presented to ensure the consistency among the discretized equations as well as to prevent the numerical oscillation of the solution, and a correction method for vorticity was also proposed to compute the vorticity field satisfying the solenoidal condition. The improved VIC simulation for the turbulent channel flow highlighted that the time evolution of the flow is fully performed and that the statistically steady turbulent flow is favorably obtained. It also demonstrated that the organized flow structures, such as the streaks and the streamwise vortices appearing in the near wall region, are successfully captured and that the turbulence statistics, such as the mean velocity and the Reynolds stress, agree well with the existing DNS results. The author [7] applied the method to simulate the motion of small air bubbles in water and the behavior of a vortex ring launched toward the bubbles. The simulation highlighted that the bubbles are entrained into the vortex ring and that the entrained bubbles are transported by the convection of the vortex ring. The objective of this study is to propose a simulation method, which is based on the improved VIC method, for incompressible gas flow laden with small solid particles. It computes the time evolution of the flow by calculating the behavior of the vortex element as well as the particle motion with the Lagrangian approach. Though Walther and Koumoutsakos [8] proposed a VIC method for particle-laden gas flow, the present method based on the improved VIC method is
T. Uchiyama / Powder Technology 235 (2013) 376–385
expected to have higher numerical accuracy. It should also be noted that the present method is considered to be usefully applied to the DNS for particle-laden turbulent flow. Because the applicability and validity of the method were demonstrated through the DNS for a single-phase turbulent channel flow [6]. This study also applies the proposed method to simulate a free fall of small spherical solid particles in air so as to demonstrate the validity of the method.
2.1. Governing equations for gas and particle When the gas-phase is incompressible, the mass and momentum conservation equations are expressed as: ð1Þ
1 1 2 þ ug ⋅∇ ug ¼ − ∇p þ ν∇ ug − F D ρg ρg ∂t
∂ug
ð2Þ
ð3Þ
where ω is the vorticity. It is postulated that the dominant forces on the particles are the drag and gravitational forces. It is also assumed that the virtual mass force, the Basset force, and the pressure gradient force are negligible. The equation of motion for a particle is written as: dup f ¼ u −up þ τp g dt
1−
! ρg g ρp
ð4Þ
where τp is the particle response time, ρpd 2/(18ρgν). f is the drag factor defined by the following equation using the drag coefficient CD and the particle Reynolds number Rep (=d|ug − up|/ν): f ¼ C D Rep =24:
ð5Þ
According to Schiller and Naumann [9], f at Rep ≤ 800 is given as: 0:687
f ¼ 1 þ 0:15Re
:
as: ð7Þ
The number of particles per unit volume, np, is written by the following equation: π 3 np ¼ α p = d 6
ug ¼ ∇ϕ þ ∇ ψ:
∇⋅ψ ¼ 0
ð10Þ
ð11Þ
Taking the curl of Eq. (10) and substituting Eq. (11) into the resultant equation, the vector Poisson equation for ψ is derived: ∇ ψ ¼ −ω:
ð12Þ
When substituting Eq. (10) into Eq. (1) and rewriting the resultant equation by using the relation ∇⋅ (∇× ψ) = 0, the Laplace equation for ϕ is obtained: 2
∇ ϕ ¼ 0:
ð13Þ
When ψ and ϕ are calculated by solving Eqs. (12) and (13) respectively, the velocity ug is computed from Eq. (10). The vorticity ω in Eq. (12) is estimated from Eq. (3). The Vortex in Cell (VIC) method discretizes the vorticity field with vortex elements to calculate the distribution of ω by tracing the convection of each vortex element. It is postulated that the position vector and vorticity for the vortex element v are xv = (xv, yv, zv) and ωv respectively. The Lagrangian form of the vorticity equation, Eq. (3), is written as: dxv ¼ uðxv Þ dt
ð14Þ
dωv 1 2 ¼ ∇⋅ðuðxv Þωðxv ÞÞ þ ν∇ ωðxv Þ− ∇ F D : ρg dt
ð15Þ
When the position and vorticity of vortex element are known at t = t, the values at t = t + Δt are computed from Eqs. (14) and (15). In the VIC method, the flow field is divided into computational grid cells to define ψ, ϕ and ω on the grids. If ω is defined at a position xk = (xk, yk, zk), the vorticity ω is given to xk, or a vortex element with ω is redistributed onto xk
ð6Þ
Considering Eq. (4), the drag force acting on a particle, Fp, is given π 3 f F p ¼ d ρp ug −up : 6 τp
According to the Helmholtz theorem, the gas velocity ug is represented as the summation of the gradient of a scalar potential ϕ and the curl of a vector potential ψ: the gas velocity ug is expressed as:
2
where FD is the force exerted by the particle acting on the gas-phase per unit volume. Taking the curl of Eq. (2) and substituting Eq. (1) into the resultant equation, the vorticity equation for the gas is derived: ∂ω 1 2 þ ∇⋅ ωug ¼ ∇⋅ ug ω þ ν∇ ω− ∇ F D ρg ∂t
2.2. Orthogonal decomposition of velocity and discretization of vorticity
The velocity calculated from Eq. (10) remains unaltered even when any gradient of scalar potential function is added to ψ. To remove this arbitrariness, a solenoidal condition is imposed on ψ:
2. Basic equations and numerical method
∇⋅ug ¼ 0
377
ωðxk Þ ¼
Nv X v
ωv W
x −x y −y z −z k v v W k v W k Δx Δy Δz
ð16Þ
where Nv is the number of vortex elements, and Δx, Δy and Δz are the grid widths. For the redistribution function W, the following equation is employed [10]: 8 2 3 jεjb1 < 1−2:5ε þ 1:5jεj W ðεÞ ¼ 0:5ð2−jεjÞ2 ð1−jεjÞ 1≤jε j≤2 : : 0 jε j > 2
ð17Þ
ð8Þ 2.3. Discretization by using staggered grid
where αp is the particle volume fraction. Since the force FD in Eqs. (2) and (3) is expressed as npFp, the following relation is derived from Eqs. (7) and (8). F D ¼ α p ρp
f ug −up τp
ð9Þ
For incompressible flow simulations, the MAC method and the SMAC method solve the Poisson equation for pressure, which is derived from the equation for pressure gradient and the continuity equation. They employ a staggered grid to ensure the consistency between the discretized equations as well as to prevent the numerical oscillation of the solution.
378
T. Uchiyama / Powder Technology 235 (2013) 376–385
In the VIC method, the Poisson equation for ψ and the Laplace equation for ϕ are derived. Therefore, the use of a staggered grid is considered to be indispensable for their discretization. But the existing VIC method has less explained the grids, and some simulations seem to employ a regular grid. The author proposed a VIC method, which uses a staggered grid for single-phase flow [6]. This simulation also employs the staggered grid illustrated in Fig. 1. The scalar potential ϕ and the particle volume fraction αp are defined at the center. The gas velocity ug is defined at the sides, while the vorticity ω and the vector potential ψ are defined on the edges.
D x y
2.4. Correction of vorticity In the VIC method, the vorticity field is discretized with vortex elements, and the field is expressed by superpositioning the vorticity distribution around each vortex element. The superposition is performed by Eq. (16). The vorticity field obtained by Eq. (16), denoted by ωr, does not always satisfy the solenoidal condition [10]. When the vorticity satisfying the condition is denoted by ωs (=∇× ug), ωr is expressed as [4]:
ωr ¼ ∇F þ ωs
z Fig. 2. Initial arrangement of particles.
ð18Þ
where F is a scalar function. Eq. (18) corresponds to the Helmholtz decomposition of ωr. Taking the divergence of Eq. (18), the Poisson equation for F is obtained: 2
∇ F ¼ ∇⋅ωr :
ð19Þ
If F is calculated from Eq. (19) and substituted into Eq. (18), the vorticity ωs is recalculated [4]. This correction for vorticity needs to solve the Poisson equation, causing the increase of the computational time. The author [6] has proposed a simplified correction method.
The vorticity before the correction, ωr, is connected with ψr through Eq. (12). When taking the curl of Eq. (12) and substituting Eq. (19) into the resultant equation, the following relations are obtained: 2
∇ ð∇⋅ψr Þ ¼ −∇⋅ωr ¼ −∇2 F:
ð20Þ
On the contrary to Eq. (11) assuming the solenoidal condition for ψ, the following equation is derived from Eq. (20) if the vorticity is not solenoidal. ∇⋅ψr ¼ −F:
ð21Þ
i-1/2
xi x i+1/2
x/D
y/D 0
3 -3 0
-3 -1
z k+1/2
0 Particles at t*=0
zk
Δz
z k-1/2
z x
Δx
y
5
Gravity
y j+1/2
yj
z/D
y j-1/2
Δy
Gas velocity Gas vorticity Vector potential Scalor potential Particle volume fraction
ug
vg
ω x , ψx
ωy , ψ y
Fig. 1. Staggered grid.
φ, αp
wg ω z , ψz 10 11 Fig. 3. Computational domain.
3
T. Uchiyama / Powder Technology 235 (2013) 376–385 Table 1 Simulation conditions.
379
t*=0
1.2 kg/m3 1.5 × 10−5 m2/s 1 mm 15 kg/m3 0.0463 s 0.216 m/s 100, 200, 400 0.0125, 0.025, 0.05 20d 0.5D/ut 1 6D × 6D × 12D 40 × 40 × 80 5.4 × 10−3
Gas density ρg Gas kinematic viscosity ν Particle diameter d Particle density ρp Particle velocity response time τp Particle terminal velocity ut Number of particles Np Initial particle volume fraction αp0 Diameter of particle-cluster D Time-scale of gas flow τg Stokes number St Computational domain Number of grids Time increment Δt/τp
2.7
x/D
1
0
-1 -1
0 y/D 1
4.59
6.88
9.18
11.6
When calculating ψr from Eq. (12) with the use of ωr and determining ϕr from Eq. (13), the curl of ur is transformed as follows: ∇ ur ¼ ∇ ð∇ϕr þ ∇ ψr Þ ¼ ∇ð∇⋅ψr Þ−∇2 ψr
ð22Þ
¼ −∇F þ ωr ¼ ωs :
-1
Fig. 5. Time evolution for particle distribution projected onto the horizontal plane (αp0 = 0.05).
t*=0 Eq. (22) demonstrates that the curl of the velocity ur calculated from ωr yields a vorticity ωs satisfying the solenoidal condition. If this recalculation of vorticity, Eq. (22), or the correction of vorticity is performed just after the calculation of Eq. (10), the discretization error of vorticity is completely removed and the time evolution of the flow is accurately simulated without solving the Poisson equation [Eq. (19)]. It should be noted that the staggered grid is indispensable to make the transformation in Eq. (22) applicable to the corresponding discretized equations.
z/D
0
2.7 2
4.59 2.5. Numerical procedure When the flow at time t is known, the flow at t + Δt is simulated by the following procedure:
6.88
4
(1) (2) (3) (4) (5) (6) (7) (8)
9.18 6
11.6
Calculate the particle motion from Eq. (4). Calculate the time variation of ω at every grid point from Eq. (15). Calculate the convection of vortex element from Eq. (14). Calculate ω from Eq. (16). Calculate ψ from Eq. (12). Calculate ϕ from Eq. (13). Calculate ug from Eq. (10). Correct the vorticity, or take the curl of ug to calculate the corrected vorticity.
3. Application to free falling of solid particles in air 3.1. Computational conditions
8 -2
-1
0
1
2
x/D Fig. 4. Time evolution for particle distribution projected onto the vertical plane (αp0 =0.05).
To confirm the validity of the proposed simulation method, the behavior of solid particles falling in an unbounded air is simulated by the method. At the initial time t ∗ (=t/τp) = 0, spherical solid particles are arranged within a spherical region of diameter D in a quiescent
380
T. Uchiyama / Powder Technology 235 (2013) 376–385
-1
air, as shown in Fig. 2. The particles are made to fall, and their fall induces the air flow around them, resulting in a particle-laden air flow. Fig. 3 shows the computational domain of 6D × 6D × 12D. It is resolved into 40× 40× 80 uniform grid cells. The particle position at the commencement of the fall (t∗ = 0) is given by using random numbers. The origin of coordinates is set at the center of the spherical region, and the z-axis is taken in the vertically downward direction. The particle diameter d and density ρp are 1 mm and 15 kg/m 3 respectively. The diameter of spherical region D is set at 20d. To investigate the effect of the number of particles Np, the simulations are performed for Np =100, 200 and 400. Accordingly, the particle volume fraction αp0 at t∗ = 0 ranges from 0.0125 to 0.05. The particle terminal velocity ut is 0.216 m/s. When the time scale of the induced air flow τg is assumed to be 0.5D/ut, the Stokes number St, defined as the ratio of the particle response time τp to τg, is 1. The time increment Δt is set at Δt/τp = 5.4 × 10−3. The computational conditions are listed in Table 1. Since the particles fall in an unbounded air, ϕ is zero and the analysis of Eq. (13) can be omitted. The second-order Adams–Bashforth method is applied to the Lagrangian calculation of the particle motion. The collision between the particles is ignored, because the particle volume fraction is low. For the calculation of the particle volume fraction, αp, the volume of the particles existing in each computational grid cell is computed, and the ratio of the volume against the cell volume is calculated. It should be mentioned that the trial simulations employing the existing vortex in cell method collapsed in the cases of Np = 200 and 400 due to the divergence of the gas velocity field.
t*=0
0
2.7 2
z/D
4.59
6.88
4
9.18 6
3.2. Results for particle volume fraction of αp0 = 0.05
11.6
Fig. 4 shows the time evolution for the particle distribution in the case of Np = 400 or αp0 = 0.05. The particle positions at six time points are projected onto a vertical plane (x–z plane). The particles arranged in the half upper part and half lower one of the spherical region at t∗ = 0 are colored by red and blue respectively. The particles fall with destroying their initial arrangement, namely they tend to disperse in the horizontal direction. Most of the particles initially arranged in the half upper part sediment relatively toward the bottom of the particlecluster, while the particles initially arranged in the half lower part rise relatively toward the top. It is found that the particle behavior depends on the particle initial position. The particle dispersion in the horizontal (x) direction becomes markedly when t∗ > 2.7. These particle motions are caused by the air flow induced by them, as explained later.
8 -2
-1
0
1
2
x/D Fig. 6. Time evolution for particle distribution near the central vertical section (αp0 =0.05).
0
z/D
4
8
t*=2.7 -3
t*=6.88 0
x/D
3
-3
0
x/D
t*=11.6 3
-3
0
x/D
Fig. 7. Time evolution for gas velocity on the central vertical section and particle distribution (αp0 = 0.05).
3
T. Uchiyama / Powder Technology 235 (2013) 376–385
381
z/D
-1
1.4
t*=2.7 3.8 6.2 8.6 11
-3
-3
y/D
0
0
x/D
3
Particle
Q
ωz
Q
ωz
Q
ωz
z/D
-1
1.4
t*=6.88 3.8 6.2 8.6 11
-3
-3
y/D
x/D
0
0 3
Particle z/D
-1
1.4
t*=11.6 3.8
6.2
8.6 11
-3
-3
y/D
0
0
x/D
3
Particle
Fig. 8. Particle distribution and gas vortical structures (αp0 = 0.05).
382
T. Uchiyama / Powder Technology 235 (2013) 376–385
-1 2
t*=0
up /ut
0
z/D
1
0
5
10
2.97
15
t*
2
Fig. 9. Time variation for falling velocity of particles (αp0 = 0.05).
5.13 The time evolution for the particle distribution projected onto a horizontal plane (x–y plane) is shown in Fig. 5, where the distributions at the same instance as Fig. 4 are plotted. Though the particles fall with dispersing in the horizontal direction, they distribute within a circular region. Most of the particles initially arranged in the half lower part, colored by blue, concentrate on the periphery of the circular region when t∗ ≥ 6.88. This is because the particles are whirled up by a vortex ring induced by them, as explained later. Fig. 6 shows the particle distribution near a vertical plane (x–z plane) passing through the center of the particle-cluster (x=y=0). The particles within a vertical slit region of |y/D|≤0.0625 are plotted. When t∗ = 2.7, the particles in the half upper part sediment relatively in the downward direction, and all particles distribute in a hemispherical region. The particle distributions of concave and thin layer types are observed at t∗ = 4.59. The distributions are attributable to the sedimentation of the particles in the half upper part. They are also caused by the fact that the particles near the periphery of the particle-cluster are relatively whirled up. Accordingly, the particles initially arranged in the half lower part of the spherical region, plotted by blue, concentrate on the periphery of the particle-cluster at t∗ ≥6.88. The particles disperse more in the horizontal direction with the passage of time. The abovementioned relative sedimentation, while up and dispersion are caused by the air flow induced by the falling particles. The time evolution for the air velocity on the vertical plane (x–z plane) passing through the center of the particle-cluster (x = y = 0) is
1.5
rp’/D
1
0.5
zp’/D
0 0.5
0.25
0
5
10
t* Fig. 10. Time variation for dispersion of particles (αp0 = 0.05).
15
4
7.69
10.1 6
12.7
8 -2
-1
0
1
2
x/D Fig. 11. Time evolution for particle distribution projected onto vertical plane (αp0 =0.025).
depicted in Fig. 7. The velocity at three time points and the corresponding particle distributions are plotted. When t∗ = 2.7, a downward air flow with higher velocity occurs just above the particle-cluster. The gravitational acceleration, which causes the particles to fall, induces the air flow through the drag force acting on the particles. The relative sedimentation of the particles initially arranged in the half upper part shown in Figs. 4 and 6 is due to this downward air flow. A pair of large-scale eddies appears around the downward air flow. Such organized vortical structure is the wake of the particle-cluster. It is a vortex ring, as explained later. At t ∗ = 6.88, the developed vortex ring whirls up the particles, and accordingly it produces the particle distribution of thin layer type. When t ∗ = 11.6, an air flow in the horizontal direction occurs just below the vortex ring, making the particles disperse in the direction. Experimental studies on a plane mixing layer [11] and a twodimensional wake [12] reported that the particles of the Stokes number St≃1 are focused on thin layers, outlining the boundaries of the largescale eddies. In this simulation, the particles distribute on the periphery of the vortex ring induced by them. Considering that St value of this simulation is 1, such particle distribution corresponds to the abovementioned preferential distribution. According to a two-dimensional simulation for the falling of particles of St=1 by Chen and Marshall [13], an air flow is induced behind the particles, and a concave distribution of the particles occurs due to the wake flow. These support the validity of the present simulation.
T. Uchiyama / Powder Technology 235 (2013) 376–385
383
0
z/D
4
8
t*=2.97 -3
0
t*=7.69 3
-3
x/D
t*=12.7
0
3
-3
x/D
0
3
x/D
Fig. 12. Time evolution for gas velocity on the central vertical section and particle distribution (αp0 = 0.025).
Fig. 8 shows the particle distribution and the air vortical structures. The second invariant of the velocity gradient tensor Q and the vertical (z) component 2of the vorticity ωz are presented with their iso-surfaces of Q = 1=τ p = 0.0016 and ωz/(1/τp) = ± 0.007. A vortex ring occurs just above the particles at t ∗ = 2.7. The scale of the vortex ring increases with the passage of time. Vortex tubes with positive and negative values of ωz entangle behind the particles. One can observe the three-dimensional vortical flow having vertical vorticities. The time evolution for the particle falling velocity ūp is shown in Fig. 9, where ūp is the velocity for the center of gravity of the particles. The velocity ūp increases rapidly just after the commencement of the particle fall. It reaches its maximum value at t ∗ = 3.1. This is attributable to the fact that the particles are accelerated by the downward air flow induced by them. The velocity ūp lessens monotonously at t ∗ > 3.1, and it tends towards the terminal velocity ut. Because the distance between the vortex ring and the particle-cluster increases with the time evolution as seen in Figs. 7 and 8, and accordingly the effect of the vortex ring on the particles reduces. A trial simulation using finer grids (48 × 48 × 96) was also performed, but the particle velocity shown in Fig. 9 was not affected by the grids. Therefore, the present grids are found to be appropriate.
To grasp the particle dispersion, the variances in the horizontal and vertical directions,r ′ p andz′ p respectively, are estimated by the following equations: 0 r
′
p
1 ¼@ Np
11=2 N p 2 2 X A X−xp þ Y−yp
0
11=2 Np 2 X 1 zp¼@ Z−zp A Np p ′
ð24Þ
where (X,Y,Z) is the position vector for the center of gravity of the particles. The time variations for r ′ p =D and z′ p =D are presented in Fig. 10. r ′ p =D increases with the passage of time. A marked increment is observed at t∗ ≥ 2.54. This is because the particles are transported in the horizontal direction by the developing vortex ring. z′ p =D takes its minimum value at t ∗ =3.24 because the particles sediment relatively within the particle-cluster, and therefore the vertical distance between the particles lessens, as found from the particle distribution at t∗ = 2.7 shown
2
up /u t
αp0 =0.05 1
0.025
0
5
ð23Þ
p
0.0125
t*
10
Fig. 13. Effect of αp0 on falling velocity of particles.
15
384
T. Uchiyama / Powder Technology 235 (2013) 376–385
1.5
1
, r p /D
αp0 =0.05 0.025
0.5
0.0125
0
, z p /D
0.5
αp0 =0.05 0.25
0.025 0.0125 0
5
10
15
t* Fig. 14. Effect of αp0 on particle dispersion.
in Fig. 4. z′ p =D increases at 3.24b t ∗ ≤9.85. This is because the particles are whirled up by the vortex ring behind them. 3.3. Effects of initial particle volume fraction αp0 In the case of Np = 200 or αp0 = 0.025, the particles distribute on the vertical plane (x–z plane) as shown in Fig. 11. When compared with the result for αp0 = 0.05 presented in Fig. 4, the particle sedimentation, whirl up and dispersion are smaller. Fig. 12 shows the distributions for particles and air velocity in the case of αp0 = 0.025. The air velocity is lower than that for αp0 = 0.05. The vortex ring is not observed. The effect of αp0 on the particle falling velocity ūp is demonstrated in Fig. 13. The velocity ūp becomes smaller with the decrement of αp0. This is because the drag force on the particles is smaller, and accordingly the downward air flow accelerating the particles in the downward direction has lower velocity. The velocity ūp for αp0 = 0.0125 is almost parallel to the particle terminal velocity ut. One can find that the air flow induced by the particles is markedly weak. The time variation for the particle dispersion is plotted in Fig. 14. The dispersion reduces with the decrement of αp0. In the case of αp0 = 0.0125, the dispersion is so small, indicating that the particles fall with maintaining their initial arrangement. 4. Concluding remarks A numerical simulation method for incompressible gas flow laden with small solid particles is proposed. It is based on a Vortex in Cell method, of which accuracy for single-phase flow simulation was successfully improved by the author in a prior study. The proposed simulation method is applied to compute a free fall of small solid particles in an unbounded air. The particles, initially arranged
within a spherical region in a quiescent air, are made to fall, and their fall induces the air flow around them, resulting in a particle-laden air flow. The particles are accelerated by the induced downward air flow just after the commencement of their fall. After the acceleration, they are whirled up by a vortex ring produced around the downward air flow. These interactions between the particle motion and the air flow are favorably compared with the existing measured and simulated results. In this study, the simulation for solid particles having a relatively low density was carried out. It should be noted that the proposed method is indeed applicable to the simulation for the higher density. This study applied the proposed method to the simulation for gas–particle two-phase flow in an unbounded space. It was confirmed in the author's previous study that the proposed method can perform the DNS for a singlephase turbulent channel flow. Therefore, the method promises to be usefully applied to wall turbulent flow laden with solid particles. Notations CD d D FD f g p Q t t⁎ u ut ūp W
drag coefficient of particle particle diameter diameter of particle cluster at t ∗ = 0 force exerted by particle acting on gas-phase per unit volume Stokes drag factor gravitational constant pressure second invariant of velocity gradient tensor time nondimensional time = t/τp velocity terminal velocity of particle falling velocity of particles redistribution function
T. Uchiyama / Powder Technology 235 (2013) 376–385
x, y, z orthogonal coordinates α volume fraction Δt time increment Δx, Δy, Δz grid width v kinematic viscosity of gas-phase ρ density τp particle response time ϕ scalar potential ψ vector potential ω vorticity = ∇ × ug Subscripts 0 g p x, y, z
initial condition gas particle components in directions of x, y, z
References [1] C.T. Crowe, M. Sommerfeld, Y. Tsuji, Multiphase Flow with Droplets and Particles, CRC Press, 1992. [2] I.P. Christiansen, Numerical simulation of hydrodynamics by the method of point vortices, Journal of Computational Physics 13 (1973) 363–379.
385
[3] G.-H. Cottet, P. Poncet, Advances in direct numerical simulations of 3D wall-bounded flows by Vortex-in-Cell methods, Journal of Computational Physics 193 (2003) 136–158. [4] R. Cocle, G. Winckelmans, G. Daeninck, Combining the vortex-in-cell and parallel fast multipole methods for efficient domain decomposition simulations, Journal of Computational Physics 227 (2008) 9091–9120. [5] P. Chatelain, A. Curioni, M. Bergdorf, D. Rossinelli, W. Andreoni, P. Koumoutsakos, Billion vortex particle direct numerical simulations of aircraft wakes, Computer Methods in Applied Mechanics and Engineering 197 (2008) 1296–1304. [6] T. Uchiyama, Y. Yoshii, H. Hamada, Direct numerical simulation of a turbulent channel flow by an improved vortex in cell method, International Journal of Numerical Methods for Heat and Fluid Flow (in press). [7] T. Uchiyama, Y. Yoshii, Numerical simulation of entrainment and transport of gas bubbles by vortex in cell method, in: Proc. ASME-JSME-KSME Joint Fluids Engineering Conference, Hamamatsu, 2011, (on USB flash drive). [8] J.H. Walther, P. Koumoutsakos, Three-dimensional vortex methods for particle-laden flows with two-way coupling, Journal of Computational Physics 167 (2001) 39–71. [9] L. Schiller, A.Z. Naumann, Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung, Zeitschrift des Vereins Deutscher Ingenieure 77 (1933) 318–321. [10] G.-H. Cottet, P. Koumoutsakos, Vortex Methods: Theory and Practice, Cambridge University Press, 2000. [11] F. Wen, N. Kamalu, J.N. Chung, C.T. Crowe, T.R. Troutt, Particle dispersion by vortex structures in plane mixing layers, Transactions of the American Society of Mechanical Engineers, Journal of Fluids Engineering 114 (1992) 657–666. [12] Y. Yang, C.T. Crowe, J.N. Chung, T.R. Troutt, Experiments on particle dispersion in a plane wake, International Journal of Multiphase Flow 26 (2000) 1583–1607. [13] H. Chen, J.S. Marshall, A Lagrangian vorticity method for two-phase particulate flows with two-way phase coupling, Journal of Computational Physics 148 (1999) 169–198.