Journal of ThermalScience andTechnology

Vol.7, No.1, 2012

Initial Stage of Nucleate Boiling:Molecular Dynamics Investigation∗

Takahiro YAMAMOTO∗∗ and Mitsuhiro MATSUMOTO∗∗,∗∗∗∗∗ Department of Mechanical Engineering and Science, Kyoto University

Yoshida-Honmachi, Sakyo-ku, Kyoto 606-8501, Japan∗∗∗ Advanced Research Institute of Fluid Science and Engineering, Kyoto University

Kyoto daigaku-Katsura, Nishikyo-ku, Kyoto 615-8530, Japan

E-mail: matsumoto@kues.kyoto-u.ac.jp

AbstractWe investigated the initial stage of nucleate boiling on ideally smooth surface with amolecular dynamics simulation technique. Lennard-Jones (LJ) model liquid was con-fined in a rectangular simulation cell, contacting with a flat smooth solid wall. Thewall consists of fcc crystal of LJ-like particles. After the system was thermally equi-librated, the temperature of wall particles was raised to transfer thermal energy to theliquid. We examined two cases, the overall heating where the surface temperature iskept constant all over the area, and the partial (spot) heating where two regions of heat-ing and cooling are placed. In both cases, when the liquid in the vicinity of the heatingsurface obtains sufficient energy, it thermally expands and its pressure decreases, lead-ing to formation of bubble nuclei of atomic size. The inception time of nucleation wasfound to be affected by surface wettability as well as the surface temperature. When thesurface is hydrophobic and the heating area is small, size oscillation of the generatedbubbles was observed.

Key words : Nucleate Boiling, Nonequilibrium Molecular Dynamics Simulation,Wettability, Microbubbles, Young-Laplace Equation

1. Introduction

Nucleate boiling in liquids is initiated by bubble formation on heating wall surface. Pre-existing nuclei (PEN) theory(1) provides some understandings for its detailed mechanism;nucleate boiling in liquids starts from trapped gas (preexisting nuclei) in pits and scratcheson the surface. Most of experimental and theoretical investigations on boiling heat transferassume activated PEN of various scale(2) – (5).

A problem about the PEN theory is that it predicts extremely large superheats (the orderof hundreds degrees in some situations) for nucleation on smooth surfaces. Several experimen-tal works have been reported about boiling on smooth surface of nanometer-scale roughness,such as vapor-deposited metal nanofilms(6), silicon dioxide films(7), and electropolished metalsurfaces(8); in all cases, nucleate boiling was observed to occur at much smaller superheats.

The PEN theory relies on the pressure balance of a small bubble nucleus. The Young-Laplace (YL) equation for a spherical bubble in equilibrium,

Pb = P0 +2γR

(1)

tells us that, for a bubble nucleus of radius R to stably exist, the pressure inside the bubble(vapor or gas) Pb should be larger than the pressure of ambient liquid P0 by 2γ/R, where γ isthe surface tension. Here raised is a question on how nano-scale bubble nuclei can appear inthe case of nucleate boiling on a smooth surface. When the bubble size decreases, the pressuredifference ΔP = 2γ/R diverges as R−1; for example, ΔP = 120 atm for a bubble of R = 10 nmin the case of water at 373 K with γ = 0.059 N/m. This suggests that bubble nuclei of 10 nm

∗Received 1 Feb., 2012 (No. 12-0080)[DOI: 10.1299/jtst.7.334]

Copyright c© 2012 by JSME

334

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 1 Properties of nano-scale bubbles in Lennard-Jones liquid; size dependence ofsurface tension and surrounding liquid (ambient) pressure. Data are taken fromRef. (9). The Young-Laplace relation, Eq. (1), is shown as curves in the bottomfigure.

can appear under the normal pressure only when the inside vapor pressure exceeds 120 atm,and therefore that initiation of the nucleate boiling on ideal surface seems almost impossible.

There are several possible scenarios: (i) the YL equation holds only for macroscopicbubbles, and it fails for nano-scale phenomena, or (ii) we can still use Eq. (1), but the surfacetension γ depends on R for micro-scale bubbles; the paradox would be dissolved if γ ∝ R, forexample.

