Tony LaddUniversity of Florida

(1)

Simulations of Particle-Fluid Suspensions with the Lattice-Boltzmann Equation

Tony Ladd

University of Florida

With thanks to MPIP-Mainz and AvH Foundation

(2)

( ) ( ) ( )

2 t

t

hydro th

f ermo

p m

r r h

r r

¶ + Ñ = -Ñ +

+

Ñ

¶ + Ñ ×

= +

+

Ñ ×

=

+ W´

σ

u uu u

u 0 u U R

F F

f R&&

Derive a variety of micromechanical models of complex fluids from same basic ingredients

Inertial Particle Motion

~Incompressible (M < 0.1) No slip

Particle-fluid Fluctuating stress

LBE is a suitable computational framework

Solutions of polymers Liquid crystals

Colloids

and biopolymers Porous Media

Inertial & Local 1. Solid particles-Newtonian Mechanics

2. Continuum fluid-Navier-Stokes equations 3. Stick boundary conditions couple particles

and fluid.Valid for particles > 30nm (Add: charge, chemical bonds, inertia) Computational framework for HI in a wide

range of materials, flows, and scales.

(3)

Outline: Applications of DNS to suspensions and particle fluid systems

 Lattice-models for fluid dynamics

 Lattice-Boltzmann method

 3 examples

 Settling of particle clusters at small Re

 Reactive flows in porous media (Stokes flow)

 Polymer solutions (with Brownian motion)

 Closing thoughts

(4)

Lattice-gas models for suspensions

 Lattice-gas models were introduced to simplify kinetic theory (Square lattice-HPP)

 FHP (‘86) showed that a hexagonal lattice gas could solve Navier-Stokes equations in 2D.

 LCF (‘88) used the FHP model to calculate viscosity and self-diffusion in a 2D colloidal suspension

 Projected 4D FCHC model for 3D simulations (Henon ’87)

 Moving boundary condition (FL ’89)

 Hydrodynamic interactions (LF ’90)

 But: LG models are too noisy; Sc ~ 1: Not Galilean invariant

 LBE (HS-with linearized collision operator)

(5)

LBE model introduces a discrete velocity

distribution: local collisions and propagation

Initial State:

x momentum + xy shear stress:

Post-Collision:

x momentum only

Propagation

( , ) ( , ) [ ( , )

_i _i^EQ

( , )]/

i i i

n t n t

n r c + D + D = t t t n r t - r - r t

(6)

Hydrodynamic fields are moments of the discrete velocity distribution n

_i

(r,t)

18

0 18

0

18

0 18

0

( , ) ( , ) ( , ) ( , ) ( , )

( , ) ( , ) ( , )

i i

i

EQ

i i i

i EQ

i i i i

i

t n t

t t n t

p t t t t n t

t n t n t

r r

r s

=

+ =

é ù

= - ë - û

å å

r r

r u r r c

r r u r u r r c c

r r r c c

Mass

Momentum

3D model has 19 velocities c_i: 000, 100 & 110 directions

Viscous Stress

Euler Stress

(7)

Macrodynamic behavior from Chapman-Enskog analysis

[ ( , ) ( , )]

( , ) ( , )

n n i iEQ n

i i i i i i

i i i

n t n t

n t t t n t

t

+ D + D = - -

å

^{r c} ^c

å

^{r c}

å

^r ^r ^c

1 2

1 1 2

; ; ;

eq

i i i

n = n +en r = er t = et t = e t

(

²

)

2 2

( , ) :

2

i i i s

EQ c i

i

s s

n a c

c c

r r

r = ^éêr + ^× + ^- ^ùú

ê ú

ë û

uu c c 1 u u c

Equilibrium distribution is chosen to give correct Euler stresses (Same low-order moments as Maxwell-Boltzmann distribution)

Define macroscopic length and time scales:

(8)

Expand space and time derivatives to 2nd order and collect terms

( )

( ) ( )

1

2 1

2 2 1

