7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Numerical simulation of 3D pool boiling using Cahn-Hilliard equation


K. Tsujimoto*, A. Nakamura*, T. Shakouchi*. and T. Ando*

Division of Mechanical Engineering ,Graduate School of Engineering, Mie University
1577 Kurimamachiy-cho, Tsu, 514-8507, Japan

tujimoto~mach.mie-u.ac~jp
Keywords: Cahn-Hilliard equation, Pool boiling, Diffuse interface method




Abstract

We propose the attractive numerical scheme in which the Cahn-Hilliard (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 one-dimensional
evaporation problem. Moreover, as a more practical problem, the results of the two- and three-dimensional 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 Cahn-Hilliard 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 above-mentioned
Cahn-Hilliard equation.
In the present study, we propose the numerical scheme
in which the Cahn-Hilliard 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
one-dimensional 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~ ~~ 1e-A -+- (4)


anbr~(x


Vp x Vp
|Vp|


xR)nlds = 0V |V|I 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 Cahn-Hilliard equation.
The plwsical quantities are defined with those of each
phase and 9, iLe.


Cahn-Hilliard equation is derived.


+ V(pu)
iL


1 P
Pe


Above-mentioned 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,(1-p)/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 second-order finite difference scheme on the stag-
gered grid. Because of the large-density 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.
Cahn-Hilliard 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))+ -k|Vw(.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 Cahn-Hilliard equa-
tion
In the computation of Cahn-Hilliard 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 Cahn-Hilliard 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 quasi-steady solution. In the
present papers, instead the Eyre's scheme, we propose
the new semi-implicit 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(]. TosolvetheCahn-Hilliard
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.
Two-dimensional 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 wall-temperature 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 wall-heat 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.
Three-dimensional 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 Cahn-Hilliard 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 semi-transparent 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 3D-boiling 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 three-dimesional 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
one-dimensional 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 Cahn-Hilliard 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.662-682, 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.66-649, pp.2405-2412,
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 Cahn-Hilliard equation, The
-2D 19th Sym. on Comutaional Fluid Dynamics B5-3,
S3D 2005 (in Japanese).

[12] H. A. van der Vorst. BiCGSTAB: A fast and
smoothly converging variant of Bi-CG for the solu-
tion of non-symmetric linear systems, SIAM J. Sci.
Stat. Comput., (2000), Vol.13, pp.631-644, 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.9-817.
and Mass Transfer, Vol.47, pp.5451-5461, 2004. [15] D.J.Eyre. An Unconditionally stable one-step

[4] S. Tanguy, T.Minard and A. Berlemont. A Level shm o rdetssespern,19
Set Method for vaporizing two-phase flows, J.
Comp. Physics, Vol.221, pp.837-857, 2007.

[5] F. Gibou, L. Chen, D. Nguyen and S. Baner-
jee. A level set based sharp interface method
for the multiphase incompressible Navier-Stokes
equations with phase change, J. Comp. Physics,
Vol.222, pp.536-555, 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.497-502, 2002.

[7] J.W. Cahn, JE. Hilliard. Free energy of a nonuni-
form system I, J. Chem. Phys., Vol.28, pp.258-267,
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.371-397,
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.335-354, 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.242-257, 2002.

[11] Y. Maesoba, K. Tsujimoto, Y Mizutani, T.
Shaokuchi and T. Ando, Numerical Simulation of