We investigated the mechanical stability of a nano-scale bubble with a molecular dynam-ics (MD) simulation technique, both for monatomic model liquid(9) and for water(10). Exam-ples of the results are shown in Fig. 1 for the monatomic liquid. The surface tension and thesaturated vapor pressure (not shown in the figure) were found to be almost independent of thebubble size. We therefore concluded that the YL equation still holds for nano-scale bubbles,which leads to a large negative pressure for the surrounding liquid at the equilibrium. Park etal.(11) reported a slight size dependence of γ (increase by up to 15%) for nano-scale bubbles,which does not dissolve the above difficulty, either.

In this paper, we try to elucidate what happens at the initial stage of nucleate boiling onideally smooth surfaces. For a heating process of model liquid, we employed an MD simula-tion technique similar to that used in the investigation of homogeneous bubble nucleation(12)

and nanobubble properties(9). The changes of local temperature, pressure, and density duringthe nucleation were monitored, which enables us to discuss the mechanism of nuclei genera-tion. Novak et al. adopted a similar MD technique to investigate nucleation from asperity orindentation of heater surface(13). Here we focus on nucleation on ideally smooth surface.

2. NomenclatureA : area of wall surface [m2]E : total energy [J]J : transition rate as 1/τi [s−1]

kB : Boltzmann constant [J/K]L : cell size [m]

m : particle mass [kg]N : number of fluid particles [–]P : pressure [Pa]

Pext : external pressure [Pa]Q : energy exchange rate [W]q : heat flux [W/m2]

335

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

R : bubble radius [m]r : interparticle distance [m]

T : temperature [K]V : volume of a bubble [m3]

ΔTs : superheat [K]ε : Lennard-Jones energy parameter [J]φll : potential function between fluid

(liquid) particlesφww : potential function between wall

particlesφwl : potential function between a wall

particle and a fluid particleγ : surface tension [J/m2]κ : specific heat ratio [–]ρ : density [kg/m3]σ : Lennard-Jones size parameter [m]τ : unit of time [s]τi : inception time [s]τo : period of bubble oscillation [s]

Subscript0 : ambientb : bubblec : cooling area

eq : equilibrium

h : heating areal : liquid or fluid

w : wall

3. Simulation Method

In the previous report,(14) we examined the nucleation process with a smooth potentialwall model as a heating surface. With this wall model, however, we cannot discuss the physicaldetails of the energy transfer between liquid and the wall because the wall has no kineticenergy. In this paper, we have developed a new model, namely a particle wall model, which isan fcc crystal of wall particles. By controlling the temperature (mean kinetic energy) of wallparticles, we can simulate heat transfer from the wall to liquid.

3.1. SystemA fixed number of particles are confined in a three dimensional rectangular cell of Lx ×

Ly × Lz, as shown in Fig. 2. Periodic boundary conditions are assumed in the horizontal direc-tions (x and y). A solid wall (described in Sec. 3.2) is situated at the bottom (z = 0), whereasthe top ceiling (z = Lz) is a simple reflective boundary. With appropriate initial conditionsand at a moderate temperature, vapor phase occupies the upper region, which enables us toexamine the nucleate boiling under an ambient pressure condition.

We assume the Lennard-Jones (LJ) 12-6 model potential for the particle interaction influid:

φll(r) = 4ε

[(σ

r

)12−

(σ

r

)6]

(2)

where r is the particle-particle distance. Two parameters ε and σ correspond to the potentialdepth and the particle diameter, respectively. Throughout this paper, we use ε, σ, and the

Fig. 2 Cross-sectional view of the simulation system for the case of (left) overallheating and (right) partial heating. The fluid (liquid and vapor) is confined in aconstant-volume rectangular cell.

336

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Table 1 Units of the simulation and values corresponding to liquid water. The symbolkB is the Boltzmann constant, 1.3807 × 10−23 J/K.

Units in simulation Units for waterlength σ 2.9 × 10−10 menergy ε 6.2 × 10−21 Jmass m 3.0 × 10−26 kgtemperature ε/kB 450 Ktime τ = σ

√m/ε 6.4 × 10−13 s

pressure ε/σ3 2.5 × 108 Pasurface tension ε/σ2 7.3 × 10−2 J/m2

heat flux ε/(σ2τ) 1.1 × 1011 W/m2

particle mass m as the units of energy, length, and mass, respectively. Although this is a modelmainly used for simple monatomic fluids, we can give a rough correspondence to water, bymatching the liquid properties (the temperature and the density at the triple point and thevapor-liquid critical point); summarized in Table 1 are the units and the values.

