7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Validating colloids agglomeration with existing theoretical Kernel

Mikael MOHAUPT*, Anne TANIERE* and Jean-Pierre MINIERt

LEMTA UMR 7563 CNRS, ESSTIN, Universit6 Henry Poincar6 Nancy 1,

2 rue Jean Lamour, F-54529 Vandoeuvre-Lbs-Nancy, France

Electricity de France, Div. R&D, MFTT, 6 Quai Watier, 78400 Chatou, France

mikael.mohaupt@esstin.uhp-nancy.fr, anne.taniere@esstin.uhp-nancy.fr and jean-pierre.minier@edf.fr

Keywords: colloid, Brownian motion, collision algorithm, kernel

Abstract

The temporal tracking of the agglomeration of very small particles up to very large-inertia particles is investigated

considering spherical particles. A deterministic collision algorithm capable of handling the multi-scale aspect

underlying this study is developed in this paper. Both the collisions due to random motion of Brownian particles and

those induced by the deterministic movement of large-inertia particle must be treated which means strict conditions on

the choice of time step for the simulation are required. The total cost of computation would be dramatically increased

which is not in accordance with the main goal of the present study. Two different limit behaviours are examinated

here to study the efficiency of the present algorithm through the calculation of collision kernels according to existing

proposals in research literature on the subject.

Nomenclature

Roman symbols

B diffusion coefficient (m s3/2)

dp diameter of particles (nm)

9 gravitational constant (m s 1)

l~ (t) Stochastic integral on position (m)

1v (t) Stochastic integral on velocitiy (m s1)

KB Boltzmann constant (m2 kg s 2 K 1)

L length of the numerical control domain (m)

Mnp mass of particles (kg)

N particles number (-)

Nc computed number of collisions (-)

P pressure (N m 2)

r radius of particles (m)

T temperature (K)

Tp particle agitation (m2 s 2)

t time (s)

to initial time (s)

tma, total simulation time (s)

v velocity (m s-1)

W(t) Wiener process (s1/2)

x position (m)

collision time step (s)

mean free path (m)

standard Gaussian random variable (-)

standard Gaussian random variable (-)

medium fluid viscosity (Pa s)

density (kg m 3)

relaxation time (s)

Subscripts

eq

max

m 1, 2,3

p

f

Superscripts

i

ij

Special notation

d- (t)

equilibrium

maximum

space direction

discrete phase (particle)

continuous phase (fluid)

particle i

particle j

relative value (j i)

n

Euclidean norm

mean value, i.e. (1/N) 1(-)

time increment

Greek symbols

At

Atini

time step (s)

initial time step (s)

Introduction

Particle deposition processes can be divided into several

steps (or elementary mechanisms): dispersion, single-

particle deposition, agglomeration and clogging. The

numerical simulation of the agglomeration mechanisms

for colloids transported in particle-laden flows is a key

issue for a wide range of engineering processes and

applications, such as nuclear safety, food technology or

aerosol concerns. The complexity of these simulations

stems from the fact that widely different scales are

involved in the overall physical phenomena. At a

microscopic level, colloids moving under a Brownian

influence can agglomerate due to electro-chemical

forces between them. At a macroscopic level, these

particles are affected by hydrodynamic forces and

transported by the fluid turbulence. These two levels

are not independent both because particles may grow in

time (due to agglomeration) which can modify the way

they are transported and because different macroscopic

agitations can result in different agglomeration rates.

To study agglomeration phenomena, it is first nec-

essary to define an algorithm for the detection of

particle collisions. Contrary to probabilistic methods

to treat particle-particle collisions (Oesterle & Petitjean

(1993)), we chose to adopt a deterministic point of

view as previously done by Chen & Kontomaris &

McLaughlin (1998) in a turbulent particle laden flow,

for instance. However, the smallest time scale encoun-

tered in this kind of study is the turbulent time scale

which is, for particles large enough, of the same order of

the particle relaxation time (denoted Tp). This imposes

therefore that the time step of their numerical simulation

is reasonable in the sense that At < Tp. In the present

study, we want to simulate the Brownian movement

of particles. The constraint At < 7p considerably

increases the total computational cost since the colloid

relaxation time is smaller than the particle relaxation

time considered in the above mentioned turbulent

studies. More precisely, our objective is to establish an

algorithm of collision capable of treating both Brownian

collisions and collisions of large-inertia particles with a

reasonable time step, thus decreasing the computation

total cost. The case of large-inertia particles is not

in itself a problem since the relaxation time scale is

