Lattice Boltzmann Simulation of Collision between 2D Circular Particles Suspension in Couette Flow

Abstract

Collision between 2D circular particles suspension in Couette flow is simulated by using multiple-relaxation-time based lattice Boltzmann and direct forcing/fictitious domain method in this paper. The patterns of particle collisions are simulated and analyzed in detail by changing the velocity of top and bottom walls in the Couette flow. It can be seen from the simulation results that, while the velocity is large enough, the number of collisions between particles will change little as this velocity varies.

1. Introduction

Multiphase flow is a very important branch of fluid mechanics, while fluid-solid two-phase flow is the main part of such flow. And fluid-solid two-phase flow is very common in nature and industry, such as raindrop formation, material science, chemical industry, aerosol deposition, fluidized beds, and injection molding machine [1, 2]. Particle collisions, which will influence the performance of product, happen occasionally in those processes. So it is very meaningful to do some research to understand particle collisions and then take control to improve the performance of product.

Even with the help of the precise and advanced particle image velocity (PIV) instrument, it is very difficult to observe the phenomenon of particle collisions in detail. Comparing the experimental measurement, numerical simulation has great advantages to investigate particle collisions in fluid-solid two-phase flow, especially for direct numerical simulation (DNS) methods. The lattice Boltzmann method (LBM), one of the best DNS methods with several remarkable advantages, was first proposed by Ladd [3] and then improved by Aidun et al. [4] to simulate particles suspended in a viscous fluid. Feng and Michaelides [5] and Tian et al. [6] united immersed boundary method (IBM) into LBM to deal with fluid-solid interface problem. In the past several decades the LBM method was proved robust and efficient for particulate flows, especially in the case of large number of particles [7–12]. Nie and Lin [13] developed a single-relaxation-time (SRT) based lattice Boltzmann-direct forcing/fictitious domain (SRT LB-DF/FD) method to simulate particle suspensions and then improved to multiple-relaxation-time (MRT LB-DF/FD).

Several papers investigated particle collisions in three-dimensional homogeneous isotropic turbulence [14–18] but not by using DNS method. In this paper, DNS method is adopted to simulate particle collisions. Because the MRT model has better ability of computing pressure and more time saving than the SRT model, the MRT LB DF/FD method is utilized. Firstly, this method is introduced in detail in Section 2. Secondly, simulation problem is described in Section 3. And finally, several simulation results and conclusions are presented in Section 4.

2. Numerical Method

2.1. The MRT LB DF/FD Method

The lattice Boltzmann method based on the multiple-relaxation-time (MRT) collision model is adopted in this paper [13]. The discrete equations can be written as follows: 𝑓_𝑖𝐱+𝐞_𝑖Δ𝑡,𝑡+Δ𝑡=𝑓_𝑖(𝐱,𝑡)−𝑀⁻¹𝑆𝐦−𝐦^(eq)+𝐅Δ𝑡, (1) where 𝑆 is the diagonal collision matrix, 𝑆=diag(0,𝑠_𝑒,𝑠_𝜀,0,𝑠_𝑞,0,𝑠_𝑞,𝑠_𝑣,𝑠_𝑣), and 𝑀 is the transform matrix, 𝑀=111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1100−20211−1−101−11−10000000001−11−1. (2)

The moment space spanned by 𝐦_𝑖 and the velocity space spanned by 𝑓_𝑖 are related by a linear mapping, 𝐦=𝑀𝑓; that is,𝑓=𝑀⁻¹𝐦. The D2Q9 model is used in 2D simulation, and the discrete velocities are listed as the follows: 𝐜_𝑖=()()()()()0,0𝑖=0,±1,0𝑐,0,±1𝑐𝑖=1−4,±1,±1𝑐,±1,±1𝑐𝑖=5−8, (3) where 𝑐=Δ𝑥/Δ𝑡 is the lattice speed, among which Δ𝑥 is the lattice spacing, and Δ𝑡 is the time step. For D2Q9 model, the corresponding nine moments are given by [19] 𝐦=𝜌,𝑒,𝜀,𝑗_𝑥,𝑞_𝑥,𝑗_𝑦,𝑞_𝑦,𝑝_𝑥𝑥,𝑝_𝑥𝑦^†, (4) where 𝜌 is the mean density and 𝑗_𝑥=𝜌𝑢_𝑥, 𝑗_𝑦=𝜌𝑢_𝑦 are the conserved moments. Moreover, the other nonconserved moments are listed as follows: 𝑒^(eq)=−2𝜌+3𝐣⋅𝐣𝜌,𝜀^(eq)=𝜌−3𝐣⋅𝐣𝜌,𝑞^(eq)_𝑥=−𝑗_𝑥,𝑞_𝑦^(eq)=−𝑗_𝑦,𝑝^(eq)_𝑥𝑥=𝑗²_𝑥−𝑗²_𝑦𝜌,𝑝^(eq)_𝑥𝑦=𝑗_𝑥𝑗_𝑦𝜌. (5)

