The effects of wall superheat and surface wettability on nucleation site interactions during boiling

International Journal of Heat and Mass Transfer 146 (2020) 118820

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

The effects of wall superheat and surface wettability on nucleation site interactions during boiling Lei Zhang a,b, Tao Wang a,b,⇑, Seolha Kim a, Yuyan Jiang a,b,* a b

Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences, Beijing 100049, China

a r t i c l e

i n f o

Article history: Received 9 June 2019 Received in revised form 17 September 2019 Accepted 29 September 2019

Keywords: Lattice Boltzmann method Wettability Nucleation site interactions Boiling

a b s t r a c t Most boiling heat transfer experimental correlations assume that the formation and growth of bubbles at adjacent nucleation sites are independent, but experimental results show that interactions do occur between adjacent nucleation sites, while, the experimental results have some divergence because of the complexity of boiling. In this paper, we conducted a numerical study by lattice Boltzmann method on the nucleation site interactions during pool boiling. The LBM is featured of a multiple-relaxationtime algorithm with hybrid thermal scheme which can carefully control the parameters that influence boiling and keep other factors unchanged. The numerical results show that temperature at the given nucleation sites has strong dependence at short separation distance, and consequently the wall superheat could change the thermal and hydrodynamic interactions between nucleation sites. Moreover, wettability also plays a critical role in boiling, so the effects of wettability on nucleation site interactions were also studied, results show that nucleation site interactions on hydrophobic surface are mainly promotive, which is different from hydrophilic surface. These findings can be used to better understand the mechanism of nucleation site interactions and explain the divergence of previous experimental results. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Nucleate boiling is one of the most effective heat transfer modes, by which energy is transported from a heating surface through the movement of bubbles. In evaluating the heat transfer performance of nucleate boiling, nucleation sites density is an important parameter. But available predictive relations or experimental correlations assumed implicitly that the formation and growth of bubbles at adjacent nucleation sites was independent, which means that the formation and growth of a bubble at one nucleation site would have no effect on the formation and growth of a bubble at adjacent nucleation sites. But recent investigation shows that interactions do occur between bubbles at adjacent nucleation sites. In the study of Eddington and Kenning [1], they found that some cavities that should have been active did not form bubbles when boiling conditions were established, whereas other cavities that should not have been active did form bubbles when boiling conditions were established. They explained that this behavior was caused by thermal interference at adjacent nucleation sites. ⇑ Corresponding authors at: Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing 100190, China. E-mail addresses: [email protected] (T. Wang), [email protected] (Y. Jiang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118820 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

Shoukri and Judd [2] also found that some cavities theoretically active could not nucleate under test conditions, and then Judd and his cooperators started to study nucleation site interactions [3–7]. They chose dimensionless separation distance S/Dd as a key parameter to describe the interference range. S represents the separation distance between nucleation sites and Dd is the average bubble departure diameter, which varied with wall superheat. Their results showed that for dimensionless separation distance S/Dd < 1, the effects of nucleation site interactions can be described as promotive, which the bubble formation at the active nucleation site can promote the bubble formation at the less active nucleation site, while for dimensionless separation distance 1 < S/ Dd < 3, the effects of nucleation site interactions can be described as inhibitive, which the bubble formation at one nucleation site would inhibit bubble formation at the adjacent nucleation sites and for dimensionless separation distance S/Dd > 3, the effects of nucleation site interactions can be described as independent, there was no interaction between adjacent nucleation sites. Mallozzi et al. [8] studied the bubble interaction using computer simulation called BOILSIM, their results were in good agreement with the experiment results conducted by Judd and Chopra [7] which also divided the interaction effects into promotion, inhibition and independent three regions.

2

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

Nomenclature a, b c cp Cov Dd e E f F Fads Fb FB Fd Fi Fm g hfg Ja l.u. LH l0 m} q s S t*

parameter in P-R equation of state lattice speed (m
s1 ) 1 specific heat (J
kg
K1 ) covariance bubble departure diameter lattice velocity vector (m
s1 ) expectation distribution function body force vector (N) fluid-solid interaction force (N) gravitational force (N) buoyancy force (N) drag force (N) added mass inertia force (N) interaction force (N) gravitational acceleration (m
s2 ) 1 latent heat (J
kg ) Jacob number lattice unit length of heater (m) characteristic length (m) steam mass flow heat flux (W
m2 ) entropy (J
K1 ) distance between nucleation sites, force (N) dimensionless time