naturally large with respect to the time step. What we

need is a model for the Brownian case which bypasses

this limitation of the time scale. The recent work

of (Peirano et al. (2006)) focusing on this problem

proposed numerical modelling and associated schemes

allowing the simulation of a Brownian trajectory in a

discrete sense without respecting the severe condition

At < Tp.

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

The present study will thus focus on the multi-scale

nature of these collisions which will be treated by a

single deterministic algorithm without any constraint

or threshold on the time step. We propose here to

study different test cases corresponding to two limits

of particle behaviour. The first is related to the motion

of small particles which move rapidly according to

Brownian motion and where trajectories are mainly

random diffusivee limit). The second limit case deals

with very large inertia particles. These particles move

in straight lines with constant velocities between colli-

sions. In order to validate these physical configurations,

we choose to check the efficiency of the algorithm

of collisions by calculating collision rates for which

the theoretical expressions are available in the previ-

ous works (Smoluchowski (1917); Friedlander (2000)).

This paper presents our work as follows. First

we shall describe the details of the present collision

algorithm. Then, the testing of the algorithm for the two

limits cases is discussed and the validation is carried out

comparing the computed results of collision rates and

those issued from theoretical expressions. An analysis

and discussion of the obtained results are given in

conclusion.

1 Particle equation movement and numerical

integration scheme.

We shall approach the Brownian movement of particles

from the Langevin point of view. For simplicity's sake,

we will consider motion in one dimension. Newton's

equation of motion for the particle (radius r, mass mp,

position xp(t), velocity vp(t), particle relaxation time

Tp, and fluid medium viscosity /) is:

( dxp(t)

d(t) F(t)

S(dt1)

where the force, F(t), due to the interaction of the Brow-

nian particle with the surrounding medium, thus con-

tains both frictional force proportional to the velocity

of the Brownian particle and random force. Therefore,

Langevin's model reads:

dxp (t)W

dvp(tW

p (t)dt

vP(t)dt + BdW(t),

Tp

B is the diffusion coefficient and dW(t) is the Wiener

process whose increments over small time steps are

stationary and independent, < dW(t) > 0, and <

dWj(t)dWj(t') = 0 > with t / t'. This system be-

ing physically valid when dt < Tp.

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

If we neglected the random force in Eq.2, the latter equa-

tion would become:

Sdxp(t)

Sdvp(t)

dt

Vp(t)dt

vP(

which has the familiar solution

xp(t) xp(O) + o vp(s)ds

Vp () = e'vp (0)

(3)

given by

< 1 >-

< IV >

< Ixl >-

(BTp)2 {At

Att]

OTp]

-A;]}t

C 1P ]

{BTp,[1 c- e

(10)

Therefore, relations for position and velocity of a

(4) Brownian particle given by the system (7) associated

with relations (9) and (10) are:

predicting that the velocity of Brownian particle decay to

zero at long times which cannot in fact be true. By con-

sidering in equilibrium where KB is Boltzmann constant

and T is the temperature of fluid where the particles are

immersed, the application of equipartition theorem stip-

ulates that

KBT

< -/) >eq

while the expression of the velocity from the system (4)

gives

2t

< I) >eq= e (6)

< v (0) >eq 0.

Therefore, the random force is necessary to obtain a

finite value of the velocity variance corresponding to

the correct equilibrium. The stochastic calculus gives

the correct solution. The recent work of Peirano et al.

(2006) helped us to apply the rules of stochastic calculus

in Ito' sense. The solution to system (2) at t to + At

is therefore given by :

{ xp(t)

vp(t)

At

x~p(to)(o)+ )Tp(1 -e 1 )+ I(t),

o At + ().

vp(to)c P 14(t).

where the stochastic integrals Ix (t) and I, (t) are defined

I(t)

1- W

BTp tt dW(s) BTpe ftt e dW(s),

Be f e dW(s).

(8)

In order to facilitate numerical applications, Peirano et

al. (2006) applied the Choleski algorithm to (I I,) as:

I1(t)

Ift(t)