The equation of motion for each particle is integrated using the leapfrog algorithm withthe time step Δt = 0.001 τ. The interaction is truncated at rc = 3.5σ and no long-tail correc-tions are made.

The number of fluid particles is N = 64, 000. The cell size is Lx = 426√2σ, Ly =

21√

36√2σ and Lz = 120σ, where Lx and Ly are determined by the size of the bottom wall

consisting of an fcc crystal.

3.2. Heating method with particle wall modelIn our new wall model, the heating surface is composed of five layers of fcc (111) plane

with 8,820 particles. Position of particles in the bottom layer is fixed in order to prevent thewall from migration. The above four layers, the top of which contacts with the liquid, arecomposed of temperature-controlled particles. The following modified LJ interaction φww isassumed only between the nearest wall particles so that the wall should not collapse:

φww(r) =

⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

4εww

[(σ

r

)12−

(σ

r

)6]

: r ≤ rw

εww

⎡⎢⎢⎢⎢⎢⎢⎣14413

(7

26

) 76 rσ− 637

169

⎤⎥⎥⎥⎥⎥⎥⎦ : r > rw

(3)

We adopted a constant retracting force beyond the inflection point of the LJ potential, rw =

(26/7)1/6 σ, to prevent large lattice deformation and sublimation at higher temperatures. Weassume no specific material for the wall; instead, we have rather arbitrarily chosen the wall-wall energy parameter εww 10 times larger than ε in order to keep the perfect crystal structure.The same mass m and the same diameter σ are used for the wall particles so that the thermalresistance at the wall-liquid interface is as small as possible(15).

The interaction between fluid particles and wall particles is assumed to be also the LJtype:

φwl(r) = 4εwl

[(σ

r

)12−

(σ

r

)6]

(4)

We can control the surface wettability by varying the wall-liquid energy parameter εwl. Fromthe contact angle evaluation for a droplet on a similar model wall(16), the wall is roughlycategorized as

εwl

ε=

⎧⎪⎨⎪⎩ < 1 hydrophobic> 1 hydrophilic

(5)

In this study, we mainly examined three cases, εwl/ε = 0.5 (hydrophobic), 1.0 (neutral), and1.5 (hydrophilic) to investigate the effects of surface wettability on boiling behavior.

337

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

We compare two ways of heating:• Overall heating, or uniform heating: All particles in the top four layers of the wall are

temperature-controlled with a simple velocity scaling technique so that the mean temperatureis kept constant at Tw. This corresponds to an ordinary setup for macroscale boiling experi-ments. We examined the temperature range 1.1 ≤ Tw ≤ 1.7. Note that the vapor-liquid criticalpoint temperature is about 1.33 for the LJ fluids.• Partial heating, or spot heating: In the overall heating case, bubble nuclei are expected

to appear simultaneously at various places on the wall. As a simpler case, we limit the nu-cleation place in a small region. The wall is divided into two regions, the temperatures ofwhich are controlled separately. The circular area, indicated by red in Fig. 2 (right), is the“heating spot” with temperature Th. Its radius was chosen to be 10σ. The surrounding area isindependently temperature-controlled (cooled) at Tc, so that nucleation is possible only on theheating spot. We mainly examined the case with Th = 3.0 and Tc = 0.7. To prevent direct heatconduction between the two wall regions, no particle interaction between them is assumed;slight lattice deformations were observed due to this heterogeneity of the wall, but the crystalstructure with the interaction, Eq. (3), was sufficiently stable. This situation is rather artificialfor the sake of easier data analysis, but has some resemblance to spot heating with short-pulselaser on metal nanoparticles(17).

After equilibrating the system at the temperature Teq = 0.7, which is slightly above thetriple point temperature of LJ fluid, we start to heat up the liquid from the bottom by keepingthe wall at higher temperature. If the wall temperature is raised instantaneously to Tw (overallheating) or Th (partial heating), we found that strong and long-lived shocks (pulse waves)are generated and propagated in liquid. Thus we gradually increase the wall temperature in50,000 steps (50 τ) as a linear function of time.

4. Results: Partial Heating Case

A typical change for the partial heating case is shown in Fig. 3, which is sequentialsnapshots (sectional views); εwl/ε = 1.0 (neutral wettability), Th = 3.0, and Tc = 0.7. Abubble nucleus appeared at the very initial stage (∼ 50 τ or 32 ps); it grew up rapidly andreached a seemingly steady state.

4.1. Local distributionTo examine how the bubble nucleus appears on the wall, we look at the change of local