In the studies of Zhang and Shoji [9,10], they manufactured twin cavities as the nucleation sites and chose laser irradiation as the heating source. Based on their results, they classified the interaction effects into three factors: hydrodynamic interaction, thermal interaction and bubble coalescence. Different with the results of Judd, their results showed that there were four different regions: two promotion regions because of hydrodynamic interaction and bubble coalescence, respectively, one inhibitive region because of thermal interaction and the independent region. While Hutter et al. [11] fabricated artificial cavities by means of microfabrication to make well defined and exactly positioned nucleation sites to study the interval influence. They changed the wall superheat to control the bubble departure diameter to switch the dimensionless separation distance S/Dd. They investigated the interactions between neighboring cavities in terms of the bubble departure frequency. Their results showed that the strong inhibitive and promotive regions of the departure frequency with decreasing space interval between cavities didn’t occur in their study. From the literatures review above, it can be seen that not only the criteria of dimensionless separation distance S/Dd for each region are different in these researches, but also the number of interaction regions and whether the interaction has effects or not on bubble departure frequency is not consistent. This is due to the complexity of boiling phenomenon. As indicated by Shoji [10], boiling is a conjugate phenomenon, and many subprocesses interact with each other. It is difficult to investigate the basic mechanism of each process separately in isolated conditions, because it is hard to change one effect factor and keep all others unchanged. That’s why we need computer simulation in studying of boiling. Due to the unknown underlying interaction mechanisms of the system, it is difficult to carry out comprehensive simulation work on boiling with conventional numerical method. In recent years, the lattice Boltzmann (LB) method, which is a mesoscopic approach based on the kinetic Boltzmann equation,

t0 T v

characteristic time (s) temperature (K) velocity

Greek symbol x acentric factor W effective mass k thermal conductivity (W
m1
K1 ) r surface tension (N
m1 ), standard deviation t kinematic viscosity (m2
s1 ) s relaxation time q density (kg
m3 ), correlation coefficient h contact angle (°) l kinetic viscosity (Pa
s), mean value Subscripts or superscripts ave average c critical eq equibilium l liquid t time v vapor sat saturation w wall

has received great attention and shows promising future in phase change phenomena [12]. In pseudopotential LB model, there is an interparticle potential to mimic the interaction between fluid molecules, so the phase separation between different phases can naturally be achieved without using the interface-tracking or interface-capturing techniques. The first pseudopotential thermal LB model was proposed by Zhang and Chen [13], they introduced temperature-dependent body force to couple the pseudopotential multiphase LB method with a scalar temperature equation and calculated liquid-vapor boiling process. Recently, Gong [14] and Li [15] proposed more accurate and more robust schemes, which were widely adopted to solve phase change problems. Gong and Cheng [16] investigated the bubble ebullition cycle under different contact angles and heat fluxes. They also investigated the cavity effects on bubble interaction and boiling performance [17]. Gong [18] and Ma [19] also used lattice Boltzmann method to investigate the boiling performance of surface with mixed wettability and Lee [20] adopted pseudopotential thermal LB model to design the radial distribution of patterned wettability. Li et al. [21] studied boiling performance on surface with mixed wettability and microstructure by means of pseudopotential LB model. Those pseudopotential models exclude empirical correlations, and the nucleation process could be simulated automatically. In this study, the pseudopotential thermal lattice Boltzmann model was adopted to carefully control the surface wettability and wall superheat to investigate their effects on nucleation site interactions.

2. Model description A multiple-relaxation-time (MRT) thermal lattice Boltzmann model is used in this work, and a brief description of the model will be presented in this part. The MRT collision operator, which can improve the stability of the numerical model, is an important extension of the relaxation LB method proposed by Higuera et al.

3

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

[22,23]. To solve boiling problem, we need to calculate the temperature field, and there are two ways to solve the thermal LB equation. The first one is using a temperature distribution function to mimic the macroscopic temperature equation and the other one is calculating the temperature field with a traditional numerical method such as finite-volume method or finite-difference method. In this work, a finite-difference method is adopted to solve the temperature equation. 2.1. MRT pseudopotential lattice Boltzmann method for multiphase flow The particle distribution functions of MRT collision operator can be expressed as follows:



 0 eq f a ðx þ ea dt ; t þ dt Þ ¼ f a ðx; t Þ  M 1 KM f b  f b þ dt F a ; ab

ð1Þ

where f a ðx; tÞ is the density distribution function and ea is the diseq crete velocity at location x, f b is the equilibrium density distribution function, M is an orthogonal transformation matrix and can be 0

found in [24],F a is the forcing term Using transformation matrix M, Eq. (1) can be rewritten as:

  K m ¼ m  Kðm  meq Þ þ dt I  S; 2 eq

moment space. m ¼ Mf ,meq ¼ mf stands for the distribution function after transformation. This is the collision process and the streaming process can be expressed as 



ð3Þ

1

where f ¼ M m . For the D2Q9 lattice, discrete velocity vectors ea can be described as

8 > < ð0; 0Þ ea ¼ ð1; 0Þc; ð1; 0Þc > : ð1; 1Þc

a¼ 0 a¼ 1-4 a ¼ 5 - 8;

and the equilibria meq ¼ mf obtained as

eq

in the moment space can be

ð5Þ where jv j2 ¼ v 2x þ v 2y . The macroscopic density and velocity are defined as

q¼

a

f a ; qv ¼

X a

dt ea f a þ F; 2