1

2

1 1 1

0 0 0 1

2

t

t s

t s s

c T

s

n

c n

c c n

c

r r

r r r

r r r s

t s r t

¶ + Ñ = =

¶ + Ñ + = =

¶ + + Ñ = =

é ù

= ë Ñ + Ñ û

u

u uu

uu uI

u u

To first order:

To 2nd order:

( ) ( ) ( ( ) ( ) )

2

2 2

1 1 1 1 1

0

/ 2

t

T

t t cs cs

r

r r r s

¶ =

¶ u + D Ñ Ñ u + Ñ u = Ñ ×

Incompressible on t₂ scale

"Lattice viscosity“-eliminates grid diffusion

(9)

Lattice-Boltzmann approximates Navier-Stokes on “large” scales

( )

2

2 2 2

2

0

; 1 ; 2 1

3

t

s s s

p

p c c x c t

t

r r

r r h

r h t r

¶ + Ñ × =

é ù

¶ + ×Ñ = -Ñ + Ñë + ÑÑ× û

= = D = - D

D u

u u u u u

Combining results from different time scales:

Navier-Stokes fluid dynamics in low velocity limit Leading order errors are M ² and Dx².

0.3 M <

(10)

Moving boundary condition by additional mass transfer-continuously varying velocity

Mass transfer prevents artificial pressure gradients Boundary conditions conserve global fluid mass

Momentum transferred into particle forces and torques

Stationary Boundary

Moving Boundary

n i

D µ ×u c

(11)

Lubrication forces important in dense suspensions; dominant in shear flows

Impractical to resolve flow in gap by any multi-

particle method: grid based, multipole, or boundary element.

Add lubrication forces pair by pair

Single patch point ~ 0.5Dx Similar results for other

components of F & T 2 additional patch points (independent of a)

0 10 20 30 40 50

0 0.1 0.2 0.3

h/R F/(6pmRU)

● LBE + lubrication

● LBE a ~ 5 Dx

(12)

Settling of a cluster of particles shows strong inertial effects even for Re ~ 1.

Cluster of 100-1000 particles

a) Re_c < 1: Cluster maintains shape Gradually sheds particles b) Re_c > 1: Forms ring structure

No shedding of particles

Breaks into smaller rings (Nicolai) Re_c 2rU R^c ^c

= m

Batchelor

& Nitsche

(13)

Computational details

1812 particles: diameter 5.4 Dx: R_c ~ 15a :  ~ 0.55: Re_c ~ 5 Periodic unit cell: 1024 x 400 x 400

~160 million grid points; 100,000 steps

16 P4 Xeons connected by Gigabit ethernet: 32 cpu’s

32 MSUPS aggregate performance: Run time ~150 hours New cluster: 192 dual-core P4’s with Gigabit ethernet

Observed good scaling up to 96 processors (~300 MSUPS) But still only limited inter-switch bandwidth (20Gbits/sec) Good scaling requires high performance switch

Extreme Networks x450a-48t ($6500)

(14)

• Sample size 15.2 ´ 9.9 cm

• Initial mean aperture mm

• dissolved until at Pe = 54 and Pe = 216

• high resolution data on fracture topography

Dissolution in a rough fracture. Modeling experiments by Detwiler et al., (GRL 2003)

water

KDP (potassium-

-dihydrogen phosphate) rough glass surface

0 0.126 h =

2 0

h = h

(15)

Velocity field calculated from implicit LBE

3D Stokes equations

• Sub-grid scale boundary conditions

• Steady-state solution determined directly, using conjugate gradients

more than 2 orders of magnitude faster than standard LBE

2

0 p Ñ × =

ìí

hÑ = Ñ î

v v

(Verberg, Ladd, 2000)

1.0 1.0 1.0 0.9

0.1

0.05 0 0

0.4

(16)

Random walk improvements

Classical random walk:

~ 10³particles per cell needed for accurate calculation of

' 1 m

m ¹ ( )0

