Presentation

Simulation of MHD Flows using the

Lattice Boltzmann Method

Kannan N. Premnath & Martin J. Pattison

MetaHeuristics LLC

Santa Barbara, CA 93105



Phase II SBIR – DOE Grant No. DE-FG02-03ER83715

Main Topics







Lattice Boltzmann Models for MHD



New Lattice Boltzmann Model for high

Hartmann number MHD



Simulation results in 2D and 3D



Summary and Conclusions

Magnetohydrodynamic (MHD) Equations



Fluid dynamical equations



   (  u)  0

t

u 1

   (uu)  p  NFLorentz    visc

t Re

Lorentz  J  B

Lorentz force F

Magnetic induction equation

B 1

   (uB - Bu)  2 B B  0

t Re m

Dimensionless numbers

u0 L  B02 L

Reynolds number Re  Stuart number N = magnetic force /inertial force

 u0

Hartmann number Ha  N Re = magnetic force /viscous force

u0 L

Magnetic Reynolds number Re m  ,



Lattice Boltzmann Model for MHD



I. Hydrodynamics

Macroscopic fields

b b

1

   f  u   f e p  c2 

 0  0 3 3D

moments of distribution function





Coincident lattices

Scalar distribution function for hydrodynamics

1 e ,   0,1,...,14

f ( x  e  t , t   t )  f ( x , t )   ( f  feq )

f   ,   0,1,..., 6

streaming collision



Equilibrium distribution functions



feq  feq   , u, B, mag 

functions of macro fields

Lattice Boltzmann Model for MHD (cont…)



II. Magnetic induction

Vector distribution function for magnetic induction

1

g ( x     t , t   t )  g ( x, t )   ( g  g )

eq



m 3D

streaming collision





Equilibrium distribution functions Coincident lattices

g  g  u, B 

eq eq

functions of macro fields

e ,   0,1,...,14



Macroscopic fields   ,   0,1,..., 6

bm

B   g

 0

 bm 

 ijk   e i g j   ui B j - Bi u j 

1 Nm 3, 2 D

 J k     B k   Nm  

mag c 2 m mag   0  4, 3D



moments of distribution function

Lattice Boltzmann Model for MHD (cont…)

Transport coefficients (Diffusivities)

c2  1 c2  1 

   f    t    m  t Prm  Magnetic Prandtl number

3 2 Nm  2 



Boundary conditions B  B a + Bi    



Special Bounce-back B i ( )  0 Insulating 

cases



Specular B i  n

reflection ( )  0 Conducting

n 

General Electromagnetic domain extending

Extrapolation method (proposed)

Case outside fluid flow domain



Fluid boundary feq

f

 

 feq u f m

f m

Post-collision f-1  2 f0  f1 Interior layer



Electromagnetic  Wall

f

Ghost layer

boundary

geq

 g ( Bm )

eq

m





Post-collision g  2 g  g

-1 0 1 



 B  B a + B i ( 0)

Advantages of LBM for MHD flows



Avoids time-consuming solution of Poisson-type pressure equation

All information obtained locally

Naturally amenable for implementation on parallel computers

Efficient calculation procedure for handling large problems



Well suited to MHD flows in complex geometries



Other advantages



Complete field formulation



Calculated fields solenoidal to machine round-off error



Current density as higher moment of the distribution function

(no finite differencing)

Comparison of LBM with Projection Method



Sample problem:



Flow through rectangular duct

Different grid sizes



12

10

For 20 timesteps Speed LB/PM 8

6

same machine 4

same problem: 2

0

20^3 30^3 40^3 50^3 60^3

Grid size

MHD Results in 2D



Orszag -Tang vortex



Time evolution









ˆ 1

J  jk  B ˆ

  k    u

 mag

Vorticity

Current density

Results comparable to other sources

MHD Results in 2D - II



Hartmann Flow









Ba Induced magnetic field profiles

Velocity profiles

y

ux ( y)

x

MHD Results in 2D -III

MHD Lid-driven Cavity





Fluid flow Domain:

128128

Electromagnetic

domain:

y 162162



Induced magnetic

fields set to zero on

the electromagnetic

boundary





x



Streamlines: Re = 100, Ha = 15.2, Rem = 100

MHD Results in 2D -III



MHD Lid-driven Cavity









u-velocity profiles v-velocity profiles

A new LB model for high Ha MHD flows - I

Adjustable magnetic Prandtl numbers (Prm) for liquid metals



High Hartmann number (Ha) flows require the resolution of

various thin viscous boundary or shear layers

1

Hartmann layers  m ~

1

Ha

1

Side layers  m2 ~

Ha

Ludford free shear layers from sharp bends



Standard LB MHD model restricted to uniform lattice grids

Standard LB MHD model uses a single relaxation time (SRT), which

restricts stability for a given resolution and variations in Prm

A new Multiple Relaxation Time (MRT) Interpolation Supplemented

Lattice Boltzmann Model (ISLBM) developed for non-uniform or stretched