ð6Þ

ð7Þ

where qav e is the averaged fluid density of the computational domain and g is the gravitational acceleration. The interaction force is the most important force in pseudopotential model which controls the phase separation and given as follows [25]

F m ¼ GWðxÞ

X

wa Wðx þ ea Þea ;

ð9Þ

h

pffiffiffiffiffiffiffiffiffiffi i2

uðT Þ ¼ 1 þ ð0:37464 þ 1:54226x  0:26992x2 Þð1  T=T c Þ , with x being the acentric factor, chosen as 0.344 for water, other variables in P-R equation of state were chosen as a = 3/49, b = 2/21 and R = 1. The fluid-solid interaction force is given by

F ads ¼ Gw WðxÞ

X

xa Sðx þ ea Þea ;

ð10Þ

a

where Gw can be used to adjust the interaction strength to obtain different contact angles, xa ¼ wa =3, Sðx þ ea Þ ¼ uðxÞsðx þ ea Þ, uðxÞ can be set to WðxÞ and sðx þ ea Þ is an indicator function that equals 0 for fluid and 1 for solid phase. In MRT LB model, the forces are incorporated via an improved forcing scheme and given by [26]

3 0 2 7 6 j 7 6 6v
F þ w2 dcjðFsm0:5 Þ 7 t e 6 7 6 7 6 cjF m j2 6 6v
F  w2 d ðs 0:5Þ 7 t f 7 6 7 6 7 6 F x 7 6 7 S¼6 7 6 F x 7 6 7 6 7 6 Fy 7 6 7 6 F 7 6 y 7 6 6 2ðv F  v F Þ 7 5 4 x x y y

ð11Þ

where F is the total force and jF m j2 ¼ F 2m;x þ F 2m;y , v is the macroscopic velocity,c is a parameter used to tune the mechanical stability condition. 2.2. Energy equation for temperature field The temperature equation is derived on the local balance law for entropy, which is given by (neglecting the viscous heat dissipation)

qT

where F ¼ ðF x ; F y Þ is the total force acting on the system, basically it contains gravitational force F b , interaction force between fluid particles F m and interaction force between fluid and solid F ads . The gravitational force is given by

  F b ¼ q  qav e g;

where

qRT auðT Þq2  ; 1  bq 1 þ 2bq  b2 q2

v xFy þ v yFx

ð4Þ


T meq ¼ q 1; 2 þ 3jv j2 ; 1  3jv j2 ; v x ; v x ; v y ; v y ; v 2x  v 2y ; v x v y ;

X

pEOS ¼

2

ð2Þ

  0 where I is the unit tensor and I  K2 S ¼ MF is the force term in

f a ðx þ ea dt ; t þ dt Þ ¼ f a ðx; tÞ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ðpEOS  qc2s Þ=Gc2 is positive, pEOS stands for pressure and obtained from the equation of state. In this work, the PengRobinson (P-R) equation of state was used, which is given by

WðxÞ ¼

ð8Þ

a

where G is a sign parameter, which can make sure the quantity under the square root sign of the ‘‘effective mass”

Ds ¼ r
ðkrT Þ; Dt

ð12Þ

where s is the entropy and k is the thermal conductivity. Substituting the following thermodynamic relation into Eq. (12)

Tds ¼ cv dT þ T

    @pEOS 1 ; d @T q q

ð13Þ

temperature equation can be rewritten as

qcv

  DT @pEOS ¼ r
ðkrT Þ  T r
v: Dt @T q

Using relation

Dð
Þ Dt

¼ @ t ð
Þ þ v
rð
Þ. Eq. (14) can be transformed

to

@ t T ¼ v
rT þ

ð14Þ

1 T r
ðkrT Þ  qc v qc v

  @pEOS r
v: @T q

ð15Þ

The forth-order Runge-Kutta scheme is adopted to solve this unsteady temperature equation. Using K ðT Þ to represent the right-hand side of Eq. (15), we can get

4

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

T tþdt ¼ T t þ

dt ðh1 þ 2h2 þ 2h3 þ h4 Þ; 6

ð16Þ

      h1 ¼ K T t ; h2 ¼ K T t þ d2t h1 ; h3 ¼ K T t þ d2t h2 ;  t  h4 ¼ K T þ dt h3 : Conventional conjugate thermal method [27] is used to get temperature and flux continuity at the solid-fluid interface. where

2.3. Code validation and computational setup 2.3.1. Code validation A thin film evaporation problem was adopted to validate our code, this benchmark problem was also used by other researchers to validate their code [28,29]. The computational domain for thin film evaporation is illustrated as Fig. 1(a), if the input heat flux at the bottom wall is small enough, no vapor bubble will emerge and liquid only evaporate on the liquid-vapor interface. The relation of heat flux and evaporate rate can be given as:

q m ¼ ; hfg 00

Fig. 2. Effect of gravity on bubble departure diameter.

ð17Þ

00

