Transcript Document

Midterm-Review
GLY-5835
Basic Fluid Dynamics
Momentum
• p = mv
• F = dp/dt = m dv/dt
Viscosity
•
•
•
•
Resistance to flow; momentum diffusion
Low viscosity: Air
High viscosity: Honey
Kinematic viscosity
Poiseuille Flow
• In a slit or pipe, the velocities at the walls are 0
(no-slip boundaries) and the velocity reaches its
maximum in the middle
• The velocity profile in a slit is parabolic and
given by:
G 2
2
u( x ) 
(a  x )
2
• G can be gravitational acceleration
or (linear) pressure gradient (Pin –
Pout)/L
u(x)
x=0
x=a
Entry Length Effects
Reynolds Number
• The Reynolds Number (Re) is a non-dimensional
number that reflects the balance between viscous and
inertial forces and hence relates to flow instability (i.e.,
the onset of turbulence)
• Re = v L/n
• L is a characteristic length in the system
• Dominance of viscous force leads to laminar flow (low
velocity, high viscosity, confined fluid)
• Dominance of inertial force leads to turbulent flow (high
velocity, low viscosity, unconfined fluid)
Re << 1 (Stokes Flow)
Eddies and Cylinder Wakes
Re = 30
Re = 40
Re = 47
Re = 55
Re = 67
Re = 100
Re = 41
Laplace Law
• There is a pressure difference between
the inside and outside of bubbles and
drops.
• The pressure is always higher on the
inside of a bubble or drop (concave side) –
just as in a balloon.
• The pressure difference depends on the
radius of curvature and the surface tension
for the fluid pair of interest: DP = g/r
Young-Laplace Law
• When solid surfaces are involved, in addition to
the fluid1/fluid2 interface – where the interaction
is given by the surface tension -- we have
interfaces between both fluids and the surface
• Often one of the fluids preferentially ‘wets’ the
surface
• This phenomenon is captured by the contact
angle
• Zero contact angle means perfect wetting;
• DP = g cos q/r
Basic Boltzmann Gas
Concepts
Kinetic Theory
• Complete set of position (r) and
momentum (p) coordinates for all
individual particles gives exact dynamical
state of system
• Together with classical mechanics, allows
exact prediction of future states
• However, this level of description is
essentially not possible
Statistical Mechanics
• Represent system by ensemble of many
copies
• Describe by distribution function
f(N)(rN,pN,t); N is number of particles
• Changes in f(N)(rN,pN,t) with time given by
Liouville equation (6N variables)
• Usually interested in low order distribution
functions (N = 1, 2)
First Order Distribution Function
• f(1)(r,p,t) gives probability of finding a particular
molecule with given position and momentum;
positions and momenta of remaining N-1
molecules unspecified
• No experiment can distinguish between
molecules, so the choice of which molecule
doesn’t matter
• f(1) adequate for describing all gas properties
that don’t depend on relative positions of
molecules (dilute gas with long mfp)
• External force Xi (small relative to
intermolecular forces)
• For each component i there is an fi(1)(r,pi,t)
such that probable number of type i
molecules with position coordinates in the
range r ±dr and momentum coordinates pi
±dpi is fi(1)(r,pi,t) dr dpi
Expected Evolution w/o Collisions
• If no collisions, then at time t + dt, the new
positions of molecules starting at r are [r +
(pi/mi)dt]
• New momenta are pi = pi +Xidt
• Thus,
fi
(1)
r, p i , t drdp i 
fi
(1)
  pi 

 r  dt , p i  X i dt , t  dt drdp i


m
i



Collisions
• But there are collisions that result in some phase
points starting at (r,pi) not arriving at (r + pi/mi dt,
pi+Xi dt) and some not starting at (r,pi) arriving
there too
• Set Gij(-)drdpidt equal to the number of molecules
that do not arrive in the expected portion of
phase space due to collisions with type j
particles during time dt
• Similarly, set Gij(+)drdpidt equal to the number of
molecules that start somewhere other than (r, pi)
and arrive in the portion of phase space due to
collisions with type j particles during time dt
fi
(1)
r, p i , t drdp i 
fi
(1)
fi
(1)
  pi 

 r  dt , p i  X i dt , t  dt drdp i