grids with improved stability

A new LB model for high Ha MHD flows - II



Scalar distribution function for hydrodynamics



 1 

f ( x  e  t , t   t )  f ( x, t )   ( f   f  )   I    S 

eq MRT

 2  Model

Forcing term representing

streaming collision Lorentz force





  Components of the MRT matrix



Vector distribution function for magnetic induction

1

g ( x     t , t   t )  g ( x, t )   ( g  g )

eq

c t

m



Interpolation step x

Non-uniform Grid



f ( x, t  t )  F  f ({ xneighbors } + e  t , t   t ) 



g ( x, t  t )  G  g ({xneighbors } + e t , t  t )  Second order Interpolation of distribution

functions

A new LB model for high Ha MHD flows - III

Hartmann Flow



0.012



Ha = 700

0.01





0.008





0.006





0.004





0.002





0

-1.01 -1 -0.99 -0.98 -0.97







Non-uniform grid with

simple step-changes

Velocity profile (Ha = 700) in grid resolutions



Boundary layer stretching transformations (e.g. Roberts transformation) can be

used to further increase Ha

3D MHD Flows - I

Bza



z

y

Developed MHD duct flow

x









Induced magnetic

Velocity profile

field profile



Hartmann walls – perfectly insulating,

Side walls - perfectly insulating (Ha = 30)

3D MHD Flows - I



Developed MHD duct flow



Side wall jets









Velocity profile Induced magnetic

field profile

Hartmann walls – conducting,

Side walls - perfectly insulating (Ha = 30)

3D MHD Flows - II



3D Developing MHD Duct Flow – Sterl problem







B0 hydrodynamic MHD effect

Bza 

1  e x / x0





z

y





x









Streamwise sharp gradient in the Pressure Variation along streamwise

applied magnetic field direction (Ha = 44)

3D MHD Flows - II



3D Developing MHD Duct Flow – Sterl problem









Velocity profile at the exit plane Induced magnetic field at the

exit plane

Summary and Conclusions







Lattice Boltzmann simulations for for 2D and 3D MHD

performed

Simulations of MHD test problems in 2D and 3D show

qualitative and quantitative agreement

A new multiple relaxation time (MRT) interpolation

supplemented lattice Boltzmann model (ISLBM) for

simulating high Ha liquid metal MHD flows

Ongoing and Future Work



Code Version 1 Capabilities Code Release - end of June, 2005



MHD flows at intermediate Hartmann numbers

Multiphase flows

Heat transfer with non-uniform thermal conductivities

Complex geometries

Parallel processing using MPI

Pre-processor: Cart3D from NASA

Post-processor: FieldView



Code will be implemented on a smaller cluster at MetaHeuristics and a larger

cluster at National Center for Supercomputing Applications (NCSA)

Ongoing and Future Work





Code Version 2 Additional Capabilities

3D MHD flows at high Hartmann

numbers with multiple relaxation time

(MRT) model

3D complex geometries

Non-uniform grids

Turbulence modeling capability using

Smagorinsky type large eddy simulation

(LES) model



Code Release - end of October, 2005

Multiphase Flow Capabilities

Example Problem - I: Drop Collisions









Head-on collision resulting in Off-center collision resulting in

reflexive separation stretching separation

Example Problem - II: Drop subjected to Magnetic

Field

Example Problem - III: Rayleigh Instability and

Satellite Droplet Formation









Liquid Cylindrical Column Liquid Cylindrical Column

perturbed by shorter wavelength perturbed by longer wavelength

surface disturbance surface disturbance

Supplementary Slides

Pre-conditioning LBM for Accelerating

Convergence to Steady State

New Pre-conditioned LB MHD Model for

Acceleration to Steady-State

Evolution equations

1 

f ( x  e  t , t   t )  f ( x, t )   ( f  feq )  1 

1   e  u   F 1

 S  t S  feq

  2   cs2 f

1

g ( x    t , t   t )  g ( x, t )   ( g  g )

eq



m

Equilibrium distribution functions Pre-conditioning parameters:

 e  u 1   e  u  2 u  u 

 

feq  w  1   2    4  2 





cs  f  2cs

 2cs 

 0   f 1

 1    

(0)

0  m 1

g  W  B 

eq

 (0)  uB  Bu



 m   

Macroscopic Fields

bm

b c2

B   g

b

   f

1 1

 u   f e  Ft p f 

 0  0 2f 3  0



Transport coefficients

c2  1 c2  1

  f   f  t   m  m    t

3  2 Nm  2

Local Grid Refinement Technique for

LB MHD model

New Local Grid Refinement Schemes for LBM with

Forcing Terms and SRT/MRT Models - I

LBE with forcing term with single relaxation time (SRT) model

xc

1  1 

f ( x  e  t , t   t )  f ( x, t )   ( f  feq )  1   S  t

  2 



where forcing term is given by S  feq 

 e  u  F {tc, c} c

 cs2

