7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Numerical simulation of 3D pool boiling using CahnHilliard equation
K. Tsujimoto*, A. Nakamura*, T. Shakouchi*. and T. Ando*
Division of Mechanical Engineering ,Graduate School of Engineering, Mie University
1577 Kurimamachiycho, Tsu, 5148507, Japan
tujimoto~mach.mieu.ac~jp
Keywords: CahnHilliard equation, Pool boiling, Diffuse interface method
Abstract
We propose the attractive numerical scheme in which the CahnHilliard (CH) equation is used to track the interface;
the temperature recovery method is introduced to represent the phase change.In the proposed scheme, phase change
can be easily taken into consideration by using CH equation involving the additional source term. As the quantitative
validation of proposed scheme, we compare the calculation result with the analytic solution for onedimensional
evaporation problem. Moreover, as a more practical problem, the results of the two and threedimensional pool
boiling are demonstrated. Although there is a numerical difficulty due to the large density difference between liquid
and vapor, we demonstrate that the stable computation is achieved using the proposed scheme.
Introduction
Numerical procedure
Governing equations and numerical scheme
As assuming incompressible flow for each phase, the
governing equations are as follows:
In various industrial applications, boiling flows gener
ally appear, however, it is difficult to accurately predict
using numerical simulations. The major reasons for nu
merical difficulty are as follows: the existence of large
density ratio between liquid and vapor and the strong
surface tension, lead to severe numerical instabilities. In
addition, since the heat transfer due to the phase change
induces the generation of mass flux at the phase bound
ary, it is difficult to strictly establish the conservation
equations. For multiphase flow computations includ
ing the phase change, the interface between the different
phases is tracked using several numerical schemes, such
as the volume of fluid (VOF) method[1,2], front tacking
method[3], level set method[4,5] and MARS method[6].
Recently as the interface tracking scheme, the phase
field method including the CahnHilliard equation[7] is
too much attention paid in the multiphase flow compu
tations[8]. However, to our knowledge, there is no ap
plication to the boiling flow using the abovementioned
CahnHilliard equation.
In the present study, we propose the numerical scheme
in which the CahnHilliard equation is used to track
the interface; the temperature recovery method is intro
duced to represent the phase change. Then we demon
strate the validity of proposed scheme based on both the
onedimensional analytical problem and two and three
dimensional pool boiling.
Af 1 1 /
1 1
Vp + V p( Vu +VuT)
P P
Du
DL
1 r _
+ I /c~i
P ]aJ x
xxc)nds
DT
DL ~
89 D
ag V(pu)
V kVT
P1 V p P1 P 1 1
V~ ~~ 1eA + (4)
anbr~(x
Vp x Vp
Vp
xR)nlds = 0V VI I
where, a is the velocity, p, p and p are the pressure, the
density and the viscosity, respectively. Af is the mass
flux due to the phase change. p, and pi are the density
of vapor and liquid, respectively. a, /c, Gand n denote the
coefficient of surface tension, the curvature of surface,
delta function and a unit vector normal to the interface,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
respectively. T, C), and k are the temperature, the spe
cific heat at constant pressure and the thermal conduc
tivity, respectively. Pe and Pc denote the Peclet Num
ber and the chemical potential. Note that the additional
terms including the mass flux, Af is introduced in both
the continuum equation and the CahnHilliard equation.
The plwsical quantities are defined with those of each
phase and 9, iLe.
CahnHilliard equation is derived.
+ V(pu)
iL
1 P
Pe
Abovementioned discussion is introduced in many lit
eratures (e.g.[9]). However, since the order parameter
relates to the density, Eq.(9) should be compatible with
the continuity equation. Thus the following is derived
by substituting the Eq.(6) into the continuity equation:
P = pi(1 + )/2 +p,(1 )/2
k = k1(1+p)/2+k,(1p)/2
1 = C1 (1+ p)/2+ C,,(1 p)/2
li Ps P i 1
Of V(u
Eq.(4) is derived for the phase change problem, so that
the both equations, Eq.(9) and(10) should be established,
Phase change model
So far, the phase change is estimated by two different
methods : One is that the heat flux difference on the in
terface is computed and relates to the phase change[1
5]. Another is the temperature recovery method[6]; in
a considering control volume including an interface, the
phase change is estimated from the excess temperature
for the saturation temperature. In the present simula
tion, in order to avoid the reconstruction of interface and
the recognition of interface locations, the temperature
recovery method is used.
The third term of ch.s. of Eq.(2) is expressed using the
continuous surface force (CSF) model[9], which means
the capillary force and is discretized using Eq.(5) known
as the surface stress (CSS) model[10]. We have already
confirmed that such as the surface stress model cor
rectly operate through the simulation of a single bubble
flow[11]. The special discretization is conducted with
the secondorder finite difference scheme on the stag
gered grid. Because of the largedensity ratio between
two phases, it is difficult to iteratively solve the pressure
Poisson equation. In the present simulation, BICGSTAB
method[12] is adopted as iterative method. It is well
known that the introduction of upwinding scheme for the
convective terms is useful to suppress the numerical in
stability. In the present simulation Adaptative QUICK
EST method[13] is employed to discretize the convec
tive terms.
CahnHilliard equation
In the phase field method[9], the interface profiles are
constructed as follows: the interface between the differ
ent phase is assumed to have finite thickness: the plws
ical quantities are continuously distributed inside the fi
nite thickness of interface: the profiles of interface are
determined based on the free energy of the considering
system.
The order parameter, 9 corresponding to the mass
density of, or the concentration of fluids is introduced,
then the Free energy, F(p) can be described as a func
tion of 9:
Flp] = ((po(.r))+ kVw(.r) ")r (7)
6F[p]
Pe(@) b(r f'(p(.r)) E V p(.r) (8)
where f'(p) is f'(p) = 3 p. E corresponds to the
interface thickness. The chemical potential, pe can be
found by minimizing the free energy, F. The inter
face diffusion fluxes are assumed to be proportional to
the gradient of chemical potential, finally the connective
pe ,AT
where L is the latent heat. When the temperature of a
cell including the interface exceeds the saturation tem
perature, T""t, the phase change takes place. According
to the excess temperature, AT, the mass flux in a cell
is estimated from eq.(12). The generated mass flux per
unit area, Af is derived as follows:
th _peAT
Discretization of convective terms
Since the discontinuity of plwsical quantities in differ
ent phases exists, in general, high accurate discretization
of convective terms is conducted in the multiphase flow
simulations. Since the central finite difference scheme
induces the numerical instability, some types of upwind
ing schemes are examined in the present simulation.
Adopting the divergence form for the discretiza
tion of convective terms, a divergence of numerical
flux,i)(uf)/Arr can be expressed using second order fi
nite difference scheme as followings;
[(a f)]i%
/ = aJr
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.4 0.5 0.6 y*
(a)Eyre'scheam
Figure 1: Schematic diagram of liquid and gas with
phase change(1D)
where f means transport quantities, [u f]+l means
the flux on the cell boundary. In the present simulations,
first order upwinding scheme, QUICK scheme[14],
Adaptative QUICKEST scheme[13] are applied to eval
uate the performance of numerical accuracy. Not shown
here, we confirm that the Adaptative QUICKEST is su
perior among the them.
Improvement of the solver for CahnHilliard equa
tion
In the computation of CahnHilliard equation, the
rapidly convergence of the computation plays a key role
to determine the interface profile. If the convergence is
not enough, the interface profile is not correctly deter
mined. In our previous research[11], we found that the
implicit scheme proposed by Eyre[15] well operates to
converge the CahnHilliard equation, compared to the
explicit scheme which requires a great number of itera
tion to converge. Eyre[15]'s scheme is as follows;
0.4 0.5 0.6
(b)Improved scheme iteration 20
Figure 2: Time evolution of interface profiles
~n 1 _* 2 (_a 4 1 +2V, l)At/Pe
V 2(5)3 3V ~)At/P e (17)
(V ( L)3 3V2 L)At/Pe (18)
The iteratively solved ~* in eq.(17) is substituted into
eq.(18).
In order to demonstrate the effectiveness of proposed
scheme, 1D phase change problem shown in Fig.1 is
examined. The computational length is 4i7r the grid
number is 128. As the initial profiles, the interface is
located at the center of computational area. To mimic
the phase change, a constant strength source is continu
ously given at the interface ( = 0 ). Figures 2 show the
time evolution of interface profiles. From Fig.2(a), since
Eyre's scheme is not able to correct the interface pro
files, the discontinuity at y* = 0.5 is remained through
di
Pe 20V4 n 1 .2 n1 +
1 (a i~!
3V (15)
However, if the phase change takes place at the inter
face, the unphysical discontinuity of order parameter
at the interface is formed. Thus the interface profiles
should quickly attain a quasisteady solution. In the
present papers, instead the Eyre's scheme, we propose
the new semiimplicit scheme in which the explicit terms
in Eyre's scheme are implicitly dealt with. The detail is
as follows;

S01
S06
S11
16
S2 1
S26
0 1
(a) temperature
0 6
2 1
2 6
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
T=TSat
Free boundary condition
Liquid
x
Vapor
T=TSat+AT Wall condition
Figure 3: Schematic diagram of liquid and vapor phase
for phase change problem
the computation. While the improved scheme is ca
pable of correcting the discontinuity of interface from
Fig.2(b), the interface is smoothly formed and continues
to move from left to right. These findings suggest that
the proposed scheme is capable of constructing the in
terface profiles under the phase change problem. In par
ticular, the new scheme demonstrates the possibility of
stable computation in boiling flows with the large den
sity ratio.
Results and discussions
One dimensional analytical problem
In order to evaluate the validation of the numerical
scheme, we consider the simple analytical model of
evaporation. As shown in Fig.3, the upper portion is
fulfilled with the liquid phase, the lower portion is the
vapor phase. The periodic boundary condition for the
horizontal directions and the free boundary at an upper
boundary and the wall and the constant temperature at a
lower boundary are assumed. The temperature of both
phases TI, T, is assumed as Ti < T,, thus, the evapo
ration occurs at the interface between two phases. The
computational length is 300, the grid number is 704.
The parameters used in the calculation are pi 125
[kg/m ], p,/pl 0.001, Op;1 4.88 [kJ/(kg K()],
Ci. 1.0 [kJ/(kg K)],k ki 19.6 [mW/(m K)],
k, kl, L 20.54 [kJ/kg], dt 10".
Figure 4(a) shows the time evolution of temperature
profiles. At an initial time (t = 0), because of no oc
currence of thermal diffusion, the temperature profile is
discontinuous at the interface position. As the computa
tion proceeds, the gas temperature around the interface
gradually decreases due to the evaporation, while the liq
uid temperature remains to be constant. Similarly Fig.
4(b) shows the time evolution of velocity profiles. At an
initial time (t 0 ), because of no occurrence of phase
change, the liquid velocity is zero. After that, the liq
(b) velocity
Figure 4: Distributions of temperature and velocity at
various time steps
uid is pushed out of computational domain by the phase
change. Because of the assumption of incompressible
flow, both the phase velocity should be constant, in
particular, the gas velocity remains to be zero. While
the liquid velocity takes place according to the phase
change. The above mentioned results demonstrated that
the discontinuous property of temperature and velocity
distribution in the analytical problem is well reproduced
using our proposed numerical scheme. In order to eval
uate the validation of the numerical scheme for the case
of high density ratio (p,/pl 0.001). The moving in
terface position is tracked. The computational length is
300, the grid number is 704. The parameters used in
the calculation are pi 125 [kg/m ], p,/pl 0.001,
Op; 4.88 [kJ/(kg K)], Ci. 1.0 [kJ/(kg K)],
ki 19.6 [mW/(m K()], k, kl, L 20.54 [kJ/kg],
di = 10". The analytical solution is demonstrated in
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
I
"
~
..?
(a> velocity
10000
Step number
erface disph
free boundary c<
Figure 5:
wall condition
Figure 6: Schematic illustration of boiling flow
(b) temperature
[1] as follows:
Figure 7: Instantaneous velocity vector plot and contour
of temperature around bubbles at t* = 60
wall at lower boundary are assumed in the computation.
The physical parameters, which are assumed to be water
and vapor, are the surface tension, a 0.07 [kg/s ],
gravity, g 9.8 [m/s2], liquid density, pi 1000
[kg/m3], the density ratio,p,/pl 0.001, the specific
heat of liquid, Op;1 4.18 [kJ/(kg K()], the heat ra
tio Cr /Op;1 0.5, the thermal conductivity of liquid,
kl 237.5 [mW/(m K()], the ratio of thermal con
ductivity, k,/kl 0.0402, the latent heat, L 2264
[kJ/kg]. At the initial instant, the water under the satu
ration temperature, T 9500, is fulfilled in the whole
region, and then the water is heating at the wall. The
saturation temperature is Tsat 10000. The wall tem
perature keeps constant, T, Tsat + AT, thus the heat
flux through the heating wall varies according to the rate
of generation of vapor bubbles. In the computations, the
temperature difference,AT T,,1 Tsat are set as
AT=10, 30, 80, 120[K(]. TosolvetheCahnHilliard
equation, the proposed scheme is used. As for whether
the phase change on the wall occurs or not, the recovery
y* (0) +2ai AT, t (19)
y* (t)
Figure 5 shows the evolution of the interface position.
In the figure, a line denotes the analytical solution, and
symbols denote the present computational results. As
you can see, the numerical result is good agreement with
the analytical solution using the proposed scheme.
Twodimensional pool boiling
Figure 6 shows the schematic illustration for the two
dimensional pool boiling. In order to simulate the pool
boiling, not shown here, the governing equations are
nondimensionalized with the physical quantities con
cerning a surface tension, a, gravity, g and a liquid den
sity, p Therefore, the reference length, velocity and
time are defined d, = J~~) as J; and
is Jd respectively. The computational domain
is taken to be a rectangular volume of 0 < x: < 4x7dds
and 0 < y < 4x7dds. The x:, y directions are hor
izontal and vertical directions, respectively. The pe
riodic boundary condition for the horizontal directions
and the free boundary at upper boundary and the heating
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(byt* = 28
(a~t* = 0
o o
O On
(e)t* = 40
O
1
(d)t* = 36
(f)t* = 42
O O
O O
O O
O
(g~t* = 45
(i~t* = 60
(1)t* = 80
(h)t* = 50
(k)t* = 70
(j~t* = 67
Figure 8: Time evolution of the generated bubbles at AT = 30
temperature method judges the phase change based on
the difference between the temperature at each grid point
and the saturation temperature, thus we do not need the
artificial model specifying the nucleation.
Figures 7 show the vector plots and the contour of
temperature. From these figures, the fluid around the
rising bubbles moves downward and the temperature of
bubbles is higher than that of fluid around bubbles, sug
gesting that the present computation is capable of repro
ducing the details of boiling phenomena.
Figures 8 show the time evolution of nucleation at
AT = 30. After the start of wall heating, the water
being excess of the saturation temperature become va
por, and after that vapor film is formed on the wall in
Fig.8(b). The instability of interface occurs in Fig.8(c),
and after that bubbles are formed on the wall. From
Fig.8(e), the generated bubbles are rising due to buoy
ancy force. The boiling phenomena occurs again at the
position where the bubble moves upward in Fig.8(g).
Figures 9 show the instantaneous picture of the flow
field under some cases of walltemperature difference,
AT. In the figures the dark colored regions denote the
bubbles. The origin of bubbles is continuously generated
on the heated wall, then the bubbles are growing and
move upward. As increasing AT, the flow regime seems
to change from nucleation boiling to film like boiling.
Figure 10 shows the wallheat flux versus AT. Be
cause of 2D simulations, the boiling curve is not agree
ment with real 3D flow quantitatively, but the fact that
the linearly increasing the heat flux at low AT region
(c)t* = 34
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
C
R.
rl~ R,
/1 n
c,
G,
cS
In,,,
c,
C3
>
(d) aT = 120"C
O
o O no r) R,
Ci,
nnQO
O
000
O O
1 nZAn
(c) AT = 8000
(b) AT = 3000
Figure 9: Instantaneous distribution of vapor bubbles at high density ratio,p,/pl = 0.001
and the peak point at high AT exist, is quantitatively
in accordance with the experimental findings. These fu
tures denote the reasonable behavior of bubble genera
tion, suggesting the validation of the proposed scheme.
Threedimensional pool boiling
Figure 11 shows the schematic illustration for the three
dimensional pool boiling. In order to simulate the pool
boiling, not shown here, the governing equations are
nondimensionalized with the physical quantities con
cerning a surface tension, a, gravity, g and a liquid den
sity, p Therefore, the reference length, velocity and
time are defined d, = J~~) as and
is Jd respectively. The computational domain
is taken to be a rectangular volume, H, x H, x Hz =
2x ~d, x 2i7dds x 2i7dds. The x:, z directions are
horizontal and the y is vertical directions. The grid num
ber is A, x N, x Nz = 64 x 64 x 64. The periodic
boundary condition for the horizontal directions and the
free boundary at upper boundary and the heating wall
at lower boundary are assumed in the computation. As
well as 2D computation, the physical parameters are as
sumed to be water and vapor. The surface tension is a
0.07 [kg/s ], gravity, g 9.8 [m/s2], liquid density,
pi 1000 [kg/m3], the density ratio,p,/pl 0.001,
the specific heat of liquid, Cpt1 4.18 [kJ/(kg K()],
the heat ratio Gr /Op;1 0.5, the thermal conductiv
ity of liquid, kl 237.5 [mW/(m K()], the ratio of
thermal conductivity, k,/kl 0.0402, the latent heat,
L 2264 [kJ/kg]. At the initial instant, the water un
der the saturation temperature, T 950 0, is fulfilled
in the whole region, and then the water is heating at
the wall. The saturation temperature is Tsat 10000.
The wall temperature keeps constant, Ta, Tsat + aT
The temperature difference,AT T,,,au Tsat is set as
AT =30 [K(]. To solve the CahnHilliard equation, our
proposed scheme is extend to a three dimensional code.
O~ O
(a) AT = 1000
alli 1 I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
101 102 d
Figure 10: Boiling curve
(a) velocity vector
Hx(Nx) \
YV
.u
.w
Figure 11: Schematic illustration of 3D Boiling flow
In the 3D computation, the artificial model specifying
the nucleation is not introduced as well as 2D simula
tion. Figures 12 show the vector plots and the contour
of temperature at a given cross section. In the figures,
the semitransparent colored regions correspond to the
vapour surfaces. From these figures, the fluid around the
rising bubbles and the temperature of bubbles demon
strate that the present computation is capable of repro
ducing the details of 3Dboiling phenomena.
Figures 13 show the time evolution of nucleation. Af
ter the start of wall heating, the water being excess of
the saturation temperature become vapor, and after that
vapor film is formed on the wall in Fig. 13(a) to (b).
The instability of interface occurs in Fig. 13(c) to (f).
Once the larger scale bubles are formed on the wall in
Fig.13(g). Afte the Fig.13(h), the bubbles having the
nearly same size, distribute intermittently in space and
time, and are rising due to buoyanct force.
Figure 14 shows the time evolution of heat flux, q"
[W/m2] at AT =30 in two and threedimesional com
putations. At initial time, the wall heat flux is the high
(b) temperature
Figure 12: Instantaneous distribution of vapor for p, /pi
0.001
est rate since the water in the whole region is lower than
the saturation temperature. After that, since the vapor
film is formed on the wall, the heat flux markedly de
creases. Then the vapor film breaks due to the buoyancy
instability and becomes the bubbles. The bubbles, which
are continuously formed, enhance the mixing over the
heating wall. As a consequence the heat flux is nearly
constant. This behavior of heat flux can be explained
based on the instantaneous flow pictures.
Conclusions
We propose the numerical scheme in which the Cahn
Hilliard equation is used to track the interface under
the phase change problem, and demonstrate the validity
and effectiveness of proposed scheme based on both the
onedimensional analytical problem and two and three
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a~t* = 0
(d)t* = 36
(e)t* = 40
(h)t* = 50
(f)t* = 42
(i~t* = 60
(1)t* = 80
(g~t* = 45
(j~t* = 67
(k)t* = 70
Figure 13: Time evolution of the generated bubbles at AT = 30
dimensional pool boiling. Conclusions are as follows;
1. It is demonstrated that the improvement of Eyre's
implicit scheme for solving the CahnHilliard equa
tion is capable of simulating the phase change prob
lem with large density ratio.
2. Compared to the analytical solution in the one
dimensional evaporation problem with the large
density ratio, it is demonstrated the proposed
scheme is capable of giving the reasonable results.
3. It is revealed that the phenomena of phase change
including vapor and water in two and three
dimensional nucleate pool boiling, is stably simu
lated using the proposed scheme and that according
to the heating rate, the pool boiling is qualitatively
reproduced under the condition with large density
ratio.
References
[1] S. Welch and J. Wilson, A Volume of Fluid based
Method for Fluid Flows with Phase Change, J.
Comp. Physics, Vol.160, pp.662682, 2000.
[2] H. Shirakawa, Y. Takata, T. Kuroki T. Ito and
S. Satonaka, Numerical Solution of Thermal and
Fluid Flow with Phase Change by VOF Method,
Trans. JSME Ser. B, Vol.66649, pp.24052412,
2003 (in Japanese).
[3] A. Esmaeeli and G. Tryggvason. Computations of
(byt* = 28
(c)t* = 34
'\r
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
two phase flow using CahnHilliard equation, The
2D 19th Sym. on Comutaional Fluid Dynamics B53,
S3D 2005 (in Japanese).
[12] H. A. van der Vorst. BiCGSTAB: A fast and
smoothly converging variant of BiCG for the solu
tion of nonsymmetric linear systems, SIAM J. Sci.
Stat. Comput., (2000), Vol.13, pp.631644, 1992.
[13] V.G. Ferreira, C.M. Oishi, F.A. Kurokawa, M.K.
Kaibara, J.A. Cuminato, A. Castelo, N. Man
giavacchi, M.F. Tome, and S. Mckee. A combi
o .. nation of implicit and adaptative upwind tools for
0 10000 20000 the numerical solution of incompressible free sur
Step number face flows, Commun. Numer. Meth. Eng., Vol.23,
Figure 14: Time evolution of surface heat flux q" at thep.4945207
heating wall [14] B.P. Leonard. A stable and accuracy convec
tion modeling procedure based on quadratic up
stream interpolation, Comp. Meth. Appl.Mech.
film boiling Part I: numerical method, Int. J. HeatEn.Vo19p.9817.
and Mass Transfer, Vol.47, pp.54515461, 2004. [15] D.J.Eyre. An Unconditionally stable onestep
[4] S. Tanguy, T.Minard and A. Berlemont. A Level shm o rdetssespern,19
Set Method for vaporizing twophase flows, J.
Comp. Physics, Vol.221, pp.837857, 2007.
[5] F. Gibou, L. Chen, D. Nguyen and S. Baner
jee. A level set based sharp interface method
for the multiphase incompressible NavierStokes
equations with phase change, J. Comp. Physics,
Vol.222, pp.536555, 2007.
[6] T. Kunugi, N Saito, Y. Fujita and A. Serizawa.
Direct Numerical Simulation of Pool and Forced
Convective Flow Boiling Phenomena, Heat Trans
fer 2002, Vol.3, pp.497502, 2002.
[7] J.W. Cahn, JE. Hilliard. Free energy of a nonuni
form system I, J. Chem. Phys., Vol.28, pp.258267,
1958.
[8] V.E. Badalassi, H.D. Ceniceros, S. Banerjee. Com
putation of multiphase systems with phase field
models, J. Comp. Physics, Vol.190, pp.371397,
2003.
[9] J.U. Brackbill, D.B. Lothe and C. Zemach. A con
tinuum method for modeling surface tension, J.
Comp. Physics, Vol.100, pp.335354, 1992.
[10] D. Lakehal, M. Meier, M. Fulgosi, Interface track
ing towards the direct simulation of heat and mass
transfer in multiphase flows, Int. J. Heat and Fluid
Flow, Vol.23, pp.242257, 2002.
[11] Y. Maesoba, K. Tsujimoto, Y Mizutani, T.
Shaokuchi and T. Ando, Numerical Simulation of