temperature and number density. Assuming a rotational symmetry around the z axis throughthe center of the heating area, we divide the system into rings of thickness (2/3)σ and height(2/3)σ, and take average of particle kinetic energy and number of particles as functions ofradial coordinate r and height z. Examples are shown in Fig. 4. At the inception of nucleateboiling (left), the liquid has high temperature (≥ 1.2) only for the region contacting with the

Fig. 3 Nucleate boiling on partial heating wall; sequential snapshots of sectional view.Example for a surface with neutral wettability (εwl/ε = 1.0) at the temperatureof the heating area Th = 3.0.

338

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 4 Local distribution of temperature and number density. The distribution wasobtained in the cylindrical coordinates with assumption of axial symmetry. Thesimulation conditions are the same as in Fig. 3.

heating spot. This high temperature region thermally expands to lead density decrease. Theinhomogeneity of temperature and density seems to be much localized in space, and thus weexpect that the system size little affects the initial dynamics of boiling, as discussed later inSec. 5.4. After a sufficiently long time (right), there exists a steady temperature gradient in theliquid. The local pressure distribution was also evaluated, but its large fluctuations concealedthe detail.

4.2. Void analysisIt is not straightforward to define and detect the bubble nuclei without ambiguity from

the simulation data. Here we adopt a simple criterion based on the instantaneous positions ofparticles; the simulation system is divided into cubic sub-cells of size 1.5σ, and the sub-cellis defined “empty” when it is not occupied by any fluid particles. A mutually connected groupof the empty sub-cells is defined as a void or a bubble. In this partial heating case, there existsonly one bubble at most.

4.3. Heat fluxThe heat flux q is relevant for boiling phenomena. Here we evaluate q with two methods;

( 1 ) The total energy El of the fluid is monitored as a function of time. Since the fluid

339

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 5 Comparison among energy exchange rates. The rate Qh ≡ qhAh for the exchangebetween the fluid and the heating region and Qc ≡ qcAc for the exchangebetween the fluid and the cooing region, where Ah and Ac are the area of eachregion, respectively, and Aw ≡ Ah + Ac is the total wall area. The flux ql isevaluated from the change of fluid total energy. Bezier spline curves are shownsince instant values of Q contain large fluctuations.

exchanges energy only with the bottom wall, the flux is estimated as

ql =El

Aw(6)

where Aw ≡ LxLy is the total area of the bottom wall.( 2 ) We also monitor the energy exchange rate Q between the fluid and the wall through

the scaling factor of the temperature control technique. Since the wall has two independentareas in the partial heating case, Qh for the heating area and Qc for the cooling area are storedseparately. The total heat flux is evaluated as

qw =Qh + Qc

Aw(7)

Both Qh and Qc generally contain large fluctuations, but Qc should be negative on average asthe heat is transferred from the warmed liquid to the cooling area.

Comparison is made in Fig. 5 for the case of εwl/ε = 1.0. Agreement between ql andqw is satisfactory in general, except for the initial transient stage with very large heat flux. Inthe following analysis, ql is mainly used for the heat flux q since fluctuations of qw are muchlarger.

4.4. DiscussionsFigure 6 summarizes the changes of bubble volume Vb, the pressure P, the temperature

T , and the heat flux q. The pressure and the temperature are spatially averaged values in theregion close to the heating area. A fluid region of 0 ≤ z ≤ 10 and r ≤ 10 is used to evaluateT . For P, a region of 5 ≤ z ≤ 10 and r ≤ 10 is adopted to avoid the “structured” liquidregion on the wall. The pressure was evaluated with the conventional virial theorem; it shouldcontain spatial anisotropy due to the generated nucleus and the temperature heterogeneity, butfluctuations are too large to analyze the detail.

In general, the bubble nucleus is generated more easily on hydrophilic surface [Fig. 6 (c)]than on hydrophobic one [(a)]. This is caused by the difference of heat flux q, which suggeststhat fluid particles are more attracted to hydrophilic wall due to the stronger interaction εwl andreceive more energy from the wall. Since we gradually raised the heating area temperature Th

in 50 τ, the heat flux takes a maximum value around t = 50 τ; the flux subsequently decreaseswith the increase of liquid temperature T .

In the cases of εwl/ε = 1.0 and 1.5, a discernible bubble nucleus appears when q takes amaximum value. Size oscillation of the bubble nucleus is interesting. A typical period of theoscillation τo is 100 τ [Fig. 6 (b)], which roughly agrees with a simple Minnaert model taking