According to Guo et al. [20], the forcing term 𝐅 is defined as follows: 𝐅=𝑀⁻¹1𝐼−2𝑆𝑀𝐅,𝐹_𝑖=𝑤_𝑖𝐜_𝑖⋅𝜶𝑐²_𝑠+𝐮𝜶∶𝐜_𝑖𝐜_𝑖𝑐⁴_𝑠−𝐮𝜶∶𝑐²_𝑠𝐼𝑐⁴_𝑠, (6) where 𝜶 is the acceleration due to the external force and 𝑤_𝑖 are weights related to the lattice model which are chosen as 𝑤₀=4/9; 𝑤_𝑖=1/9, 𝑖=1∼4; 𝑤_𝑖=1/36, 𝑖=5∼8. The speed of sound 𝑐_𝑠 is equal to 𝑐/3^1/2. By applying the Taylor expansion techniques and the Chapman-Enskog analysis, (1) leads to the hydrodynamic equations that are shown as follows: ∇⋅𝐮=0,𝜌_𝑓𝜕𝐮+()𝐮𝜕𝑡𝐮⋅∇=−∇𝑝+𝜇∇²𝐮+𝐚. (7)

The shear viscosity and the bulk viscosity can be defined as follows [21]: 1𝜈=31𝑠_𝜈−121𝑐Δ𝑥,𝜁=31𝑠_𝑒−12𝑐Δ𝑥. (8)

2.2. Collisions Model

A collision model is needed to avoid particles overlapping. Then the short-range repulsive force developed by Wan and Turek [22] is utilized in this paper because it is easily carried out in the code. For particle-particle collisions, the short-range repulsive force is calculated by 𝐹^𝑃_𝑖,𝑗=0𝑑_𝑖,𝑗>𝑅_𝑖+𝑅_𝑗1+𝜉,𝜀^′_𝑃𝐗_𝑖−𝐗_𝑗𝑅_𝑖+𝑅_𝑗−𝑑_𝑖,𝑗𝑑_𝑖,𝑗≤𝑅_𝑖+𝑅_𝑗,1𝜀_𝑃𝐗_𝑖−𝐗_𝑗𝑅_𝑖+𝑅_𝑗+𝜉−𝑑_𝑖,𝑗²𝑅_𝑖+𝑅_𝑗≤𝑑_𝑖,𝑗≤𝑅_𝑖+𝑅_𝑗+𝜉, (9) where 𝐗_𝑖 and 𝐗_𝑗 are the mass center coordinates of the 𝑖th and 𝑗th particle, 𝑅_𝑖 and 𝑅_𝑗 are their radius, 𝑑_𝑖,𝑗 is the distance between their mass centers which equals to |𝐗_𝑖−𝐗_𝑗|, 𝜉 is the force range which is usually set to be one or two lattice spacing, and 𝜀^′_𝑝 and 𝜀_𝑝 are two small positive stiffness parameters for particle-particle collisions.

But for particle-wall collisions, the corresponding short-range repulsive force is determined as follows: 𝐹^𝑊_𝑖,𝑗=0𝑑^′_𝑖>2𝑅_𝑖1+𝜉,𝜀^′_𝑊𝐗_𝑖−𝐗^′_𝑖2𝑅_𝑖−𝑑^′_𝑖𝑑^′_𝑖≤2𝑅_𝑖,1𝜀_𝑊𝐗_𝑖−𝐗^′_𝑖2𝑅_𝑖+𝜉−𝑑^′_𝑖²2𝑅_𝑖≤𝑑^′_𝑖≤2𝑅_𝑖+𝜉, (10) where 𝐗^′_𝑖 is the coordinate of the nearest imaginary particle located on the boundary nearby the 𝑖th particle, 𝑑^′_𝑖 is the distance between the imaginary particle and the real one which equals to |𝐗_𝑖−𝐗^′_𝑖|, and 𝜀^′_𝑊 and 𝜀_𝑊 are the other two small positive stiffness parameters for particle-wall collisions. In this paper, 𝜀^′_𝑝, 𝜀_𝑝, 𝜀^′_𝑊, and 𝜀_𝑊 are all set to be 1e-7.

