Simulation of flow and soot particle distribution in wall-flow DPF based on lattice Boltzmann method

Accepted Manuscript Simulation of Flow and Soot Particle Distribution in Wall-Flow DPF Based on Lattice Boltzmann Method Xiangjin Kong, Zhijun Li, Boxi Shen, Yue Wu, Yanke Zhang, Dong Cai PII: DOI: Reference:

S0009-2509(19)30254-4 https://doi.org/10.1016/j.ces.2019.03.039 CES 14862

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

30 August 2018 5 March 2019 17 March 2019

Please cite this article as: X. Kong, Z. Li, B. Shen, Y. Wu, Y. Zhang, D. Cai, Simulation of Flow and Soot Particle Distribution in Wall-Flow DPF Based on Lattice Boltzmann Method, Chemical Engineering Science (2019), doi: https://doi.org/10.1016/j.ces.2019.03.039

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Simulation of Flow and Soot Particle Distribution in Wall-Flow DPF Based on Lattice Boltzmann Method Xiangjin Kong, Zhijun Li, Boxi Shen, Yue Wu, Yanke Zhang, Dong Cai State Key Laboratory of Engines, Tianjin University, Tianjin 300072, China Abstract A two-dimensional mesoscopic gas-solid two-phase flow model has been developed to investigate the flow and soot loading in the micro-channel of diesel particulate filters. Soot particle size examined is in the range of 10nm to 1µm. The flow is solved by an incompressible lattice Boltzmann model and the transport of solid particles is described by the cell automation probabilistic model. The lattice Boltzmann-cell automation probabilistic model (LB-CA model) is validated with the results of previous studies. The effects of different upstream velocities on the flow field in channels are investigated. The distribution and deposition of soot particles with different sizes in clean channels are simulated based on the LB-CA method and the LBLagrangian method respectively. The effects of deposited soot particles on flow field are evaluated in real soot particle capture process. The results show that the distributions of velocity field and pressure field in the channel are significantly affected by the upstream velocity. Compared with the effect of the particle size, the upstream velocity is more influential on the particle deposition distributions. The profiles of deposition distribution from the LB-CA method are in close agreement with those from the LB-Lagrangian method. The deposition distributions of particles with different diameters at the top of the porous wall are similar to the distributions of wall velocity along the channel length. Generally, the deposited soot particles increase the axial pressure and decrease the axial velocity in the inlet channel. The evolution trend of the areas where wall velocity undergoes changes is consistent with that of the solid nodes made of the captured soot particles.



Corresponding author at: State Key Laboratory of Engines, Tianjin University, 92 weijin Road, 300072 Tianjin, China. E-mail address: [email protected] (Zhijun Li).

1

Keywords: lattice Boltzmann method, cell automation probabilistic model, Lagrangian model, diesel particulate filter, soot distribution, porous media. 1. Introduction Due to the favorable heat and mass transfer characteristics, porous media is extensively used in scientific and industrial applications, such as fuel cells, heat exchangers and particulate filters (Yuan and Sundén, 2014; Guan et al., 2015). In particulate filter engineering applications, gas-solid two-phase flow through thin layers of porous media has been of continuing interest. As an example, Diesel Particulate Filters (DPFs), which were developed in the 1980s for controlling particulate emission from diesel engines, have become an indispensable aftertreatment device for diesel engines up to now, to satisfy the increasingly stringent emission regulations on particulate matter (PM) and particulate number (PN) (Guan et al., 2015; Torregrosa et al., 2017; Orihuela et al., 2018). Wall-flow DPFs are the dominant type of DPFs for PM (PM also known as soot) reduction with a quite impressive efficiency, which is in excess of 95% in mass and >99% in number over a wide range of engine operating conditions (Guan et al., 2015; Mohankumar and Senthilkumar, 2017). Wall-flow DPF is ceramic monolith of honeycomb multi-channel structure with alternatively blocked inlet and outlet channels. The exhaust gas is forced to flow through the porous walls of channels and soot particles are captured in blocked inlet channels, shown in Fig. 1. As the filter operates, the collected soot particles increase the backpressure of the diesel engine and need to be periodically removed by regeneration process. It becomes evident that the optimization of wall-flow DPFs requires detailed information on filtration process and regeneration process. The filtration process is particularly important since it provides the initial condition to the regeneration process. Therefore, profound understanding of fundamental filtration dynamics is required to assure the optimal design and operation of DPF systems.

2

Figure 1. Structure scheme and working principle of wall-flow DPF. The filtration process of the DPF has been extensively investigated by different numerical approaches. Using the Eulerian-Eulerian approach to investigate the filtration process, which is first proposed by Konstandopoulos et al. (1989, 2000), has grown in popularity and acceptance for lower computational costs. In this approach, the porous wall is virtual and represented by spherical unit cells with a mean pore size and a mean porosity. These spherical unit cells are evenly distributed on porous wall grids, which is unable to reflect the true structure of porous wall with a pore size distribution. Due to this virtual structure of porous wall, some important phenomena in the filtration process cannot be simulated, such as bridge formation inside the porous channels (Sanui and Hanamura, 2016) and hills-like structures formation around the border region of a superficial pore (Choi and Lee, 2013). The pressure drop caused by the deposited soot particles is also virtual, which must be coupled with the filtration process. The filtration process regimes (deep bed filtration regime, transition regime and cake filtration regime) are also divided artificially by some empirical parameters for the same reason, and how to divide the filtration process regimes has been the focus of the Eulerian-Eulerian approach optimization. Konstandopoulos et al. (1989, 2000) define a partition coefficient of the first slab to control the transition from deep bed to cake filtration regime, which means that the cake layer is completely controlled by the porous wall properties instead of the growth dynamic of the soot deposits collected on it. Mohammed et al. (2006) consider the effects of the soot deposits collected in the first slab on cake layer 3

formation by introducing the filtration efficiency of a cake layer. However, the maximum value of the filtration efficiency must be set from empirical data. According to the experimental observations, Bensaid et al. (2009, 2010) set soot packing limit to control the transition from deep bed regime to cake filtration regime without consideration of the effects of the soot deposits collected in the first slab on the cake layer formation. Tandon et al. (2010) divide the deep bed filtration regime into two stages (in fact they can be viewed as deep bed filtration regime and transition regime) by setting the transition permeability and a critical value of collectors. Although this approach is shown to calculate filtration efficiency accurately, it is not coupled with a pressure drop model. Bollerhoff et al. (2012) develop the Eulerian-Eulerian model for nonhomogeneous porous wall structures. In their model, three filtration process regimes are divided by a critical limit of porosity and the formation of a dense cake layer structure. Wall saturation index is used to control soot collection in deep bed filtration regime. Gong and Rutland (2015) develop a heterogeneous multiscale filtration (HMF) model, in which a pore size probability density function is introduced to represent heterogeneity and multiscale characteristics of the porous wall. In their model, partition coefficient (Gong et al., 2017), similar to that in Konstandopoulos' model (Konstandopoulos and Johnson, 1989; Konstandopoulos et al., 2000), is used to divide the filtration process regimes. In order to keep low soot penetration thickness inside the porous wall, Serrano et al. (2016) divide the porous wall into two layers. The upper layer is used for filtration of soot particles and the other is kept clean in the filtration process. In their model, three filtration process regimes are divided by saturation coefficient with different values. Some empirical parameters are set in the simulation, such as the shape factor and the surface correcting factor, for obtaining the same phenomenon appeared in experimental findings. From these developments in Eulerian-Eulerian models, it is found that due to the virtual porous structure, more and more improvements on porous wall structures and the filtration process are required to approximate the actual filtration process. Recent studies indicate that the interaction between particulates and the porous substrate has a significant impact on the filtration performance (Bollerhoff et al., 2012; Swanson et al., 2013; Yamamoto and Sakai, 2016; Meng et al., 2016). Accordingly, the non-virtual porous structure of the filter substrate is critical for profound understanding of fundamental filtration dynamics. 4

Moreover, the Eulerian-Eulerian approach is based on Navier-Stokes equations for gas phase simulation. For soot particles smaller than 1µm, the non-continuous molecular structure of the flow should be considered (Sbrizzai et al., 2005). The Knudsen (Kn) number, which is the ratio between the mean free path of the gas molecules and the particle diameter, is a characteristic parameter to measure whether the continuum assumption of Navier-Stokes equations is reasonable or not. Fig. 2 describes different regimes of fluid flow in channels depending on the Knudsen number (Roy et al., 2003). When the temperature is 600K, the mean free path value (for air) is 0.179µm (Konstandopoulos et al., 2000). The value of the Knudsen number is in the range of 0.17917.9 for particles from 1µm to 10nm, which is beyond the validity regime of the Navier-Stokes equations. The same issues of the virtual porous structure and the continuum assumption of Navier-Stokes equations also exist in the Eulerian-Lagrangian approach (Sbrizzai et al., 2005; Liu et al., 2009; Wu et al., 2011). In the EulerianLagrangian approach, it is assumed that soot particles will stick to the surface of the porous wall once they reach the surface (Sbrizzai et al., 2005; Liu et al., 2009; Wu et al., 2011). Particular difficulties arise from the irregular surface caused by the non-virtual porous structure and the deposited particles. The EulerianLagrangian approach thus has to use high-resolution and adaptive grids near the non-virtual porous structure and the deposit, on the cost of computational efficiency and precision. Besides, the Eulerian-Lagrangian approach is also not easily applicable to the simulation of cake formation in the filtration process (Bensaid et al., 2009), although the Lagrangian approach is a suitable way to model particle behavior with a higher precision, especially in simulating the particle trajectories (Wang et al., 2013). Quite recently, the lattice Boltzmann method (LBM) has been widely used for alternative fluid simulation method in DPF (Yamamoto and Sakai, 2016; Yamamoto et al., 2006, 2007, 2009; Yamamoto and Yamauchi, 2013; Yamamoto and Sakai, 2015). There is no Knudsen number limit to lattice Boltzmann equation for channel flow shown in Fig. 2. Main advantages of LBM are the simplicity of the algorithm and the flexibility for complex geometries of the porous wall. However, the pressure variables, which satisfy