Grid refinement factors

   c 1 c {tf, f} f

m  xc  tc  f   m  c  

1 1 q

  p

 xf  tf 2  2  f 1 f



Transformation Relations

 

f( c )  f( eq , f )  mp f( f )  f( eq , f ) 

1

2

1  p  S  tc xf



f( f )  f( eq ,c ) 

m



1 1 ( c )

p f  f( eq ,c )   1  p 1  S  tf

1

2

Here, tilde refers to post-collision value







f( c )  f( eq )  mq  f( f )  f( eq )  

1

 q  1 S tc Similar transformation relations can be

2 developed for the vector distribution

 f( eq )  q 1  f( c )  f( eq )    q 1  1 S  tf

1 1 function representing magnetic induction

f( f )

m 2

New Local Grid Refinement Schemes for LBM with

Forcing Terms and SRT/MRT Models - II

LBE with forcing term with multi relaxation time (MRT) model

 1 

f ( x  e  t , t   t )  f ( x, t )   ( f   f eq )   I    S  t

 2 



where forcing term is given by S  f eq  e  u  F xc

 cs2

ˆ

  T 1T  diag (s0 , s1, s2 ,..., s8 )

{tc, c} c

Grid refinement factors

 

m  xc  tc {tf, f} f

 xf  tf

1 1  1 1

  m    , i  0,1, 2,..., b

si f 2  sic 2 

xf

P    c 1  I   1  I 

 1

f





Q  c 1 f

New Local Grid Refinement Schemes for LBM with

Forcing Terms and SRT/MRT Models - III



LBE with forcing term with multi relaxation time (MRT) model (cont…)

Transformation Relations



 

f ( c )  f ( eq , f )  mP f ( f )  f ( eq , f ) 

1

2

 I  P  S tc



f ( f )  f ( eq ,c ) 

m



1 1 ( c )

P f  f ( eq,c )    I  P 1  S tf

1

2





f ( c )  f ( eq )  mQ  f ( f )  f ( eq )  

1

 Q  1 S tc

2



Q  f  f ( eq )    Q 1  1 S tf

1 1 ( c ) 1

f ( f )  f ( eq ) 

m 2



Here, tilde refers to post-collision value

Curved Boundary Treatment for MHD Flows

using LBM

New Curved Boundary Treatment - I

Scalar LBE with forcing term

1  1 

f ( x  e  t , t   t )  f ( x, t )  f ( x, t )  ( f  feq )  1   S  t

  2 

Here, tilde refers to post-collision value



Reconstructed distribution function

from the solid side wall

ff

f ( xb , t )  (1   ) f ( x f , t )   f(*) ( xb , t )

2  w  x  x f

 e  uw w

cs2

 1 b

     S  (1   ) S  ( x ,t )  t

where  2 f





 e u  e  u f  u f  u f 

2



f(*) ( xb , t )  w  ( x f , t ) 1  



 

bf



 2

cs 2cs 4

2cs 

2

x f  xw

   e  e

 1 1  x f  xb

ubf  1   u f  uw  ubf  u ff 

    1  1 cs 

1

  2  1   c

2  1  2   2 3

  2 

 



New Curved Boundary Treatment - II

Vector LBE

1

g ( x    t , t   t )  g ( x, t )  g ( x, t )  ( g  g )

eq



m

Here, tilde refers to post-collision value



Reconstructed distribution g ( x , t )  (1   ) g ( x , t )   g (*) ( x , t )  2W   (0)

  m 

function from the solid side

b m f b

  w

 W  m ( Bw  B f )  s( Bbf  B f ) 

 

   (0) 

g ( xb , t )  W  Bbf 

(*) bf

 (0)  uB  Bu

where



   

 1 1 

Bbf  1   B f  Bw 

   Bbf  B ff 

 

 1 1   bf   ff

(0) (0)



(0)  1   (0)  (0) 

bf

 

f



w

 1 2  1 

 1

  m   

2  1  2 m  2  2

m 



m

  m  (1   m )( m  1)   m  s  0



 m  (1   m )( m  1)   m  m 

s  is a free parameter



m  

s 

Sub Grid Scale (SGS) Turbulence Modeling

for MHD Flows using LBM

Sub Grid Scale (SGS) Modeling of MHD Turbulent

Flows For LES using LBM

Evolution equation of “coarse-grained” LBE

1

f ( x  e  t , t   t )  f ( x, t )   ( f  feq )



Total relaxation time    0 t



c2  1

Laminar kinematic viscosity  o   o    t

3 2



Smagorinsky SGS eddy

 Smag   Cs   S , S  Sij Sij , Cs ~ 0.09

2

viscosity



Effective Eddy  Ba

2



 eddy   Smag  exp    , Cm ~ 0.2

viscosity due to    Cm    Smag 

2



magnetic field  

Magnetic damping factor

(Shimomura, Phys. Fluids., 3: 3098 (1991))



c2  1

Total kinematic “viscosity”  total   0  eddy      t

3 2