3. Problem Description

The main purpose of this paper is to study collisions between 2D circular particles suspension in Couette flow. The schematic diagram of geometry is shown in Figure 1 which is constructed by four solid walls, 𝑙𝑥=𝑙𝑦=8 cm, and mesh element is 640000. The left and right walls are fixed, while the top wall moves in the right direction with a velocity of 𝑈₀ and the bottom wall moves in the left direction at the same velocity. 128 circular particles (8 rows, 16 columns) are placed in the middle of this domain, and the distance of center of two neighborhood particles in horizontal and vertical direction is two times of the diameter of particle. In order to ignore the influence of gravity, the density of particles is set to be same as the suspension fluid, 𝜌_𝑝=𝜌_𝑓=1.0 g/cm³. The kinematic viscosity of the fluid is set to be 0.01 cm² s⁻¹. For each particle, the diameter is 0.2 cm and the arrangement of Lagrangian points is shown in Figure 2.

Figure 1 

Schematic diagram of geometry.

Figure 2 

Arrangements of Lagrangian points.

4. Simulation Results and Conclusions

In this paper, when 𝑡^′=𝑡/10⁴s≤1.5, the circular particles are fixed to generate a fully developed Couette flow which can be seen from Figures 3(a), 4(a), and 5(a). And when 𝑡^′>1.5, those particles begin to move in the suspension fluid. At first, 𝑈₀ is set to be very small which equals to 0.1 m/s and the instantaneous vorticity contours can be seen from Figures 3(a)–3(i). Obviously, the vorticity is not strong enough which makes the number of collisions very big.

(a) 𝑡′=1.5

(b) 𝑡′=3.0

(c) 𝑡′=4.0

(d) 𝑡′=5.5

(e) 𝑡′=6.5

(f) 𝑡′=7.0

(g) 𝑡′=8.0

(h) 𝑡′=9.0

(i) 𝑡′=10.0

Figure 3 

Vorticity contours when 𝑈0=0.1 m/s.

(a) 𝑡′=1.5

(b) 𝑡′=3.0

(c) 𝑡′=4.0

(d) 𝑡′=5.3

(e) 𝑡′=6.6

(f) 𝑡′=7.0

(g) 𝑡′=7.4

(h) 𝑡′=7.8

(i) 𝑡′=8.2

Figure 4 

Vorticity contours when 𝑈0=0.2 m/s.

(a) 𝑡′=1.5

(b) 𝑡′=3.0

(c) 𝑡′=4.0

(d) 𝑡′=5.3

(e) 𝑡′=6.6

(f) 𝑡′=7.0

(g) 𝑡′=7.4

(h) 𝑡′=7.8

(i) 𝑡′=8.2

Figure 5 

Vorticity contours when 𝑈0=0.3 m/s.

Instantaneous vorticity contours for 𝑈₀=0.2 m/s at several interval times are shown in Figures 4(a)–4(i). It can be seen from those figures that the particles are moving by the influence of the two main symmetry vortices. Even when 𝑡^′=8.2, the vortex contour is rotational symmetry from the center of the simulation domain.

But when 𝑈₀=0.3 m/s, Figures 5(a)–5(i) show that the circular particles are agglomerated near the four sidewalls by the influence of the two main vortices. And when 𝑡^′>7.8, the particles are no longer rotational symmetry from the center of the simulation domain, distributed at completely random. The vortex patterns are totally different when 𝑈₀ varies from 0.1 m/s to 0.3 m/s; however, the number of particle collisions varies little when 𝑈₀ is large enough which can be seen from Figure 6.

Figure 6 

Comparison of collision number with different 𝑈0.

Acknowledgment

This work is supported by the National Natural Science Foundation of China with Grant 11132008.