m
i





pi 
  r  dt , pi  X i dt , t  dt drdpi


m
i





 f i (1) r, pi , t drdpi   Gij   Gij  drdpi dt
j
Maxwell-Boltzmann Distribution
0.004
98 K
198 K
298 K
398 K
498 K
598 K
Fraction of Molecules
0.0035
0.003
0.0025
0.002
0.0015
0.001
0.0005
0
0
500
1000
Speed (m/s)
1500
2000
Cellular Automata and
Lattice Gas Cellular
Automata
CA Components
Transition Rule:
State of Neighbors/Self Last Time Step
New State
Wolfram’s Binary Rule Numbers
• 3 neighbor, 2 state CAs
7
6
5
4
3
2
1
0
27
26
25
24
23
22
21
20
Block diagram from
http://mathworld.wolfram.com/ElementaryCellularAutomaton.html
FHP Lattice Gas Cellular Automaton
Bit Value
128
64
32
16
8
4
2
1
A
0
0
0
0
0
0
0
1
B
0
0
0
0
0
0
1
0
C
0
0
0
0
0
1
0
0
D
0
0
0
0
1
0
0
0
E
0
0
0
1
0
0
0
0
F
0
0
1
0
0
0
0
0
S
0
1
0
0
0
0
0
0
R
1
0
0
0
0
0
0
0
S = Solid, R = Random
Boolean variables n = (n1, n2, …, n8)
8 bits  256 possibilities
• Fundamental basis is mass
and momentum conservation
• All particles have the same
mass and speed so that
momentum conservation
reduces to conservation of
the vector sum of the
velocities
C
B
D
A
E
F
Maximum of 1 particle per direction
Zero net momentum, head-on, 2and 3-particle collisions
Pre-collision
Possible post-collision configurations
Choose based on random bit R
Pre-collision
Post-collision
Unchangeable configurations
Pre-collision
Post-collision
No configurations other than the original conserve mass and momentum
All 5 and 6 particle collisions are similar; no configurations other than the original
conserve momentum
Collision ‘Look up’ Table
• Configurations 64 through 127 and
192 through 255 (01000000 through
01111111 and 11000000 through
11111111) are on solids
• Bounce back
C
B
D
A
E
In-State Bit Value
Out-State Bit Value
128
R
64
S
32
F
16
E
8
D
4
C
2
B
1
A
64
0
1
0
0
0
0
0
0
65
0
1
0
0
0
0
0
66
0
1
0
0
0
0
1
128
R
64
S
32
F
16
E
8
D
4
C
2
B
1
A
64
0
1
0
0
0
0
0
0
1
72
0
1
0
0
1
0
0
0
0
80
0
1
0
1
0
0
0
0
0
1
1
1
1
1
1
1
…
127
F
…
0
1
1
1
1
1
1
1
127
Collision ‘Look up’ Table
C
B
• Two-particle head-on collisions
D
A
E
In-State Bit Value
128
R
F
Out-State Bit Value
64 32
S F
16
E
8
D
4
C
2
B
1
A
128
R
64
S
32
F
16
E
8
D
4
C
2
B
1
A
9
(AD)
0
0
0
0
1
0
0
1
18
(BE)
0
0
0
1
0
0
1
0
137
(AD)
1
0
0
0
1
0
0
1
164
(CF)
1
0
1
0
0
1
0
0
Collision ‘Look up’ Table
C
B
• Three-particle head-on collisions
D
A
E
In-State Bit Value
F
Out-State Bit Value
128
R
64
S
32
F
16
E
8
D
4
C
2
B
1
A
42
(BDF)
0
0
1
0
1
0
1
0
0
21
(ACE)
0
0
0
1
0
1
0
1
0
1
170
(BDF)
1
0
1
0
1
0
1
0
1
0
149
(ACE)
1
0
0
1
0
1
0
1
128
R
64
S
32
F
16
E
8
D
4
C
2
B
1
A
21
(ACE)
0
0
0
1
0
1
0
1
42
(BDF)
0
0
1
0
1
0
1
149
(ACE)
1
0
0
1
0
1
170
(BDF)
1
0
1
0
1
0
Remapings
Noise
• Averaging
• Real thermodynamics
Lattice Boltzmann Model
Unit Vectors ea
e1
e2
e6
e3
e4
FHP
e5
D2Q9
Velocities
0,1
1,1
-1,1
-1,0
0,0
1,0
1,-1
-1,-1
0,-1
D2Q9
Lattice Boltzmann Model
Unit Vectors ea
f2
Direction-specific particle
densities fa
e1
e2
f1
f6
f3
e6
f7 (rest)
f4
f5
e3
e4
e5
Macroscopic flows
   fa
