- Email: [email protected]
Boiling heat transfer on hydrophilic-hydrophobic mixed surfaces: A 3D lattice Boltzmann study
Applied Thermal Engineering 142 (2018) 846–854
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Boiling heat transfer on hydrophilic-hydrophobic mixed surfaces: A 3D lattice Boltzmann study
T
⁎
Y. Yu, Z.X. Wen, Q. Li , P. Zhou, H.J. Yan School of Energy Science and Engineering, Central South University, Changsha 410083, China
H I GH L IG H T S
boiling performance on a type of mixed surfaces is studied using a 3D LB model. • The optimal contact angle of hydrophilic region decreases with increasing the wall superheat. • The hydrophilic region plays an increasingly important role when the wall superheat increases. • The • The orthogonal array tests are performed to identify the most important influential factor.
A R T I C LE I N FO
A B S T R A C T
Keywords: Boiling Hydrophilic-hydrophobic Mixed wettability Lattice Boltzmann model
In this paper, the boiling heat transfer performance on a type of hydrophilic-hydrophobic mixed surfaces is numerically investigated using a three-dimensional (3D) thermal multiphase lattice Boltzmann model with liquid-vapor phase change. The mixed surface is square-pillar textured with the pillar being composed of hydrophilic side walls and a hydrophobic top surface. Numerical simulations are carried out to investigate the influences of the contact angle of the hydrophilic region of the mixed surface, the pillar width, and the pillar height. It is found that an appropriate increase of the contact angle of the hydrophilic region, θphi, can promote the bubble nucleation on the bottom substrate and hence enhances the nucleate boiling on the hydrophilichydrophobic mixed surface, but a larger θphi may make the boiling enter the transition or film boiling regime. Moreover, it is shown that the heat flux initially increases with the increase of the pillar width and shows a declining trend after reaching its peak value. The optimal θphi and the optimal pillar width are found to decrease with the increase of the wall superheat, revealing that the hydrophilic region of the mixed surface should play an increasingly important role when the wall superheat increases. The orthogonal array tests are performed and it is found that the contact angle of the hydrophilic region is the most important influential factor among the three investigated factors.
1. Introduction Boiling is widely utilized in industry and is one of the key technologies in thermal power plants, nuclear reactors, and steel manufacturing [1] as it is the best way to cool a hot body and to retrieve heat from a heat source effectively. Numerous research efforts have been made regarding boiling heat transfer and a lot of enhanced heat transfer techniques have been developed to promote the performance of boiling heat transfer [1–3], among which the surface modification approach such as controlling the surface wettability [4,5] and applying micro/ nano-scale structures [6] has attracted significant attention. At the very dawn of boiling studies, the importance of surface wettability was not recognized since the wettability effect was included
⁎
in the effects of surface characteristics such as roughness and material properties [1]. Through many experimental and numerical studies, the importance of surface wettability on boiling heat transfer has been revealed. Generally, the advantage of a hydrophobic surface lies in that the onset of boiling on a hydrophobic surface requires a lower wall superheat than that on a hydrophilic surface [7,8]. Conversely, the advantage of a hydrophilic surface is that its critical heat flux is much higher than that of a hydrophobic surface [9]. An ideal heating surface may make use of the hydrophobicity to promote boiling onset and the hydrophilicity to achieve a high critical heat flux. Therefore a recent trend in surface modification is to combine the advantages of hydrophobicity and hydrophilicity to develop a variety of mixed-wettability surfaces for enhancing boiling heat transfer
Corresponding author. E-mail address: [email protected] (Q. Li).
https://doi.org/10.1016/j.applthermaleng.2018.07.059 Received 11 June 2018; Received in revised form 2 July 2018; Accepted 3 July 2018 1359-4311/ © 2018 Elsevier Ltd. All rights reserved.
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
[4,10–17]. Betz et al. [10] performed a pioneering study in this field. They prepared two types of flat surfaces with mixed wettability. One is a surface with hydrophobic hexagonal islands on a hydrophilic network and the other is that with hydrophilic hexagonal islands on a hydrophobic network. Both patterns were shown to be capable of enhancing nucleate boiling heat transfer. Jo et al. [4,13] fabricated a heterogeneous wettability surface composed of a hydrophilic substrate with hydrophobic dots and showed that the surface with heterogeneous wettability can provide better nucleate boiling heat transfer than that with homogeneous wettability. Dai et al. [11] proposed a type of hydrophilic-hydrophobic mixed surfaces, which was synthesized from functionalized multiwall carbon nanotubes by introducing hydrophilic functional groups on the pristine multiwall carbon nanotube surfaces, and demonstrated that the mixed surface is superior in enhancing nucleate boiling. Jo et al. [14] suggested a mixed surface comprising a wetting pattern located at the head of surface microstructures and showed that hydrophilic self-assembled monolayer patterns and spatially fabricated micro-post structures can enhance nucleate boiling heat transfer. Kumar et al. [15] have investigated the effects of mixed wettability on the boiling performance of cylindrical copper surfaces. Recently, Li et al. [16] numerically investigated the boiling heat transfer performance on a type of hydrophilic-hydrophobic mixed surfaces, which was textured with pillars consisting of hydrophilic side walls and hydrophobic tops. They showed that the hydrophobicity of the tops of pillars promotes bubble nucleation and reduces the required wall superheat for boiling onset. Moreover, it was found that increasing the contact angle of the tops of pillars leads to a leftward shift of the boiling curve and a leftward and upward shift of the heat transfer coefficient curve. Meanwhile, Shen et al. [17] studied the boiling performance on mixed surfaces with square-pillars in millimeters at low heat fluxes. They showed that, under the same geometric size and heating power conditions, the heat transfer coefficients of mixed surfaces are higher than those of spatially uniform wetting surfaces, regardless of mixed modes. In recent years, the rapid development of computer technology has brought a drastic advancement in numerical simulations of boiling phenomena for the understanding of boiling fundamentals [3]. The molecular dynamics (MD) method is a promising way to grasp an understanding of the fundamental physics of boiling phenomena. However, the MD method is usually limited to systems of very small size [1]. The lattice Boltzmann (LB) method, which is a popular mesoscopic numerical approach for simulating fluid flow and heat transfer [18–22], has been applied to investigate phase-change heat transfer in the past decade [16,23–29]. The LB method can be viewed as a special discrete solver for mimicking the kinetic Boltzmann equation. An important advantage of the LB method for modeling interfacial phenomena is that the interface between different phases can arise, deform, and migrate naturally without using any technique to track or capture the interface [20]. Following a recent study of Li et al. [16], in the present work we aim to investigate the boiling performance on a type of hydrophilic-hydrophobic mixed surfaces using a three-dimensional (3D) thermal multiphase LB model with liquid-vapor phase change. The mixed surface is square-pillar textured with the pillar being composed of hydrophilic side walls and a hydrophobic top surface. The rest of the present paper is organized as follows. A 3D thermal multiphase LB model with liquid-vapor phase change is introduced in Section 2. In Section 3, numerical simulations of boiling heat transfer on the hydrophilic-hydrophobic mixed surface are carried out to investigate the effects of the contact angle of the hydrophilic region of the mixed surface, the pillar width, and the pillar height. Some discussions of the related physical mechanisms are also provided there. A brief summary is given in Section 4.
2. Numerical model In the LB method, the fluid flow is simulated by tracking the evolution of the density distribution function and then the macroscopic averaged properties are obtained by accumulating the density distribution function. The viscous effect is modeled through a linearized collision operator such as the Bhatnagar-Gross-Krook (BGK) collision operator [30,31] and the multiple-relaxation-time (MRT) collision operator [16,26,32–36]. Using the MRT collision operator, the LB equation, which governs the evolution of the density distribution function, can be written as follows:
¯ αβ (fβ −f eq )|(x , t ) + δt (Gα−0.5Λ ¯ αβ Gβ )|(x , t ) , fα (x + eα δt , t + δt ) = fα (x, t )−Λ β (1) where fα is the density distribution function, f αeq is the equilibrium density distribution function, x is the spatial position, eα is the discrete velocity in the α th direction, t is the time, δt is the time step, Gα is the ¯ αβ = (M−1ΛM)αβ is the forcing term in the discrete velocity space, and Λ collision operator, in which M is the transformation matrix and Λ is a diagonal matrix. With the aid of the transformation matrix M, the right-hand side of Eq. (1), i.e., the collision process, can be implemented in the moment space:
Λ m∗ = m−Λ(m−meq) + δt ⎛I− ⎞ S , ⎝ 2⎠
(2)
meq
Mf eq
= where I is the unit matrix, m = Mf , , and S = MG is the forcing term in the moment space. The streaming process is then implemented as follows:
fα (x + eα δt , t + δt ) = f α∗ (x, t ), f∗
M−1m∗
∑ fα ,
ρu =
(3)
M−1
and is the inverse matrix of the transformation where = matrix. The macroscopic density and velocity are calculated through the following relationships:
ρ=
α
∑ eα fα α
+
δt F, 2
(4)
where F is the total force exerted on the system. In the classical MRT-LB method, the transformation matrix M is an orthogonal matrix. Recently, some studies [37–40] have shown that the transformation matrix of an MRT-LB model is not necessary to be an orthogonal one. In the present study, a non-orthogonal 3D MRT-LB model proposed in Ref. [40] is employed, which is constructed based on the D3Q19 lattice model. In comparison with the usual orthogonal 3D MRT-LB model, the non-orthogonal MRT-LB model can retain the numerical accuracy while considerably simplifying the transformation matrix M and its inverse matrix M−1 [40]. The lattice velocities {eα} of the D3Q19 lattice are given by 0 0 ⎤ ⎡0 1 − 1 0 0 0 0 1 − 1 1 − 1 1 − 1 1 − 1 0 0 eα = ⎢ 0 0 0 1 − 1 0 0 1 − 1 − 1 1 0 0 0 0 1 − 1 1 − 1⎥ . 0 0 1 −1 −1 1 1 −1 −1 1 ⎦ ⎣0 0 0 0 0 1 − 1 0 0
(5) The equilibria meq = Mf eq in Eq. (2) are given by
m0eq = ρ , m1eq = ρu x , m2eq = ρu y , m3eq = ρuz , m4eq = ρ + ρ |u|2 , m5eq = ρ (2u x2−u y2−uz2), m6eq = ρ (u y2−uz2), m7eq = ρu x u y , m8eq = ρu x uz , m9eq = ρu y uz , eq eq eq eq eq = ρcs2 u y , m11 = ρcs2 u x , m12 = ρcs2 uz , m13 = ρcs2 u x , m14 = ρcs2 uz , m10 eq = ρcs2 u y , m15 eq eq = φ + ρcs2 (u x2 + u y2), m17 = φ + ρcs2 (u x2 + uz2), m16 eq = φ + ρcs2 (u y2 + uz2), m18
(6) where 847
cs2
= 1 3 and φ =
ρcs4 (1−1.5
|u|2 ) .
The relaxation matrix Λ (the
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
a quantity ϕ , the spatial gradient of ϕ and the Laplacian of ϕ can be calculated by the following second-order isotropic difference schemes, respectively
matrix for relaxation rates) in Eq. (2) is defined as follows:
Λ = diag(1, 1, 1, 1, se, s v, s ν, s ν, s ν, s ν, sq, sq, sq, sq, sq, sq, sπ , sπ , sπ ), (7)
∂i ϕ (x) ≈
where se and s ν determine the bulk and shear viscosities, respectively, and sq and sπ are related to the high-order non-hydrodynamic moments. The forcing term S in Eq. (2) is given by
∇2 ϕ (x) ≈
0 Fx Fy Fz
⎤ ⎡ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 6σ | Fm|2 + 2 F · u ⎢ ψ2δt (se−1 − 0.5) ⎥ ⎥ ⎢ ⎢ 2(2Fx u x −Fy u y−Fz uz ) ⎥ ⎥ ⎢ 2(Fy u y−Fz uz ) ⎥ ⎢ Fx u y + Fy u x ⎥ ⎢ ⎥ ⎢ Fx uz + Fz u x ⎥ ⎢ Fy uz + Fz u y ⎥, S=⎢ ⎥ ⎢ 2 c F y s ⎥ ⎢ 2 ⎥ ⎢ cs Fx ⎥ ⎢ cs2 Fz ⎥ ⎢ ⎥ ⎢ 2 c F s x ⎥ ⎢ ⎥ ⎢ cs2 Fz ⎥ ⎢ 2 cs Fy ⎥ ⎢ ⎥ ⎢ 2 ⎢ 2cs (u x Fx + u y Fy ) ⎥ ⎢ 2c 2 (u F + u F ) ⎥ x x z z s ⎥ ⎢ ⎢ 2cs2 (u y Fy + uz Fz ) ⎥ ⎦ ⎣
2 cs2 δt2
α
∑ ωα [ϕ (x + eα δt )−ϕ (x)], α
(11)
(12)
3. Numerical results and discussions 3.1. The simulation setup and the boiling curve The computational domain is taken as L x × L y × Lz = 150 l.u.×150 l.u.×260 l.u. (l.u. represents lattice units) and a square pillar is located at the center of the bottom substrate, as shown in Fig. 1. The bottom surface and the side walls of the square pillar are hydrophilic (the contact angle is denoted by θphi ) whereas the top surface of the pillar is hydrophobic (the contact angle is denoted by θpho ). The pillar width and height are represented by W and H, respectively. The initial setting of the computational domain is a liquid (0 ⩽ z ⩽ 150 l.u.) below its vapor state and the initial temperature of the computational domain is taken as Ts = 0.8Tc except that the temperature of the heating surface is set to Tb , which leads to the wall superheat ΔT = Tb−Ts . The temperature of the top boundary is maintained at Ts . The Zou-He boundary scheme [48] is employed at the heating surface and the periodic boundary is applied in the x and y directions. The specific heat at constant volume is taken as cV = 10 and the thermal conductivity λ = ρcV χ is chosen to be proportional to the density with cV χ = 0.07 . The gravitational acceleration is set to g = 3 × 10−5 . Our simulations are performed without the temperature solver during the first 1000 time steps. The influence of the contact angle of the top surface of the pillar has been studied in Ref. [16] and we employ θpho ≈ 103∘ throughout the present study. In this section, the pillar width and height are taken as W = 40 l.u. and H = 20 l.u., respectively, and the contact angle θphi is chosen as θphi ≈ 45. 3∘ . Fig. 2 displays some snapshots of the boiling processes on the hydrophilic-hydrophobic mixed surface at ΔT = 0.015, 0.017 , and 0.032 . From the figure we can see that in the case of ΔT = 0.015 the bubbles are just nucleated on the top surface of the pillar, i.e., the hydrophobic region of the mixed surface, and there are no bubbles on the bottom substrate. With the increase of the wall superheat from 0.015 to 0.017, some bubbles are also nucleated on the bottom substrate, as seen in
(8)
∑ wα ψ (x + eα δt ) eα ,
(9)
α
where G is the interaction strength, ψ (x) is the pseudopotential, and wα are the weights. For the D3Q19 lattice, the weights wα in Eq. (9) are given by w1 − 6 = 1 6 and w7 − 18 = 1 12 . The constant σ in the fifth moment of the forcing term S given by Eq. (8) is employed to adjust the mechanical stability condition of the pseudopotential multiphase LB approach so as to achieve the thermodynamic consistency [41–44]. The total force F includes the pseudopotential interaction force Fm and the buoyant force given by Fb = (ρ−ρave ) g , where ρave is the average density in the computational domain and g = (0, 0, −g ) is the gravitational acceleration. The pseudopotential ψ (x) is taken as ψ (x) = 2(pEOS −ρcs2) Gc 2 [41], in which pEOS is a non-ideal equation of state and c is the lattice constant. In the present study, the Peng-Robinson equation of state is utilized [26,45]. By neglecting the viscous heat dissipation, the governing equation for the temperature field is given by [26,46]
∂t T = −u ·∇T +
∑ ωα ϕ (x + eα δt ) eαi,
where cs = c 3 is the lattice sound speed and ωα are given by ω1 − 6 = 1 18 and ω7 − 18 = 1 36. In the present study, the saturation temperature is set to Ts = 0.8Tc (Tc is the critical temperature), which corresponds to the liquid density ρL ≈ 7.2 and the vapor density ρV ≈ 0.197 . All the quantities are taken in lattice units, i.e., the units in the LB method, which are based on the lattice constant c = δx δt = 1.
where F is the total force exerted on the system. The interaction force in the pseudopotential multiphase LB approach is given by [20]:
Fm = −G ψ (x)
1 cs2 δt
1 T ⎛ ∂pEOS ⎞ ∇ ·(λ∇T )− ∇·u, ρcV ρcV ⎝ ∂T ⎠ ρ ⎜
⎟
(10)
where λ is the thermal conductivity and cV is the specific heat at constant volume. In a recent paper, Gong et al. [47] claimed that they established a modified temperature equation for liquid-vapor phase change. Here we emphasize that the study of Gong et al. [47] was based on a misunderstanding of the derivation of local entropy balance law and they employed the relationship of ideal gases (e = c v T ) to derive the temperature equation for non-ideal gases, which is incorrect. The issues associated with the temperature equation for liquid-vapor phase change have been clarified in Ref. [46]. The temperature equation given by Eq. (10) is solved with the classical fourth-order Runge-Kutta scheme for time discretization [26]. The isotropic central schemes are applied for spatial discretization. For
Fig. 1. Schematic illustration of the hydrophilic-hydrophobic mixed surface. The bottom substrate and the sides of the square pillar are hydrophilic, whereas the top surface of the pillar is hydrophobic. 848
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
0.0021
heat flux
0.0018 0.0015 0.0012 0.0009
(a) ΔT = 0.015
0.0006
0.01
0.02
0.03
0.04
wall superheat Fig. 3. Boiling curve of the hydrophilic-hydrophobic mixed surface with θphi ≈ 45. 3∘ and θpho ≈ 103∘ . The pillar width and pillar height are taken as W = 40 l.u. and H = 20 l.u., respectively.
and H = 20 l.u., respectively. The variations of the heat flux with the contact angle θphi in the cases of ΔT = 0.015, 0.018, and 0.021 are presented in Fig. 4, where the heat flux is the time average of the transient heat flux during t ∈ [103δt , 2 × 10 4δt ]. From Fig. 4 we can see that in the case of ΔT = 0.015 the heat flux experiences a slight change when θphi is increased from 31. 9∘ to 50. 8∘. After that, the heat flux rises rapidly, reaching its peak at θphi ≈ 66∘ ∼ 73∘, and then the heat flux decreases. In the cases of ΔT = 0.018 and 0.021, the heat flux initially increases with the increase of θphi . After reaching its peak value, the heat flux basically shows a declining trend. Particularly, it can be found that the peak values of the heat fluxes in the three cases are achieved with different θphi , i.e., θphi ≈ 66∘ ∼ 73∘, 57∘, and 40∘ , respectively, which means that a smaller θphi or a larger wettability difference Δθ = θpho−θphi is required to optimize the heat flux when the wall superheat increases. We take the case of ΔT = 0.018 as an example to illustrate the variations of the heat flux with the contact angle θphi . Figs. 5, 6, 7 and 8 display some snapshots of the boiling processes on the hydrophilichydrophobic mixed surface at ΔT = 0.018 with θphi ≈ 35. 6∘ , 40∘ , 50. 8∘, and 73. 2∘ , respectively. The nucleate boiling can be observed in Fig. 5, where θphi is set to 35. 6∘. However, it can be found that the bubble nucleation in Fig. 5 is completely attributed to the hydrophobicity of the top surface of the pillar and there are no bubbles nucleated on the
(b) ΔT = 0.017
(c) ΔT = 0.032 Fig. 2. Snapshots of the boiling processes at (a) ΔT = 0.015, (b) 0.017 , and (c) 0.032 . From left to right: t = 7000δt , 12000δt , and 18000δt . The contact angles θphi ≈ 45. 3∘ and θpho ≈ 103∘ . The pillar width and pillar height are taken as W = 40 l.u. and H = 20 l.u., respectively.
T = 0.015 T = 0.018 T = 0.021
0.0020
heat flux
Fig. 2(b), which enhance the nucleate boiling. When the wall superheat increases to 0.032, a large portion of the heating surface is covered by vapor and the boiling process enters the transition boiling regime, as shown in Fig. 2(c). The boiling curve is shown in Fig. 3, where the heat flux is the time average of the transient heat flux during t ∈ [103δt , 2 × 10 4δt ], which is in turn the spatial average of the local heat flux over the whole heating surface. Initially, the heat flux increases smoothly as the wall superheat is increased from 0.013 to 0.015. Later, it experiences a sharp rise and reaches the critical heat flux around ΔT = 0.025 ∼ 0.028. After that, the heat flux gradually decreases due to the appearance of transition boiling. Overall, the boiling curve displayed in Fig. 3 is qualitatively in good agreement with the typical boiling curve for water [49].
0.0015
0.0010
0.0005 30
40
50
contact angle 3.2. The effect of the contact angle θphi
60 phi
70
80
(deg.)
Fig. 4. Variations of the heat flux with the contact angle θphi in the cases of ΔT = 0.015, 0.018, and 0.021. The pillar width and height are fixed at W = 40 l.u. and H = 20 l.u., respectively. The contact angle of the top surface of the pillar is taken as θpho ≈ 103∘ .
In this section, numerical simulations are carried out to investigate the influence of the contact angle of the hydrophilic region of the mixed surface, i.e., θphi . The pillar width and height are fixed at W = 40 l.u. 849
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
Fig. 5. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 35. 6∘ . Cross-sectional view at y = 0.5L y . (a) t = 2000δt , (b) 4000δt , (c) 6000δt , and (d) 10000δt .
bottom substrate. When θphi increases from 35. 6∘ to 40∘, some bubbles are nucleated on the bottom substrate, as shown in Fig. 6(b). Subsequently, the bubbles nucleated on the bottom substrate interact with the bubble generated on the top surface of the pillar, leading to a larger bubble, as can be seen in Fig. 6(c). Owing to the bubble nucleation on the bottom substrate as well as the interaction between the bubbles nucleated on the bottom substrate and the bubble on the top of the pillar, the boiling heat transfer is enhanced and therefore the heat flux obtained with θphi ≈ 40∘ is larger than that obtained with θphi ≈ 35. 6∘ in the case of ΔT = 0.018, which can be observed in Fig. 4. Furthermore, by comparing Fig. 7 with Fig. 6, we can find that the nucleate boiling is further enhanced as the contact angle θphi increases from 40∘ to 50. 8∘ since the bubbles nucleated on the bottom substrate appear earlier in Fig. 7. The snapshots shown in Figs. 7(a) and 6(a) are
both taken at t = 2500δt . Nevertheless, it can be seen that small bubbles have been nucleated on the bottom substrate in Fig. 7(a) while there are no bubbles on the bottom substrate in Fig. 6(a). Accordingly, the heat flux obtained with θphi ≈ 50. 8∘ is larger than that obtained with θphi ≈ 40∘ , as seen in Fig. 4. However, when the contact angle θphi is increased to 73. 2∘ , the heat flux decreases considerably because a large amount of vapor is generated on the bottom surface, resulting in the transition boiling, as shown in Fig. 8. In summary, the numerical results shown in Figs. 5–8 reveal that an appropriate increase of the contact angle θphi can promote the bubble nucleation on the bottom substrate and hence enhances the nucleate boiling on the hydrophilic-hydrophobic mixed surface, but a larger θphi may make the boiling enter the transition or film boiling regime. In Fig. 4 it has been shown that the optimal θphi is different for different
Fig. 6. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 40∘ . (Top) Three-dimensional view and (bottom) cross-sectional view at y = 0.5L y . (a) t = 2500δt , (b) 4500δt , (c) 6000δt , and (d) 13000δt . 850
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
Fig. 7. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 50. 8∘ . Cross-sectional view at y = 0.5L y . (a) t = 2500δt , (b) 4500δt , (c) 7000δt , and (d) 11000δt .
0.0018
wall superheats, and decreases with the increase of the wall superheat, which is reasonable since a smaller contact angle can prevent the formation of a continuous vapor film on the heating surface when the vapor generation is enhanced by a higher wall superheat [26].
T = 0.018 T = 0.0165 T = 0.015
heat flux
0.0016
3.3. The effects of H and W and the orthogonal array tests In this section, the effects of the pillar height and width are studied. The contact angles θphi and θpho are chosen as θphi ≈ 50. 8∘ and θpho ≈ 103∘ , respectively. First, we investigate the effect of the pillar width and the pillar height is fixed at H = 20 l.u. Fig. 9 depicts the variations of the heat flux with the pillar width in the cases of ΔT = 0.015, 0.0165, and 0.018. From the figure we can see that in the case of ΔT = 0.015 the heat flux initially increases with the increase of the pillar width and reaches its peak at W = 80 l.u. After that, the heat flux gradually declines. A similar trend can be observed in the case of ΔT = 0.0165 but the peak value of the heat flux is achieved at W = 70 l.u. In the case of ΔT = 0.018, the heat flux rises rapidly as the pillar width increases from W = 20 l.u. to 40 l.u., and then the heat flux shows a downward trend, although there is a slight rise of the heat flux when the pillar width is increased from W = 50 l.u. to 60 l.u. Similar to Fig. 4, Fig. 9 shows that the optimal pillar width varies with the wall superheat. More specifically, it decreases with the increase of the wall superheat. In Ref. [16], two-dimensional (2D) LB simulations have been performed to investigate the effect of the pillar width on the boiling performance of the hydrophilic-hydrophobic mixed surface and it has been shown that a large pillar width may result in a substantial decrease of the heat flux. Such a phenomenon can also
0.0014 0.0012 0.0010 0.0008 20
40
60
80
100
120
pillar width (l.u.) Fig. 9. Variations of the heat flux with the pillar width in the cases of ΔT = 0.015, 0.0165 , and 0.018. The pillar height is fixed at H = 20 l.u.
be observed in Fig. 9 in the cases of ΔT = 0.0165 and 0.018, which is attributed to the fact that, when a very large pillar width is applied, the boiling process may enter the transition or film boiling regime because of the hydrophobicity of the top surface of the pillar. The results that the optimal θphi and the optimal pillar width decrease with the increase of the wall superheat reveal that the hydrophilic region of the mixed surface should play an increasingly important role when the wall superheat increases. As previously mentioned, the
Fig. 8. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 73. 2∘ . Cross-sectional view at y = 0.5L y . (a) t = 2000δt , (b) 4000δt , (c) 7000δt , and (d) 11000δt . 851
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
advantage of hydrophobicity is promoting boiling onset at a low wall superheat, while the advantage of hydrophilicity lies in that it provides a relatively high critical heat flux. Hence, with the increase of the wall superheat, the hydrophilic region gradually becomes important and decreasing the pillar width can reduce the influence of the hydrophobic region (the top surface of the pillar). Similarly, a smaller θphi can enhance the influence of the hydrophilic region when the wall superheat increases. The differences between the present 3D simulations and the previous 2D simulations [16] should be noticed. In 2D simulations, increasing the pillar width does not affect the bubble nucleation at the two corners formed by the bottom substrate and the pillar. However, in 3D simulations, increasing the pillar width will influence the bubble nucleation on the four edges formed by the bottom substrate and the square pillar, as shown in Fig. 10, where some snapshots of the boiling processes at ΔT = 0.0165 are displayed. By comparing Fig. 10(b) with Fig. 10(a), we can see that the bubble nucleation on the bottom substrate is affected when the pillar width changes from 40 l.u. to 70 l.u. For example, some small bubbles nucleated on the bottom substrate can be observed in Fig. 10(b) at t = 8000δt and 11000δt , whereas there are no similar bubbles on the bottom substrate in Fig. 10(a). This is also the reason why the heat flux obtained with W = 70 l.u. is larger than that obtained with W = 40 l.u. in the case of ΔT = 0.0165. The effect of the pillar height is shown in Fig. 11, where the variations of the heat flux with the pillar height in the cases of ΔT = 0.015, 0.0165, and 0.018 are presented. The pillar width is fixed at W = 40 l.u. The figure illustrates that in the case of ΔT = 0.015 the heat flux rises smoothly as the pillar height increases from H = 20 l.u. to 50 l.u. In the other two cases, the heat flux fluctuates with the increase of the pillar height, but the fluctuation is very small. Overall, the influence of the pillar height is relatively small in comparison with the effects of the
0.0018
T = 0.018 T = 0.0165 T = 0.015
heat flux
0.0016 0.0014 0.0012 0.0010 0.0008 20
30
40
50
pillar height (l.u.) Fig. 11. Variations of the heat flux with the pillar height in the cases of ΔT = 0.015, 0.0165 , and 0.018. The pillar width is fixed at W = 40 l.u.
contact angle θphi and the pillar width. Increasing the pillar height does not change the area of the top surface of the pillar, but it reduces the ratio of the hydrophobic region to the whole heating surface, i.e., Spho Swhole , which varies from 6.23% to 5.25% when the pillar height is increased from H = 20 l.u. to 50 l.u. It can be readily found that such a change is much smaller than that in the previous 2D simulations [16], in which the ratio Spho Swhole degenerates to L pho L whole . This may be the reason why the heat flux reported in the 2D simulations considerably changes with the pillar height [16] whereas in the present 3D simulations the variations of the heat flux with the pillar height are relatively small.
Fig. 10. Snapshots of the boiling process at ΔT = 0.0165 with θphi ≈ 50. 8∘ and θpho ≈ 103∘ . (a) W = 40 l.u. and (b) 70 l.u. From left to right: t = 3500δt , 6000δt , 8000δt , and 11000δt . 852
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
Table 1 The orthogonal array tests and the result analsyis. Test serial number
Heat flux
Factors A/θphi
B/pillar width
A× B
C/pillar height
A× C
Kj1 × 103
1 (40°) 1 (40°) 1 (40°) 1 (40°) 2 (57°) 2 (57°) 2 (57°) 2 (57°) 3 (66°) 3 (66°) 3 (66°) 3 (66°) 4 (73.2°) 4 (73.2°) 4 (73.2°) 4 (73.2°) 4.9361
1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 4.6191
1 2 3 4 2 1 4 3 3 4 1 2 4 3 2 1 4.8702
1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (50 l.u.) 3 (40 l.u.) 4 (50 l.u.) 1 (20 l.u.) 2 (30 l.u.) 4 (50 l.u.) 3 (40 l.u.) 2 (30 l.u.) 1 (20 l.u.) 2 (30 l.u.) 1 (20 l.u.) 4 (50 l.u.) 3 (40 l.u.) 4.4447
1 2 3 4 4 3 2 1 2 1 4 3 3 4 1 2 5.0041
Kj2 × 103
5.6491
4.7152
4.4585
4.9961
4.9611
Kj3 × 103
4.7576
5.2888
4.4322
4.8724
4.7815
Kj4 × 103
3.8650
4.5847
5.4469
4.8946
4.4611
K¯ j1 × 103
1.2340
1.1548
1.2176
1.1112
1.2510
K¯ j2 × 103
1.4123
1.1788
1.1146
1.2490
1.2403
K¯ j3 × 103
1.1894
1.3222
1.1081
1.2181
1.1954
K¯ j4 × 103
0.9663
1.1462
1.3617
1.2237
1.1153
Range × 103 Rank Optimization
0.4460
0.1760
0.2536
0.1379
0.1358
1 A2
3 B3
2 (A × B)4
4 C2
5 (A × C)1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0.0011666 0.0012148 0.0012758 0.0012789 0.0012131 0.0014257 0.0016457 0.0013646 0.0011142 0.0013971 0.0012915 0.0009548 0.0011252 0.0006776 0.0010758 0.0009864
contact angle θphi and the pillar width has an important influence. Similarly, the “Range” of the interaction A× C is almost the same as that of factor C. Hence, the interaction effects cannot be neglected. As a result, the optimal combination cannot be directly determined by the largest K¯2i ( j = 2 corresponds to factor B) and the largest K¯ 4i ( j = 4 corresponds to factor C). According to the analysis of orthogonal array tests with two-factor interfactions, when the interaction effects have an important influence, a comparison of the results (the heat fluxes in the present study) should be made so as to find the optimal combination. According to the results displayed in Table 1, the optimal combination is A2B3C1.
Finally, the orthogonal array tests are performed to identify the most important factor among the aforementioned three influential factors: the contact angle θphi , the pillar width, and the pillar height. The wall superheat is taken as ΔT = 0.018 and we consider four levels for each factor, namely, 40∘, 57∘, 66∘, and 73. 2∘ for θphi , 20 l.u., 30 l.u., 40 l.u., and 80 l.u. for the pillar width, and 20 l.u., 30 l.u., 40 l.u., and 50 l.u. for the pillar height. Besides, two two-factor interaction effects are also considered, i.e., the interaction effect between θphi and the pillar width, and the interaction effect between θphi and the pillar height. Hence the orthogonal array tests include five factors and four levels shown in Table 1, i.e., L16(45), containing 16 tests in total. In the table, A× B and A× C represent the interaction effect between θphi and the pillar width, and the interaction effect between θphi and the pillar height, respectively. The results of the orthogonal array tests are shown in Table 1, where Kji denotes the sum of the heat fluxes corresponding to the i-th level of the j-th factor, K¯ ji is the average value of Kji , and the “Range” is the difference between the maximum and the minimum of the four average values (K¯ j1, K¯ j2 , K¯ j3 , and K¯ j4 ) in each group. From the table we can see that the “Range” of factor A, i.e., the contact angle θphi , is much larger than the others, which means that θphi is the most important influential factor among the three investigated factors. Meanwhile, it can be found that, the “Range” of factor C (the pillar height) and the “Range” of the interaction A× C are much smaller than the others, which indicates that the pillar height is the least influential one. These findings are consistent with those revealed by Figs. 4, 9, and 11. In addition, it can be found that K¯12 ( j = 1 corresponds to factor A) is the largest of the four K¯1i (i = 1, 2, 3, and 4). Hence, the level A2, i.e., θphi ≈ 57∘ , is the best of the four levels for factor A when the wall superheat is taken as ΔT = 0.018. According to Table 1, if the interaction effects are not taken into account, the optimal combination is A2B3C2. However, the table shows that the “Range” of the interaction A× B is larger than that of factor B (the pillar width), which implies that the interaction effect between the
4. Conclusions In this paper, we have investigated the boiling performance on a type of hydrophilic-hydrophobic mixed surfaces using a 3D thermal LB model with liquid-vapor phase change, which consists of a 3D pseudopotential multiphase LB model for simulating the flow field and a finite-difference solver for simulating the temperature field. The mixed surface is square-pillar textured with the pillar being composed of hydrophilic side walls and a hydrophobic top surface. The effects of the contact angle θphi , the pillar width, and the pillar height have been investigated. The main findings and conclusions are summarized as follows. (i) An appropriate increase of the contact angle θphi can promote the bubble nucleation on the bottom substrate and enhances the nucleate boiling on the hydrophilic-hydrophobic mixed surface, but a larger θphi may make the boiling enter the transition or film boiling regime. (ii) Different from the previous 2D simulations, the present 3D simulations showed that the heat flux initially increases with the increase of the pillar width and basically shows a declining trend after reaching its peak value. 853
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
[19] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. [20] Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combust. Sci. 52 (2016) 62–105. [21] A. Xu, W. Shyy, T. Zhao, Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries, Acta Mech. Sin. 33 (2017) 555–574. [22] Y. Shi, G.H. Tang, H.H. Xia, Lattice Boltzmann simulation of droplet formation in Tjunction and flow focusing devices, Comput. Fluids 90 (2014) 155–163. [23] G. Hazi, A. Markus, On the bubble departure diameter and release frequency based on numerical simulation results, Int. J. Heat Mass Transf. 52 (2009) 1472–1480. [24] A. Márkus, G. Házi, On pool boiling at microscale level: the effect of a cavity and heat conduction in the heated wall, Nucl. Eng. Des. 248 (2012) 238–247. [25] L. Biferale, P. Perlekar, M. Sbragaglia, F. Toschi, Convection in multiphase fluid flows using lattice Boltzmann methods, Phys. Rev. Lett. 108 (2012) 104502. [26] 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. [27] Q. Li, P. Zhou, H.J. Yan, Pinning-depinning mechanism of the contact line during evaporation on chemically patterned surfaces: a lattice Boltzmann study, Langmuir 32 (2016) 9389–9396. [28] Y. Yu, Q. Li, C.Q. Zhou, P. Zhou, H.J. Yan, Y. Yu, Q. Li, C.Q. Zhou, P. Zhou, H.J. Yan, Investigation of droplet evaporation on heterogeneous surfaces using a three-dimensional thermal multiphase lattice Boltzmann model, Appl. Therm. Eng. 127 (2017) 1346–1354. [29] Q. Li, Q.J. Kang, M.M. Francois, A.J. Hu, Lattice Boltzmann modeling of self-propelled Leidenfrost droplets on ratchet surfaces, Soft Matter 12 (2016) 302–312. [30] Y.H. Qian, D. d'Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479–484. [31] S. Chen, H. Chen, D. Martnez, W. Matthaeus, Lattice Boltzmann model for simulation of magnetohydrodynamics, Phys. Rev. Lett. 67 (1991) 3776–3779. [32] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546–6562. [33] D. d'Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, L.S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. Ser. A 360 (2002) 437–451. [34] K.N. Premnath, J. Abraham, Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow, J. Comput. Phys. 224 (2007) 539–559. [35] D. Sun, M. Zhu, J. Wang, B. Sun, Lattice Boltzmann modeling of bubble formation and dendritic growth in solidification of binary alloys, Int. J. Heat Mass Transf. 94 (2016) 474–487. [36] A. Xu, L. Shi, T.S. Zhao, Accelerated lattice Boltzmann simulation using GPU and OpenACC with data management, Int. J. Heat Mass Transf. 109 (2017) 577–588. [37] K.N. Premnath, S. Banerjee, Incorporating forcing terms in cascaded lattice Boltzmann approach by method of central moments, Phys. Rev. E 80 (2009) 036702. [38] D. Lycett-Brown, K.H. Luo, Multiphase cascaded lattice Boltzmann method, Comput. Math. Appl. 67 (2014) 350–362. [39] Q. Liu, Y.L. He, D. Li, Q. Li, Non-orthogonal multiple-relaxation-time lattice Boltzmann method for incompressible thermal flows, Int. J. Heat Mass Transf. 102 (2016) 1334–1344. [40] Q. Li, D.H. Du, L.L. Fei, K.H. Luo, Y. Yu, Three-dimensional non-orthogonal multiple-relaxation-time lattice Boltzmann model for multiphase flows, arXiv preprint, (2018) arXiv:1805.08643. [41] 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 87 (2013) 053301. [42] D. Zhang, K. Papadikis, S. Gu, Three-dimensional multi-relaxation time latticeBoltzmann model for the drop impact on a dry surface at large density ratio, Int. J. Multiph. Flow 64 (2014) 11–18. [43] A. Xu, T.S. Zhao, L. An, L. Shi, A three-dimensional pseudo-potential-based lattice Boltzmann model for multiphase flows with large density ratio and variable surface tension, Int. J. Heat Fluid Flow 56 (2015) 261–271. [44] Y. Shi, G.H. Tang, H.H. Xia, Investigation of coalescence-induced droplet jumping on superhydrophobic surfaces and liquid condensate adhesion on slit and plain fins, Int. J. Heat Mass Transf. 88 (2015) 445–455. [45] P. Yuan, L. Schaefer, Equations of state in a lattice Boltzmann model, Phys. Fluids 18 (2006) 042101. [46] Q. Li, P. Zhou, H.J. Yan, Improved thermal lattice Boltzmann model for simulation of liquid-vapor phase change, Phys. Rev. E 96 (2017) 063303. [47] W. Gong, Y.Y. Yan, S. Chen, E. Wright, A modified phase change pseudopotential lattice Boltzmann model, Int. J. Heat Mass Transf. 125 (2018) 323–329. [48] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591–1598. [49] T.L. Bergman, A.S. Lavine, F.P. Incropera, D.P. Dewitt, Fundamentals of Heat and Mass Transfer, JohnWiley & Sons, New York, 2011.
(iii) The optimal θphi and the optimal pillar width change with the wall superheat. Particularly, both of them decrease with the increase of the wall superheat, which reveals that the hydrophilic region of the mixed surface should play an increasingly important role when the wall superheat increases. (iv) The orthogonal array tests showed that the contact angle θphi is the most important influential factor among the three investigated factors and the pillar height is the least influential one. Moreover, it is found that the interaction effect between the contact angle θphi and the pillar width has an important influence. In future work, we will consider the mixed surfaces with an array of square pillars, through which the effect of the width between two pillars can be investigated. Such a study requires a much larger computational domain than that employed in the present work and also much more computational cost. Acknowledgments This work was supported by the National Natural Science Foundation of China (No. 51506227) and the Foundation for the Author of National Excellent Doctoral Dissertation of China (No. 201439). References [1] Y. Koizumi, M. Shoji, M. Monde, Y. Takata, N. Nagai, Boiling: Research and Advances, Elsevier, 2017. [2] V.K. Dhir, Boiling heat transfer, Annu. Rev. Fluid Mech. 30 (1998) 365–401. [3] V.K. Dhir, G.R. Warrier, E. Aktinol, Numerical simulation of pool boiling: a review, J. Heat Transfer 135 (2013) 061502. [4] H.J. Jo, H.S. Ahn, S.H. Kang, M.H. Kim, A study of nucleate boiling heat transfer on hydrophilic, hydrophobic and heterogeneous wetting surfaces, Int. J. Heat Mass Transf. 54 (2011) 5643–5652. [5] G. Lu, X.D. Wang, Y.Y. Duan, A critical review of dynamic wetting by complex fluids: from Newtonian fluids to non-Newtonian fluids and nanofluids, Adv. Colloid Interface Sci. 236 (2016) 43–62. [6] D.E. Kim, D.I. Yu, D.W. Jerng, M.H. Kim, H.S. Ahn, Review of boiling heat transfer enhancement on micro/nanostructured surfaces, Exp. Therm Fluid Sci. 66 (2015) 173–196. [7] B. Bourdon, R. Rioboo, M. Marengo, E. Gosselin, J.D. Coninck, Influence of the wettability on the boiling onset, Langmuir 28 (2012) 1618–1624. [8] H.J. Jo, M. Kaviany, S.H. Kim, M.H. Kim, Heterogeneous bubble nucleation on ideally-smooth horizontal heated surface, Int. J. Heat Mass Transf. 71 (2014) 149–157. [9] S.G. Kandlikar, A theoretical model to predict pool boiling CHF incorporating effects of contact angle and orientation, J. Heat Transfer 123 (2000) 1071–1079. [10] A.R. Betz, J. Xu, H. Qiu, D. Attinger, Do surfaces with mixed hydrophilic and hydrophobic areas enhance pool boiling? Appl. Phys. Lett. 97 (2010) 103115. [11] X. Dai, X. Huang, F. Yang, X. Li, J. Sightler, Y. Yang, C. Li, Enhanced nucleate boiling on horizontal hydrophobic-hydrophilic carbon nanotube coatings, Appl. Phys. Lett. 102 (2013) 161605. [12] B.J. Suroto, Effects of hydrophobic-spot periphery and subcooling on nucleate pool boiling from a mixed-wettability surface, J. Therm. Sci. Technol. 8 (2013) 294–308. [13] H.J. Jo, S.H. Kim, H.S. Park, M.H. Kim, Critical heat flux and nucleate boiling on several heterogeneous wetting surfaces: controlled hydrophobic patterns on a hydrophilic substrate, Int. J. Multiph. Flow 62 (2014) 101–109. [14] H.J. Jo, D.I. Yu, H. Noh, H.S. Park, M.H. Kim, Boiling on spatially controlled heterogeneous surfaces: wettability patterns on microstructures, Appl. Phys. Lett. 106 (2015) 181602. [15] C.S.S. Kumar, Y.W. Chang, P.H. Chen, Effect of heterogeneous wettable structures on pool boiling performance of cylindrical copper surfaces, Appl. Therm. Eng. 127 (2017) 1184–1193. [16] 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. [17] C. Shen, C. Zhang, Y. Bao, X. Wang, Y. Liu, L. Ren, Experimental investigation on enhancement of nucleate pool boiling heat transfer using hybrid wetting pillar surface at low heat fluxes, Int. J. Therm. Sci. 130 (2018) 47–58. [18] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364.
854
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Boiling heat transfer on hydrophilic-hydrophobic mixed surfaces: A 3D lattice Boltzmann study
T
⁎
Y. Yu, Z.X. Wen, Q. Li , P. Zhou, H.J. Yan School of Energy Science and Engineering, Central South University, Changsha 410083, China
H I GH L IG H T S
boiling performance on a type of mixed surfaces is studied using a 3D LB model. • The optimal contact angle of hydrophilic region decreases with increasing the wall superheat. • The hydrophilic region plays an increasingly important role when the wall superheat increases. • The • The orthogonal array tests are performed to identify the most important influential factor.
A R T I C LE I N FO
A B S T R A C T
Keywords: Boiling Hydrophilic-hydrophobic Mixed wettability Lattice Boltzmann model
In this paper, the boiling heat transfer performance on a type of hydrophilic-hydrophobic mixed surfaces is numerically investigated using a three-dimensional (3D) thermal multiphase lattice Boltzmann model with liquid-vapor phase change. The mixed surface is square-pillar textured with the pillar being composed of hydrophilic side walls and a hydrophobic top surface. Numerical simulations are carried out to investigate the influences of the contact angle of the hydrophilic region of the mixed surface, the pillar width, and the pillar height. It is found that an appropriate increase of the contact angle of the hydrophilic region, θphi, can promote the bubble nucleation on the bottom substrate and hence enhances the nucleate boiling on the hydrophilichydrophobic mixed surface, but a larger θphi may make the boiling enter the transition or film boiling regime. Moreover, it is shown that the heat flux initially increases with the increase of the pillar width and shows a declining trend after reaching its peak value. The optimal θphi and the optimal pillar width are found to decrease with the increase of the wall superheat, revealing that the hydrophilic region of the mixed surface should play an increasingly important role when the wall superheat increases. The orthogonal array tests are performed and it is found that the contact angle of the hydrophilic region is the most important influential factor among the three investigated factors.
1. Introduction Boiling is widely utilized in industry and is one of the key technologies in thermal power plants, nuclear reactors, and steel manufacturing [1] as it is the best way to cool a hot body and to retrieve heat from a heat source effectively. Numerous research efforts have been made regarding boiling heat transfer and a lot of enhanced heat transfer techniques have been developed to promote the performance of boiling heat transfer [1–3], among which the surface modification approach such as controlling the surface wettability [4,5] and applying micro/ nano-scale structures [6] has attracted significant attention. At the very dawn of boiling studies, the importance of surface wettability was not recognized since the wettability effect was included
⁎
in the effects of surface characteristics such as roughness and material properties [1]. Through many experimental and numerical studies, the importance of surface wettability on boiling heat transfer has been revealed. Generally, the advantage of a hydrophobic surface lies in that the onset of boiling on a hydrophobic surface requires a lower wall superheat than that on a hydrophilic surface [7,8]. Conversely, the advantage of a hydrophilic surface is that its critical heat flux is much higher than that of a hydrophobic surface [9]. An ideal heating surface may make use of the hydrophobicity to promote boiling onset and the hydrophilicity to achieve a high critical heat flux. Therefore a recent trend in surface modification is to combine the advantages of hydrophobicity and hydrophilicity to develop a variety of mixed-wettability surfaces for enhancing boiling heat transfer
Corresponding author. E-mail address: [email protected] (Q. Li).
https://doi.org/10.1016/j.applthermaleng.2018.07.059 Received 11 June 2018; Received in revised form 2 July 2018; Accepted 3 July 2018 1359-4311/ © 2018 Elsevier Ltd. All rights reserved.
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
[4,10–17]. Betz et al. [10] performed a pioneering study in this field. They prepared two types of flat surfaces with mixed wettability. One is a surface with hydrophobic hexagonal islands on a hydrophilic network and the other is that with hydrophilic hexagonal islands on a hydrophobic network. Both patterns were shown to be capable of enhancing nucleate boiling heat transfer. Jo et al. [4,13] fabricated a heterogeneous wettability surface composed of a hydrophilic substrate with hydrophobic dots and showed that the surface with heterogeneous wettability can provide better nucleate boiling heat transfer than that with homogeneous wettability. Dai et al. [11] proposed a type of hydrophilic-hydrophobic mixed surfaces, which was synthesized from functionalized multiwall carbon nanotubes by introducing hydrophilic functional groups on the pristine multiwall carbon nanotube surfaces, and demonstrated that the mixed surface is superior in enhancing nucleate boiling. Jo et al. [14] suggested a mixed surface comprising a wetting pattern located at the head of surface microstructures and showed that hydrophilic self-assembled monolayer patterns and spatially fabricated micro-post structures can enhance nucleate boiling heat transfer. Kumar et al. [15] have investigated the effects of mixed wettability on the boiling performance of cylindrical copper surfaces. Recently, Li et al. [16] numerically investigated the boiling heat transfer performance on a type of hydrophilic-hydrophobic mixed surfaces, which was textured with pillars consisting of hydrophilic side walls and hydrophobic tops. They showed that the hydrophobicity of the tops of pillars promotes bubble nucleation and reduces the required wall superheat for boiling onset. Moreover, it was found that increasing the contact angle of the tops of pillars leads to a leftward shift of the boiling curve and a leftward and upward shift of the heat transfer coefficient curve. Meanwhile, Shen et al. [17] studied the boiling performance on mixed surfaces with square-pillars in millimeters at low heat fluxes. They showed that, under the same geometric size and heating power conditions, the heat transfer coefficients of mixed surfaces are higher than those of spatially uniform wetting surfaces, regardless of mixed modes. In recent years, the rapid development of computer technology has brought a drastic advancement in numerical simulations of boiling phenomena for the understanding of boiling fundamentals [3]. The molecular dynamics (MD) method is a promising way to grasp an understanding of the fundamental physics of boiling phenomena. However, the MD method is usually limited to systems of very small size [1]. The lattice Boltzmann (LB) method, which is a popular mesoscopic numerical approach for simulating fluid flow and heat transfer [18–22], has been applied to investigate phase-change heat transfer in the past decade [16,23–29]. The LB method can be viewed as a special discrete solver for mimicking the kinetic Boltzmann equation. An important advantage of the LB method for modeling interfacial phenomena is that the interface between different phases can arise, deform, and migrate naturally without using any technique to track or capture the interface [20]. Following a recent study of Li et al. [16], in the present work we aim to investigate the boiling performance on a type of hydrophilic-hydrophobic mixed surfaces using a three-dimensional (3D) thermal multiphase LB model with liquid-vapor phase change. The mixed surface is square-pillar textured with the pillar being composed of hydrophilic side walls and a hydrophobic top surface. The rest of the present paper is organized as follows. A 3D thermal multiphase LB model with liquid-vapor phase change is introduced in Section 2. In Section 3, numerical simulations of boiling heat transfer on the hydrophilic-hydrophobic mixed surface are carried out to investigate the effects of the contact angle of the hydrophilic region of the mixed surface, the pillar width, and the pillar height. Some discussions of the related physical mechanisms are also provided there. A brief summary is given in Section 4.
2. Numerical model In the LB method, the fluid flow is simulated by tracking the evolution of the density distribution function and then the macroscopic averaged properties are obtained by accumulating the density distribution function. The viscous effect is modeled through a linearized collision operator such as the Bhatnagar-Gross-Krook (BGK) collision operator [30,31] and the multiple-relaxation-time (MRT) collision operator [16,26,32–36]. Using the MRT collision operator, the LB equation, which governs the evolution of the density distribution function, can be written as follows:
¯ αβ (fβ −f eq )|(x , t ) + δt (Gα−0.5Λ ¯ αβ Gβ )|(x , t ) , fα (x + eα δt , t + δt ) = fα (x, t )−Λ β (1) where fα is the density distribution function, f αeq is the equilibrium density distribution function, x is the spatial position, eα is the discrete velocity in the α th direction, t is the time, δt is the time step, Gα is the ¯ αβ = (M−1ΛM)αβ is the forcing term in the discrete velocity space, and Λ collision operator, in which M is the transformation matrix and Λ is a diagonal matrix. With the aid of the transformation matrix M, the right-hand side of Eq. (1), i.e., the collision process, can be implemented in the moment space:
Λ m∗ = m−Λ(m−meq) + δt ⎛I− ⎞ S , ⎝ 2⎠
(2)
meq
Mf eq
= where I is the unit matrix, m = Mf , , and S = MG is the forcing term in the moment space. The streaming process is then implemented as follows:
fα (x + eα δt , t + δt ) = f α∗ (x, t ), f∗
M−1m∗
∑ fα ,
ρu =
(3)
M−1
and is the inverse matrix of the transformation where = matrix. The macroscopic density and velocity are calculated through the following relationships:
ρ=
α
∑ eα fα α
+
δt F, 2
(4)
where F is the total force exerted on the system. In the classical MRT-LB method, the transformation matrix M is an orthogonal matrix. Recently, some studies [37–40] have shown that the transformation matrix of an MRT-LB model is not necessary to be an orthogonal one. In the present study, a non-orthogonal 3D MRT-LB model proposed in Ref. [40] is employed, which is constructed based on the D3Q19 lattice model. In comparison with the usual orthogonal 3D MRT-LB model, the non-orthogonal MRT-LB model can retain the numerical accuracy while considerably simplifying the transformation matrix M and its inverse matrix M−1 [40]. The lattice velocities {eα} of the D3Q19 lattice are given by 0 0 ⎤ ⎡0 1 − 1 0 0 0 0 1 − 1 1 − 1 1 − 1 1 − 1 0 0 eα = ⎢ 0 0 0 1 − 1 0 0 1 − 1 − 1 1 0 0 0 0 1 − 1 1 − 1⎥ . 0 0 1 −1 −1 1 1 −1 −1 1 ⎦ ⎣0 0 0 0 0 1 − 1 0 0
(5) The equilibria meq = Mf eq in Eq. (2) are given by
m0eq = ρ , m1eq = ρu x , m2eq = ρu y , m3eq = ρuz , m4eq = ρ + ρ |u|2 , m5eq = ρ (2u x2−u y2−uz2), m6eq = ρ (u y2−uz2), m7eq = ρu x u y , m8eq = ρu x uz , m9eq = ρu y uz , eq eq eq eq eq = ρcs2 u y , m11 = ρcs2 u x , m12 = ρcs2 uz , m13 = ρcs2 u x , m14 = ρcs2 uz , m10 eq = ρcs2 u y , m15 eq eq = φ + ρcs2 (u x2 + u y2), m17 = φ + ρcs2 (u x2 + uz2), m16 eq = φ + ρcs2 (u y2 + uz2), m18
(6) where 847
cs2
= 1 3 and φ =
ρcs4 (1−1.5
|u|2 ) .
The relaxation matrix Λ (the
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
a quantity ϕ , the spatial gradient of ϕ and the Laplacian of ϕ can be calculated by the following second-order isotropic difference schemes, respectively
matrix for relaxation rates) in Eq. (2) is defined as follows:
Λ = diag(1, 1, 1, 1, se, s v, s ν, s ν, s ν, s ν, sq, sq, sq, sq, sq, sq, sπ , sπ , sπ ), (7)
∂i ϕ (x) ≈
where se and s ν determine the bulk and shear viscosities, respectively, and sq and sπ are related to the high-order non-hydrodynamic moments. The forcing term S in Eq. (2) is given by
∇2 ϕ (x) ≈
0 Fx Fy Fz
⎤ ⎡ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 6σ | Fm|2 + 2 F · u ⎢ ψ2δt (se−1 − 0.5) ⎥ ⎥ ⎢ ⎢ 2(2Fx u x −Fy u y−Fz uz ) ⎥ ⎥ ⎢ 2(Fy u y−Fz uz ) ⎥ ⎢ Fx u y + Fy u x ⎥ ⎢ ⎥ ⎢ Fx uz + Fz u x ⎥ ⎢ Fy uz + Fz u y ⎥, S=⎢ ⎥ ⎢ 2 c F y s ⎥ ⎢ 2 ⎥ ⎢ cs Fx ⎥ ⎢ cs2 Fz ⎥ ⎢ ⎥ ⎢ 2 c F s x ⎥ ⎢ ⎥ ⎢ cs2 Fz ⎥ ⎢ 2 cs Fy ⎥ ⎢ ⎥ ⎢ 2 ⎢ 2cs (u x Fx + u y Fy ) ⎥ ⎢ 2c 2 (u F + u F ) ⎥ x x z z s ⎥ ⎢ ⎢ 2cs2 (u y Fy + uz Fz ) ⎥ ⎦ ⎣
2 cs2 δt2
α
∑ ωα [ϕ (x + eα δt )−ϕ (x)], α
(11)
(12)
3. Numerical results and discussions 3.1. The simulation setup and the boiling curve The computational domain is taken as L x × L y × Lz = 150 l.u.×150 l.u.×260 l.u. (l.u. represents lattice units) and a square pillar is located at the center of the bottom substrate, as shown in Fig. 1. The bottom surface and the side walls of the square pillar are hydrophilic (the contact angle is denoted by θphi ) whereas the top surface of the pillar is hydrophobic (the contact angle is denoted by θpho ). The pillar width and height are represented by W and H, respectively. The initial setting of the computational domain is a liquid (0 ⩽ z ⩽ 150 l.u.) below its vapor state and the initial temperature of the computational domain is taken as Ts = 0.8Tc except that the temperature of the heating surface is set to Tb , which leads to the wall superheat ΔT = Tb−Ts . The temperature of the top boundary is maintained at Ts . The Zou-He boundary scheme [48] is employed at the heating surface and the periodic boundary is applied in the x and y directions. The specific heat at constant volume is taken as cV = 10 and the thermal conductivity λ = ρcV χ is chosen to be proportional to the density with cV χ = 0.07 . The gravitational acceleration is set to g = 3 × 10−5 . Our simulations are performed without the temperature solver during the first 1000 time steps. The influence of the contact angle of the top surface of the pillar has been studied in Ref. [16] and we employ θpho ≈ 103∘ throughout the present study. In this section, the pillar width and height are taken as W = 40 l.u. and H = 20 l.u., respectively, and the contact angle θphi is chosen as θphi ≈ 45. 3∘ . Fig. 2 displays some snapshots of the boiling processes on the hydrophilic-hydrophobic mixed surface at ΔT = 0.015, 0.017 , and 0.032 . From the figure we can see that in the case of ΔT = 0.015 the bubbles are just nucleated on the top surface of the pillar, i.e., the hydrophobic region of the mixed surface, and there are no bubbles on the bottom substrate. With the increase of the wall superheat from 0.015 to 0.017, some bubbles are also nucleated on the bottom substrate, as seen in
(8)
∑ wα ψ (x + eα δt ) eα ,
(9)
α
where G is the interaction strength, ψ (x) is the pseudopotential, and wα are the weights. For the D3Q19 lattice, the weights wα in Eq. (9) are given by w1 − 6 = 1 6 and w7 − 18 = 1 12 . The constant σ in the fifth moment of the forcing term S given by Eq. (8) is employed to adjust the mechanical stability condition of the pseudopotential multiphase LB approach so as to achieve the thermodynamic consistency [41–44]. The total force F includes the pseudopotential interaction force Fm and the buoyant force given by Fb = (ρ−ρave ) g , where ρave is the average density in the computational domain and g = (0, 0, −g ) is the gravitational acceleration. The pseudopotential ψ (x) is taken as ψ (x) = 2(pEOS −ρcs2) Gc 2 [41], in which pEOS is a non-ideal equation of state and c is the lattice constant. In the present study, the Peng-Robinson equation of state is utilized [26,45]. By neglecting the viscous heat dissipation, the governing equation for the temperature field is given by [26,46]
∂t T = −u ·∇T +
∑ ωα ϕ (x + eα δt ) eαi,
where cs = c 3 is the lattice sound speed and ωα are given by ω1 − 6 = 1 18 and ω7 − 18 = 1 36. In the present study, the saturation temperature is set to Ts = 0.8Tc (Tc is the critical temperature), which corresponds to the liquid density ρL ≈ 7.2 and the vapor density ρV ≈ 0.197 . All the quantities are taken in lattice units, i.e., the units in the LB method, which are based on the lattice constant c = δx δt = 1.
where F is the total force exerted on the system. The interaction force in the pseudopotential multiphase LB approach is given by [20]:
Fm = −G ψ (x)
1 cs2 δt
1 T ⎛ ∂pEOS ⎞ ∇ ·(λ∇T )− ∇·u, ρcV ρcV ⎝ ∂T ⎠ ρ ⎜
⎟
(10)
where λ is the thermal conductivity and cV is the specific heat at constant volume. In a recent paper, Gong et al. [47] claimed that they established a modified temperature equation for liquid-vapor phase change. Here we emphasize that the study of Gong et al. [47] was based on a misunderstanding of the derivation of local entropy balance law and they employed the relationship of ideal gases (e = c v T ) to derive the temperature equation for non-ideal gases, which is incorrect. The issues associated with the temperature equation for liquid-vapor phase change have been clarified in Ref. [46]. The temperature equation given by Eq. (10) is solved with the classical fourth-order Runge-Kutta scheme for time discretization [26]. The isotropic central schemes are applied for spatial discretization. For
Fig. 1. Schematic illustration of the hydrophilic-hydrophobic mixed surface. The bottom substrate and the sides of the square pillar are hydrophilic, whereas the top surface of the pillar is hydrophobic. 848
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
0.0021
heat flux
0.0018 0.0015 0.0012 0.0009
(a) ΔT = 0.015
0.0006
0.01
0.02
0.03
0.04
wall superheat Fig. 3. Boiling curve of the hydrophilic-hydrophobic mixed surface with θphi ≈ 45. 3∘ and θpho ≈ 103∘ . The pillar width and pillar height are taken as W = 40 l.u. and H = 20 l.u., respectively.
and H = 20 l.u., respectively. The variations of the heat flux with the contact angle θphi in the cases of ΔT = 0.015, 0.018, and 0.021 are presented in Fig. 4, where the heat flux is the time average of the transient heat flux during t ∈ [103δt , 2 × 10 4δt ]. From Fig. 4 we can see that in the case of ΔT = 0.015 the heat flux experiences a slight change when θphi is increased from 31. 9∘ to 50. 8∘. After that, the heat flux rises rapidly, reaching its peak at θphi ≈ 66∘ ∼ 73∘, and then the heat flux decreases. In the cases of ΔT = 0.018 and 0.021, the heat flux initially increases with the increase of θphi . After reaching its peak value, the heat flux basically shows a declining trend. Particularly, it can be found that the peak values of the heat fluxes in the three cases are achieved with different θphi , i.e., θphi ≈ 66∘ ∼ 73∘, 57∘, and 40∘ , respectively, which means that a smaller θphi or a larger wettability difference Δθ = θpho−θphi is required to optimize the heat flux when the wall superheat increases. We take the case of ΔT = 0.018 as an example to illustrate the variations of the heat flux with the contact angle θphi . Figs. 5, 6, 7 and 8 display some snapshots of the boiling processes on the hydrophilichydrophobic mixed surface at ΔT = 0.018 with θphi ≈ 35. 6∘ , 40∘ , 50. 8∘, and 73. 2∘ , respectively. The nucleate boiling can be observed in Fig. 5, where θphi is set to 35. 6∘. However, it can be found that the bubble nucleation in Fig. 5 is completely attributed to the hydrophobicity of the top surface of the pillar and there are no bubbles nucleated on the
(b) ΔT = 0.017
(c) ΔT = 0.032 Fig. 2. Snapshots of the boiling processes at (a) ΔT = 0.015, (b) 0.017 , and (c) 0.032 . From left to right: t = 7000δt , 12000δt , and 18000δt . The contact angles θphi ≈ 45. 3∘ and θpho ≈ 103∘ . The pillar width and pillar height are taken as W = 40 l.u. and H = 20 l.u., respectively.
T = 0.015 T = 0.018 T = 0.021
0.0020
heat flux
Fig. 2(b), which enhance the nucleate boiling. When the wall superheat increases to 0.032, a large portion of the heating surface is covered by vapor and the boiling process enters the transition boiling regime, as shown in Fig. 2(c). The boiling curve is shown in Fig. 3, where the heat flux is the time average of the transient heat flux during t ∈ [103δt , 2 × 10 4δt ], which is in turn the spatial average of the local heat flux over the whole heating surface. Initially, the heat flux increases smoothly as the wall superheat is increased from 0.013 to 0.015. Later, it experiences a sharp rise and reaches the critical heat flux around ΔT = 0.025 ∼ 0.028. After that, the heat flux gradually decreases due to the appearance of transition boiling. Overall, the boiling curve displayed in Fig. 3 is qualitatively in good agreement with the typical boiling curve for water [49].
0.0015
0.0010
0.0005 30
40
50
contact angle 3.2. The effect of the contact angle θphi
60 phi
70
80
(deg.)
Fig. 4. Variations of the heat flux with the contact angle θphi in the cases of ΔT = 0.015, 0.018, and 0.021. The pillar width and height are fixed at W = 40 l.u. and H = 20 l.u., respectively. The contact angle of the top surface of the pillar is taken as θpho ≈ 103∘ .
In this section, numerical simulations are carried out to investigate the influence of the contact angle of the hydrophilic region of the mixed surface, i.e., θphi . The pillar width and height are fixed at W = 40 l.u. 849
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
Fig. 5. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 35. 6∘ . Cross-sectional view at y = 0.5L y . (a) t = 2000δt , (b) 4000δt , (c) 6000δt , and (d) 10000δt .
bottom substrate. When θphi increases from 35. 6∘ to 40∘, some bubbles are nucleated on the bottom substrate, as shown in Fig. 6(b). Subsequently, the bubbles nucleated on the bottom substrate interact with the bubble generated on the top surface of the pillar, leading to a larger bubble, as can be seen in Fig. 6(c). Owing to the bubble nucleation on the bottom substrate as well as the interaction between the bubbles nucleated on the bottom substrate and the bubble on the top of the pillar, the boiling heat transfer is enhanced and therefore the heat flux obtained with θphi ≈ 40∘ is larger than that obtained with θphi ≈ 35. 6∘ in the case of ΔT = 0.018, which can be observed in Fig. 4. Furthermore, by comparing Fig. 7 with Fig. 6, we can find that the nucleate boiling is further enhanced as the contact angle θphi increases from 40∘ to 50. 8∘ since the bubbles nucleated on the bottom substrate appear earlier in Fig. 7. The snapshots shown in Figs. 7(a) and 6(a) are
both taken at t = 2500δt . Nevertheless, it can be seen that small bubbles have been nucleated on the bottom substrate in Fig. 7(a) while there are no bubbles on the bottom substrate in Fig. 6(a). Accordingly, the heat flux obtained with θphi ≈ 50. 8∘ is larger than that obtained with θphi ≈ 40∘ , as seen in Fig. 4. However, when the contact angle θphi is increased to 73. 2∘ , the heat flux decreases considerably because a large amount of vapor is generated on the bottom surface, resulting in the transition boiling, as shown in Fig. 8. In summary, the numerical results shown in Figs. 5–8 reveal that an appropriate increase of the contact angle θphi can promote the bubble nucleation on the bottom substrate and hence enhances the nucleate boiling on the hydrophilic-hydrophobic mixed surface, but a larger θphi may make the boiling enter the transition or film boiling regime. In Fig. 4 it has been shown that the optimal θphi is different for different
Fig. 6. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 40∘ . (Top) Three-dimensional view and (bottom) cross-sectional view at y = 0.5L y . (a) t = 2500δt , (b) 4500δt , (c) 6000δt , and (d) 13000δt . 850
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
Fig. 7. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 50. 8∘ . Cross-sectional view at y = 0.5L y . (a) t = 2500δt , (b) 4500δt , (c) 7000δt , and (d) 11000δt .
0.0018
wall superheats, and decreases with the increase of the wall superheat, which is reasonable since a smaller contact angle can prevent the formation of a continuous vapor film on the heating surface when the vapor generation is enhanced by a higher wall superheat [26].
T = 0.018 T = 0.0165 T = 0.015
heat flux
0.0016
3.3. The effects of H and W and the orthogonal array tests In this section, the effects of the pillar height and width are studied. The contact angles θphi and θpho are chosen as θphi ≈ 50. 8∘ and θpho ≈ 103∘ , respectively. First, we investigate the effect of the pillar width and the pillar height is fixed at H = 20 l.u. Fig. 9 depicts the variations of the heat flux with the pillar width in the cases of ΔT = 0.015, 0.0165, and 0.018. From the figure we can see that in the case of ΔT = 0.015 the heat flux initially increases with the increase of the pillar width and reaches its peak at W = 80 l.u. After that, the heat flux gradually declines. A similar trend can be observed in the case of ΔT = 0.0165 but the peak value of the heat flux is achieved at W = 70 l.u. In the case of ΔT = 0.018, the heat flux rises rapidly as the pillar width increases from W = 20 l.u. to 40 l.u., and then the heat flux shows a downward trend, although there is a slight rise of the heat flux when the pillar width is increased from W = 50 l.u. to 60 l.u. Similar to Fig. 4, Fig. 9 shows that the optimal pillar width varies with the wall superheat. More specifically, it decreases with the increase of the wall superheat. In Ref. [16], two-dimensional (2D) LB simulations have been performed to investigate the effect of the pillar width on the boiling performance of the hydrophilic-hydrophobic mixed surface and it has been shown that a large pillar width may result in a substantial decrease of the heat flux. Such a phenomenon can also
0.0014 0.0012 0.0010 0.0008 20
40
60
80
100
120
pillar width (l.u.) Fig. 9. Variations of the heat flux with the pillar width in the cases of ΔT = 0.015, 0.0165 , and 0.018. The pillar height is fixed at H = 20 l.u.
be observed in Fig. 9 in the cases of ΔT = 0.0165 and 0.018, which is attributed to the fact that, when a very large pillar width is applied, the boiling process may enter the transition or film boiling regime because of the hydrophobicity of the top surface of the pillar. The results that the optimal θphi and the optimal pillar width decrease with the increase of the wall superheat reveal that the hydrophilic region of the mixed surface should play an increasingly important role when the wall superheat increases. As previously mentioned, the
Fig. 8. Snapshots of the boiling process at ΔT = 0.018 with θphi ≈ 73. 2∘ . Cross-sectional view at y = 0.5L y . (a) t = 2000δt , (b) 4000δt , (c) 7000δt , and (d) 11000δt . 851
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
advantage of hydrophobicity is promoting boiling onset at a low wall superheat, while the advantage of hydrophilicity lies in that it provides a relatively high critical heat flux. Hence, with the increase of the wall superheat, the hydrophilic region gradually becomes important and decreasing the pillar width can reduce the influence of the hydrophobic region (the top surface of the pillar). Similarly, a smaller θphi can enhance the influence of the hydrophilic region when the wall superheat increases. The differences between the present 3D simulations and the previous 2D simulations [16] should be noticed. In 2D simulations, increasing the pillar width does not affect the bubble nucleation at the two corners formed by the bottom substrate and the pillar. However, in 3D simulations, increasing the pillar width will influence the bubble nucleation on the four edges formed by the bottom substrate and the square pillar, as shown in Fig. 10, where some snapshots of the boiling processes at ΔT = 0.0165 are displayed. By comparing Fig. 10(b) with Fig. 10(a), we can see that the bubble nucleation on the bottom substrate is affected when the pillar width changes from 40 l.u. to 70 l.u. For example, some small bubbles nucleated on the bottom substrate can be observed in Fig. 10(b) at t = 8000δt and 11000δt , whereas there are no similar bubbles on the bottom substrate in Fig. 10(a). This is also the reason why the heat flux obtained with W = 70 l.u. is larger than that obtained with W = 40 l.u. in the case of ΔT = 0.0165. The effect of the pillar height is shown in Fig. 11, where the variations of the heat flux with the pillar height in the cases of ΔT = 0.015, 0.0165, and 0.018 are presented. The pillar width is fixed at W = 40 l.u. The figure illustrates that in the case of ΔT = 0.015 the heat flux rises smoothly as the pillar height increases from H = 20 l.u. to 50 l.u. In the other two cases, the heat flux fluctuates with the increase of the pillar height, but the fluctuation is very small. Overall, the influence of the pillar height is relatively small in comparison with the effects of the
0.0018
T = 0.018 T = 0.0165 T = 0.015
heat flux
0.0016 0.0014 0.0012 0.0010 0.0008 20
30
40
50
pillar height (l.u.) Fig. 11. Variations of the heat flux with the pillar height in the cases of ΔT = 0.015, 0.0165 , and 0.018. The pillar width is fixed at W = 40 l.u.
contact angle θphi and the pillar width. Increasing the pillar height does not change the area of the top surface of the pillar, but it reduces the ratio of the hydrophobic region to the whole heating surface, i.e., Spho Swhole , which varies from 6.23% to 5.25% when the pillar height is increased from H = 20 l.u. to 50 l.u. It can be readily found that such a change is much smaller than that in the previous 2D simulations [16], in which the ratio Spho Swhole degenerates to L pho L whole . This may be the reason why the heat flux reported in the 2D simulations considerably changes with the pillar height [16] whereas in the present 3D simulations the variations of the heat flux with the pillar height are relatively small.
Fig. 10. Snapshots of the boiling process at ΔT = 0.0165 with θphi ≈ 50. 8∘ and θpho ≈ 103∘ . (a) W = 40 l.u. and (b) 70 l.u. From left to right: t = 3500δt , 6000δt , 8000δt , and 11000δt . 852
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
Table 1 The orthogonal array tests and the result analsyis. Test serial number
Heat flux
Factors A/θphi
B/pillar width
A× B
C/pillar height
A× C
Kj1 × 103
1 (40°) 1 (40°) 1 (40°) 1 (40°) 2 (57°) 2 (57°) 2 (57°) 2 (57°) 3 (66°) 3 (66°) 3 (66°) 3 (66°) 4 (73.2°) 4 (73.2°) 4 (73.2°) 4 (73.2°) 4.9361
1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (80 l.u.) 4.6191
1 2 3 4 2 1 4 3 3 4 1 2 4 3 2 1 4.8702
1 (20 l.u.) 2 (30 l.u.) 3 (40 l.u.) 4 (50 l.u.) 3 (40 l.u.) 4 (50 l.u.) 1 (20 l.u.) 2 (30 l.u.) 4 (50 l.u.) 3 (40 l.u.) 2 (30 l.u.) 1 (20 l.u.) 2 (30 l.u.) 1 (20 l.u.) 4 (50 l.u.) 3 (40 l.u.) 4.4447
1 2 3 4 4 3 2 1 2 1 4 3 3 4 1 2 5.0041
Kj2 × 103
5.6491
4.7152
4.4585
4.9961
4.9611
Kj3 × 103
4.7576
5.2888
4.4322
4.8724
4.7815
Kj4 × 103
3.8650
4.5847
5.4469
4.8946
4.4611
K¯ j1 × 103
1.2340
1.1548
1.2176
1.1112
1.2510
K¯ j2 × 103
1.4123
1.1788
1.1146
1.2490
1.2403
K¯ j3 × 103
1.1894
1.3222
1.1081
1.2181
1.1954
K¯ j4 × 103
0.9663
1.1462
1.3617
1.2237
1.1153
Range × 103 Rank Optimization
0.4460
0.1760
0.2536
0.1379
0.1358
1 A2
3 B3
2 (A × B)4
4 C2
5 (A × C)1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0.0011666 0.0012148 0.0012758 0.0012789 0.0012131 0.0014257 0.0016457 0.0013646 0.0011142 0.0013971 0.0012915 0.0009548 0.0011252 0.0006776 0.0010758 0.0009864
contact angle θphi and the pillar width has an important influence. Similarly, the “Range” of the interaction A× C is almost the same as that of factor C. Hence, the interaction effects cannot be neglected. As a result, the optimal combination cannot be directly determined by the largest K¯2i ( j = 2 corresponds to factor B) and the largest K¯ 4i ( j = 4 corresponds to factor C). According to the analysis of orthogonal array tests with two-factor interfactions, when the interaction effects have an important influence, a comparison of the results (the heat fluxes in the present study) should be made so as to find the optimal combination. According to the results displayed in Table 1, the optimal combination is A2B3C1.
Finally, the orthogonal array tests are performed to identify the most important factor among the aforementioned three influential factors: the contact angle θphi , the pillar width, and the pillar height. The wall superheat is taken as ΔT = 0.018 and we consider four levels for each factor, namely, 40∘, 57∘, 66∘, and 73. 2∘ for θphi , 20 l.u., 30 l.u., 40 l.u., and 80 l.u. for the pillar width, and 20 l.u., 30 l.u., 40 l.u., and 50 l.u. for the pillar height. Besides, two two-factor interaction effects are also considered, i.e., the interaction effect between θphi and the pillar width, and the interaction effect between θphi and the pillar height. Hence the orthogonal array tests include five factors and four levels shown in Table 1, i.e., L16(45), containing 16 tests in total. In the table, A× B and A× C represent the interaction effect between θphi and the pillar width, and the interaction effect between θphi and the pillar height, respectively. The results of the orthogonal array tests are shown in Table 1, where Kji denotes the sum of the heat fluxes corresponding to the i-th level of the j-th factor, K¯ ji is the average value of Kji , and the “Range” is the difference between the maximum and the minimum of the four average values (K¯ j1, K¯ j2 , K¯ j3 , and K¯ j4 ) in each group. From the table we can see that the “Range” of factor A, i.e., the contact angle θphi , is much larger than the others, which means that θphi is the most important influential factor among the three investigated factors. Meanwhile, it can be found that, the “Range” of factor C (the pillar height) and the “Range” of the interaction A× C are much smaller than the others, which indicates that the pillar height is the least influential one. These findings are consistent with those revealed by Figs. 4, 9, and 11. In addition, it can be found that K¯12 ( j = 1 corresponds to factor A) is the largest of the four K¯1i (i = 1, 2, 3, and 4). Hence, the level A2, i.e., θphi ≈ 57∘ , is the best of the four levels for factor A when the wall superheat is taken as ΔT = 0.018. According to Table 1, if the interaction effects are not taken into account, the optimal combination is A2B3C2. However, the table shows that the “Range” of the interaction A× B is larger than that of factor B (the pillar width), which implies that the interaction effect between the
4. Conclusions In this paper, we have investigated the boiling performance on a type of hydrophilic-hydrophobic mixed surfaces using a 3D thermal LB model with liquid-vapor phase change, which consists of a 3D pseudopotential multiphase LB model for simulating the flow field and a finite-difference solver for simulating the temperature field. The mixed surface is square-pillar textured with the pillar being composed of hydrophilic side walls and a hydrophobic top surface. The effects of the contact angle θphi , the pillar width, and the pillar height have been investigated. The main findings and conclusions are summarized as follows. (i) An appropriate increase of the contact angle θphi can promote the bubble nucleation on the bottom substrate and enhances the nucleate boiling on the hydrophilic-hydrophobic mixed surface, but a larger θphi may make the boiling enter the transition or film boiling regime. (ii) Different from the previous 2D simulations, the present 3D simulations showed that the heat flux initially increases with the increase of the pillar width and basically shows a declining trend after reaching its peak value. 853
Applied Thermal Engineering 142 (2018) 846–854
Y. Yu et al.
[19] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. [20] Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combust. Sci. 52 (2016) 62–105. [21] A. Xu, W. Shyy, T. Zhao, Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries, Acta Mech. Sin. 33 (2017) 555–574. [22] Y. Shi, G.H. Tang, H.H. Xia, Lattice Boltzmann simulation of droplet formation in Tjunction and flow focusing devices, Comput. Fluids 90 (2014) 155–163. [23] G. Hazi, A. Markus, On the bubble departure diameter and release frequency based on numerical simulation results, Int. J. Heat Mass Transf. 52 (2009) 1472–1480. [24] A. Márkus, G. Házi, On pool boiling at microscale level: the effect of a cavity and heat conduction in the heated wall, Nucl. Eng. Des. 248 (2012) 238–247. [25] L. Biferale, P. Perlekar, M. Sbragaglia, F. Toschi, Convection in multiphase fluid flows using lattice Boltzmann methods, Phys. Rev. Lett. 108 (2012) 104502. [26] 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. [27] Q. Li, P. Zhou, H.J. Yan, Pinning-depinning mechanism of the contact line during evaporation on chemically patterned surfaces: a lattice Boltzmann study, Langmuir 32 (2016) 9389–9396. [28] Y. Yu, Q. Li, C.Q. Zhou, P. Zhou, H.J. Yan, Y. Yu, Q. Li, C.Q. Zhou, P. Zhou, H.J. Yan, Investigation of droplet evaporation on heterogeneous surfaces using a three-dimensional thermal multiphase lattice Boltzmann model, Appl. Therm. Eng. 127 (2017) 1346–1354. [29] Q. Li, Q.J. Kang, M.M. Francois, A.J. Hu, Lattice Boltzmann modeling of self-propelled Leidenfrost droplets on ratchet surfaces, Soft Matter 12 (2016) 302–312. [30] Y.H. Qian, D. d'Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479–484. [31] S. Chen, H. Chen, D. Martnez, W. Matthaeus, Lattice Boltzmann model for simulation of magnetohydrodynamics, Phys. Rev. Lett. 67 (1991) 3776–3779. [32] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546–6562. [33] D. d'Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, L.S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. Ser. A 360 (2002) 437–451. [34] K.N. Premnath, J. Abraham, Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow, J. Comput. Phys. 224 (2007) 539–559. [35] D. Sun, M. Zhu, J. Wang, B. Sun, Lattice Boltzmann modeling of bubble formation and dendritic growth in solidification of binary alloys, Int. J. Heat Mass Transf. 94 (2016) 474–487. [36] A. Xu, L. Shi, T.S. Zhao, Accelerated lattice Boltzmann simulation using GPU and OpenACC with data management, Int. J. Heat Mass Transf. 109 (2017) 577–588. [37] K.N. Premnath, S. Banerjee, Incorporating forcing terms in cascaded lattice Boltzmann approach by method of central moments, Phys. Rev. E 80 (2009) 036702. [38] D. Lycett-Brown, K.H. Luo, Multiphase cascaded lattice Boltzmann method, Comput. Math. Appl. 67 (2014) 350–362. [39] Q. Liu, Y.L. He, D. Li, Q. Li, Non-orthogonal multiple-relaxation-time lattice Boltzmann method for incompressible thermal flows, Int. J. Heat Mass Transf. 102 (2016) 1334–1344. [40] Q. Li, D.H. Du, L.L. Fei, K.H. Luo, Y. Yu, Three-dimensional non-orthogonal multiple-relaxation-time lattice Boltzmann model for multiphase flows, arXiv preprint, (2018) arXiv:1805.08643. [41] 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 87 (2013) 053301. [42] D. Zhang, K. Papadikis, S. Gu, Three-dimensional multi-relaxation time latticeBoltzmann model for the drop impact on a dry surface at large density ratio, Int. J. Multiph. Flow 64 (2014) 11–18. [43] A. Xu, T.S. Zhao, L. An, L. Shi, A three-dimensional pseudo-potential-based lattice Boltzmann model for multiphase flows with large density ratio and variable surface tension, Int. J. Heat Fluid Flow 56 (2015) 261–271. [44] Y. Shi, G.H. Tang, H.H. Xia, Investigation of coalescence-induced droplet jumping on superhydrophobic surfaces and liquid condensate adhesion on slit and plain fins, Int. J. Heat Mass Transf. 88 (2015) 445–455. [45] P. Yuan, L. Schaefer, Equations of state in a lattice Boltzmann model, Phys. Fluids 18 (2006) 042101. [46] Q. Li, P. Zhou, H.J. Yan, Improved thermal lattice Boltzmann model for simulation of liquid-vapor phase change, Phys. Rev. E 96 (2017) 063303. [47] W. Gong, Y.Y. Yan, S. Chen, E. Wright, A modified phase change pseudopotential lattice Boltzmann model, Int. J. Heat Mass Transf. 125 (2018) 323–329. [48] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591–1598. [49] T.L. Bergman, A.S. Lavine, F.P. Incropera, D.P. Dewitt, Fundamentals of Heat and Mass Transfer, JohnWiley & Sons, New York, 2011.
(iii) The optimal θphi and the optimal pillar width change with the wall superheat. Particularly, both of them decrease with the increase of the wall superheat, which reveals that the hydrophilic region of the mixed surface should play an increasingly important role when the wall superheat increases. (iv) The orthogonal array tests showed that the contact angle θphi is the most important influential factor among the three investigated factors and the pillar height is the least influential one. Moreover, it is found that the interaction effect between the contact angle θphi and the pillar width has an important influence. In future work, we will consider the mixed surfaces with an array of square pillars, through which the effect of the width between two pillars can be investigated. Such a study requires a much larger computational domain than that employed in the present work and also much more computational cost. Acknowledgments This work was supported by the National Natural Science Foundation of China (No. 51506227) and the Foundation for the Author of National Excellent Doctoral Dissertation of China (No. 201439). References [1] Y. Koizumi, M. Shoji, M. Monde, Y. Takata, N. Nagai, Boiling: Research and Advances, Elsevier, 2017. [2] V.K. Dhir, Boiling heat transfer, Annu. Rev. Fluid Mech. 30 (1998) 365–401. [3] V.K. Dhir, G.R. Warrier, E. Aktinol, Numerical simulation of pool boiling: a review, J. Heat Transfer 135 (2013) 061502. [4] H.J. Jo, H.S. Ahn, S.H. Kang, M.H. Kim, A study of nucleate boiling heat transfer on hydrophilic, hydrophobic and heterogeneous wetting surfaces, Int. J. Heat Mass Transf. 54 (2011) 5643–5652. [5] G. Lu, X.D. Wang, Y.Y. Duan, A critical review of dynamic wetting by complex fluids: from Newtonian fluids to non-Newtonian fluids and nanofluids, Adv. Colloid Interface Sci. 236 (2016) 43–62. [6] D.E. Kim, D.I. Yu, D.W. Jerng, M.H. Kim, H.S. Ahn, Review of boiling heat transfer enhancement on micro/nanostructured surfaces, Exp. Therm Fluid Sci. 66 (2015) 173–196. [7] B. Bourdon, R. Rioboo, M. Marengo, E. Gosselin, J.D. Coninck, Influence of the wettability on the boiling onset, Langmuir 28 (2012) 1618–1624. [8] H.J. Jo, M. Kaviany, S.H. Kim, M.H. Kim, Heterogeneous bubble nucleation on ideally-smooth horizontal heated surface, Int. J. Heat Mass Transf. 71 (2014) 149–157. [9] S.G. Kandlikar, A theoretical model to predict pool boiling CHF incorporating effects of contact angle and orientation, J. Heat Transfer 123 (2000) 1071–1079. [10] A.R. Betz, J. Xu, H. Qiu, D. Attinger, Do surfaces with mixed hydrophilic and hydrophobic areas enhance pool boiling? Appl. Phys. Lett. 97 (2010) 103115. [11] X. Dai, X. Huang, F. Yang, X. Li, J. Sightler, Y. Yang, C. Li, Enhanced nucleate boiling on horizontal hydrophobic-hydrophilic carbon nanotube coatings, Appl. Phys. Lett. 102 (2013) 161605. [12] B.J. Suroto, Effects of hydrophobic-spot periphery and subcooling on nucleate pool boiling from a mixed-wettability surface, J. Therm. Sci. Technol. 8 (2013) 294–308. [13] H.J. Jo, S.H. Kim, H.S. Park, M.H. Kim, Critical heat flux and nucleate boiling on several heterogeneous wetting surfaces: controlled hydrophobic patterns on a hydrophilic substrate, Int. J. Multiph. Flow 62 (2014) 101–109. [14] H.J. Jo, D.I. Yu, H. Noh, H.S. Park, M.H. Kim, Boiling on spatially controlled heterogeneous surfaces: wettability patterns on microstructures, Appl. Phys. Lett. 106 (2015) 181602. [15] C.S.S. Kumar, Y.W. Chang, P.H. Chen, Effect of heterogeneous wettable structures on pool boiling performance of cylindrical copper surfaces, Appl. Therm. Eng. 127 (2017) 1184–1193. [16] 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. [17] C. Shen, C. Zhang, Y. Bao, X. Wang, Y. Liu, L. Ren, Experimental investigation on enhancement of nucleate pool boiling heat transfer using hybrid wetting pillar surface at low heat fluxes, Int. J. Therm. Sci. 130 (2018) 47–58. [18] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364.
854