p  c2  / 3

p

and the density



are not two independent

(Yamamoto and Sakai, 2016; Yamamoto et al., 2006, 2007, 2009; Yamamoto

and Yamauchi, 2013; Yamamoto and Sakai, 2015). The particle speed 5

c

is usually a fixed value in LBM. When

pressure drop of DPF is less than 15kPa, the gas can be viewed as incompressible flow (Konstandopoulos, 2003). It requires that the pressure instead of

p

should be an independent variable. Moreover, the soot concentration,

the individual soot particle, is considered in filtration process (Yamamoto and Sakai, 2016;

Yamamoto et al., 2006, 2007, 2009; Yamamoto and Yamauchi, 2013; Yamamoto and Sakai, 2015). Therefore, the motion information of soot particles with different particle sizes is not available. The particle deposition model (the cellular automata model) proposed by Chopard et al. (2000) is used in their simulation. Although the cellular automata model is convenient to implement dynamic boundary conditions and complicated microscopic mechanisms of solid particles (e.g., deposition, toppling, and erosion) (Chopard et al., 2000; Masselot and Chopard, 1998; Chopard and Masselot, 1999) by some simple and intuitive rules, it is only capable to describe particle behaviors on the level of qualitative simulation, for some empirical or tentative formulations or parameters are introduced in those rules (Wang et al., 2012). Lattice Boltzmann Equation

Euler Eqns.

0

Navier-Stokes Equations No-slip

Slip-conditions

10-3

Kn Continuum Flow

Burnett Equation

10-2

10-1

100

101

Kn

∞

Transition Flow Free-molecule Flow

Slip Flow

Figure 2. Knudsen number regimes. The objective of this study is to investigate particulate filtration process with consideration of the nonvirtual porous structure by using a fundamental computational model. In order to deal with the complex boundary of the non-virtual porous walls and dynamic boundary conditions of soot deposits in the filtration process, an incompressible lattice Boltzmann model is introduced to simulate gas phase (Guo et al., 2000). In this gas phase model, pressure variable

p

is a truly independent variable, which is convenient to reflect the

effect of the deposited soot particles on the pressure field. An improved cell automation probabilistic model (the CA model) is used to simulate solid phase, which is able to quantitatively describe particle transport under combined effects of drag forces and Brownian diffusion (Wang et al., 2012). To acquire the motion information of soot particles with different particle sizes, soot particles, rather than soot concentration, are considered in the 6

CA model. The LB-CA model is first validated by the results of recent numerical studies, in which the drag, Saffman lift, gravity and Brownian forces are considered (Jafari et al., 2010; Brandon and Aggarwal, 2001; Afrouzi et al., 2015; Tong et al., 2017). The porous wall is generated by quartet structure generation set (QSGS) (Wang et al., 2007). Then, the effects of different upstream velocities on the pressure field and velocity field are investigated. The distribution and deposition of particles with different sizes in clean channels are simulated based on the LB-CA method and the LB-Lagrangian method respectively. Lastly, the effects of the deposited soot particles on flow field are evaluated in real soot particle capture process. 2. Model description 2.1 The lattice Boltzmann method for fluid flow simulation LBM is based on lattice Boltzmann equation (LBE), which can be directly derived from the discretized Boltzmann equation in both time and phase space, shown in Eq. (1) (He and Luo, 1997). f  x  e t , t   t   f  x, t   

(1) In Eq. (1),

f  x, t 

is the single-particle distribution function at computing node

x

at time

t

,

microscopic velocity,  t is time step. f  x  e   t, t   t  is of a state after advection and changes due to

e

is the

. 

is

collision function and can be calculated by using various approximations, such as the linearization of the collision operator and the Bhatnagar-Gross-Krook (BGK) approximation. In the following analysis, the lattice Boltzmann equation with the BGK approximation is used as the evolution equation of distribution function, shown in Eq. (2). f  x  e t , t   t   f  x, t   

1



 f  x, t   f  x, t   eq

(2) Here  is the relaxation time due to collision, and

f eq  x, t 

is the Maxwell equilibrium distribution function. 7

The lattice Boltzmann method consists of three ingredients: a lattice mode (a discrete velocity model), the evolution equation and equilibrium distribution functions. The lattice mode with b velocities on a simple cubic lattice of dimension d (denoted DdQb) (Qian et al., 1992) is widely accepted as basic lattice of LBM. The present work employs the two-dimensional nine-velocity lattice model (D2Q9) shown in Fig. 3. e2

e6

e5

e3

e0

e7

e1

e8

e4

Figure 3. D2Q9 model of lattice Boltzmann method. And its lattice velocities are given by  0,0                                                      0   e  c cos   1  / 2 ,sin   1  / 2                 1  4   2c cos  2  9   / 4  ,sin  2  9   / 4           5  8









(3) where subscript  indicates the velocity direction, the particle speed

c

is defined as

c   x / t

and  x is the

lattice spacing. After the discrete procedure of velocity, the evolution equation of the lattice Boltzmann equation with the BGK collision approximation is expressed as Eq. (4). f  x  eα tt , t   t   f  x, t   

1



 f  x, t   f  x, t    0,...,8 

eq



(4) And the equilibrium distribution function

feq  x, t 

is defined as (Guo et al., 2000):

8

p   0  4d 0 c 2  0 s0  u          =0   p feq = d1 2  0 s  u                 1  4  c  p d 2 2  0 s  u                 5  8  c

(5) where

0

is a constant satisfying 0   f eq , 

s  u    3 eα  u  / c 2  4.5  eα  u  / c 4  1.5u2 / c 4    2

, with the weight

coefficient: 4 / 9        0   =1 / 9        1  4 , d 0 , d1 , d 2 1 / 36      5  8 

are parameters satisfying

d1  d2  d0 , d1  d2  1 / 2 .

The fluid density, velocity, and

pressure are obtained from the distribution function:    f 

(6) u   ceα f 

(7) p  0

c2   f  s0  u   4d0   0

(8) The lattice Boltzmann model described prior has two advantages: firstly, it completely eliminates the compressibility effect that lies in other existing lattice Boltzmann models and can be used for both steady and unsteady flows. And secondly,

p

is a truly independent variable and



almost remains constant. Through

multiscaling expansion, the incompressible Navier-Stokes equations can be derived from this incompressible model (Guo et al., 2000).

9

2.2 Cell automation probabilistic approach for particle transport The cell automation probabilistic approach is proposed by Masselot and Chopard (Chopard et al., 2000; Masselot and Chopard, 1998; Chopard and Masselot, 1999), in which, the discrete point-like particles are constrained to exist only at the same regular lattice nodes as the fluid particles of the LBM. Generally speaking, the realizing of the cell automation probabilistic approach involves three steps: (1) Calculating displacement vector

X p

of particles during a time-step ( t ~ t   t ).

(2) Determining the value of the discrete motion probabilities

p , a = 1,2,3,4 .

(3) Generating random numbers and obtaining new location of particles at t   t . In Wang's improved CA model (Wang et al., 2012), the displacement vector

X p

is calculated by the

Lagrangian approach. The motion of each particle is directly calculated from the Newtonian equation under consideration of different external forces. In this paper, drag, Saffman lift, gravity and Brownian force are all considered. The Newtonian equation of a particle is given as (Ansari et al., 2015):  pVp dup / dt   f (up  u)  (  p   )Vp g  FB  FL

(9) in which,

u p  dX p / dt

,

f  3 d p / Cc , Cc  1  (1.257  0.4e

vector of the particle respectively.

u

( 1.1d p / 2 )

is velocity vector of the gas.

is the vector of external force acceleration.

FB

Cc

f

,

up

are density, volume and velocity

is Stokes-Cunningham drag coefficient.

is Brownian force from Brownian diffusion.

to the velocity distribution change around particles. diameter and

)× 2 / d p ,  p , V p

Xp

FL

g

is lift force due

is displacement vector of the particle.

dp

is particle

is an empirical correction factor called Cunningham slip correction factor.  is mean free path

of gas, which is expressed as gas temperature respectively.

    M / 2RT R