a
u
fe
a a
a

Density
Velocity
Fundamental LBM Equations
and their Implementation
Collide and Stream

f x, t   f x, t 
f x  e Dt , t  Dt   f x, t  
a
a
a
a
Streaming

eq
a
Collision (i.e., relaxation towards local equilibrium)
Collision and streaming steps must be separated if solid
boundaries are present (bounce back boundary is a separate
collision)
Computation of macroscopic
density and velocity
   fa
a
u
fe
a a
a

Computation of feq
2
2

e a  u 9 e a  u  3 u 
eq
f a x   wa  (x) 1  3 2 

4
2 
2 c
2c 
c

Building an LBM Model
Define velocity vectors
D2Q9
%define es
6
ex(0)= 0; ey(0)= 0
ex(1)= 1; ey(1)= 0
ex(2)= 0; ey(2)= 1
ex(3)=-1; ey(3)= 0
ex(4)= 0; ey(4)=-1
ex(5)= 1; ey(5)= 1
ex(6)=-1; ey(6)= 1
ex(7)=-1; ey(7)=-1
ex(8)= 1; ey(8)=-1
2
e2
e6
3
5
e3
0
e5
e1
1
e4
e8
e7
7
4
• In MATLAB, can write ex(0+1)=0, ex(1+1)=1, etc.
8
Problem definition
LY=10
LX=20
tau = 1
g=0.00001
%set solid nodes
is_solid_node=zeros(LY,LX)
for i=1:LX
is_solid_node(1,i)=1
is_solid_node(LY,i)=1
end
%
if ~is_interior_solid_node(j,i)
LY
LX
Initialize density and fs (assuming
zero velocity)
%define initial density and fs
rho=ones(LY,LX);
f(:,:,1) = (4./9. ).*rho;
f(:,:,2) = (1./9. ).*rho;
f(:,:,3) = (1./9. ).*rho;
f(:,:,4) = (1./9. ).*rho;
f(:,:,5) = (1./9. ).*rho;
f(:,:,6) = (1./36.).*rho;
f(:,:,7) = (1./36.).*rho;
f(:,:,8) = (1./36.).*rho;
f(:,:,9) = (1./36.).*rho;
2
2




e

u
9
e

u
3
u
eq
a
a
f a x   wa  (x ) 1  3 2 

4
2
c
2
c
2
c


// Computing macroscopic density, rho, and velocity, u=(ux,uy).
for( j=0; j<LY; j++)
{
for( i=0; i<LX; i++)
{
u_x[j][i] = 0.0;
u_y[j][i] = 0.0;
rho[j][i] = 0.0;
if( !is_solid_node[j][i])
{
a
for( a=0; a<9; a++)
a
{
rho[j][i] +=
f[j][i][a];
u_x[j][i] += ex[a]*f[j][i][a];
u_y[j][i] += ey[a]*f[j][i][a];
a a
}
a
u_x[j][i] /= rho[j][i];
u_y[j][i] /= rho[j][i];
}
}
}
  f