(

/ <2 < W<2 >

<12 > C

J~p,,

(9)

where (C and &C are two independent standard Gaus-

sian random variables (zero mean and variance equals

to unity). The components of the covariance matrix are

At

axp(t) = p(to) + p(to)Tp(l1 c )

SAt 112

A+BTt [At 2 (1 j

a(t1e e )

] B r>(1 ) [27 (1 2 n 1)] ,

v(t) -vp(to)e 1 +[ (-e 2 A .

(11)

This system of equations can be used directly in discrete

time which amounts to treating the model and its numer-

ical scheme without distinction. In addition, to be con-

sistent with analytical solutions of the system (7), the

numerical schemes established by Peirano et al. (2006)

can provide the correct system limit of EDS. They are

linked with the physical behaviour of particles although

the separation of scale, 0 < dt < Tp, imposed by the

system (7) was not respected. In other words, these

numerical schemes allow to perform simulations with

a time step At such as At > Tp but also in the case

At < Tp as shown by Fig. displaying the particle vari-

ance position in the case where At < T7 and in Fig.2

when At > Tp. Taking the following value of the dif-

fusion coefficient (12) obtained from the equipartition

theorem for translational kinetic energy, i.e.

B KBT (12)

V rnpTp

the variance of the particle position as well as the vari-

ance of the velocity can be obtained. As can be seen

in Figs. 1 and 2, the calculated values and those given

by the exact solution correspond well. Moreover, the

physics of a Brownian particle is respected since the val-

ues of particle variance positions reach a state of equi-

librium (Fig.2). We also observed that the mean parti-

cle position evolves around zero (results which are not

shown here). Therefore, the above results are consis-

tent with Einstein's observation, i.e. for long diffusion

times, one has < x(t) >~ (BT)2t. Concerning the re-

sults on the particle velocity, same remark can be made

regarding to the consistency of the computed results and

those which came from exact solutions (see Figs.3 and

4) whatever the value of At with respect to Tp.

i

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

times) times)

Figure 1: Brownian particle position variance. N =

10000 particles, case where At < Tp, At 1 x 10 0 s,

t,,m = 1 x 105s, Tp 3 10 6.

6E-08

5E-08 computed

5E-08 exact solution

thermal equilibrium

4E-08

A

V%'3E-08

V

2E-08

1E-08

0 200 400 600 800 1001

times)

Figure 2: Brownian particle position variance. N

10000 particles, case where At > Tp, At 1 x 10 2s,

t,,m = 1000s, Tp 3 x 106s.

To sum up, the use of relations 2 and the associated

numerical model allows us to study the diffusive limit

case corresponding to Brownian motion whatever the

value of At with respect to Tp. This kind of model

includes the drag force as well as the Brownian force

which both acting at the same level on colloids. For the

Figure 3: Brownian particle velocity variance. N

10000 particles, case where At < Tp, At 1 x 10 ls,

tm = 1 x 10 5s, Tp 3 x 10 6.

8E-05 6

6E-06

A

V

v

4E-06

2E-06

0

-- computed

exact solution

thermal equilibrium

200 400 600 800

times)

Figure 4: Brownian particle velocity variance. N

10000 particles, case where At > Tp, At 1 x 102s,

t,,m = 1000s, Tp 3 x 106.

smallest dimension in the range of diameter of colloids,

the particle inertia effect is negligible with respect to the

Brownian one. Therefore, stochastic integrals (9) pre-

vail. The system of equations (11) is thus approximated

6E-16

5E-16

4E-16

A

%1-3E-16

v

2E-16

1E-16

by neglecting the drag force:

xp (t) = xp(to + BTp At 2 (1

(1 L (1 e

+BT(1 e )2 [2Tp(1 -e

Vp(t)= Y- [1e ]

At) 1/2

1P )/

-1/2,

(13)

In the other limit case corresponding to the large-inertia

particle, the stochastic effect disappears and the parti-

cle displacement is in a straight line during a time step.

Under this assumption, the system (11) yields:

A,

{ Xp(t) = X(to) + Vp(to)Tp(l (14)

t (14)

VP (t) = Vp(to)e .

Thanks to the system of equations (11), the diffusive

behaviour of particles can be simulated whatever may

be values of the particle relaxation time with respect to

the time step. In the Brownian case (the most restrictive

because it corresponds to the smaller value of Tp) this ad-

vantage can help us our aim was to devise an algorithm

that can simulate collisions of Brownian motion without

respecting the constraint (At < Tp). For large-inertia

particle corresponding to the ballistic limit, this condi-

tion is by no means a restriction (the extreme case being

that Tp -- oc). Between these two limit cases, which

may reflect a complex particle laden flow problem, the

upper limit of the time step of the simulation will be im-

posed by the physical phenomena involved (turbulence,

for instance). Based on the idea developed by Peirano

Table 1: Considered test cases.

Test Type of particle motion Time step condition

case 1 Ballistic limit At < 7T

case 2 Diffusive limit At < 7p

case 3 Diffusive limit At > Tp

et al. (2006),our aim is to construct a collision algo-

rithm including both the effects of Brownian simulation

and effects associated with large inertia of the particles.

Some test cases (described in Table. 1) will be examined

in a later section after a presentation of the details of the

collision algorithm.

Algorithmic details

A brief review.

A large number of collision algorithms have already

proved their consistency in one physical case or an-

other, corresponding specifically to diffusive behaviours

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

or ballistic behaviours of particles. The large multi-scale

character (from diffusive limit to ballistic one) seems

not to be taken into account in existing deterministic

algorithms of collision. The only way to do this is to

adapt the time step simulation to keep it much smaller

than the particle relaxation time in order to obtain parti-

cle displacement in a straight line. Therefore, whatever

the diameter of particles and forces acting upon them,

the free-molecular regime up to the ballistic one can be

simulated. However, the major drawback lies in the im-

portant computational cost. This can be reduced how-

ever if a probabilistic collision algorithm is used. For

instance, Oesterle & Petitjean (1993) proposed a La-

grangian simulation technique of a non-dilute gas-solid

suspension flow where a probabilistic method is used to

simulate particle-to-particle collisions during the trajec-

tory calculation of a particle. The principle consists of

introducing artificial inter-particle collisions during the

trajectory calculation, the probability of these collisions

depending on the local concentration and velocity of the

solid phase. The main advantage of this method lies in

the fact that dense or moderately dense flows can be sim-

ulated by computing a reasonable number of simultane-

ous trajectories. This simulation technique is not appro-

priate for our study since our main goal is to analyse

the physical phenomenon of the agglomeration process

and even to develop new agglomeration kernels. Thus

a collision algorithm of deterministic nature is essential

to achieve our goal. As stated earlier, this type of al-

gorithms is often time-consuming because the time step

is usually very short compared to the relaxation time of

the particle. Moreover, these algorithms are often based

on the overlapping detection. It consists of regarding

each pair of particles at each time step in the compu-

tation domain and testing if the distance between two

particles is smaller than the sum of their radii. An exam-

ple of schematic of overlapping algorithm is presented

in Fig.5. An expanded study on this simple kind of al-

gorithm (Wunsch & Fede & Simonin (2008)) stipulates

that to obtain good results on this particular case, severe

conditions have to be respected, the most important be-

ing based on the ratio between 6- avoiding thus to miss

overlapping. The relation is:

61 3 At

-dp= 2 Tp

-l fT~

1 1p

0.13,

where 61 is the mean free path of particles during one

time step (in m) and Tp 1/3(v1 + vp2 + v3) is

the particle agitation. This is a severe limit for the

maximum time step allowed when dp is very small, as

in the case of a colloid particle. Unfortunately, there is

no generalization on the time step which has to be used

for a given physical factor. Algorithms appearing in the

literature are elaborated for a given physical case. For

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Figure 5: A simple collision algorithm based on over-

lapping detection.

example, when studying turbulent flow in a dispersed

phase with larger particles whose relaxation times are

greater than unity, we can give more flexibility to the

value of time step, thereby reducing the total numerical

time cost. Another way to reduce this computational

cost while keeping At < 7, is to use the kind of

geometrical tests proposed by Lavieville & Deutsch &

Simonin (1995). They carried out geometrical tests

on the mesh in order to limit the number of possible

neighbours of a given particle, restricting therefore the

number of numeric tests from N2 to N log(N).

The present collision algorithm.

An overview of the elaborated algorithm is given

in Fig. 6. The key issue in a collision algorithm

considering spherical solid particles is to focus on the

behaviour of each couple of particles in the flow. This

is the most inefficient process since no collisions are

missed but it remains the simplest way of making sure

all collisions are covered. We are aware that better

numerical costs could be obtained by focusing on the

definition of the vicinity and by using a new criterion

which define possible neighbours of a given particle.

However, this is not our main goal as we seek to reduce

computational costs by working on numerical schemes

associated with the physics of colloids rather than on a

technical definition of the vicinity of a given particle.

When considering perfect elastic rebounds, two dif-

Figure 6: Limit case collision algorithm based on

straight line displacement.

ferent ways can be used in order to study the collision

detection for a given pair of particles. The first is to de-

tect and to treat an overlap as shown in Fig. 5. The sec-

ond one lies in detection of a collision during a time step

which is based on criteria derived from the molecular

dynamic theory. These criteria are explained below. The

first criterion defines the space domain covered during

one time step by the relative movement of both particles.

This domain is given by a cylinder defined as (Eq. 16)

and presented in Fig. 7 which is defined for a particle

couple (i, j).

base = 7(ri + r3)2

length

direction

11 ,i, v 1 '

:-(Vi Vi

volume = 7(r' + r3- || ;.|| t.. (16)

One necessary condition for detecting a collision at a

given time between a pair of particles (i, j) is that an-

other particle (denoted j) is located in the cylinder de-

fined above (Eq. 16). This condition is not sufficient

because it does not suggest that particle trajectories are

crossed. According to that, the possible collision can

be given using the second criterion defined by relation

(Eq. 17)

wX3 V3 < 0, (17)

which consists of observing whether the two particles

yo x2

zorx3

II V" At

Figure 7: Collision algorithm. First detection collision

criterion.

are moving towards each other (condition is true) or if

they move away from each other (condition is false).

Figure 8: Collision algorithm. Second detection colli-

sion criterion.

Once these two stages are successfully passed we

need to calculate the particle collision time step (At')

which defines the necessary time for the particle denoted

i to enter in collision with another (denoted j here). The

calculus of such a time step is well-known and a good

explanation is available in the work of Sigurgeirsson et

al. (2001). Consequently, we have to solve (Eq. 18)

defined as:

||xJ(t + At,) x(t + At,)| =r' + r, (18)

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

which is equivalent to solve (19) :

(v' (t))2At +2(v (t) x (t))Atc+(x' (t))2 (r'+r')2.

(19)

Atc will be deduced by taking the minimum of positive

solution of 19. Finally, the end of the process will take

the minimum of Atc among those calculated for each

particle i. This is a key point in our algorithm and

could be seen as a highly strict. However, this path is

the best way to approach processes of agglomeration

that may occur at different scales (from electrochemical

forces for the smallest scales up to hydrodynamic

forces for the largest ones). The choice of a minimum

collision time step allows us to bring the numerical

representation of the flow to a state approaching nearer

to a collision event from any state of a particle laden

flow. The agglomeration process can thus be treated

with the utmost precision contrary to the most common

treatment based on a probability of agglomeration.

Testing the algorithm : ballistic regime

This kind of behaviour corresponds to a case when par-

ticles move in straight lines with constant velocities dur-

ing one time step. This is typically the case of kinetic

theory of gases in which the only way to change particle

direction is the collision with another.

Kinetic gas theory considers what is happened in a per-

fect gas, focusing on how molecules hit each other to

define the pressure in the gas. This is dependent on

the number of collisions encountered, and thus depen-

dent on molecule concentration. The objective of this

section is to test the collision algorithm in the case of

gases by calculating the well-known collision frequency

of molecules according to the kinetic theory of gases

given by the following relation (Reif (1965))

N 2 16KBTN

Kth = d

L w -rmp 2

The idea is to compare this kind of theoretical frequency

with the measured frequency which is simply deducted

from the the number of collisions Nc per unit volume

per time during the time simulation, i.e.

2NcL3

N(N 1)t

where L is the dimension of the numerical domain.

Simulation overview.

The numerical domain is a cubic control volume with

periodic conditions in each direction. The length of the

cubic volume is taken by 10nm. Particles i,-'"'I are

uniformly distributed in the numerical domain with an

initial distribution of molecules velocities obeying to a

Maxwellian distribution in a directional sense (Fig. 9) :

F(vpm) rn xp

2KBT

2KBT

written in a normalized sense (Fig. 9), it yields :

F(iivp) I 47r (2MPT) 3 Vpl 2exp

2K T

0.0008

0.0006

S0.0004

0- 2

0.0002

2

2KBT )

(23)

Figure 9: Initial distribution of molecules velocity. 250

helium molecules.

Since particles have a constant velocity and are dis-

placed in a straight line, particle positions are advanced

in time following this numerical scheme derived from

relations (14) using mathematical limited development

at A around zero.

5,

xp(t) = x(to) + V(to)At. (24)

The physical conditions used herein correspond to

ideal gases under a 4 bar pressure which allows us to

obtain a sufficient concentration of molecules of gases

and therefore to expect more particles collisions in the

numerical domain. Theses gases are listed in Table.2

where details of their molar mass are given along with

details of their Van Der Waals's diameter, including

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Table 2: Tested ideal gases and used properties.

Gas Molar mass (g) VdW radius (pm)

Helium 4 140

Neon 20 154

Argon 40 188

Krypton 83.3 202

Xenon 131.3 216

totally the spatial geometry of molecules according

to the hard sphere modelling of molecules. This was

thus a perfect case to validate the collision detection

process. The time step used is At 1014s and the

total time of simulation is taken equal to tm,,a 10 s,

corresponding to a number of encountered collision

between 500000 and 1500000 depending on the ideal

gas considered.

Results.

Results are presented for some of the gases listed

(Table. 2) and are displayed in Fig.10 as a function

of the current time. The ratio between the theoretical

collision frequency and the computed one from sim-

ulation evolves around unity whatever the gas tested

attesting that the collision process is validated. In order

to complement this validation, we have carried out

simulations by varying simulation's parameters such as

the time step At. Results are reported in Fig.11 where

only the Helium gas is tested. Good results were also

obtained.

We can thus conclude that the collision detection is

validated since no dependency on particle diameters or

time steps is observed. The next step of validating col-

lision algorithm is to test the detection collision process

in the other limit case, i.e. for a diffusive movement of

particles, typically the Brownian case.

Testing the algorithm: diffusive regime

The two limit diffusive cases corresponding to Cases 2

and 3 are used here in order to calculate the Brownian

kernel with the present collision algorithm. For simplic-

ity's sake and to keep an assumption of monodispersed

flow, perfect elastic collisions are considered. The goal

is to compute Brownian kernels issued from numeri-

cal simulation (21) and to compare them with those ob-

tained from literature (Smoluchowski (1917)):

s 8KBT

3p

14 1 1 1

1 2'

E

'A

.A

08

C 06 -

06

Helium

2 A krypton

02- Neon

Xenon

a 5E-08 1E-C

time (s)

Figure 10: Dimensionless collision kernel in kinetic gas

theory 250 molecules.

S4

1 2

.2

?08

o06

.2

S04 dt=105 s

S

E dt=103 s

02 -- dt10 14s

02

0 5E-09 1E-0

time (s)

Figure 11: Dimensionless collision kernel in case of

Helium: time independency. 250 molecules.

Recently, Trzeciak et al. (2005) raises the issue of treat-

ment of elastic collisions when simulates a random mo-

tion of particles. Indeed, after the collision, particles can

remain close to each other and overestimating the re-

collision delay depending on the pair distribution func-

tion. If particles stay in each other's vicinity after the

collision, a completely different pair distribution func-

tion can be obtained, hence different coagulation rate.

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

The simple way to overcome this problem, in our case,

is to treat collision events by disappearing one of the

collision partners and to redistribute it in the numerical

domain thanks to the initial distribution of particle posi-

tions (uniform distribution). Indeed, the pair distribution

function is not affected since Brownian motion admits a

zero mean value for the particle position. Details of sim-

ulation parameters are available in Table.3.

Whatever the cases (2 and 3), the system of equations

Table 3: Summary of gas and particle properties used in

computations.

Variable Value Units

KB 1.3806503. 10 23 m3 kg K- s-2

7 T 296.15 K

p 1.83245 10-5 Pa s

pp 1000.00 kg m3

pf 1.00 kg m-3

tmax 100.00 s

(11 is used in order to obtain positions and velocities. In

Case 2, we integrate the system (11) in time on a time

step which is considered very small compared with the

physical phenomenon response. Therefore, the particle

displacement could be seen as a straight line during a

time step. Such an assumption is not valid in Case 3

where At > Tp. Nonetheless, according to Figs.2 and

4, the statistical Brownian motion is well reproduced.

The results of Case 2 and Case 3 are reported respec-

tively in Fig.12 and in Fig.13 versus the current time.

Figures 12 and 13 show the ratio between the calculated

Brownian Kernel and the theoretical kernel as a func-

tion of the current time. 50 particles are used to pro-

ceed to the computation of the ratio in Case 1 whereas

in Case 2 1000 particles are considered. Therefore these

results are preliminary. These results required relatively

8 short computation time (about a week) and nonetheless

show whether the present collision algorithm is capa-

ble of calculating the frequency of Brownian collisions.

Two types of results are given. The first type of results,

as shown in Fig. 12, are carried out when At < Tp

whereas in Fig. 13, the calculations were made for the

case where At > Tp. These two opposite cases are

not incidental to our approach. It should be noted of

course that our ultimate goal is to follow the formation

of a cluster in time and for a reasonable cost of compu-

tation time. For a lower cost, it can be simulated with a

time step such that At > Tp while respecting the physi-

cal phenomena, i.e. the diffusive behaviour of Brownian

particles. This is what we attempted to achieve by pre-

senting a case where At < 7p where we are sure that

physics is respected. The results (Fig. 12) show that our

1 8

16

14

" 12

E

E 08

m 06

04

02

0

time (s)

Figure 12: Brownian kernel, case where At < Tp.

N = 50 particles in a control cubic volume (edge length:

0.5pim ). At 3 x 10 s, 7p 3.03 x 10 6s for

dp = 10 and 7- 7.58 x 107s ford = 10 m.

E

Figure 13: Brownian kernel, case where At > Tp. N

1000 particles in a control cubic volume (edge length:

1mm ), At 102s, Tp 3.03 x 106s for dp

10 ',, and Tp 7.58 x 107s for dp 10 m.

algorithm gives good results, since the ratio of kernels

evolves around the unit. Results seem not to depend on

particle diameter. However, results should improve by

taking more than 50 particles and currently we are work-

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

ing along these lines. The results presented for the other

case (Fig. 13) defined as At > Tp are inconclusive. The

ratio can reach values 30 times greater than a unit and

depends on the value of particle diameter.

As exposed above, the used model given by relation (11)

can predict the position and the velocity of a particle

whatever the value of the time step with respect to Tp.

Nevertheless, when the time step is chosen much higher

than the particle relaxation time (At > Tp), the given

results must be interpreted with care since during a time

step the particle displacement is not a straight line. The

probabilistic nature of Brownian motion has to be taken

into account. A simple illustration of this problem is de-

scribed in Fig. 14. Moreover, for a given pair of particle

(i, j), exposed criteria do not necessarily lead to a col-

lision detection in the case where At > Tp (Fig. 15(b))

rather than in the case where At < Tp (Fig. 15(a)) for

which straight line approximation is acceptable inducing

thus a collision.

-----Q

(xp(l), Vp()) (xp(t + Al), vp( + At))

r \ t / I

\ I I I

(XW(t). vpt)) ( (X,(t + At), V,(t + At))

Figure 14: Opposition between particle real displace-

ment and numerical displacement in the two studied

cases for Brownian motion.

At first sight, our algorithm and more precisely the

conditions of collisions (Fig. 7 and Fig. 8) highlight the

problem, i.e. the detection of collision is elaborated

only for a straight displacement during one time step.

The collision cylinder (Eq. 16) must be redefined to

enlarge the vicinity of the pair of considered particles

taking their random trajectories into account. Moreover,

the calculus of Atc (Eq. 19) is also distorted for the

same reasons. These two steps have to be re-examined

from another perspective. The idea is to include the ran-

dom nature of displacement between t and t + At in the

collision algorithm. The calculus of the collision time

step for a rectilinear motion applied to Brownian ran-

dom walk, without any particular adaptation, conducts

to an excessive number of collisions in mismatch with

the simulated physical phenomenon.

S |At < Tp (a)

-4-

(*'(t + At)

S< AtI (b)

,-, / -'

S I '

(t) r -

t- .,(t + At)

Figure 15: Possible relative displacement of a pair of

particles.

Further study could be devoted to a special aspect of

Brownian theory, i.e. the Brownian bridge which could

be seen as a valid possible solution for this problem.

This probabilistic tool offers the possibility to restitute

a possible random trajectory between t and t + At in ac-

cordance with statistical properties of Brownian motion.

Conclusion

In a context of a study devoted to the particle agglom-

eration process, it is necessary to first establish a colli-

sion algorithm which could be applied to a large range of

physical phenomena such as Brownian motion of parti-

cles as well as the motion of very large-inertia particles.

This kind of multi-scale approach to treating collisions

is present in an agglomeration process since the particle

diameters can grow in time due to the agglomeration.

It seems therefore imperative to include this multi-scale

dimension directly in the collision algorithm. This algo-

rithm must be able to handle both collisions at very small

particle relaxation scales (like Brownian collisions) and

collisions between very large-inertia particles. The main

goal of the present work was therefore to present prelim-

inary steps in constructing of such a collision algorithm

of particles while respecting the following constraint of

a reasonable numerical time step in order to not to in-

crease too much the total computational cost. To vali-

date this approach, some tests were carried out for the

two limit physical cases. The first corresponds to the

motion of large-inertia particles which have a straight

line displacement during one time step and the second

is devoted to Brownian particles where trajectories are

mainly random between t and t + At. The efficiency

of the algorithm was tested on the calculation of colli-

sion kernels by comparing results issued from numeri-

cal computations to the well-known theoretical relations

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

which are available in the literature.

The problem seems to be easier when a time step is cho-

sen to be smaller compared to the particle relaxation

time. In such a case, whatever Brownian motion or the

movement of heavy particles, the time step chosen was

so small that we can consider that particles move in a

straight line. However, this kind of simulation is ex-

tremely expensive for the Brownian case. It is not dif-

ficult to reduce the total computational cost in the case

of a hard sphere model in the kinetic theory of gases,

the particle relaxation time Tp being considered as infi-

nite. Choosing a reasonable time step is therefore easier

in this limit physical case. The collision algorithm was

tested for this case and a good level of compliance with

theoretical results was obtained. The collision kernel ra-

tio evolves around unity whatever the diameter of tested

particles and also whatever the chosen time step. This

case is not really problematic in contrary to Brownian

motion since the latter varies extremely rapidly between

two time steps and the well-know statistics, as the vari-

ance of the position of particles at long times (the result

of Einstein), could be not reached if it is observed at

At > Tp. However, recently Peirano et al. (2006) pro-

vided an opportunity to solve Brownian motion even if

At > Tp adapting the analytical solution of Langevin

equation in a discrete sense. We used their numerical

schemes and good results on the physics of Brownian

motion were obtained whatever the value of the time step

with respect to the particle relaxation time. Nonethe-

less, when the Brownian collision kernels are calculated

in case of At < Tp, the ratio evolves around unity con-

trary to the case for which At > Tp. For this case, we

think that the Brownian character has to be directly in-

cluded in the collision algorithm. Nevertheless, the dif-

ficulty is that between two time steps, the trajectories of

a pair of Brownian particles vary rapidly. The more the

time step is large, the more the difficulty of predicting

the displacement of the particle pair in a cylindrical do-

main, as it has been defined, seems arduous. However,

the literature on the Brownian bridge theory could be a

mean to get other results. That is the next step of the

present study.

Acknowledgements

We thank Pascal Fede and Tomasz Trzeciak for their nu-

merous useful remarks and discussions about this work.

References

Chen M. and Kontomaris K. and McLaughlin J.B., Di-

rect numerical simulation of droplet collisions in a tur-

bulent channel flow. Part I: collision algorithm, Int. J.

Multiphase Flow, Volume 24, pp. 1079-1103, 1998

7th International Conference on Multiphase Flow,

ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Chen M. and Kontomaris K. and McLaughlin J.B., Di-

rect numerical simulation of droplet collisions in a tur-

bulent channel flow. Part II: collision rates, Int. J. Multi-

phase Flow, Volume 24, pp. 1105-1138, 1998

Friedlander S.K., Smoke, Dust, and Haze, Fundamentals

of Aerosol Dynamics Second Edition, Oxford Univer-

sity Press, 2000

Lavieville J. and Deutsch E. and Simonin O., Large

eddy simulation of interactions between colliding parti-

cles and a homogeneous isotropic turbulent fields, Gas-

Particle Flows, Volume 228, pp. 347-357, 1995

Minier J.P. and Peirano E., The PDF approach to turbu-

lent polydispersed two-phase flows, Phys. Reports, Vol-

ume 352, p. 1-214, 2001

Oesterl6 B. and Petitjean A., Simulation of particle-to-

particle interactions in gas-solid flows, Int. J. Multiphase

Flow, Volume 19, No. 1, pp. 199-211, 1993

Peirano E. and Chibbaro S. and Pozorski J. and Minier

J. P., Mean field/PDF numerical approach for polydis-

persed turbulent two-phase flows, Progress in Energy

and Combustion Science, Volume 32, pp. 315-371, 2006

Reif F., Fundamentals of statistical and thermal physics,

McGraw-Hill, ISBN 07-051800-9, chap. 7, 1965

Sigurgeirsson H. and Stuart A. and Wan W. L., Algo-

rithms for particle-field simulation with collisions, Jour-

nal Of Computational Physics, 172, pp. 766-807, 2001

M. von Smoluchowski, Z. Phys., 92, 129, 1917

Trzeciak T. and Podgorski A. and and Marijnissen, J.

Stochastic calculation of collision kernels: Brownian co-

agulation in concentrated systems, Proc.WCPT5, 117e.,

2006

Wunsch D. and Fede P. and Simonin O., Development

and validation of a binary collision detection algorithm

for a polydispersed particle mixture, ASME Fluids En-

gineering Division Summer Conference, Jacksonville,

Florida USA, August 10-14, 2008