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

Authors: Li-Zhong Huang; De-Ming Nie
Journal: Abstract and Applied Analysis
 (2013)

---

## 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.


---

## Body

## 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:
(1)fi(x+eiΔt,t+Δt)=fi(x,t)-M-1S[m-m(eq)]+FΔt,
where S is the diagonal collision matrix, S=diag(0,se,sε,0,sq,0,sq,sv,sv), and M is the transform matrix,
(2)M=[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].The moment space spanned bymi and the velocity space spanned by fi are related by a linear mapping, m=Mf; that is,f=M-1m. The D2Q9 model is used in 2D simulation, and the discrete velocities are listed as the follows:
(3)ci={(0,0)i=0,(±1,0)c,(0,±1)ci=1-4,(±1,±1)c,(±1,±1)ci=5-8,
where c=Δx/Δt is the lattice speed, among which Δx is the lattice spacing, and Δtisthe time step. For D2Q9 model, the corresponding nine moments are given by [19]
(4)m=(ρ,e,ε,jx,qx,jy,qy,pxx,pxy)†,
where ρ is the mean density and jx=ρux, jy=ρuy are the conserved moments. Moreover, the other nonconserved moments are listed as follows:
(5)e(eq)=-2ρ+3j·jρ,ε(eq)=ρ-3j·jρ,qx(eq)=-jx,qy(eq)=-jy,pxx(eq)=(jx2-jy2)ρ,pxy(eq)=jxjyρ.According to Guo et al. [20], the forcing term F is defined as follows:
(6)F=M-1(I-12S)MF-,Fi¯=wi[ci·αcs2+uα:cicics4-uα:cs2Ics4],
where α is the acceleration due to the external force and wi are weights related to the lattice model which are chosen as w0=4/9; wi=1/9, i=1~4; wi=1/36, i=5~8. The speed of sound cs is equal to c/31/2. By applying the Taylor expansion techniques and the Chapman-Enskog analysis, (1) leads to the hydrodynamic equations that are shown as follows:
(7)∇·u=0,ρf(∂u∂t+(u·∇)u)=-∇p+μ∇2u+a.The shear viscosity and the bulk viscosity can be defined as follows [21]:
(8)ν=13(1sν-12)cΔx,ζ=13(1se-12)cΔx.


### 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
(9)Fi,jP={0di,j>Ri+Rj+ξ,1εP′(Xi-Xj)(Ri+Rj-di,j)di,j≤Ri+Rj,1εP(Xi-Xj)(Ri+Rj+ξ-di,j)2Ri+Rj≤di,j≤Ri+Rj+ξ,
where Xi and Xj are the mass center coordinates of the ith and jth particle, Ri and Rj are their radius, di,j is the distance between their mass centers which equals to |Xi-Xj|, ξ is the force range which is usually set to be one or two lattice spacing, and εp′ and εp 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:(10)Fi,jW={0di′>2Ri+ξ,1εW′(Xi-Xi′)(2Ri-di′)di′≤2Ri,1εW(Xi-Xi′)(2Ri+ξ-di′)22Ri≤di′≤2Ri+ξ,
where Xi′ is the coordinate of the nearest imaginary particle located on the boundary nearby the ith particle, di′ is the distance between the imaginary particle and the real one which equals to |Xi-Xi′|, and εW′ and εW are the other two small positive stiffness parameters for particle-wall collisions. In this paper, εp′, εp, εW′, and εW are all set to be 1e-7.


## 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:
(1)fi(x+eiΔt,t+Δt)=fi(x,t)-M-1S[m-m(eq)]+FΔt,
where S is the diagonal collision matrix, S=diag(0,se,sε,0,sq,0,sq,sv,sv), and M is the transform matrix,
(2)M=[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].The moment space spanned bymi and the velocity space spanned by fi are related by a linear mapping, m=Mf; that is,f=M-1m. The D2Q9 model is used in 2D simulation, and the discrete velocities are listed as the follows:
(3)ci={(0,0)i=0,(±1,0)c,(0,±1)ci=1-4,(±1,±1)c,(±1,±1)ci=5-8,
where c=Δx/Δt is the lattice speed, among which Δx is the lattice spacing, and Δtisthe time step. For D2Q9 model, the corresponding nine moments are given by [19]
(4)m=(ρ,e,ε,jx,qx,jy,qy,pxx,pxy)†,
where ρ is the mean density and jx=ρux, jy=ρuy are the conserved moments. Moreover, the other nonconserved moments are listed as follows:
(5)e(eq)=-2ρ+3j·jρ,ε(eq)=ρ-3j·jρ,qx(eq)=-jx,qy(eq)=-jy,pxx(eq)=(jx2-jy2)ρ,pxy(eq)=jxjyρ.According to Guo et al. [20], the forcing term F is defined as follows:
(6)F=M-1(I-12S)MF-,Fi¯=wi[ci·αcs2+uα:cicics4-uα:cs2Ics4],
where α is the acceleration due to the external force and wi are weights related to the lattice model which are chosen as w0=4/9; wi=1/9, i=1~4; wi=1/36, i=5~8. The speed of sound cs is equal to c/31/2. By applying the Taylor expansion techniques and the Chapman-Enskog analysis, (1) leads to the hydrodynamic equations that are shown as follows:
(7)∇·u=0,ρf(∂u∂t+(u·∇)u)=-∇p+μ∇2u+a.The shear viscosity and the bulk viscosity can be defined as follows [21]:
(8)ν=13(1sν-12)cΔx,ζ=13(1se-12)cΔx.


## 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
(9)Fi,jP={0di,j>Ri+Rj+ξ,1εP′(Xi-Xj)(Ri+Rj-di,j)di,j≤Ri+Rj,1εP(Xi-Xj)(Ri+Rj+ξ-di,j)2Ri+Rj≤di,j≤Ri+Rj+ξ,
where Xi and Xj are the mass center coordinates of the ith and jth particle, Ri and Rj are their radius, di,j is the distance between their mass centers which equals to |Xi-Xj|, ξ is the force range which is usually set to be one or two lattice spacing, and εp′ and εp 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:(10)Fi,jW={0di′>2Ri+ξ,1εW′(Xi-Xi′)(2Ri-di′)di′≤2Ri,1εW(Xi-Xi′)(2Ri+ξ-di′)22Ri≤di′≤2Ri+ξ,
where Xi′ is the coordinate of the nearest imaginary particle located on the boundary nearby the ith particle, di′ is the distance between the imaginary particle and the real one which equals to |Xi-Xi′|, and εW′ and εW are the other two small positive stiffness parameters for particle-wall collisions. In this paper, εp′, εp, εW′, and εW 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 Figure1 which is constructed by four solid walls, lx=ly=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 U0 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, ρp=ρf=1.0 g/cm3. The kinematic viscosity of the fluid is set to be 0.01 cm2 s−1. 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, whent′=t/104s≤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 t′>1.5, those particles begin to move in the suspension fluid. At first, U0 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.Vorticity contours whenU0=0.1 m/s.
(a)
t
′
=
1.5
(b)
t
′
=
3.0
(c)
t
′
=
4.0
(d)
t
′
=
5.5
(e)
t
′
=
6.5
(f)
t
′
=
7.0
(g)
t
′
=
8.0
(h)
t
′
=
9.0
(i)
t
′
=
10.0Vorticity contours whenU0=0.2 m/s.
(a)
t
′
=
1.5
(b)
t
′
=
3.0
(c)
t
′
=
4.0
(d)
t
′
=
5.3
(e)
t
′
=
6.6
(f)
t
′
=
7.0
(g)
t
′
=
7.4
(h)
t
′
=
7.8
(i)
t
′
=
8.2Vorticity contours whenU0=0.3 m/s.
(a)
t
′
=
1.5
(b)
t
′
=
3.0
(c)
t
′
=
4.0
(d)
t
′
=
5.3
(e)
t
′
=
6.6
(f)
t
′
=
7.0
(g)
t
′
=
7.4
(h)
t
′
=
7.8
(i)
t
′
=
8.2Instantaneous vorticity contours forU0=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 t′=8.2, the vortex contour is rotational symmetry from the center of the simulation domain.But whenU0=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 t′>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 U0 varies from 0.1 m/s to 0.3 m/s; however, the number of particle collisions varies little when U0 is large enough which can be seen from Figure 6.Figure 6
Comparison of collision number with differentU0.



---

*Source: 131025-2013-11-25.xml*