u
fe

Periodic Boundaries
• On boundary,
‘neighboring’ point is
on opposite boundary
ip = ( i<LX-1)?( i+1):( 0 );
in = ( i>0 )?( i-1):( LX-1);
jp = ( j<LY-1)?( j+1):( 0 );
jn = ( j>0 )?( j-1):( LY-1);
where
LHS = (COND)?(TRUE_RHS):(FALSE_RHS);
means
if( COND) { LHS=TRUE_RHS;} else{ LHS=FALSE_RHS;}
// Compute the equilibrium distribution function, feq.
f1=3.;
f2=9./2.;
f3=3./2.;

ea  u
eq
f a x  wa  (x) 1  3 2 
for( j=0; j<LY; j++)
c
{

for( i=0; i<LX; i++)
{
if( !is_solid_node[j][i])
{
rt0 = (4./9. )*rhoij;
rt1 = (1./9. )*rhoij;
rt2 = (1./36.)*rhoij;
ueqxij = uxij; //Can add forcing here; see MATLAB code
ueqyij = uyij;
uxsq = ueqxij * ueqxij;
uysq = ueqyij * ueqyij;
uxuy5 = ueqxij + ueqyij;
uxuy6 = -ueqxij + ueqyij;
uxuy7 = -ueqxij + -ueqyij;
uxuy8 = ueqxij + -ueqyij;
usq = uxsq + uysq;

2
9 e a  u  3 u 2 


4
2 c
2 c2 
feqij[0] = rt0*( 1.
- f3*usq);
feqij[1] = rt1*( 1. + f1*ueqxij + f2*uxsq
- f3*usq);
feqij[2] = rt1*( 1. + f1*ueqyij + f2*uysq
- f3*usq);
feqij[3] = rt1*( 1. - f1*ueqxij + f2*uxsq
- f3*usq);
feqij[4] = rt1*( 1. - f1*ueqyij + f2*uysq
- f3*usq);
feqij[5] = rt2*( 1. + f1*uxuy5 + f2*uxuy5*uxuy5 - f3*usq);
feqij[6] = rt2*( 1. + f1*uxuy6 + f2*uxuy6*uxuy6 - f3*usq);
feqij[7] = rt2*( 1. + f1*uxuy7 + f2*uxuy7*uxuy7 - f3*usq);
feqij[8] = rt2*( 1. + f1*uxuy8 + f2*uxuy8*uxuy8 - f3*usq);
}
}
}
2
2




e

u
e

u
9
3
u
eq
a
a
f a x   wa  (x) 1  3 2 

4
2 
2
2
c
c
c


Collide
// Collision step.
for( j=0; j<LY; j++)
for( i=0; i<LX; i++)
if( !is_solid_node[j][i])
for( a=0; a<9; a++)
{
fij[a] = fij[a] + ( fij[a] - feqij[a])/tau;
}
f a x, t  Dt   f a x, t  

Dt f a x, t   f