where m and hfg are steam mass flow and latent heat of phase change at 0.86T c , respectively. It can be seen from Fig. 1(b) that the LB simulation results are consistent with analytical results. To validate our code in terms of bubble departure, the effect of gravity on bubble departure diameter was also tested. Fritz [31] proposed the following correlation for bubble departure diameter from a horizontal heated surface:

"

r  D ¼ 0:0208h  g ql  qv

#0:5

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l0 ¼



r

g ql  qv

; u 0 ¼

pffiffiffiffiffiffi l0 gl0 ; t 0 ¼ ; u0

where r is the surface tension, ql and qv are densities of saturated liquid and vapor, respectively. The dimensionless length, dimensionless velocity and dimensionless time can be calculated by 

:

ð18Þ

Eq. (18) suggests that bubble departure diameter is proportional to g 0:5 , Fig. 2 shows the simulation results and its fitting curve, it can be seen that the results agree quite well with the empirical correlation. 2.3.2. Computational setup In this work, the capillary length l0 is chosen as the reference length. Capillary length is a typical characteristic variable for describing liquid-vapor two phase flow, and it is widely used in previous boiling literatures [32–34]. The corresponding reference velocity u0 and reference time t0 are given by:

ð19Þ

l ¼ l=l0 ; u ¼ u=u0 ; t  ¼ t=t 0 :

ð20Þ

Fig. 3 illustrate the computational domain, as can be seen in the figure, a heating plate located at the central part of the bottom wall and two cavities, as the nucleation sites, are located on the top surface of the heating plate. The size of heated plate is LH ¼ 200 lattices in length and H = 20 lattices in thickness. The width of the two cavities is different with each other, the right one is 4 lattices wider than the right one to nucleate more easily, so we can study the interference effects of one active nucleation site on its adjacent nucleation sites. After grid dependency checking, we chose a 600  800 lattice size for the rectangular computational domain which could ensure accuracy and minimize computation time. Initially, the computation domain is filled with saturated water at saturated temperature T sat ¼ 0:86T c and the top side of the

Fig. 1. Thin film evaporation: (a) Schematic of computational domain; (b) Comparison of LBM results and analytical results [30].

5

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

3. Results and discussions 3.1. The effects of wall superheat

Fig. 3. Computational setup.

computation domain is outlet boundary with constant pressure, while periodic boundary condition is imposed on the left and right boundaries. The corresponding liquid properties are represented in lattice units: saturated liquid and vapor density ql ¼ 6:5, qv ¼ 0:38, kinetic viscosity tl ¼ 0:04, tv ¼ 0:1732. The specific latent heat is calculated using the theoretical model proposed by Gong and Cheng [16], and the value is related with the parameters in equation of state and temperature. In our study, the value is hfg ¼ 0:5720; and thermal conductivity is chosen as kl ¼ 0:48. For solid part, cp ¼ 1:5, ks ¼ 1. The temperature of the bottom side of the heating plate is maintained at T w , (expressed in terms of dimensionless number Ja ¼ cp;l ðT w  T sat Þ=hfg , where Tw and Tsat are wall temperature and saturation temperature, cp;l and hfg are the specific heat of water and latent heat, respectively.) and the left and right sides of the heater are adiabatic.

To evaluate the effects of thermal interaction, the temperature fluctuation with time at two nucleation sites was recorded. Fig. 4 shows the temperature relation of the two adjacent nucleation sites with different separation distance. Overall, the temperature relation between these two points is weakening as the separation distance increasing. The temperature distribution at Point A and Point B is closely spaced in a narrow area when dimensionless separation distance S/Dd = 0.65, as shown in Fig. 4(a). As the dimensionless separation distance increases, the temperature distributes in a larger area but the temperature still shows strong relation at Point A and Point B, as shown in Fig. 4(b). However, when the dimensionless separation distance exceeds a certain value, the temperature values at Point A and Point B are irregularly dispersed, as shown in Fig. 4(c) and Fig. 4(d), indicating that the thermal interaction of these two nucleation sites is gradually weakened. For the evaluation of the temperature dependency of two nucleation sites quantitatively, the correlation coefficient between the temperature fluctuation was adopted. According to the definition of the correlation coefficient, the correlation coefficient between two variables can be defined as the ratio of its covariance to standard deviation, namely:

qX;Y ¼

Cov ðX; Y Þ

rX rY

¼

  

E X  lX Y  lY

rX rY

;

ð21Þ

where q is correlation coefficient and Cov is covariance. l and r represent the mean value and standard deviation, respectively. E is expectation. To calculate the correlation coefficient of a sample, the correlation coefficient can be calculated via:

Fig. 4. The temperature history of two nucleation sites (Point A and Point B) (a) S/Dd = 0.65; (b) S/Dd = 1.6; (c) S/Dd = 2.3; (d) S/Dd = 2.9.

6

qX;Y

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820  
 Pn
Yi  Y i¼1 X i  X rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2 Pn
Pn
X  X Y  Y i i i¼1 i¼1