J c

( _s 0) J = k c -c

Variable mass random walk:

• Tracking one particle at a time

• Works for linear kinetics only

(17)

Aperture growth at Pe = 54

experiment simulation

7 h0

0

2 0

h = h

• Channels form, grow, and compete for the flow

• Only a few channels survive at the end

• Strongly non-linear process

(18)

Key problem in simulating polymer solutions is the very long time scales.

Characteristic polymer relaxation time

For 100 unit chain, 10² steps per monomer diffusion time

~10⁶ steps per Zimm time

Need a short cycle time (< 10^-3s) to permit useful simulations of long-chains.

Brownian dynamics restricted to chains < 100 monomers since cycle time is proportional to

Use point particles to obtain a polymer simulation method Inertial equivalent of Brownian dynamics.

( )

³ ^1.8 ²

~ / ; /

Z R bG M N M M b DM

t t = t t =

N3

(19)

Brownian motion can be added to LBE via fluctuations in fluid stress (controlled)

Add Gaussian white noise at each node

So that the fluctuation dissipation relation is satisfied

( , ) ( , ) [ ( , ) ^EQ( , )]/ _i^f

i i i i i

n r c+ D + D =t t t n r t - n r t - n r t t + n

( )

² ¹⁸ ^, ^,

0

2 ; ⁱ

f f f

xy B xy i i x i y

i

k T n c c

s m s ⁼

=

= =

å

Velocity correlation function of a suspended particle agrees

quantitatively with dissipative decay of velocity and with Boussinesq equation

(20)

Collision operators for MRT, M10, BGK

0 1 2

3 2

4 2 2 2

5 2 2

6 7 8 9

2

10 2

11 2

12 2 2

13 2 2

14 2 2

15 4 2

16 2

17

1

1 2

______________________

(3 5) (3 5) (3 5)

( )

3 6 1

(2

x y z

y z

y z z x x y

x y z

y z x

z x y

x y z

ee c

e c

e c c c

e c c

e c c e c c e c c

e c c

e c c c

e c c

e c

==

= -

= - -

= -

==

=

= -

= - +

= - ² ² ²

2 2 2

18

3)(2 )

(2 3)( _y ^x _z )^y ^z

c c c

e c c c

- -

= - -

, 0

2

0 1 3 4 9

10 18 10 18 10 18

,

(1 )

; ; 0; ;

0; 0; ~ 2

0; ~ 2

( , ) ( , )

Nb

k k j j

j

eq neq e f

k k k k k k

eq eq eq

s

e e e

f f f

eq e f

k i

i i i k

k

m e n

m m m m m

m m m c

m m m

m m m T

n t t t w m t e

l

r r r r

h h

=

- -

- - -

=

¢ = + + + +

= = = +

= =

+ D + D = ¢

×

å

u uu

f uf fu

r c r

0 e

Nb

k= k

å

_e

(21)

Collision operators and hydrodynamic size

l_k = 0 for k = 0, 1, 2, 3: conservation laws

MRT: six independent, non-zero l_k (by symmetry)

Adjust location of hydrodynamic boundary via l_10-18. M10: three l_k; l_k = -1 for k > 9.

BGK: one l_k; all l_k equal (for k > 3).

2.67 1.73

0.43 10

2.69 2.09

1.04 5

2.70 2.45

1.90 2

2.72 2.67

2.58 1

2.71 2.77

2.69 0.7

2.71 2.83

2.75 0.6

2.72 2.90

2.73 0.55

2.73 2.94

2.77 0.53

MRT M10

t BGK

• MRT: t-independent radius (a₀ = 2.7).

• Decreased computational time, since large viscosity now accessible

• Insignificant differences in speed 1250 ticks/site (P4)

(22)

Fluctuations

Fluctuations in stress (Landau):

l₇ corrects for discrete time FDT

Improved agreement with FDT by including fluctuations in m_10-18 (Adhikari et al., 2004).

( )