eq
a
x, t 
Bounceback Boundaries
• Bounceback boundaries are particularly simple
• Played a major role in making LBM popular
among modelers interested in simulating fluids
in domains characterized by complex
geometries such as those found in porous media
• Their beauty is that one simply needs to
designate a particular node as a solid obstacle
and no special programming treatment is
required
Two type of solids:
Bounceback
Boundaries
Collision (solid nodes):
Bounceback
// Standard bounceback.
temp = fij[1]; fij[1] = fij[3]; fij[3] = temp;
temp = fij[2]; fij[2] = fij[4]; fij[4] = temp;
temp = fij[5]; fij[5] = fij[7]; fij[7] = temp;
temp = fij[6]; fij[6] = fij[8]; fij[8] = temp;
// Streaming step.
for( j=0; j<LY; j++)
{
jn = (j>0 )?(j-1):(LY-1)
jp = (j<LY-1)?(j+1):(0 )
f x  ea Dt, t  Dt   f a x, t  Dt 
for( i=0; i<LX; i++)
a
{
if( !is_interior_solid_node[j][i])
{
in = (i>0 )?(i-1):(LX-1)
ip = (i<LX-1)?(i+1):(0 )
ftemp[j ][i ][0] = f[j][i][0];
ftemp[j ][ip][1] = f[j][i][1];
ftemp[jp][i ][2] = f[j][i][2];
ftemp[j ][in][3] = f[j][i][3];
ftemp[jn][i ][4] = f[j][i][4];
ftemp[jp][ip][5] = f[j][i][5];
ftemp[jp][in][6] = f[j][i][6];
ftemp[jn][in][7] = f[j][i][7];
ftemp[jn][ip][8] = f[j][i][8];
}
}
}
Incorporating Forcing (gravity)
• Gravitational acceleration is incorporated
in a velocity term
du
F  ma  m
dt
ueqxij = u_x(j,i)+tau*g; %add forcing
F
Du 

u eq  u  Du  u 
F

More LBM Boundary
Conditions
Key paper:
• Zou, Q. and X. He, 1997, On pressure and
velocity boundary conditions for the lattice
Boltzmann BGK model, Phys. Fluids 9,
1591-1598.
Choices
• Specify density (i.e., pressure via EOS)
– Velocity computed
– Dirichlet
• Specify velocity
– Density/pressure computed
– Von Neumann
• Lots of temporal/spatial flexibility
D2Q9 BCs
• For example:
Out
In
– f(4,7,8) = function of
f(1,2,3,5,6) and BC
type
Complex Geometries and
Higher Reynolds Numbers
Scaling to real world
• Choose a Reynolds number
• Choose a combination of slit width and viscosity to give
maximum velocity < ≈0.1 lu ts-1; smaller is better
• Maximum velocity 3/2 of average for Poiseuille flow in slit
• Best to use  = 1 for simple bounceback boundaries:
yields kinematic viscosity of 1/6 lu2 ts-1
• Solve for the gravitational acceleration needed to drive
the flow by rearranging Poiseuille equation
Entry Length Effects
•
•
•
r is the 'radius' -- the distance from
the center to the point of interest
a is the half-width (d/2)
u is the velocity at the point of
interest and uavg is the average
velocity
Tritton DJ (1988) Physical Fluid Dynamics, 2nd Ed. Oxford
University Press, Oxford New York
x /(d Re) is dimensionless distance down the pipe. When this distance is infinite, Poiseuille flow
is fully developed. If Re = 1, Poiseuille flow is well-developed in a short distance down the pipe:
x/d = 1 → x = d = 2a (just 2 half-widths down the pipe). As the Reynolds number increases, this
distance can become quite large. If x = 1 m with Re = 103 in a 10 cm pipe, x/(d Re) = 1 m /(10-1
m 103) = 10-2 and Poiseuille flow will not be fully developed even 1 m from the inlet.
Entry Length Effects
Tritton DJ (1988) Physical Fluid Dynamics, 2nd Ed. Oxford
University Press, Oxford New York
Velocity BCs
• Simple compressible fluid model
• Identical constant velocities at each end of the domain:
– mass of fluid will change with time
• Flow must be accompanied by a pressure gradient and hence the
pressure must be lower at the outlet
• The pressure and density are related through an ideal gas law of the
form P = /3 in this model
• Densities at the input and output must be different
• If the velocity boundaries on each end of the domain are equal,
mass will accumulate in the system because the mass flux of fluid in
(vin in) will exceed the out flux (vout out)
• Problem increases in severity as the pressure difference increases
• Incompressible model of Zou and He (1997), pressure boundaries,
or gravity-driven flow can be used to avoid this complication
Flow Past a Cylinder
• The drag force FD is defined in terms of
the drag coefficient CD as
FD  r u CD
2
0
• Gravitational Force
F  g LW  r
2