ð22Þ

where X and Y stand for the temperature at the two nucleation sites, respectively. Fig. 5 shows the correlation coefficient of two nucleation sites calculated by Eq. (22). As can be seen in the graph, the correlation coefficient decreases as the dimensionless distance increases. At dimensionless distance S/Dd = 0.65, the correlation coefficient is 0.82, implies that strong thermal interaction happens between the two nucleation sites. At S/Dd = 1.6, the correlation coefficient could still be 0.69, while, as the dimensionless distance keep increasing, the correlation coefficient decreases sharply, at S/Dd = 2.9. The correlation coefficient decreases to 0.004, close to zero, which indicates that there is no interaction between these nucleation sites. As discussed above, the thermal interaction is strong when two nucleation sites are close to each other. But the effects of wall superheat on nucleation site interactions were not investigated by previous researches to the limit of the researcher’s knowledge. Although Hutter et al. [11] changed the wall superheat to study nucleation site interactions, they just used it as way to obtain different dimensionless separation distance, rather than to analyze the effects of wall superheat. Fig. 6 shows the averaged bubble departure frequency varied with dimensionless separation distance under various wall superheat. The averaged bubble departure frequency is usually used as an indicator to evaluate the nucleation site interactions [9,11,35]. Generally, the bubble departure frequency increases as the wall superheat increases. At wall superheat Ja = 0.218, the average bubble release frequency can get 0.060 (lattice unit), while at Ja = 0.178, the average bubble release frequency is around 0.048, 20% lower than that of the case of Ja = 0.218. Moreover, at lower wall superheat (Ja = 0.178 & Ja = 0.184), the inhibitive interaction of the adjacent nucleation sites is clear, as can be seen in Fig. 6, there is a sharp decrease of average bubble departure frequency at dimensionless separation distance S/Dd ~ 1.6. While, for higher wall superheat (Ja = 0.218 & Ja = 0.195), this inhibitive interaction is not obvious. At low wall superheat (Ja = 0.178 & Ja = 0.184), the interaction can be clearly divided into promotive, inhibitive, promotive and independent four regions, which is consistent with Zhang and Shoji’s results [9], while, at high wall superheat (Ja = 0.218), the boundaries of those regions are not clear and only promotive effect is observed.

Fig. 5. Temperature correlation coefficients between two nucleation sites at different separation distance.

Fig. 6. Averaged bubble release frequency at diffident separation distance ant the effects of wall superheat.

This implies that wall superheat also plays a critical role in nucleation site interactions. Generally, the thermal interaction has inhibitive effect on adjacent nucleation site, since heat is removed from the heating surface because of liquid evaporation during bubble growth and the cooling effect would generate a negative effect on the adjacent bubble activity. For the hydrodynamics interaction, it can be analysis with force balance of bubble departure [9]. When the buoyancy force is bigger than the forces that holding bubble, the bubble will depart, namely:

FB  Fr þ Fi þ Fd;

ð23Þ

where F B is buoyancy force, F r is surface tension force, F i is added mass inertia force, F d is drag force. F d is related to hydrodynamic interaction among all these forces, which given by:

Fd ¼

    ds 1 4 2 ds  tw   tw : ql D C D 2 p dt dt

ð24Þ

In Eq. (24), the larger the liquid wake velocity tw , the smaller the drag force F d is. This indicates that the hydrodynamic interaction is promotive for the bubble departure. Based on the above analysis, the hydrodynamic interaction and thermal interaction can be illustrated as Fig. 7. At small dimensionless separation distance S/Dd, the intensity of hydrodynamic interaction is stronger

Fig. 7. Schematic illustration of thermal and hydrodynamic interaction.

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

than the intensity of thermal interaction because of bubble coalescence and in this region, nucleation site interactions are promotive. As the dimensionless separation distance S/Dd increase, no bubble coalescence happens and the intensity of thermal interaction is dominate, so the nucleation site interactions are inhibitive. As the separation distance becomes larger, both hydrodynamics interaction and thermal interaction intensity decrease, but hydrodynamics interaction still stronger than thermal interaction because of bubble and bulk liquid interaction, and the nucleation site interactions are slightly promotive in this region. The hydrodynamics interaction and thermal interaction keep weakening as the dimensionless separation distance increasing and their influence in nucleation site interactions becomes weak, so there are no nucleation site interactions happen, i.e. the independent region. As wall superheat increases, the bubble growth and departure becomes more vigorous and the hydrodynamic interaction between nucleation sites is enhanced [36,37]. Generally, there are two limiting cases of the bubble growth and departure process. One is inertia-controlled growth and the other one is heat transfer controlled growth [38,39]. While, as the wall superheat increase, the liquid superheat layer could be easily rebuilt and the cooling effect of the microlayer evaporation could be weakened. So the intensity of thermal interaction is decreasing as the wall superheat rising and the hydrodynamics interaction becomes domain gradually. This is consistent with the simulation results showed by Fig. 6. 3.2. The effects of wettability Previous researches are all based on hydrophilic surface (h = 45) and the above analysis is also on hydrophilic surface. The nucleation site interactions on hydrophobic surface (h = 110) are missing. Therefore, the nucleation site interactions on hydrophobic surface are analyzed here. Fig. 8 shows the bubble dynamics on hydrophobic surface at wall superheat Ja = 0.178 with a dimensionless separation distance S/Dd = 0.8. As can be seen in the figure, because the separation distance is shorter than base diameter of vapor bubble, the two nucleation sites meagered into one nucleation site and only one single bubble emerged on the surface (t* = 1.15  t* = 1.35). The two nucleation sites acted like one nucleation site and only one bubble growth and departure from the surface.