^m⁷^f ² ⁼ ^s ^yz² ⁼ ²^T^hl⁷²

2

10 18

~ 0.6 0.8; stress fluctuations only

1; including fluctuations in

x

j T

j m

T r

r ^-

-

=

(23)

Point forces couple polymer and fluid (Ahlers and Duenweg ~2000)

For a bead-rod or bead-spring chain + fluctuating LBE fluid:

2) Calculate velocity field at each bead by interpolation

3) Calculate force on bead based on velocity relative to fluid 4) Redistribute force to LBE nodes

5) Add fluctuating force to beads to balance frictional losses Single particle correlation matrix

Long-range correlations in random force come from fluid dynamics of fluctuating LBE model.

Studied dynamical scaling laws in long chains (N ~ 10³) but for relatively short times.

(24)

Hydrodynamic Interactions

between point particles

_F

r/∆x

∆x

U

On lattice Off lattice

(25)

Self-diffusion of an isolated chain

CPU TIME ~

Computational effort can be greatly reduced for longer chains.

Fixed

Independent of N.

1.8 3

~

~ /

Z

V N b

t N b h T

3.6 6 / N b T

g ~ 5

R Dx

b

(26)

Weak Confinement Strong Confinement

Self-diffusion of an confined polymer

[1]

[2]

[1]Brochard F, deGennes P G, J. Chem. Phys., 67,52

[2]Jendrejack et al.,

J. Chem. Phys., 119,1165

R_g H

(27)

Confined polymers in flow

N = 15 H ~ 8R_g

y/b

2

4u R0 _g

Pe = DH

u₀

(28)

Particle methods fall into two categories

Forces: MD Conservative

DPD Conservative, Dissipative, Fluctuating SPH Conservative, pressure (from EOS) Computationally intensive neighbor search ~ 1000 FLOP Collisions: DSMC Boltzmann

LG Discrete RCLG Rotational

Local collision process is faster but only applicable to gases.

Spatial resolution limited by cell size

(29)

Particle methods are not competitive with CFD or LBE for hydrodynamic problems

Statistics: u_max< 0.3 c_s to maintain incompressibility 10% accuracy requires ~ 1000 particles Maximum resolution ~ 10 particles

Computational effort 10⁴-10⁵ larger than CFD Time averaging means reducing u_max

Time scales:

Cannot enforce proper time scale separation unless (of the order of 1000 in colloids) SPH and DSMC used for large-scale, high-speed flows

~ ~

D

H B

Sc a

D k T a

t h r h r

t h s

= = ~ mk T₂^B

h s

æ ö

ç ÷

è ø

a _?s

(30)

Even DPD does not work well for HI

Dissipative forces can increase But needs very large friction, and density,

Depletion forces perturb

thermodynamics and short- range structure

No hydrodynamics at small scales

(Whittle & Dickinson JCIS 2001)

5 2 3 1

~ 1.6 10 _c

Sc ´ ^- g% %n ar^-

~ 100 g_%

~ 10 n

(31)

Some advantages of LBE

External boundaries: arbitrary shape, no added cost Simplicity of random forces; potentially very fast

Simplicity (<5000 lines) and speed (10¹² grid points/day)

Superior accuracy for relative motion between solid and fluid Permeability of random arrays

of spheres (N = 16)

LBE: solid circles (a = 2Dx) Multipole: solid line (M=200) Brinkman: dotted line

(32)

Closing thoughts

Discrete kinetic theory (LGA/LBE) developed from intuitive, physically based, models: HPP-FHP-HS

Led to numerically important constraints being built in Exact conservation laws

Isotropic momentum diffusion (weighted diagonals) Dispersion free

Models developed from physically motivated guesses: e.g.

Moving boundary condition from Monte Carlo

Past 10 years LBE has become increasingly mathematical

Improved accuracy via unphysical equilibrium distribution Improved numerics: adaptive grids, elliptic solvers, etc.

But I believe there are still opportunities for physical insight.