http://scienceworld.wolfram.com/physics/CylinderDrag.html
Tritton DJ (1988) Physical Fluid Dynamics, 2nd Ed. Oxford
University Press, Oxford New York
Re = 0.16
Re = 41
Taneda S (1956)
von Kármán Vortex Street, Re = 105
h / l  (1 /  ) sinh 1 (l )  0.281
h
l
Photo by S. Taneda. S. Taneda and the Society for Science on Form, Japan
Multiphase LB Models
Single
Component
Multiphase
Single Phase
(No Interaction)
Number of
Components
Attractive
Nature of
Interaction
Interaction
Strength
Multi- Component
Multiphase
Miscible
Fluids/Diffusion
Immiscible
Fluids
(No Interaction)
High
Repulsive
Inherent
Parallelism
Low
Equations of State
•
•
•
•
•
Perfect (Ideal) gas law
van der Waals gas law
Molar volumes
Temperature dependence and Critical Points
Liquid-vapor coexistence and the Maxwell
Construction
• Water and the non-quantitative nature of the
van der Waals equation
• Alternative presentation: P()
• Modern EOS for water
Perfect or ideal gas law (EOS)
PV  nRT 
•P is pressure (ATM)
•V is volume (L)
nRT
P
V
•n is number of mols
•R is gas constant (0.0821 L atm mol-1 K-1)
•T is temperature (K)
Molar Volume
nRT
P
V
Vm = V/n is the volume occupied by one mole of substance.
The gas laws can be re-written to eliminate the number of
moles n.
Perfect or ideal gas law
(EOS):
RT
P
Vm
Single Component D2Q9 LBM
RT
P
Vm

Pc  
3
2
s
So, if 1/Vm (mol L-1) is density, cs2 = RT
van der Waals EOS
Perfect or ideal gas law (EOS):
van der Waals EOS:
nRT
P
V
nRT
n
P
 a 
V  nb
V 
•‘a’ term due to attractive forces between molecules [atm L2 mol-2]
•‘b’ term due to finite volume of molecules [L mol-1]
2
Molar Volume
Vm is the volume occupied by one mole of substance. The gas
laws can be re-written to eliminate the number of moles n.
van der Waals gas law
(EOS):
nRT
n
P
 a 
V  nb  V 
2
 1 
RT
P
 a 
Vm  b
 Vm 
2
CO2: P-V Space
300
Ideal, 298K
250
Supercritical, 373K
P (atm)
200
Critical, 304K
150
Subcritical, 200K
100
50
0
0
0.1
0.2
Vm (L)
0.3
0.4
0.5
Critical T
van der Waals gas law (EOS):
Mathematica?
 1 
RT
P
 a 
Vm  b
 Vm 
2
2
 RT
 1  
d
 a   
Vm  b
 Vm  
dP
0
dVm
dVm
dP
T
2a
0
 3
2
dVm
Vm  b  Vm
Liquid-Vapor Coexistence:
P (atm)
CO2, P-V Space
70
69
68
67
66
65
64
63
62
61
60
Subcritical, 293K
Vapor Pressure
a = 3.592 L6 atm mol-2
b = 0.04267 L3 mol-1
Vapor Pressure
Vapor
Molar
Volume
Liquid
Molar
Volume
0
0.1
Vm (L)
0.2
0.3
Maxwell Construction:
P (atm)
CO2, P-V Space
70
69
68
67
66
65
64
63
62
61
60
Subcritical, 293K
Vapor Pressure
Maxwell Construction:
Area A = Area B
Vm , g

Vm ,l
PdVm  P(Vm, g  Vm,l )
Vapor Pressure
Liquid
Molar
Volume
0
B
Gives:
--vapor pressure
--densities of
coexisting liquid and
vapor
Vapor
Molar
Volume
A
0.1
Vm (L)
0.2
0.3
Flat interfaces only!
Water, 298K: P-V Space
1000
800
van der Waals
P (atm)
600
400
Ideal
200
0
-200
-400
-600
-800
-1000
0
0.5
1
Vm (L)
1.5
2
van der Waals Non-Quantitative:
Water, 298K: P-  Space
2000
1500
P (atm)
1000
500
Ideal Gas
Pvap
0
-500
van der Waals
-1000
Isothermal
Compressibility
of Water
-1500
-2000
0
200
400
600
800
 (kg/m3)