7

However, when the separation distance is larger than the base diameter of the vapor bubble, bubble coalescence still happened, but the three-phase contact line would retract during the departure process and the two nucleation sites did not merge into one nucleation site (t* = 0.06 & t* = 0.10). Different with hydrophilic surface, a seed bubble would leave on the heating surface after preceding bubble departure [40], so one more bubble emerged on the heating surface, as can be seen in Fig. 9 (t* = 0.17). This seed bubble grew up and coalesced with adjacent vapor bubble (t* = 0.23 & t* = 0.56). This indicates that although there is thermal interaction and locally cooling effect because of liquid evaporation, its inhibitive effect is limited on hydrophobic surface and hydrophobic surface tend to has more nucleation sites. As the separation distance increase, the coalescence mechanism becomes different. Fig. 10 shows the nucleation site interactions at S/Dd = 2.4 (upside figures) and S/Dd = 3.1 (downside figures). Although there is still one extra bubble emerged on the center of the heating surface (t* = 0.58, t* = 0.43), its mechanism is different with the case of S/Dd = 1.3. At the beginning of the bubble ebullition, new bubbles emerged at the given nucleation sites, and the larger one departed first, the wake effect dragged the bubble at the adjacent nucleation to slide on the heating surface (t* = 0.52, t* = 0.41), since hydrophobic surfaces have the property of repelling water and favoring vapor, the vapor bubble can easily slide on the surface, which is different with hydrophilic surfaces. The bubble completely separated from the original nucleation site during the sliding process and departed on the center of the two nucleation sites (t* = 0.58, t* = 0.43). The left seed bubble merged with the adjacent bubble and accelerated bubble release. This implies that the hydrodynamic interaction is stronger on hydrophobic surface than that on the hydrophilic surface. To investigate the effect of wall superheat on nucleation site interactions of hydrophobic surface, the case of Ja = 0.218 was studied numerically, as can be seen in Fig. 11. The interaction model is different with low wall superheat case, more bubble emerged on the heating surface because the wall superheat was high and bubble could easily nucleate on hydrophobic surface. At S/Dd = 0.6, the two given nucleation sites merged into one and two bubbles nucleated on the left and right side of the given nucleation sites (t* = 0.17). These bubbles coalesced with the bubble at the given nucleation sites (t* = 0.29), and a big bubble formed and departed

Fig. 8. Nucleation site interactions at S/Dd = 0.8 and Ja = 0.178.

8

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

Fig. 9. Nucleation site interactions at S/Dd = 1.3 and Ja = 0.178.

Fig. 10. Nucleation site interactions at S/Dd = 2.4 (up), S/Dd = 3.1 (down) and Ja = 0.178.

quickly (t* = 0.32  t* = 0.34). Then next cycle began. As the separation distance increases, this phenomenon also happens. As discussed above, it is difficult to use bubble departure frequency as an indicator to evaluate the effect of nucleation site interactions on hydrophobic surface, because bubbles not only emerge at the given nucleation sites, but also on the smooth heating surface. So, the averaged vapor evaporation rate is adopted to evaluate the intensity of nucleation site interactions. As can be seen in Fig. 12, basically, the average volume of evaporated vapor with two nucleation sites is bigger than that with single bubble case, because of enhanced hydrodynamic interaction has promotive effect on bubble release. Looking the graph in more details, it can be seen that the average volume of evaporated vapor with two given nucleation sites at low wall superheat (Ja = 0.178) is

lower than the case of single bubble case on hydrophobic surface. This may due to the merging of the two nucleation sites that makes it not take the merits of two nucleation sites, as discussed with Fig. 8. While, as the typical sliding phenomenon happening on hydrophobic surface, more bubbles emerge on the heating surface and boiling is enhanced by hydrodynamic interaction. At higher wall superheat, bubbles could easily nucleate on the smooth heating surface and the average volume of vapor is also higher than single bubble case. Moreover, compared with hydrophilic surface, the nucleation site interactions on hydrophobic surface in terms of vapor generated have no obvious relation with the dimensionless separation distance S/Dd, as can be seen in Fig. 12, except for the case of low wall superheat and separation distance shorter than bubble base diameter.

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