340

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 6 Partial heating; comparison among surfaces of different wettability. Bubblevolume Vb, local pressure P and local temperature T near the heating area, andheat flux q entering the fluid.

no account of surface tension and viscosity,(18)

τo = 2πR

√ρ

3κPl 95 τ (8)

with the maximum bubble radius R 3

√3Vb

2π∼ 7σ (assuming a hemispherical bubble), the

liquid density ρ 0.7 mσ−3, the specific heat ratio κ = 5/3 (monatomic fluid), and the ambientliquid pressure Pl 0.03 εσ−3 as the saturated vapor pressure at T = 1.0.(9)

In contrast to our previous wall model(14), where we observed a large negative pressure−1.0 ∼ −0.5, the pressure P in the vicinity of heating area in this simulation fluctuates very

341

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 7 Pressure fluctuations; comparison among equilibrium MD simulations of bulkliquid with various densities at T = 0.90 and the boiling simulations.

much around zero. Still it often reaches negative values as large as −0.1, which apparentlyleads to spontaneous nucleation on the wall, which looks similar to homogeneous nucleationin bulk liquid(12), (20).

To look into the pressure fluctuations in some detail, the histogram of instant pressure isshown for each εwl/ε condition in Fig. 7, in comparison with those of equilibrium MD resultsat several number densities. The equilibrium MD simulations were executed for “bulk” (i.e.,with periodic boundaries for all three directions) LJ liquid of 1,000 particles at T = 0.9,which is roughly equal to the size and the condition for the pressure evaluation in the boilingsimulations. The pressure fluctuations for the εwl/ε = 1.5 (hydrophilic) and 1.0 (neutral) casesare almost the same as those at equilibrium with number density ρ = 0.75 at T = 0.9; the meanpressure is almost zero, but instant values reaches as negatively large as −0.2, which seemssufficient for liquid-vapor instability, or spinodal transition. When considering the fact thatρ = 0.75 is slightly smaller than the density ρ 0.755 of saturated liquid(19) at T = 0.9, theliquid near the heating region is in a stretched state, as expected.

4.5. Scenario of bubble nucleus generationBased on these results, we can speculate how a nucleus appears on the heating area. We

note that the YL equation requires the ambient pressure P0 to be about −0.05 for the bubblenucleus of radius R = 5 ∼ 10 to emerge at T = 1.0 (Fig. 1). These conditions are similar tothose in the partial heating simulation; the nucleus first appears at T = 0.9 ∼ 1.0 and P oftenreaches −0.05. Thus, we come up with a scenario for the nucleate boiling to start on an ideallysmooth surface.

( 1 ) The liquid temperature near the heating area increases.( 2 ) The liquid thermally expands, and the local liquid pressure fluctuates and sometimes

reaches large negative values.( 3 ) When the pressure reaches P0 determined from the YL equation, tiny bubble nuclei

appear on the surface.( 4 ) When the bubbles are small, they may repeat generation and collapse.( 5 ) Once they exceed a certain size, they stabilize and continuous growth starts.

5. Results: Overall Heating Case

An example for the overall heating simulation is shown in Fig. 8 for the case of neutralwettability (εwl/ε = 1.0) at the wall temperature Tw = 1.3. Since the heating wall is uniform,several bubble nuclei were generated almost simultaneously on the wall.

The local temperature and the local number density were evaluated as functions of thedistance from the wall; two examples are shown in Fig. 9. As the liquid temperature risesin the vicinity of the wall, the local density gradually decreases, and finally nuclei appear.

342

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 8 Nucleate boiling for the overall heating case; sequential snapshots (top) ofsectional view and (bottom) in the vicinity (z < 1.0) of the heating wall.Examples for a surface with neutral wettability (εwl/ε = 1.0) at the walltemperature Tw = 1.3.

Similarly to the partial heating case, nucleation occurs easier on hydrophilic surface sincemore energy is supplied from the wall through the stronger interaction. Carey et al. haveproposed a model(21), which suggests that nucleation occurs slightly (a few nanometers) awayfrom the surface due to the wall–fluid interaction. The change of density profile in Fig. 9 (b)for hydrophilic surface seems in accordance with the model; the first liquid layer is so stronglyattracted to the heating wall that the density decrease takes place first in the region 4 ∼ 8σaway from the surface.