1000
1200
1400
Including Directional H-Bonds:
Water, 298K: P-  Space
2,000
1,500
P (atm)
1,000
500
Ideal Gas
0
Truskett et al
(1999)
-500
-1,000
Isothermal
Compressibility
of Water
-1,500
-2,000
0
200
400
600
800
 (kg/m3)
1000
1200
1400
Recap
• van der Waals equation gives a simple
qualitative explanation of phase separation
based on molecular attraction and finite
molecular size
• Maxwell construction gives the vapor pressure
and the densities of coexisting liquid and gas at
equilibrium, FOR FLAT INTERFACES.
• van der Waals equation fails to quantitatively
reproduce the EOS of water
LBM Fluid Cohesion

An attractive force F between nearest
neighbor fluid particles is induced as
follows:
6
F(x, t )  G (x, t ) (x  e a Dt )e a
a 1
•G is the interaction strength
• Y is the interaction potential:
 0 

     0 exp 



Other forms
possible
• Y0 and 0 are arbitrary constants
Shan, X. and H. Chen, 1994. PRE 49, 2941-2948.
D2Q9
8
F(x, t )  G (x, t ) wa (x  e a Dt , t )e a
a 1
wa is 2 for a = 1, 2, 3, 4 and is 1 for a = 5, 6, 7, 8
G <0 for attraction between particles
force is stronger when the density is higher
Incorporating Forces
F = ma = m du/dt
F
Du 

U = u + F1/ + F2/ + ...
LBM Non-ideal Equation of State
GRT
  2
P  RT 
2
Non-ideal Component
     0e

Shan and Chen, 1994. PRE 49, 2941-2948
Realistic EOS for water: Follows ideal gas law
at low density, compressibility of water at
high density and spinodal at high tension
No repulsive
potential in
LB model
Liquid/vapor coexistence at equilibrium (and flat
interface) determined by Maxwell construction
 0
Water EOS after Truskett et al., 1999. J. Chem. Phys. 111,
2647-2656.
Simplified EOS
GRT
2
  
P  RT 
2
     0e
1
RT 
3

 0
G
2
P     
3 6

Young-Laplace Eqn (1805)
1 1
DP     
 r1 r2 
∞
r
2
DP 
r
DP  0
Young-Laplace Eqn (1805)
1 1
DP     
 r1 r2 

