VOLUME 86, NUMBER 16 PHYSICAL REVIEW LETTERS 16 APRIL 2001
Lattice Boltzmann Model for Anisotropic Liquid-Solid Phase Transition
W. Miller*
Institut für Kristallzüchtung (IKZ), Max-Born-Strasse 2, 12489 Berlin, Germany
S. Succi and D. Mansutti
Istituto per le Applicazioni del Calcolo (C.N.R.), Viale del Policlinico 137, Rome 00161, Italy
(Received 18 July 2000)
We develop a simple reaction model for the liquid-solid phase transition in the context of the lattice
Boltzmann method with enhanced collisions. Calculations for a two-dimensional test problem of Ga
melting and for a two-dimensional anisotropic growth of dendrites are presented and commented on.
DOI: 10.1103/PhysRevLett.86.3578 PACS numbers: 64.70.Dv, 44.25.+f, 81.10.Aj
During the last decade lattice gas automata and lat-
tice Boltzmann (LB) methods evolved as powerful tools
for the simulation of complex fluid motion and pattern
formation. These methods have been applied to many
physical problems such as, e.g., phase separation of two
immiscible fluids [1]. Only very recently a thermal lat-
tice Boltzmann method has been adapted to the problem
of liquid-solid phase transition [2]. De Fabritiis et al. used
a thermal model with two types of quasiparticles, for liquid
and solid phases, respectively. In one dimension, this leads
to a 5 3 5 system of constraints for the reaction terms.
Besides being restricted to the one-dimensional case, this
model contained a number of empirical parameters and as-
sumptions which cast some doubts on its viability for more
realistic studies.
In this Letter, we simplify and extend at a time this
previous model by using only one type of quasiparticles
and a phase field f, which defines the solid and liquid
fraction (21 #f# 1, where 21 means totally solid and
11 totally liquid). The quasiparticles evolve in discrete
time t
on a face-centered hypercubic lattice in which ev-
ery node r
is surrounded by 24 neighbors i at a distance
j c
i
j
p
2. These 24 neighbors correspond to the full set
of permutations c
i
61, 61, 61, 61. The distribution
N
i
of particles with speed c
i
obeys the following discrete
Boltzmann equation [3,4]:
N
i
r
1 c
i
, t
1 1 N
i
r
, t
1 F
B
i
1 F
fs
i
1 F
L
i
1
X
j
A
ij
N
j
r
, t
2 N
eq
j
r
, t
. (1)
Here N
eq
j
is the local equilibrium, which is calculated
via an expansion in the lattice gas velocity u r
, t
1
r r
,t
P
i
c
i
N
i
r
, t
, where r r
, t
P
i
N
i
r
, t
is
the lattice gas density. A
ij
is the collision matrix which
controls the viscosity n
of the lattice fluid via the inverse
of its leading nonzero eigenvalue [4]. F
B
i
, F
fs
i
, F
L
i
represent the components of the buoyancy force, the drag
force, and the heat source, respectively. The buoyancy
force is defined via j
F
B
j bT , where b is related to
the expansion coefficient times gravity acceleration. The
drag force and heat source are discussed below. The tem-
perature T is treated as a scalar and identified as the
velocity component along the fourth dimension of the
FCHC lattice [5,6]. We use the collision matrix given by
Massaioli et al., which allows nonunit Prandtl numbers
[7,8], extended by a reservoir of rest particles [9]. With a
suitable choice of the local equilibrium, the hydrodynamic
limit of Eq. (1) is known to reproduce anomaly free
Navier-Stokes equations for fluid flow plus the standard
transport equation for the temperature field within the
Boussinesq approximation.
The phase-field f evolves according to a rate equation:
f r
, t
f r
, t
2 1 1 R r
, t
. (2)
Since melting/solidification is an activated process, we
postulate the following chemical expression for the reac-
tion term R r
, t
:
R r
, t
f
1
K
1
r
, t
1 2f r
, t
2 1
2 f
2
K
2
r
, t
1 1f r
, t
2 1 , (3)
where f
1
, f
2
are frequency factors, basically the inverse
time scale for solidification/melting, respectively, and K
1
,
K
2
are switch functions controlling the onset of melting/
solidification around the critical temperature T
c
. The spe-
cific form is [T
0
r
, t
T r
, t
2 T
c
]:
K
6
r
, t
1
2
1 6 tanhT
0
r
, t
7 T
s
T
w
, (4)
where T
w
controls the energy range of the transition and
T
s
defines two distinct activation energies for melting and
solidification so as to reduce thermal fluctuations at the
interface. In order to allow the phase transition only at
the interface we modify the switch functions in the follow-
ing way:
K
6
r
, t
! K
6
r
, t
3
1
8
X
Ω
1
2
1 6f r
1 c
i
, t
æ
2
. (5)
Assuming f
1
f
2
f , T
s
0, and T
0
T
w
, 1,a
first order scale analysis leads to the partial differential
3578 0031-9007 01 86(16) 3578(4)$15.00 © 2001 The American Physical Society