5.1. DiscussionsChange of the bubble volume Vb, the local pressure P and the local temperature T for the

region z < 10, and the heat flux q are shown in Fig. 10 for the case of Tw = 1.3. The volumeVb here is the sum of all “bubbles” (or vacancies) since each bubble has large fluctuations insize and shape. When the periodic boundary conditions for horizontal directions are takeninto consideration, Vb diverges to infinity at the percolation point, but the threshold is out ofthis range. The flux q is obtained from the energy change rate of the fluid, Eq. (6).

In either case, sudden transition of Vb is observed from gradual growth to one-way in-crease (or film-like boiling), from which we define the inception time τi; τi ∼ 600 τ forεwl/ε = 1.0 and ∼ 400 τ for εwl/ε = 1.5, for example; as expected, τi is shorter for walls withlarger εwl . The local pressure fluctuates very much during the nuclei generation. The behaviorof q is also similar to that of partial heating cases; it takes a peak value qmax at t ∼ 50 τ, whichis the time of the wall temperature reaching Tw.

5.2. Wall wettabilityTo show how the wall wettability affects the boiling, the peak value of heat flux qmax

is plotted for various conditions in Fig. 11, where the superheat ΔTs is defined as Tw − Teq

(Teq = 0.7 is the initial temperature of liquid). As expected, more heat is transferred onhydrophilic walls. Note that this qmax is different from the critical heat flux (CHF) of boilingheat transfer; the CHF is the limit value during steady nucleate boiling while qmax here is themaximum flux before transition to nucleate boiling. Thus qmax = 2 × 10−4 (εσ−2τ−1), forexample, roughly corresponds to the heat flux of 22 MW/m2 for water (Table 1), which is oneorder of magnitude larger than typical values of CHF.

The inception time τi is also evaluated, which is the waiting time before the transition to“film-like boiling,” as define in Sec. 5.1. The inverse of τi gives the transition rate J similar to

343

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 9 Time development of the temperature profile T (z) and the density profile ρ(z)for the overall heating simulation with wall temperature Tw = 1.3; (top) neutralwettability wall and (bottom) hydrophilic wall.

the nucleation rate(13). They are plotted against the superheat in Fig. 12. Again the transitionoccurs more easily on the wall with higher wettability. Note that J here is not per unit areanor per unit volume, in contrast to the usual definition of nucleation rate; we suppose thatτi is almost independent of the system size (area or volume). If expressed per volume asin Ref. (13), J = 10−3 τ−1 corresponds to 8.3 × 1032 m−3s−1, recovering a similar order ofmagnitude in Ref. (13).

The dependence on surface wettability may need some comments. Although many worksindicate(5) that CHF increases with improved surface wettability, it was suggested(22) that theheat flux during the steady nucleate boiling is larger for surfaces of larger contact angle (i.e.,less wettability); this is attributed to the change in the fraction of “active” nucleation sites. Incontrast, our simulation gives that both qmax and J increase on surface of larger wettability,because larger wettability is equivalent to stronger wall–fluid interaction and leads to largerheat flux under all conditions.

5.3. Patterned SurfaceTo show the surface wettability effect more clearly, we carried out a series of overall

344

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 10 Overall heating; comparison between neutral (εwl/ε = 1.0) and hydrophilic(εwl/ε = 1.5) walls at Tw = 1.3. Bubble volume Vb, local pressure P and localtemperature T near the heating wall, and heat flux q. Arrows indicate the timeof inception, or transition to film boiling.

Fig. 11 Maximum heat flux on wall of various wettability; overall heating simulations.

heating simulations on “patterned surface.” An example is shown in Fig. 13 for the caseof patterning with εwl/ε = 1.5 and εwl/ε = 1.0 regions. When the wall temperature wasmaintained uniformly at Tw = 1.3 (i.e., ΔTs = 0.6), large decrease of liquid density wasobserved first on the hydrophobic surface (at t = 420 τ) due to the less wall-liquid attraction.However, this decrease did not lead to nucleation on this area, and finally nuclei appeared onthe hydrophilic surface due to larger energy transfer from the surface, which agrees with ourexpectation.

5.4. Size effectAlthough simulations with larger system size would be preferable to see how the system

size affects the results, we have not done simulations with different system sizes yet. Insteadwe here give brief comparison to our former model(14), where liquid of 7,600 LJ particles washeated with a constant rate using a uniform potential wall.