DP 
r
Measuring Surface Tension
1
0.9
0.8
0.7
-2
DP (mu ts )
y = 14.332x - 0.0016
2
R = 0.995
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.01
0.02
0.03
0.04
1/r (lu-1)
0.05
0.06
0.07
Surface Adhesion (Adsorption)
in LBM
Wetting
http://www.hdm-stuttgart.de/projekte/printing-inks/b_sel42.jpg
Wetting
http://psii.kist.re.kr/Teams/psii/research/Con_4.jpg
Geometrically-controlled
Superhydrophobic surfaces
http://www.nature.com/nmat/journal/v1/n1/images/nmat715-f1.jpg
LBM Adhesive Force Formula
Fads (x, t )  GadsY(x, t ) wa s(x  e a Dt )e a
a
• s is a ‘switch’ that takes on value 1 if the site
at x + eaDt is a solid and is 0 otherwise
• We seem to have flexibility in the choice of
the pre-sum factor; the papers cited use  or
Y
Computation of 
• // Compute psi, Eq. (61).
• for( j=0; j<LY; j++)
•
for( i=0; i<LX; i++)
•
if( !is_solid_node[j][i])
•
{
•
psi[j][i] = 4.*exp( -200. / ( rho[j][i]));
•
}
Sforce
•
•
•
•
•
•
•
•
•
•
•
•
•
// Compute interaction force,
Eq. (66).
for( j=0; j<LY; j++)
{
jp = ( j<LY-1)?( j+1):( 0 );
jn = ( j>0 )?( j-1):( LY-1);
for( i=0; i<LX; i++)
{
ip = ( i<LX-1)?( i+1):( 0 );
in = ( i>0 )?( i-1):( LX-1);
if( !is_solid_node[j][i])
{
sum_x=0.;
sum_y=0.;
•
•
•
if( is_solid_node[j ][ip]) // neighbor 1
{ sum_x = sum_x + WM*ex[1];
sum_y = sum_y + WM*ey[1]; }
•
•
•
if( is_solid_node[jp][i ]) // neighbor 2
{ sum_x = sum_x + WM*ex[2];
sum_y = sum_y + WM*ey[2]; }
•
•
•
if( is_solid_node[j ][in]) // neighbor 3
{ sum_x = sum_x + WM*ex[3];
sum_y = sum_y + WM*ey[3]; }
Sforce
•
•
•
if( is_solid_node[jn][i ]) // neighbor 4
{ sum_x = sum_x + WM*ex[4];
sum_y = sum_y + WM*ey[4]; }
•
•
•
if( is_solid_node[jp][ip]) // neighbor 5
{ sum_x = sum_x + WD*ex[5];
sum_y = sum_y + WD*ey[5]; }
•
•
•
if( is_solid_node[jp][in]) // neighbor 6
{ sum_x = sum_x + WD*ex[6];
sum_y = sum_y + WD*ey[6]; }
•
•
•
•
if( is_solid_node[jn][in]) // neighbor 7
{ sum_x = sum_x + WD*ex[7];
sum_y = sum_y + WD*ey[7]; }
if( is_solid_node[jn][ip]) // neighbor 8
•
•
{ sum_x = sum_x + WD*ex[8];
sum_y = sum_y + WD*ey[8]; }
•
•
•
•
•
sforce_x[j][i] = -Gads * psi[j][i] * sum_x;
sforce_y[j][i] = -Gads * psi[j][i] * sum_y;
}
}
}
Contact Angles in SCMP LBM
• Cohesive force:
F(x, t )  G (x, t ) wa (x  e a Dt )e a
a
• Adhesive force:
Fads (x, t )  GadsY(x, t ) wa s(x  e a Dt )e a
a
Interplay between these forces will determine wetting
Contact Angles in SCMP LBM
F(x, t )  G (x, t ) wa (x  e a Dt )e a
a
Assume uniform liquid or vapor surroundings:
F  G 2  wa e a
a
Contact Angles in LBM
• Assume uniform surroundings:
F  G 2  wa e a
a
Liquid
F  G l
l
2
Vapor
w e
a a
a
F  G v
v
2
w e
a a
a
Contact Angles in LBM
• Assume uniform surroundings:
F  G  wa e a
a
Liquid surrounded by solid
Vapor surrounded by solid
Fads  Gads l  wa e a
Fads  Gads v  wa e a
l
a
v
a
Contact Angles in LBM
• Zero degree contact angle:
– Adhesive force equal to cohesive force for
liquid
Liquid
F  G l
l
2
Liquid surrounded by solid
w e
a a
Fads  Gads l  wa e a
l
a
 Gads  G l
a
Contact Angles in LBM
• 180 degree contact angle:
– Adhesive force on vapor equal to cohesive
force for vapor
Vapor
F  G v
v
2
Vapor surrounded by solid
w e
a a
Fads  Gads v  wa e a
v
a
 Gads  G v
a
Contact Angles in LBM
• 90 degree contact angle:
– Adhesive force on vapor equal to cohesive
force for ‘interface’  (= [l + v/2)
Interface
F  G 2  wa e a
Interface surrounded by solid
Fads  Gads  wa e a
a
 Gads  G
a