9

Fig. 11. Nucleation site interactions at S/Dd = 0.6 and Ja = 0.218.

experiment. This is because Hutter tuned the dimensionless separation distance by changing wall superheat, while wall superheat also influences nucleation sties interaction.  There is a big difference between nucleation site interactions on hydrophobic surface and hydrophilic surface. Vapor bubble can easily slide on the hydrophobic surface because it favors vapor and repel water, and this makes the hydrodynamic interaction dominate in the nucleation site interactions. Therefore, the nucleation site interactions on hydrophobic surface mainly has promotive effect.

Declaration of Competing Interest The authors declare no conflict of interest. Acknowledgements Fig. 12. Average volume of vapor generated at different nucleation sites separation distance on hydrophobic surface and hydrophilic surface.

This work is supported by National Natural Science Foundation of China (No. 51876203) and China-Japan Research Cooperative Program (No. 2016YFE0118100).

4. Concluding remarks

References

In this study, the effects of wall superheat and surface wettability on nucleation site interactions were investigated numerically and the results can be used to explain the divergence of previous experiment findings. Our conclusions are the followings.

[1] R.I. Eddington, D.B.R. Kenning, A.I. Korneichev, Comparison of gas and vapour bubble nucleation on a brass surface in water, Int. J. Heat Mass Transf. 21 (7) (1978) 855–862. [2] M. Shoukri, R.L. Judd, Nucleation site activation in saturated boiling, J. Heat Transf. 97 (1) (1975) 93–98. [3] R.L. Judd, C.H. Lavdas, The nature of nucleation site interaction, J. Heat Transf. 102 (3) (1980) 461–464. [4] M. Sultan, R.L. Judd, Interaction of the nucleation phenomena at adjacent sites in nucleate boiling, J. Heat Transf. 105 (1) (1983) 3–9. [5] A. Calka, R.L. Judd, Some aspects of the interaction among nucleation sites during saturated nucleate boiling, Int. J. Heat Mass Transf. 28 (12) (1985) 2331–2342. [6] R.L. Judd, On nucleation site interaction, J. Heat Transf. 110 (2) (1988) 475–478. [7] R.L. Judd, A. Chopra, Interaction of the nucleation processes occurring at adjacent nucleation sites, J. Heat Transf. 115 (4) (1993) 955–962. [8] R. Mallozzi, R.L. Judd, N. Balakrishnan, Investigation of randomness, overlap and the interaction of bubbles forming at adjacent nucleation sites in pool boiling, Int. J. Heat Mass Transf. 43 (18) (2000) 3317–3330. [9] L. Zhang, M. Shoji, Nucleation site interaction in pool boiling on the artificial surface, Int. J. Heat Mass Transf. 46 (3) (2003) 513–522. [10] M. Shoji, Studies of boiling chaos: a review, Int. J. Heat Mass Transf. 47 (6) (2004) 1105–1128.

 The temperature under two given nucleation sites has strong interaction at small separation distance, and the temperature correlation coefficient is around 0.8 at S/Dd = 0.65, which means they have strong dependency. While, as the separation distance increase, the correlation coefficient decreases sharply, and the thermal interaction weakening.  The wall superheat could influence the nucleation site interactions by changing the relatively intensity of thermal interaction and hydrodynamic interaction. This can be used to explain the divergence of previous experiment findings, which Zhang and Shoji divided the interaction region into four parts, while, Hutter found no obvious nucleation site interactions in their

10

L. Zhang et al. / International Journal of Heat and Mass Transfer 146 (2020) 118820