In a case of overall heating on hydrophilic wall (εwl/ε = 1.5), for example, numerousnuclei were generated at a similar time t ∼ 300 − 500 τ. The liquid temperature at the time of

345

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 12 Inception time τi and transition rate J = 1/τi for wall of various wettability;overall heating simulations.

Fig. 13 Example of nucleate boiling on patterned surface; the wall surface of left halfarea is hydrophilic (εwl/ε = 1.5) while the right half is neutral (εwl/ε = 1.0).The whole wall is heated with Tw = 1.3.

nucleus generation was about 0.90 − 0.95, which is also similar to the results in Fig. 9. Fur-thermore, the maximum heat flux qmax for the smaller system was found to be ∼ 20 MW/m2,very close to the results here. The only apparent difference is the pressure; the pressure ob-served in the smaller simulation fluctuated around a large negative value, ∼ −0.4 before thenucleus generation, while that in the simulations here is close to zero on average (Fig. 10).We thus speculate that size effects are not very large for this type of simulations at atomisticscales; the nucleation dynamics is governed very locally under these extreme conditions, asalready mentioned in Sec. 4.1. Size dependence will be more evident under milder conditions(smaller superheat ΔTs), where stochastic characters of nuclei generation will be relevant andalso the heat conduction will reach much further in space.

5.5. Constant pressure simulationWe have so far executed the boiling simulations with a constant-volume two-phase sys-

346

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

Fig. 14 Overall heating simulation under constant pressure Pext. Change of the bubblevolume and the local temperature for the wall of neutral wettability εwl/ε = 1.0at the wall temperature Tw = 1.3.

tem, as shown in Fig. 2; in this setup, we can observe the boiling under a natural conditionwith saturated vapor pressure. However, the liquid temperature increases during long runs(see Fig. 9, for example), followed by increase of vapor pressure.

Another possible type of boiling simulation is a constant-pressure MD method similar tothat used in Ref. (13). Only for the overall heating case, we carried out the constant-pressuresimulation with a movable ceiling; the conditions are Tw = 1.3 and εwl/ε = 1.0. Results(change of bubble volume Vb and the local temperature T ) are shown in Fig. 14 for severalvalues of the external pressure Pext. When we apply negative external pressure, boiling takesplace at earlier time. The results with Pext = 0 are qualitatively similar to those of the previoustwo-phase simulation with the same conditions [Fig. 10 (a)] in respect of the volume changeand the temperature change; generated nuclei started to grow fast when the local temperatureincreased up to T 1.1 in both cases. Note that the local pressure was fluctuating around zeroduring the nucleation stage for the two-phase (constant volume) simulation.

Thus we suggest that the two-phase MD simulation and the constant-pressure one essen-tially give the same results. Each method has its own disadvantage; the temperature and thevapor pressure gradually increases in the two-phase simulation, while the moving ceiling mayaffect the nucleation and bubble growth dynamics in the constant-pressure calculation. Thesewill be less influential, however, when we use a sufficiently large system, especially with avery thick liquid layer.

6. Conclusion

We investigated the initial stage of nucleate boiling on ideally smooth surfaces with MDsimulation. Mainly with a constant-volume two-phase simulation cell, boiling behavior wascompared for partial (spot) heating and overall heating situations.

Although the conventional PEN theory suggests extreme difficulty of nucleate boiling onatomically smooth surface, we observed that bubble nuclei are easily generated, especiallyon wettable (hydrophilic) surfaces. Our simulation suggests the following scenario for theincipience of nucleate boiling without preexisting nuclei: (1) The liquid near the surface isheated, leading to thermal expansion and large pressure fluctuations. (2) When the fluctuatedlocal pressure reaches a sufficiently large negative value, bubble nuclei (or voids) of atomic

347

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

scale are spontaneously generated. This may occur near the spinodal line(20), (21). (3) Once thevoid size exceeds a certain nulceation threshold, bubble growth starts, as a usual nucleationprocess.

The surface wettability largely affects the boiling inception, mainly through the amountof wall–liquid energy transfer; liquid receives more heat and starts nucleation faster on wallsof more wettability.

Finally we confirm that the constant-volume two-phase simulation and the constant-pressure simulation give essentially the same results on nucleate boiling inception.

Acknowledgement

We acknowledge the contribution by our former student, Mr. Shunsuke Kawabe, in pro-gramming and data analysis.

References

( 1 ) Bankoff, S. G., Entrapment of gas in the spreading of a liquid over a rough surface,AIChE J., Vol. 4 (1958) 24–26.