,,

M

and

T

are gas kinematic viscosity, gas molecular weight and

is the universal gas constant. Brownian force

calculated using the Eq. (10,11,12,13) for two-dimensional problems. and Saffman lift force in the x and y directions respectively. 10

FB , x , FB , y

FB

,

and Saffman lift force FL , x

and

FL , y

FL

are

are Brownian force

6kB d pT

FB , x   x

Cc

(10) 6kB d pT

FB , y   y

Cc

(11) FL , x  1.615 0.5d p2 (u x  u p , x )

du y dx

sign(

du y dx

)

(12) FL , y  1.615 0.5d p2 (u y  u p , y )

dux du sign( x ) dy dy

(13) x

,

and u p, x

 y are

and

u p, y

zero mean, unit-variance independent Gaussian random numbers.

kB

is Boltzmann constant.

ux , u y

are gas and particle velocities in the x and y directions, respectively. The particle velocity in each

moment can be calculated by Eq. (14), Eq. (15) (Ansari et al., 2015). {t }

 (    )Vp g x  FB , x  (  p   )V p g x  FB , x {t } u{pt,x t }   u x  p )] e  +[u p , x  (u x   f  FL , x f  FL , x  

 ( f  FL , x )

 pV p

t

(14) {t }

 (    )Vp g y  FB , y  (  p   )Vp g y  FB , y {t } u{pt, y t }   u y  p )] e   [u p , y  (u y   f  FL , y f  FL , y  

 ( f  FL , y )

 pV p

t

(15) gx

and

gy

are external force accelerations in the

x

and

y

directions. The particle location in each moment can

be calculated by Eq. (16), Eq. (17) (Ansari et al., 2015).  V  (    )Vp g x  FB , x   x{t  t }  x{t}   u x  p  t p p    f  F f  FL , x L, x    

{t }

  (  p   )V p g x  FB , x   u p , x   u x    f  FL , x    

11

 ( f  FL , x )   t   1  e  pV p      

(16)

y

{t   t }

y

{t }

{t } {t }  ( f  FL , y )  V   t    (  p   )Vp g y  FB , y  (  p   )Vp g y  FB , y     p p  V    uy  t   u p , y   u y   1  e p p          f  FL, y f  FL, y  f  FL, y         

p2 A XA

p3 p4

p1

△Xp

(17)

XA+e1δt B

XA+e8δt

XA+e4δt

Figure 4. Rules of particle motion in cell automation probabilistic model. Fig. 4 illustrates the rules of particle motion in cell automation probabilistic model. During a time-step, a solid particle at a node may still stay at the original node or migrate to its nearest-neighbor node, depending on the motion probabilities of the particle ( p ,  1,2,3,4 , shown in Fig. 4): X A

t + δt

 X A   1  e1   t  2  e2   t  3  e3   t  4  e4   t t

(18) where



e1  e3

, it means that 1 and

p

is a Boolean variable which is equal to 1 with probability 3

p

or equal to 0 with probability

1  p .

Since

are not equal to 1 at the same time. The magnitude of the motion probability

is proportional to the projection of the actual displacement of the particle in the direction

eα

.

 ΔX p  t   p  max  0,  u p  eα    max  0,  eα  ,  1,2,3,4 x   x 

(19) Giving the example as shown in Fig. 4, e4  e2 ,

it is clear that

with probability

p2  0 , p3  0 .

p1 1  p4  ,

p1  e1  X p /  x , p4  e4  X p /  x

. Since

p1  0

,

p4  0

and e1  e3 ,

The particle may stay at rest with probability 1  p1 1  p4  , or jumps to east

to south with probability 1  p1  p4 , to south-east with probability

p1 p4 .

The Monte

Carlo method is used to determine the new position of the particle after a time-step. First, two independent random numbers

r1

, r2 are generated through random number generator, which are obedient to a uniform

distribution in the interval [0,1]. It is proven that if the number of simulation particles is large enough binomial 12

scattering can be approximated by a Gaussian distribution (Chopard et al., 1994). Then, through comparing the values of two groups  r1, p1  ,  r2 , p4  respectively, the new site of the particle is determined, which follows Eq. (20). if   if  if   if

r1  p1 and r2  p4 , X p

t + δt 

=X A 

r1  p1 and r2  p4 , X p

 X A   e1 t

r1  p1 and r2  p4 , X p

 X A   e4 t

r1  p1 and r2  p4 , X p

 X A   e1 t  e4 t  X p   e8 t

t + δt  t + δt 

t + δt 

t

t t

t

t

(20) 2.3 Calculation domain and boundary conditions

Figure 5. A schematic showing computational domain for two-dimensional calculations. Given the symmetric structure of the wall-flow DPF, an inlet channel and two halves of adjacent outlet channels are chosen as the computational domain, which are shown in Fig. 5.

Lc , H c , Ws

and W p represent the

channel length, the channel width, the channel wall thickness and the plug thickness, respectively. The porous structure units are generated by QSGS. The main geometrical features of the filters used in the simulation are shown in Table 1. In order to ensure the purpose of fully developed flow for simulation stability, two domains of 0.2 times channel length are added to the upstream of inlet channel and the downstream of outlet channel respectively. Table 1: Main geometrical features of the filters and operative conditions used in the simulation. Filter type Channel length Lc Channel width H c Channel wall Plug thickness thickness Ws WW s p

Case 1

13

Case 2 10.8mm 1.8mm 0.2mm 0.2mm

Case 3

38.16% Porosity  0.54kg/m3 Gas density  Gas temperature T 600K Gas viscosity  3.2×10-5Pa·s 2000kg/m3 Particle density  p Upstream velocity (1m/s, (3m/s, (6m/s, 5 u Downstream 100)Pa 0) 0) 0 pressure p0 The boundary conditions of the calculation domain are specified as follows. The inlet velocity of fluid is constant:

u  u0 ,

and it is considered that the initial velocity of solid particles is the same with the fluid velocity.

At the outflow boundary, a fixed pressure p0 is imposed. The distribution functions

f

at both inlet boundary

and outlet boundary are computed by the method proposed by Guo et al. (2002). The distribution function

f

is

decomposed into equilibrium and nonequilibrium parts approximated with a first-order extrapolation of the nonequilibrium part of the distribution function at the neighboring fluid nodes. Periodic boundary condition is used for the upper and lower boundaries, that is, a fluid particle or a solid particle moving out of the domain from the upper or lower boundary will re-enter into the domain through the opposite boundary immediately. As for the surface of solid boundaries, such as the wall of porous media or the plug of the channels, no-slip boundary condition is applied and the distribution function is calculated by bounce back boundary scheme. Soot particles will stick to the surface of the porous wall once they reach the surface. 1201×201 nodes are chosen to represent the computational domain. The grid size is 20

in the simulation and the grid independence study is

provided in section 3. A program was developed in the C++ language to simulate fluid and particle motion. The flowchart of LB-CA model algorithm is shown in Fig. 6.

14

Start Initialization

Fluid Evolution

f  x  e t , t   t   f  x , t   

1  f  x, t   f eq  x, t  

Boundary Conditions Calculating Fluid Variables

Stop Criterion Reached ?

No

Yes

Releasing Particles Particles Evolution (Cell automation probabilistic approach)

Calculating Particles Variables

No

Stop Criterion Reached ? Yes End

Figure 6. Flowchart of LB-CA model algorithm. 3. Model validation In order to validate the accuracy of the gas phase model, a channel, where a laminar flow past a square block, is used. The velocity values at different positions calculated by present model are compared with the literature values. The computational domain schematic is depicted in Fig. 7, where the dimension of the block is D*D (D=1.8mm) and the flow configuration is defined by L=50D, H=8D, and L1=12.5D. Fig. 8 presents various velocity profiles at Re=100, Re=50 and Re=30 (Reynolds number of the flow is based on maximum inlet velocity

umax

and cylinder width D, Re  Dumax /  ), respectively, with the same coarse mesh of 501×81. R-

squared is a widely used goodness-of-fit measure. The R-squared values of Re=50 and Re=30 in Fig. 8(a) are 0.976 and 0.994 respectively. The corresponding R-squared values are 1.0 and 1.0 in Fig. 8(b), 0.943 and 0.932 15

in Fig. 8(c), 0.981 and 0.985 in Fig. 8(d). All the R-squared values exceed 0.93 (R-squared closer to 1) indicating that the velocity profiles at Re=50 and Re=30 are in excellent agreement with the computational data of previous studies (Guo et al., 2008). Overall the velocity profiles at Re=100 and the literature values (Breuer et al., 2000) are in agreement with each other. The R-squared values are 0.996, 0.911, 0.935 and 0.980 in Fig. 8(a)-(d) respectively. However, small discrepancies are observed in the velocity profiles. That is because flow around a square cylinder in a channel at D/H=0.125 shows a steady behavior for Re<60 and becomes unsteady when 60
y

H x D

H/2

L1

L

Figure 7. Schematic diagram of the flow past a square block in a 2D channel. (a)

0.5

(b)

Present work Re=100 Present work Re=50

1.0

Present work Re=100 Present work Re=50 Present work Re=30 M. Breuer et al. Re=100 Guo et al. Re=50 Guo et al. Re=30

0.5

0.0

-10

(c)

-5

0

5

10 15 x/D [-]

20

25

30

0.0

M. Breuer et al. Re=100 Guo et al. Re=30

-0.5

Guo et al. Re=30

35

-10

M. Breuer et al. Re=100 Guo et al. Re=50 Guo et al. Re=30

-5

0

5

10 15 x/D [-]

0.5

-2

0 y/D [-]

30

35

0.0

Present work Re=100 Present work Re=50 Present work Re=30 M. Breuer et al. Re=100 Guo et al. Re=50 Guo et al. Re=30

-0.2

Present work Re=100 Present work Re=50 Present work Re=30 -4

25

(d)

1.0

0.0

20

0.2

v (0,y)/Umax [-]

u (0,y)/Umax [-]

1.5

v (x,0)/Umax [-]

u (x,0)/Umax [-]

Present work Re=30

2

-0.4

4

-4

16

-2

0 y/D [-]

2

4

Deposition efficiency (front side) (%)

Figure 8. Velocity profiles at Re=100, Re=50 and Re=30, (a) u (x,0), (b) v (x,0), (c) u (0,y), (d) v (0,y). 80 Total number of grids 900×120 Total number of grids 1800×240 Total number of grids 2700×360 Total number of grids 3600×480

60

40

20

0 10-4

10-3

10-2 10-1 Stokes number

100

Figure 9. Mesh independent study for particle deposition. To examine the validity of the CA model, the motion of solid particles is simulated in a channel with a square cylinder obstruction and the particle deposition efficiencies on the front side of the block are calculated. The schematic diagram is the same as described in Fig. 7 whereas the flow configuration is different, L=30D (D=1.8mm), H=4D, and L1=9.5D. 10 groups of tests are conducted at Re=200, and in every group 620 particles for a given Stokes number ( Stk  Cc  p d p2um / 18 D ) are injected uniformly in the projected area of the block at a distance of 9D upstream of the block. In the CA model solid particles are constrained to exist only at the same grid nodes as the fluid particles of LBM. Thus, the motion of solid particles is affected by grid density. Before performing computations on particle transport and deposition, a mesh independence study is performed. The grid resolution of 900×120 has been proved to be good enough to obtain reliable flow fields (Jafari et al., 2010). Based on that, four different uniform mesh grids, 900×120, 1800×240, 2700×360, 3600×480, are examined. The block diameter in the lattice units is corresponding to 30, 60, 90, and 120, respectively. The impact of the mesh on the particle deposition efficiency is shown in Fig. 9. The comparison shows that the particle deposition is fairly sensitive to the grid density. The differences of particle deposition efficiency at different Stokes numbers are calculated when the number of grids increases from 900×120 to 1800×240 and from 1800×240 to 2700×360. The average percentages of those differences to the corresponding deposition efficiency of the grids 900×120 are 39.334% and 15.382% respectively. Although the accuracy of particle deposition efficiency increases with increasing number of grids, the calculation time becomes much longer. It is also found that difference between 2700×360 and 3600×480 is not significant and the average percentage of this difference has 17

dropped to 5.430%. Therefore, the grid of 2700 × 360 can be used with short calculation time and acceptable results. The grid size is 20

. The computed particle deposition efficiencies for different Stokes numbers are

compared with recent numerical studies (Jafari et al., 2010; Brandon and Aggarwal, 2001; Afrouzi et al., 2015; Tong et al., 2017) in Fig. 10. The simulation result of present work agrees with the existing results. The deposition efficiency remains roughly constant for particles with Stokes number less than 0.1 and increases

Deposition efficiency (front side) (%)

rapidly when Stokes number is greater than 0.1. 100 90 80 70 60

Brandon & Aggrawal (2001) S.Jafari et al. (2010) H. Hassanzadeh Afrouzi et al. (2015) Tong et al. (2017) Present work based on LB-CA model

50 40 30 20 10 0 10-4

10-3

10-2

10-1

100

101

Stokes number

Figure 10. Validation of particle deposition efficiency with Stokes number. 4. Results and discussion 4.1 Flow field In conventional numerical approaches (Bensaid et al., 2010; Sbrizzai et al., 2005), the flow field through porous walls is modelled by the Darcy equation without consideration of the microstructure of porous walls in DPF. In this simulation, the flow field through the porous walls is modelled directly by LBM considering the microstructure of porous walls. The porous structure is generated by QSGS. The simulated parameters are shown in Table 1. Three different upstream velocities are used in this study. In this section, the flow is simulated without soot accumulation.

18

Figure 11. Dimensionless velocity field (u/u0) inside the inlet and outlet channels at three different upstream velocities: (a) u0=1m/s, (b) u0=3m/s, (c) u0=6m/s. Fig. 11 shows velocity field inside the inlet and outlet channels at different upstream velocities. For each case, since the exit of inlet channel is closed, the velocity decreases along the inlet channel gradually. The flow is forced to pass through the porous filter wall and the velocity profile passing through the wall will be discussed later. In the outlet channel, the flow rate is increased toward the end of outlet channel. When the upstream velocity increases, the inertia of the flow becomes more important. In the inlet channel, the location where velocity decreases to the same dimensionless speed u/u0 are gradually moving farther away from the entrance with the increase of upstream velocity, which means that a larger portion of the flow chooses to cross the porous wall at the rear part. The similar phenomenon can be found in outlet channel: the location of velocity increasing to the same dimensionless speed u/u0 are gradually closer to the exit of outlet channel with the increase of upstream velocity. Fig. 12 shows pressure field inside inlet and outlet channels at different upstream velocities. For each case, the pressure inside the inlet channel is higher than that inside the outlet channel. The pressure difference is the driving force for exhaust gas flow movement across the porous wall from an inlet channel to an outlet channel. In the outlet channel, the pressure gradually decreases along the flow path. And the pressure at the same position gradually increases with the increase of upstream velocity. In the inlet channel, different distribution characteristics of the pressure are presented with the increase of upstream velocity. At low upstream velocity 19

values (u0=1m/s; 3m/s), the pressure field has a uniform distribution; At high upstream velocity values (u0=6m/s), the pressure field distribution turns to be non-uniform: the pressure values at the rear part of the channel are obviously higher than that of the rest part of the inlet channel. In order to analyze the variation of pressure field and velocity field along the channel quantitatively, the profiles of axial pressure, axial velocity and wall velocity along the channel length at different upstream velocities are shown in Fig. 13. The dimensionless channel length 0
Figure 12. Dimensionless pressure field (pre/p0) inside the inlet and outlet channels at three different upstream velocities: (a) u0=1m/s, (b) u0=3m/s, (c) u0=6m/s.

20

3.5

1 m/s 3 m/s 6 m/s

Dimensionless axial pressure p/p0 [-]

(a)

1.0020

3.0

1 m/s 3 m/s 6 m/s

1.0015 1.0010

2.5 2.0

1 m/s 3 m/s 6 m/s

1.5

1.0005

1 m/s 3 m/s 6 m/s

(b)

Dimensionless axial velocity u/u0 [-]

1.0025

1.0

1.0000

0.5

0.9995

0.0 0.0

0.2

0.4

0.6

0.8

1.0

0.0

Dimensionless channel length x/Lc [-]

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

1.0

(c)

Dimensionless wall velocity uw/u0 [-]

0.9

1 m/s 3 m/s 6 m/s

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0.0

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

Figure13. Variation of axial pressure and velocity along the channel length at different upstream velocities: (a) axial pressure profiles along the inlet channel (solid line) and outlet channel (short dot line), (b) axial velocity variation along the inlet channel (solid line) and outlet channel (short dot line), (c) wall velocity along the channel. The general trends of variation of axial velocity profiles along the channel length under different upstream velocities are basically the same in Fig. 13(b). The value of inlet channel velocity increases first and then decreases. In the outlet channel, the velocity shows an increasing profile. In Fig. 13(c) the wall flow velocity profiles on the porous surface are shown. In some areas, the wall flow velocity is zero. This is because the flow only passes through the pore of the porous surface. Although some variations are observed due to the nonuniformity of the porous structure, two peaks (one big and one small) of local wall flow velocities emerge at the fore part and rear part of the inlet channel at u0=6m/s, 3m/s. And the wall flow velocity profile is almost uniform along the channel at u0=1m/s. With the increase of upstream velocity, the linearity of the axial velocity curves in Fig. 13(b) and the uniformity of the wall velocity in Fig. 13(c) along the channel both become worse gradually. Those behaviors can be explained by the fact that at larger upstream velocity, a larger portion of the flow chooses to cross the porous wall at the rear part of the inlet channel. 21

4.2 Particle deposition and distribution of clean porous wall Solid particle capture process of clear porous wall is simulated in this section. In the following cases, solid particles will “disappear” if they deposit on the porous wall, so the shape of porous wall does not change and the porous wall keeps “clean”. This assumption ensures that deposition of every particle is an independent event free of the effect of deposited particles. With this method, the statistic information (such as probability density function) about particle deposition and distribution can be obtained. In the simulation, 10030 particles with a given particle diameter are injected uniformly for 118 times at the entrance of the inlet channel. Since the time interval between two releases is long enough (2000 time-steps in lattice units), the interaction between the gas and the soot particle is treated as one-way coupling. Table 2: Number of the deposited particles at different upstream velocities based on the LB-CA method and the LB-Lagrangian method respectively. Velocity

Position TOP

1m/s BOTTOM TOP 3m/s BOTTOM TOP 6m/s BOTTOM

Method LB-CA LB-LA LB-CA LB-LA LB-CA LB-LA LB-CA LB-LA LB-CA LB-LA LB-CA LB-LA

Diameter 1μm 100nm 4978 4981 5006 5010 4992 5001 5016 5015 4928 4975 5004 5002 4903 5004 5013 5011 4892 4950 4999 4995 4887 4949 5005 5006

10nm 4978 5002 4992 5014 5006 5015 4991 5003 4939 5004 4954 4995

Since the Lagrangian method has a higher precision in simulating the particle trajectories than the CA method (Wang et al., 2013), particles motion in the DPF channels (not including area inside the porous walls) and particles deposition on the surfaces of porous wall have been simulated (Sbrizzai et al., 2005; Liu et al., 2009). To better verify the accuracy of LB-CA model, particle depositions and distributions on the surfaces of porous wall (top surface and bottom surface) are also simulated by LB-Lagrangian model. Due to the inconvenience of the Lagrangian method in dealing with the complex boundary of the non-virtual porous walls and dynamic boundary conditions of soot deposits, the porous wall is viewed as a rectangle area in particles 22

motion simulation, though the non-virtual porous structure is considered in flow simulation. Soot particles will stick to the surface of the rectangle area once they reach the surface. Therefore, both the Lagrangian method and the CA method can be used for particles motion simulation. The number of the deposited particles at different upstream velocities based on both methods are listed in Table 2. For each method at the same upstream velocity, the number of the deposited particles on the top surface is almost the same as the number of the deposited particles on the bottom surface. For each method at the same surface of porous wall, there is little change in the number of the deposited particles at different upstream velocities. It is found that the difference of the numbers of the deposited particles between the LB-CA method and the LB-Lagrangian method respectively is not significant, which means that the LB-CA method has good ability to predict the number of the deposited particles. In order to describe particle depositions along the channel quantitatively, particles deposited on the bottom porous wall surface are counted. Fig. 14-Fig. 16 show the particle deposition efficiencies along the inlet channel at different upstream velocities respectively. Since particle depositions are independent events, the particle deposition efficiency distribution also represents the probability density function of the deposited particles. Fig. 14 illustrates the profiles of deposition efficiency of particles with different particle sizes along the channel under u0=6m/s. The results from the LB-CA method show that the deposition distributions of particles with different diameters are similar to the distributions of wall velocity at the middle part and the rear part of the inlet channel. With the increase of particle diameter, the profiles of deposition distribution gradually move toward the downstream of the inlet channel, which demonstrates that particle inertia plays a key role in determining the particle deposition distribution. It is obvious that the profiles of deposition distribution from the LB-CA method are in close agreement with those from the LB-Lagrangian method. The results from the LB-Lagrangian method show that particles do not deposit at the fore part of the inlet channel most likely due to the vena contracta condition. This phenomenon was also recorded by Sbrizzai et al. (2005) and Liu et al. (2009) based on the conventional numerical Eulerian-Lagrangian approach. With the increase of particle diameter, the initial deposition locations are gradually moving farther away from the entrance of inlet channel. In the conventional 23

numerical Eulerian-Eulerian approach (Bensaid et al., 2010), soot particles are found to deposit at the fore part of the inlet channel and the profiles of the soot layer thickness along the filter agree well with the experimental data. From this point on consideration, the LB-CA method is more consistent with the reality than the LBLagrangian method. 12

12

(b) dp=100 nm

10

Deposition efficiency (%)

Deposition efficiency (%)

(a) dp=10 nm LB-CA LB-LA

8 6 4 2 0

10

LB-CA LB-LA

8 6 4 2 0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

Dimensionless channel length x/Lc [-]

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

12

Deposition efficiency (%)

(c) dp=1 μm 10

LB-CA LB-LA

8 6 4 2 0 0.0

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

Figure 14. Deposition efficiency of particles deposited on the bottom porous wall surface along the inlet channel for different soot particle sizes at u0=6m/s, (a) =10 nm, (b) =100 nm, (c) =1µm. 12

12

10

(b) dp=100 nm

Deposition efficiency (%)

Deposition efficiency (%)

(a) dp=10 nm LB-CA LB-LA

8 6 4 2

10

0

LB-CA LB-LA

8 6 4 2 0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

Dimensionless channel length x/Lc [-]

0.2

0.4

0.6

0.8

Dimensionless channel length x/Lc [-]

24

1.0

12

Deposition efficiency (%)

(c) dp=1 μm 10

LB-CA LB-LA

8 6 4 2 0 0.0

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

Figure 15. Deposition efficiency of particles deposited on the bottom porous wall surface along the inlet channel for different soot particle sizes at u0=3m/s, (a) =10 nm, (b) =100 nm, (c) =1µm. 12

12

10

(b) dp=100 nm

Deposition efficiency (%)

Deposition efficiency (%)

(a) dp=10 nm LB-CA LB-LA

8 6 4 2

10

0

LB-CA LB-LA

8 6 4 2 0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

Dimensionless channel length x/Lc [-]

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

12

Deposition efficiency (%)

(c) dp=1 μm 10

LB-CA LB-LA

8 6 4 2 0 0.0

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

Figure 16. Deposition efficiency of particles deposited on the bottom porous wall surface along the inlet channel for different soot particle sizes at u0=1m/s, (a) =10 nm, (b) =100 nm, (c) =1µm.

25

60 50 40 30

(a) 6 m/s CA dp=10 nm CA dp=100 nm CA dp=1μm LA dp=10 nm LA dp=100 nm LA dp=1μm

70

Deposition efficiency (%)

Deposition efficiency (%)

70

20 10 0

Deposition efficiency (%)

60 50 40 30

50 40 30 20 10 0

The fore part 0
60

(b) 3 m/s CA dp=10 nm CA dp=100 nm CA dp=1μm LA dp=10 nm LA dp=100 nm LA dp=1μm

The middle part 0.15
The rear part 0.75
The middle part 0.15
The rear part 0.75
The fore part 0
The middle part 0.15
The rear part 0.75
(c) 1 m/s CA dp=10 nm CA dp=100 nm CA dp=1μm LA dp=10 nm LA dp=100 nm LA dp=1μm

20 10 0

The fore part 0
The relative standard deviation (%)

Figure 17. Deposition efficiencies of particles deposited on the different parts of bottom porous wall surface along the inlet channel, (a) u0=6m/s, (b) u0=3m/s, (c) u0=1m/s. CA dp=10 nm CA dp=100 nm CA dp=1μm LA dp=10 nm LA dp=100 nm LA dp=1μm

120

100

80

60

40

20

u0= 6 m/s

u0= 3 m/s

u0= 1 m/s

Figure 18. Variation of the relative standard deviations at different upstream velocities. The deposition distributions of particles with different diameters of both methods are gradually similar to the distributions of wall velocity with a decrease in upstream velocity in Fig. 15 and Fig. 16. The deposition efficiencies of particles deposited on the different parts of bottom surface of porous wall along the inlet channel are calculated, which is shown in Fig. 17. It is found that for particles with a given diameter, the deposition number of particles significantly decreases at the rear part of the inlet channel, and significantly increases at the fore part of the inlet channel as the upstream velocity decreases. At the middle part of the inlet channel, the 26

deposition amount increases slightly with a decrease in upstream velocity. The relative standard deviation (RSD) is calculated to assess the uniformity of particle deposition distribution. The smaller the RSD value is, the better the uniformity is. In Fig. 18 the value of RSD for particles with a given diameter gradually decreases with the decrease of the upstream velocity, which means that the uniformity of particle deposition distribution becomes better gradually as the upstream velocity decreases. From Fig. 14-Fig. 16, it is also seen that compared with the effect of the particle size, the upstream velocity influences more on the particle deposition distributions, which demonstrates that upstream velocity plays a key role in determining the particle deposition distribution.

Figure 19. Local particulate number distribution in the bottom porous wall for different soot particle sizes at u0=3m/s, (a) =10 nm, (b) =100 nm, (c) =1µm. In conventional numerical approaches, the “unit-cell” filtration model has been widely used in modeling particle distribution inside the porous wall (Konstandopoulos and Johnson, 1989; Konstandopoulos et al., 2000; Gong and Rutland, 2015). In the “unit-cell” model, the porous wall is discretized into a number of slabs along the wall direction. The porous structure of each slab is virtual, which is characterized by a mean pore size and a mean porosity (Konstandopoulos and Johnson, 1989; Konstandopoulos et al., 2000) or by a probability density 27

function (Gong and Rutland, 2015). It was found that particulate concentration gradually decreases in the wall thickness direction (Konstandopoulos and Johnson, 1989; Gong and Rutland, 2015). However, little effort has been carried out to study particle distribution in the porous wall in the axial direction. In this part only the CA method is used to simulate the particle motion in the DPF channels and inside the porous wall without assumption that the porous wall is viewed as a rectangle area. Fig. 19 shows the local particulate number distribution in the bottom porous wall for different soot particle sizes. For the sake of brevity, among the investigated velocities, only the case with upstream velocity of 3m/s is presented and discussed here. The results show there exists a gradient of particulate number across the porous wall and most of particles are deposited at the area nearest the inlet channel (the dimensionless porous wall thickness 0.7≤ y/Ws ≤1.0). They are consistent with the results in the literatures (Konstandopoulos and Johnson, 1989; Gong and Rutland, 2015). The results also show that the deposition distributions of particles with different diameters at the top of the porous wall are similar to the distributions of wall velocity along the channel length. With the increase of porous wall depth, this phenomenon is no longer obvious. By comparing the particle number distribution of particles with different sizes, it is found that the amount of particle with 100nm inside porous wall (0.7≤ y/Ws ≤1.0) is higher than that of particle with 10nm or 1µm, which means the particles of 100nm are easier to penetrate through the porous media. The same result was found by Gong and Rutland (2015). 4.3 Dynamic evolution in soot particle capture process In real situation, soot particles will deposit on the surface of the porous wall without occurrence of “clean” process. The deposited soot particles will change the porosity and the pore size of the porous wall and thus affect the flow field. It is necessary to investigate the effect of the deposited soot particles on the pressure field and velocity field in particle capture process. In soot particle capture process, soot particles will deposit on the boundary nodes of porous wall. As time goes by, more and more soot particles deposit on the boundary nodes. When the number of the deposited particles in one boundary node reaches the threshold value

, the

node is altered from boundary node to a “solid” node (Wang et al., 2012; Tong et al., 2017), just like a solid 28

point of porous wall. Particles of 100nm are used in this simulation. The upstream velocity u0 is 3 m/s. The number of particles released in this simulation is 100300 and the mode of particles released in this simulation is the same as that mentioned in the previous section 4.2.

Figure 20. Dynamic evolution of solid boundary node in soot particle capture process for 100 nm soot particles at u0=3m/s. Fig. 20 shows the dynamic evolution of solid boundary in soot particle capture process. The solid boundary with different colors represents the newly generated solid nodes within different time periods, the red for period from 6.7 ms to 66.7 ms, the green for period from 66.7 ms to 126.7 ms, the blue and the pink for period from 126.7 ms to 186.7 ms and 186.7 ms to 246.7 ms. It is found that the solid boundary nodes made of the captured soot particles appear first at the rear part of the inlet channel in Fig. 20(a) and then gradually move to the middle part of the inlet channel in Fig. 20(b). As time goes by, the solid boundary nodes are distributed along the whole channel length in Fig. 20(c) and Fig. 20(d). Through statistics and analysis, more than 50 percent of solid boundary nodes are formed in the rear part of the inlet channel. This phenomenon coincides with the obvious tendency of particles with

=100nm to deposit at the end of the inlet channel in the previous 29

section 4.2. Fig. 21 shows distributions of the solid nodes quantitatively. The solid nodes with different colors represent the newly generated solid nodes within different time periods, which are the same as those in Fig. 20. It can be seen that a thick and dense deposit layer appears from 0.75 to 1.0 dimensionless channel length. For this reason, the actual channel width is reduced at this part of channel. Thick deposit layers can also be found at some positions in the channel, like around x/Lc=0.15, 0.6, which also makes the channel width narrower.

Figure 21. Distributions of the deposited soot particles of 100 nm at 246.7ms. Fig. 22 depicts the effects of the deposited soot particles on the flow field along the channel length. In order to better reflect the effects of the deposited soot particles on the flow field, relative differences between the simulation results of “loaded” porous wall and that of “clean” porous wall in the previous section 4.2 are calculated. In Fig. 22(a), the axial relative pressure difference gradually increases as time goes by, which means that the resistance to the flow gradually increases with the deposition of soot particles. It is also found that the maximum axial relative pressure difference appears around x/Lc=0.75. This results from the vena contracta condition from 0.75 to 1.0 dimensionless channel length, which makes the axial pressure decrease. It is obvious that the trend of relative pressure difference curves varies at the fore part of the channel. The relative pressure difference shows a decreasing profile at 246.7ms, because with the deposition of soot particles a tapered channel is formed at the fore part of the channel. The tapered channel widens the scope of vena contracta condition. Fig. 22(b) shows the variation of relative axial velocity difference along the channel length at different time. On the whole, the relative axial velocity difference is less than 0, which means that the axial velocity 30

decreases in the depositing process of soot particles. However, at 246.7ms the opposite results are found at the front of the channel, around x/Lc=0.15. The tapered channel at this part makes the axial velocity increase. At x/Lc=0.15, the width of the channel reaches its minimum shown in Fig. 21 and the relative axial velocity difference reaches a maximum. A similar situation occurs at the end of the channel. A tapered channel has been formed with the deposition of soot particles and the relative axial velocity difference shows an increasing profile from 0.75 to 1.0 dimensionless channel length. At the middle part of channel, the relative axial velocity difference gradually decreases along the channel length except the area around x/Lc=0.6. The general trend of relative axial velocity difference curve at 186.7ms is basically the same with that at 246.7ms. The general trend of relative axial velocity difference curve at 126.7ms is similar with that at 66.7ms. The difference between them is mainly reflected in the variation trend at the front of the channel. In the initial stage of particle deposition (66.7ms, 126.7ms), fewer particles choose to deposit and the tapered channel at the fore part of

7x10-5 6x10

66.7 ms 186.7 ms

-5

126.7ms 246.7 ms

Dimensionless axial relative velocity difference Δu/u0 [-]

Dimensionless axial relative pressure difference Δp/p0 [-]

channel has not formed yet. Thus, relative axial velocity difference shows a decreasing profile. (a)

5x10-5 4x10-5

0.02

126.7ms 246.7 ms

(b)

0.00

-0.04 -0.06

-5

-0.08

1x10-5

-0.10

0

-0.12

-1x10-5

-0.14 0.0

0.15 0.2

0.4

0.6

0.75 0.8

1.0

0.0

Dimensionless channel length x/Lc [-]

Dimensionless relative wall velocity difference Δuw/u0 [-]

66.7 ms 186.7 ms

0.04

-0.02

3x10-5 2x10

0.06

0.15 0.2

0.4

0.6

0.75 0.8

1.0

Dimensionless channel length x/Lc [-]

0.2

(c) 0.1 0.0 -0.1 -0.2

66.7 ms 126.7ms 186.7 ms 246.7 ms

-0.3 -0.4 -0.5 0.0

0.2

0.4

0.6

0.8

1.0

Dimensionless channel length x/Lc [-]

Figure 22. Effects of the deposited soot particles of 100 nm on flow field along the inlet channel length at u0=3m/s, (a) axial relative pressure difference, (b) axial relative velocity difference, (c) relative wall velocity difference. 31

Fig. 22(c) shows the variation of relative wall velocity difference along the channel length at different time. It can be seen that the areas where wall velocity undergoes obvious changes appear first at the rear part of the inlet channel (66.7ms) and then gradually move to the middle part of the inlet channel (126.7ms). As time goes by, the areas where wall velocity undergoes changes gradually expand to the whole inlet channel (186.7ms, 246.7ms). This phenomenon is consistent with the depositing process of soot particles in Fig. 20. Generally, the absolute value of local wall velocity difference is increasing over time. 5. Conclusions A two-dimensional mesoscopic model has been developed in the present paper using the incompressible lattice Boltzmann method for fluid flow and the cell automation probabilistic approach for particle motion. The LB-CA model is validated with the results of previous studies. The porous wall of DPF is generated through QSGS. The effects of different upstream velocities on the pressure field and the velocity field are investigated. The distribution and deposition of particles with different sizes in clean channels are simulated based on the LB-CA method and the LB-Lagrangian method respectively. The effects of deposited soot particles on the flow field are evaluated in real soot particle capture process. The following conclusions can be drawn from present results: The inlet channel’s axial pressure has a decrease first at u0=6m/s, 3m/s and then the pressure increases due to the impact of inertia. At u0=1m/s the axial pressure value is quite constant along the inlet channel. The outlet pressure at different upstream velocities shows a decreasing profile. The general trends of variation of axial velocity profiles along the channel length are basically the same under different upstream velocities. The value of inlet channel velocity increases first and then decreases. In the outlet channel, the velocity shows an increasing profile. Two peaks of local wall flow velocities emerge at the fore part and rear part of the inlet channel at u0=6m/s, 3m/s. And the wall flow velocity profile is almost uniform along the channel when u0=1m/s. The linearity of the axial velocity curves and the uniformity of the wall velocity along the channel

32

both become worse gradually with the increase of upstream velocity. The distribution of velocity field and pressure field in the channel are significantly affected by the upstream velocity. The LB-CA method is capable to predict the number of the deposited particles on the surfaces of porous wall. In the LB-CA method, the deposition distributions of particles with different diameters at u0=6m/s are similar to the distributions of wall velocity at the middle part and the rear part of the inlet channel. The deposition distributions of particles with different diameters are similar to the distributions of wall velocity at u0=3 m/s, 1 m/s. With the increase of particle diameter at the same upstream velocity, the profiles of deposition distribution move toward the downstream of the inlet channel gradually. For particles with a given diameter, the deposition number of particles significantly decreases at the rear part of the inlet channel, and remarkably increases at the fore part of the inlet channel as the upstream velocity decreases. At the middle part of the inlet channel, the deposition amount increases slightly with a decrease in upstream velocity. Compared with the effect of the particle size, the upstream velocity influences more on the particle deposition distributions. The uniformity of particle deposition distribution becomes better gradually as the upstream velocity decreases. The profiles of deposition distribution from the LB-CA method are in close agreement with those from the LBLagrangian method. The LB-CA method is more consistent with the reality than the LB-Lagrangian method. A gradient of particulate number across the porous wall is found and most of particles are deposited at the area nearest the inlet channel (0.7≤ y/Ws ≤1.0). The deposition distributions of particles with different diameters at the top of the porous wall are similar to the distributions of wall velocity along the channel length. The particles of 100nm are easier to penetrate through the porous media. For 100n

particles at u0=3 m/s, the solid boundary made of the captured soot particles appears first in the

rear part of the inlet channel and then gradually moves to the middle part of the inlet channel. Finally, the solid boundary nodes are distributed along the whole channel length. On the whole, the deposited soot particles increase the axial pressure and decrease the axial velocity. The axial pressure gradually increases in the soot particle capture process. The axial velocity gradually decreases in the initial stage of particle deposition and then shows a rising tendency in some local areas of the channel. The relative axial pressure difference shows a 33

decreasing profile and the relative axial velocity difference shows an increasing profile from 0.75 to 1.0 dimensionless channel length. At the middle part of channel, the relative axial pressure difference increases and the relative axial velocity difference decreases along the channel length gradually. The variations of relative axial pressure difference and the relative axial velocity difference at the fore part of the inlet channel are not the same in different stages during the particle deposition. The evolution trend of the areas where wall velocity undergoes changes is consistent with the evolution trend of the solid nodes made of the captured soot particles. Generally, the absolute value of local wall velocity difference is increasing over time. Conflict of interest No conflict of interest. Acknowledgments This work was supported by the National Natural Science Foundation of China [grant numbers 51576140, 51276128]; the Special Fund for Development of Small and Medium Enterprises [grant number SQ2013ZOA100012]. References Ansari, V., Goharrizi, A.S., Jafari, S., Abolpour, B., 2015. Numerical study of solid particles motion and deposition in a filter with regular and irregular arrangement of blocks with using lattice Boltzmann method. Computers & Fluids 108, 170-178. https://doi.org/10.1016/j.compfluid.2014.11.022 Afrouzi, H.H., Sedighi, K., Farhadi, M., Moshfegh, A., 2015. Lattice Boltzmann analysis of micro-particles transport in pulsating obstructed channel flow. Computers & Mathematics with Applications 70(5), 11361151. https://doi.org/10.1016/j.camwa.2015.07.008 Brandon, D.J., Aggarwal, S.K., 2001. A numerical investigation of particle deposition on a square cylinder placed in a channel flow. Aerosol Science and Technology 34(4), 340-352. https://doi.org/10.1080/02786820121279 34

Bensaid, S., Marchisio, D.L., Fino, D., Saracco, G., Specchia, V., 2009. Modelling of diesel particulate filtration in wall-flow traps. Chemical Engineering Journal 154(1), 211-218. https://doi.org/10.1016/j.cej.2009.03.043 Bensaid, S., Marchisio, D.L., Fino, D., 2010. Numerical simulation of soot filtration and combustion within diesel particulate filters. Chemical Engineering Science 65(1), 357-363. https://doi.org/10.1016/j.ces.2009.06.051 Bollerhoff, T., Markomanolakis, I., Koltsakis, G., 2012. Filtration and regeneration modeling for particulate filters with inhomogeneous wall structure. Catalysis today 188(1), 24-31. https://doi.org/10.1016/j.cattod.2011.12.017 Breuer, M., Bernsdorf, J., Zeiser, T., Durst, F., 2000. Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume. International journal of heat and fluid flow 21(2), 186-196. https://doi.org/10.1016/S0142-727X(99)00081-8 Choi, S., Lee, K., 2013. Detailed investigation of soot deposition and oxidation characteristics in a diesel particulate filter using optical visualization. SAE Technical Paper 2013-01-0528. DOI: https://doi.org/10.4271/2013-01-0528 Chopard, B., Frachebourg, L., Droz, M., 1994. Multiparticle lattice gas automata for reaction diffusion systems. International Journal of Modern Physics C 5(1), 47-63. https://doi.org/10.1142/S0129183194000052 Chopard, B., Masselot, A., 1999. Cellular automata and lattice Boltzmann methods: a new approach to computational fluid dynamics and particle transport. Future Generation Computer Systems 16(2), 249-257. https://doi.org/10.1016/S0167-739X(99)00050-3 Chopard, B., Masselot, A., Dupuis, A., 2000. A lattice gas model for erosion and particles transport in a fluid. Computer Physics Communications 129(1-3), 167-176. https://doi.org/10.1016/S0010-4655(00)00104-1 Guan, B., Zhan, R., Lin, H., Huang, Z., 2015. Review of the state-of-the-art of exhaust particulate filter technology in internal combustion engines. Journal of environmental management 154, 225-258. https://doi.org/10.1016/j.jenvman.2015.02.027 35

Gong, J., Rutland, C.J., 2015. PDF-Based heterogeneous multiscale filtration model. Environ. Sci. Technol. 49(8), 4963-4970. DOI:10.1021/acs.est.5b00329 Gong, J., Viswanathan, S., Rothamer, D.A., Foster, D.E., Rutland, C.J., 2017. Dynamic heterogeneous multiscale filtration model: probing micro and macroscopic filtration characteristics of gasoline particulate filters. Environmental science & technology 51(19), 11196-11204. DOI: 10.1021/acs.est.7b02535 Orihuela, M.P., Gómez-Martín, A., Miceli, P., Becerra, J.A., Chacartegui, R., Fino, D. 2018. Experimental measurement of the filtration efficiency and pressure drop of wall-flow diesel particulate filters (DPF) made of biomorphic Silicon Carbide using laboratory generated particles. Applied Thermal Engineering 131, 41-53. https://doi.org/10.1016/j.applthermaleng.2017.11.149 Torregrosa, A.J., Serrano, J.R., Piqueras, P., García-Afonso, Ó., 2017. Experimental and computational approach to the transient behaviour of wall-flow diesel particulate filters. Energy 119, 887-900. https://doi.org/10.1016/j.energy.2016.11.051 Konstandopoulos, A.G., Johnson, J.H., 1989. Wall-Flow diesel particulate filters-their pressure drop and collection efficiency. SAE Technical Paper 890405. https://doi.org/10.4271/890405 Konstandopoulos, A.G., Kostoglou, M., Skaperdas, E., Papaioannou, E., Zarvalis, D., Kladopoulou, E., 2000. Fundamental studies of diesel particulate filters: transient loading, regeneration and aging. SAE Technical Paper 2000-01-1016. https://doi.org/10.4271/2000-01-1016 Konstandopoulos, A.G., 2003. Flow resistance descriptors for diesel particulate filters: definitions, measurements and testing. SAE Technical Paper 2003-01-0846. https://doi.org/10.4271/2003-01-0846 Roy, S., Raju, R., Chuang, H.F., Cruden, B.A., Meyyappan, M., 2003. Modeling gas flow through microchannels and nanopores. Journal of applied physics 93(8), 4870-4879. https://doi.org/10.1063/1.1559936 Serrano, J.R., Climent, H., Piqueras, P., Angiolini, E., 2016. Filtration modelling in wall-flow particulate filters of low soot penetration thickness. Energy 112, 883-898. https://doi.org/10.1016/j.energy.2016.06.121

36

Swanson, J., Watts, W., Kittelson, D., Newman, R., Ziebarth, R., 2013. Filtration efficiency and pressure drop of miniature diesel particulate filters. Aerosol Science and Technology 47(4), 452-461. http://doi.org/10.1080/02786826.2012.763087 Sbrizzai, F., Faraldi, P., Soldati, A., 2005. Appraisal of three-dimensional numerical simulation for sub-micron particle deposition in a micro-porous ceramic filter. Chemical Engineering Science 60(23), 6551-6563. https://doi.org/10.1016/j.ces.2005.05.038 Sanui, R., Hanamura, K., 2016. Scanning electron microscopic visualization of bridge formation inside the porous channels of diesel particulate filters. SAE International Journal of Fuels and Lubricants 9(3), 725733. https://doi.org/10.4271/2016-01-9079 Tandon, P., Heibel, A., Whitmore, J., Kekre, N., Chithapragada, K., 2010. Measurement and prediction of filtration efficiency evolution of soot loaded diesel particulate filters. Chemical Engineering Science, 65(16), 4751-4760. https://doi.org/10.1016/j.ces.2010.05.020 Liu, Y., Gong, J., Fu, J., Cai, H., Long, G., 2009. Nanoparticle motion trajectories and deposition in an inlet channel of wall-flow diesel particulate filter. Journal of Aerosol Science 40(4), 307-323. https://doi.org/10.1016/j.jaerosci.2008.12.001 Mohammed, H., Triana, A.P., Yang, S.L., Johnson, J.H., 2006. An advanced 1D 2-layer catalyzed diesel particulate filter model to simulate: filtration by the wall and particulate cake, oxidation in the wall and particulate cake by NO2 and O2, and regeneration by heat addition. SAE Technical Paper 2006-01-0467. DOI: https://doi.org/10.4271/2006-01-0467 Meng, Z., Fang, J., Pu, Y., Yan, Y., Wu, Y., Wang, Y., Song, Q., 2016. Experimental study on the influence of DPF micropore structure and particle property on its Filtration process. Journal of Combustion 2016. http://dx.doi.org/10.1155/2016/9612856 Mohankumar, S., Senthilkumar, P., 2017. Particulate matter formation and its control methodologies for diesel engine: A comprehensive review. Renewable and Sustainable Energy Reviews 80, 1227-1238. https://doi.org/10.1016/j.rser.2017.05.133 37

Masselot, A., Chopard, B., 1998. A lattice Boltzmann model for particle transport and deposition. EPL (Europhysics Letters) 42(3), 259. https://doi.org/10.1209/epl/i1998-00239-3 Wu, G., Kuznetsov, A.V., Jasper, W.J., 2011. Distribution characteristics of exhaust gases and soot particles in a wall-flow ceramics filter. Journal of Aerosol Science 42(7), 447-461. https://doi.org/10.1016/j.jaerosci.2011.04.003 Wang, H., Zhao, H., Guo, Z., Zheng, C., 2012. Numerical simulation of particle capture process of fibrous filters using Lattice Boltzmann two-phase flow model. Powder technology 227, 111-122. https://doi.org/10.1016/j.powtec.2011.12.057 Guo, Z., Shi, B., Wang, N., 2000. Lattice BGK model for incompressible Navier-Stokes equation. Journal of Computational Physics 165(1), 288-306. https://doi.org/10.1006/jcph.2000.6616 Jafari, S., Salmanzadeh, M., Rahnama, M., Ahmadi, G., 2010. Investigation of particle dispersion and deposition in a channel with a square cylinder obstruction using the lattice Boltzmann method. Journal of Aerosol Science 41(2), 198-206. https://doi.org/10.1016/j.jaerosci.2009.10.005 Tong, Z.X., Li, M.J., He, Y.L., Tan, H.Z., 2017. Simulation of real time particle deposition and removal processes on tubes by coupled numerical method. Applied Energy 185(2), 2181-2193. https://doi.org/10.1016/j.apenergy.2016.01.043 He, X., Luo, L.S., 1997. Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Physical Review E 56(6), 6811. DOI:https://doi.org/10.1103/PhysRevE.56.6811 Qian, Y.H., d'Humières, D., Lallemand, P., 1992. Lattice BGK models for Navier-Stokes equation. EPL ( Europhysics Letters) 17(6), 479. https://doi.org/10.1209/0295-5075/17/6/001 Guo, Z., Zheng, C., Shi, B., 2002. Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chinese Physics 11(4), 366. https://doi.org/10.1088/10091963/11/4/310

38

Guo, Z., Liu, H., Luo, L.S., Xu, K., 2008. A comparative study of the LBE and GKS methods for 2D near incompressible laminar flows. Journal of Computational Physics 227(10), 4955-4976. https://doi.org/10.1016/j.jcp.2008.01.024 Wang, H., Zhao, H., Guo, Z., He, Y., Zheng, C. 2013. Lattice Boltzmann method for simulations of gas-particle flows over a backward-facing step. Journal of Computational Physics 239, 57-71. https://doi.org/10.1016/j.jcp.2012.12.032 Wang, M., Pan, N., Wang, J., Chen, S., 2007. Mesoscopic simulations of phase distribution effects on the effective thermal conductivity of microgranular porous media. Journal of Colloid and Interface Science 311(2), 562-570. https://doi.org/10.1016/j.jcis.2007.03.038 Wurzenberger, J.C., Kutschi, S., 2007. Advanced simulation technologies for diesel particulate filters, a fundamental study on asymmetric channel geometries. SAE Technical Paper 2007-01-1137. https://doi.org/10.4271/2007-01-1137 Yamamoto, K., Satake, S., Yamashita, H., Takada, N., Misawa, M., 2006. Lattice Boltzmann simulation on porous structure and soot accumulation. Mathematics and Computers in Simulation 72(2-6), 257-263. https://doi.org/10.1016/j.matcom.2006.05.021 Yamamoto, K., Satake, S., Yamashita, H., Takada, N., Misawa, M., 2007. Lattice Boltzmann simulation on flow with soot accumulation in diesel particulate filter. International Journal of Modern Physics C 18(04), 528-535. https://doi.org/10.1142/S0129183107010760 Yamamoto, K., Satake, S., Yamashita, H., Takada, N., Misawa, M., 2009. Fluid simulation and X-ray CT images for soot deposition in a diesel filter. The European Physical Journal Special Topics 171(1), 205212. doi:10.1140/epjst/e2009-01030-x Yamamoto, K., Yamauchi, K., 2013. Numerical simulation of continuously regenerating diesel particulate filter. Proceedings of the Combustion Institute 34(2), 3083-3090. https://doi.org/10.1016/j.proci.2012.06.117 Yamamoto, K., Sakai, T., 2015. Simulation of continuously regenerating trap with catalyzed DPF. Catalysis Today 242, 357-362. https://doi.org/10.1016/j.cattod.2014.07.022 39

Yamamoto, K., Sakai T., 2016. Effect of pore structure on soot deposition in diesel particulate filter. Computation 4(4), 46. https://doi.org/10.3390/computation4040046 Yuan, J., Sundén, B., 2014. On mechanisms and models of multi-component gas diffusion in porous structures of fuel cell electrodes, International Journal of Heat and Mass Transfer 69, 358-374. https://doi.org/10.1016/j.ijheatmasstransfer.2013.10.032 Definitions/Abbreviations LB-CA

Lattice Boltzmann-cell automation dp

Particle diameter

model

probabilistic model

QSGS

Quartet structure generation set

Cc

Cunningham slip correction factor

DPF

Diesel particulate filter



Mean free path of gas

PM

Particulate matter



Gas kinematic viscosity

Kn

Knudsen number

M

Gas molecular weight

LBM

Lattice Boltzmann method

T

Gas temperature

CA model

Cell automation probabilistic model

R

Universal gas constant

p

The pressure of gas

FB , x

Brownian force in the x direction



The density of gas

FB , y

Brownian force in the y direction

c

Particle speed in lattice unit

FL , x

Saffman lift force in the x direction

LBE

Lattice Boltzmann Equation

FL , y

Saffman lift force in the y direction

x ,y

zero mean, unit-variance independent

BGK

Bhatnagar-Gross-Krook Gaussian random numbers kB

CA

Boltzmann constant

Cell automation probabilistic approach approach Gas velocity in the x direction and y LBE

ux , u y

Lattice Boltzmann Equation

direction 40

e

Microscopic velocity in lattice unit

u p, x

Particle velocity in the x direction

t

Time step in lattice unit

u p, y

Particle velocity in the y direction



Collision function

gx

External force acceleration in the x direction External force acceleration in the y τ

gy

Relaxation time due to collision

direction The lattice mode with b velocities on a DdQb



A Boolean variable

eα

Lattice velocity in  direction

simple cubic lattice of dimension d The two-dimensional nine-velocity lattice D2Q9 model 

Direction subscript

r1 , r2

two independent random numbers

x

Lattice spacing

Lc

Channel length

X p

Displacement vector

Hc

Channel width

p

Discrete motion probabilities

Ws

Channel wall thickness

p

Density of the particle

Wp

Plug thickness

Vp

Volume of the particle

f

The distribution function

up

Velocity vector of the particle

D

The width of the square block

u

Velocity vector of the gas

H

The height of the channel

f

Stokes-Cunningham drag coefficient

L

The length of the channel

FB

Brownian force

Re

Reynolds number

FL

Lift force

Stk

Stokes number

Xp

Displacement vector of the particle

The threshold value

41

• The lattice Boltzmann-cell automation probabilistic model was developed to investigate the flow and soot loading in the micro-channels of DPF. • Non-virtual porous structure wall was considered in the filtration process of soot particles. • The effects of upstream velocities on flow field and soot distribution in clear porous wall were studied. • The distribution of soot particles in clean porous wall was also simulated by lattice Boltzmann-Lagrangian model and the results of both two models were compared with each other. • The deposited soot particles dynamic evolution and its effect on flow field were studied.

Declarations of interest: No conflict of interest.