[11] C. Hutter, K. Sefiane, T.G. Karayiannis, A.J. Walton, R.A. Nelson, D.B.R. Kenning, Nucleation site interaction between artificial cavities during nucleate pool boiling on silicon with integrated micro-heater and temperature microsensors, Int. J. Heat Mass Transf. 55 (11–12) (2012) 2769–2778. [12] Y.Y. Yan, Y.Q. Zu, B. Dong, LBM, a useful tool for mesoscale modelling of single-phase and multiphase flow, Appl. Therm. Eng. 31 (5) (2011) 649–655. [13] R. Zhang, H. Chen, Lattice Boltzmann method for simulations of liquid-vapor thermal flows, Phys. Rev. E 67 (6) (2003) 066711. [14] S. Gong, P. Cheng, A lattice Boltzmann method for simulation of liquid–vapor phase-change heat transfer, Int. J. Heat Mass Transf. 55 (17–18) (2012) 4923– 4927. [15] Q. Li, Q.J. Kang, M.M. Francois, Y.L. He, K.H. Luo, Lattice Boltzmann modeling of boiling heat transfer: The boiling curve and the effects of wettability, Int. J. Heat Mass Transf. 85 (2015) 787–796. [16] S. Gong, P. Cheng, Lattice Boltzmann simulation of periodic bubble nucleation, growth and departure from a heated surface in pool boiling, Int. J. Heat Mass Transf. 64 (2013) 122–132. [17] S. Gong, P. Cheng, Two-dimensional mesoscale simulations of saturated pool boiling from rough surfaces. Part II: Bubble interactions above multi-cavities, Int. J. Heat Mass Transf. 100 (2016) 938–948. [18] S. Gong, P. Cheng, Numerical simulation of pool boiling heat transfer on smooth surfaces with mixed wettability by lattice Boltzmann method, Int. J. Heat Mass Transf. 80 (2015) 206–216. [19] X. Ma, P. Cheng, X. Quan, Simulations of saturated boiling heat transfer on bioinspired two-phase heat sinks by a phase-change lattice Boltzmann method, Int. J. Heat Mass Transf. 127 (2018) 1013–1024. [20] J.S. Lee, J.S. Lee, Conjugate heat transfer analysis for the effect of the eccentricity of hydrophobic dot arrays on pool boiling, Appl. Therm. Eng. 110 (2017) 844–854. [21] Q. Li, Y. Yu, P. Zhou, H.J. Yan, Enhancement of boiling heat transfer using hydrophilic-hydrophobic mixed surfaces: A lattice Boltzmann study, Appl. Therm. Eng. 132 (2018) 490–499. [22] F.J. Higuera, J. Jiménez, Boltzmann approach to lattice gas simulations, EPL (Europhys. Lett.) 9 (7) (1989) 663. [23] F.J. Higuera, S. Succi, R. Benzi, Lattice Gas Dynamics with Enhanced Collisions, EPL (Europhys. Lett.) 9(4) (1989) 345. [24] P. Lallemand, L.S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (6) (2000) 6546–6562. [25] X. Shan, Analysis and reduction of the spurious current in a class of multiphase lattice Boltzmann models, Phys. Rev. E 73 (4) (2006) 047701.

[26] Q. Li, K.H. Luo, X.J. Li, Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model, Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 87 (5) (2013) 053301. [27] E.M. Sparrow, M. Faghri, Fluid-to-fluid conjugate heat transfer for a vertical pipe—internal forced convection and external natural convection, J. Heat Transf. 102 (3) (1980) 402–407. [28] W.Z. Fang, L. Chen, Q.J. Kang, W.-Q. Tao, Lattice Boltzmann modeling of pool boiling with large liquid-gas density ratio, Int. J. Therm. Sci. 114 (2017) 172– 183. [29] C. Zhang, F. Hong, P. Cheng, Simulation of liquid thin film evaporation and boiling on a heated hydrophilic microstructured surface by Lattice Boltzmann method, Int. J. Heat Mass Transf. 86 (2015) 629–638. [30] L. Zhang, T. Wang, Y. Jiang, S. Kim, C. Guo, A study of boiling on surfaces with temperature-dependent wettability by lattice Boltzmann method, Int. J. Heat Mass Transf. 122 (2018) 775–784. [31] W. Fritz, Maximum volume of vapor bubbles, Physik Zeitschr 36 (1935) 379– 384. [32] S. Gong, P. Cheng, Lattice Boltzmann simulations for surface wettability effects in saturated pool boiling heat transfer, Int. J. Heat Mass Transf. 85 (2015) 635– 646. [33] C. Zhang, P. Cheng, F. Hong, Mesoscale simulation of heater size and subcooling effects on pool boiling under controlled wall heat flux conditions, Int. J. Heat Mass Transf. 101 (2016) 1331–1342. [34] C. Zhang, P. Cheng, Mesoscale simulations of boiling curves and boiling hysteresis under constant wall temperature and constant heat flux conditions, Int. J. Heat Mass Transf. 110 (2017) 319–329. [35] L.H. Chai, X.F. Peng, B.X. Wang, Nucleation site interaction during boiling, Int. J. Heat Mass Transf. 43 (23) (2000) 4249–4258. [36] Y.H. Chang, Y.M. Ferng, Experimental investigation on bubble dynamics and boiling heat transfer for saturated pool boiling and comparison data with previous works, Appl. Therm. Eng. 154 (2019) 284–293. [37] E. Teodori, T. Valente, I. Malavasi, A.S. Moita, M. Marengo, A.L.N. Moreira, Effect of extreme wetting scenarios on pool boiling conditions, Appl. Therm. Eng. 115 (2017) 1424–1437. [38] V.P. Carey, Liquid-vapor phase-change phenomena: an introduction to the thermophysics of vaporization and condensation processes in heat transfer equipment, Hemisphere Pub. Corp., Washington, D.C., 1992. [39] L.J.P.M.S. Rayleigh, VIII. On the pressure developed in a liquid during the collapse of a spherical cavity, 34(200) (1917) 94–98. [40] H. Jo, H.S. Ahn, S. Kang, M.H. Kim, A study of nucleate boiling heat transfer on hydrophilic, hydrophobic and heterogeneous wetting surfaces, Int. J. Heat Mass Transf. 54 (25–26) (2011) 5643–5652.