( 2 ) Skripov, V. P., Metastable Liquids, (1974) Wiley, New York.( 3 ) Brennen, C. E., Cavitation and Bubble Dynamics, (1995) Oxford University Press, Ox-

ford.( 4 ) Debenedetti, P. G., Metastable Liquids, (1996) Princeton University Press, Princeton.( 5 ) Dhir, V. K., Boiling heat transfer, Annu. Rev. Fluid Mech., Vol. 30 (1998) 365–401.( 6 ) Theofanous, T. G., Tu, J. P., Dinh, A. T., and Dinh, T. N., The boiling crisis phenomenon

Part I: nucleation and nucleate boiling heat transfer, Exp. Therm. Fluid Sci., Vol. 26(2002) 775–792.

( 7 ) Chen, T., Klausner, J. F., Garimella, S. V., and Chung, J. N., Subcooled boiling incip-ience on a highly smooth microheater, Int. J. Heat Mass Transf., Vol. 49 (2006) 4399–4406.

( 8 ) Bon, B., Guan, C-K., and Klausner, J. F., Heterogeneous nucleation on ultra smoothsurfaces, Exp. Therm. Fluid Sci., Vol. 35 (2011) 746–752.

( 9 ) Matsumoto, M. and Tanaka, K., Nano bubble –Size dependence of surface tension andinside pressure, Fluid Dynamics Res., Vol. 40 (2008) 546–553.

(10) Matsumoto, M., Surface tension and stability of a nano bubble in water: Molecularsimulation, J. Fluid Sci. Tech., Vol. 3 (2008) 922–929.

(11) Park, S. H., Weng, J. G., and Tien, C. L., A molecular dynamics study on surface tensionof micro bubbles, Int. J. Heat Mass Transf., Vol. 44 (2001) 1849–1856.

(12) Kinjo, T. and Matsumoto, M., Cavitation processes and negative pressure, Fluid PhaseEquilibria, Vol. 144 (1998) 343–350.

(13) Novak, B. R., Maginn, E. J., and McCready, M. J., An atomistic simulation study of therole of asperities and indentations on heterogeneous bubble nucleation, J. Heat Transf.,Vol. 130 (2008) 042411:1–9.

(14) Matsumoto, M. and Yamamoto, T., Initial stage of nucleate boiling at atomic scale, Proc.8th ASME–JSME Thermal Eng. Joint Conf. (2011) AJTEC2011-44435.

(15) Matsumoto, M., Wakabayashi, H., and Makino, T., Thermal resistance of crystal inter-face: Molecular dynamics simulation Trans. JSME Ser. B, Vol. 68 (2002) 1919-1925;English translation, Heat Transfer – Asian Research, Vol. 34 (2005) 135–146.

(16) Taura, G. and Matsumoto, M., Molecular dynamics simulation of micro droplet im-pingement on solid surface, J. Fluid Sci. Tech., Vol. 5 (2010) 207–218.

(17) Kotaidis, V., Dahmen, C., von Plessen, G., Springer, F., and Plech, A., Excitation ofnanoscale vapor bubbles at the surface of gold nanoparticles in water, J. Chem. Phys.,Vol. 124 (2006) 184702:1–7.

(18) Minnaert, M., On musical air-bubbles and the sounds of running water, Phil. Mag. Ser. 7,Vol. 16 (1933) 235–248.

348

Journal of ThermalScience and Technology

Vol.7, No.1, 2012

(19) Lotfi, A., Vrabec, J., Fischer, J., Vapor liquid equilibria of the Lennard-Jones fluid fromthe NpT plus test particle method, Molecular Physics, Vol. 76 (1992) 1319–1333.

(20) Wu, Y. W. and Pan, C., A molecular dynamics simulation of bubble nucleation in homo-geneous liquid under heating with constant mean negative pressure, Microscale Ther-mophys. Eng., Vol. 7 (2003) 137–151.

(21) Carey, V. P. and Wemhoff, A. P., Thermodynamic analysis of near-wall effects on phasestability and homogeneous nucleation during rapid surface heating, Int. J. Heat MassTransf., Vol. 48 (2005) 5431–5445.

(22) Wang, C. H. and Dhir, V. K., Effect of surface wettability on active nucleation site den-sity during pool boiling of water on a vertical surface, J. Heat Transf., Vol. 115 (1993)659–669.

349