- Email: [email protected]
A multicomponent multiphase enthalpy-based lattice Boltzmann method for droplet solidification on cold surface with different wettability
International Journal of Heat and Mass Transfer 127 (2018) 136–140
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Technical Note
A multicomponent multiphase enthalpy-based lattice Boltzmann method for droplet solidification on cold surface with different wettability Peng Xu, Sichuan Xu ⇑, Yuan Gao, Pengcheng Liu Tongji University, School of Automotive Studies, Caoan Road 4800, Shanghai 201804, PR China
a r t i c l e
i n f o
Article history: Received 8 May 2018 Received in revised form 2 July 2018 Accepted 4 July 2018
Keywords: Lattice Boltzmann method Droplet Solidification Multicomponent multiphase
a b s t r a c t In this paper, a liquid droplet solidification contacting cold solid surface is modeled by a novel multicomponent multiphase enthalpy-based lattice Boltzmann method (LBM). The freezing process of liquid droplet coupled with thermal and mass transfer on cold flat surface with different wettability (h = 160°, 135°, 60°, 20°) is investigated in detail. The solid volume fraction of droplet solidification with Fourier number are shown under different contact angles. The simulation results show that the solidification time for droplet contacting on a hydrophilic surface is much smaller than that of hydrophobic surface when liquid droplet shares the same initial radius. Ó 2018 Published by Elsevier Ltd.
1. Introduction The solidification process is a common phenomenon in many physical fields, such as air conditioning, phase change energy storage, defrosting, and etc. It is a complicated multiphase flow problem coupled with heat and mass transfer and phase change, especially for multicomponent multiphase system. Over the past several decades, many scientists have studied the liquid-solid phase change process for water droplet either by experimental observation or by macroscopic numerical analysis [1,2]. The LBM has recently gained much attention as a powerful tool to simulate complex physical problem due to advantages of capacity for investigating complicated geometries, simple implementation, high computation efficiency, easy implementation of parallel-processing. It is difficult to determine the phase change interface by conventional computational fluid dynamics (CFD) because of complicated heat and mass transfer mechanism with nonlinear solid-liquid interface characteristics. Three types of liquid-solid phase change LBM models have been developed: immersed boundary method [3], phase-field method [4] and enthalpy-based method [5–13]. Among them, the enthalpy-based LBM model is employed widely due to its simplicity and effectiveness for solid-liquid phase change problem with assumption of no solid-phase movement. Jiaung et al. [5] proposed the enthalpy-based LBM model to simulate phase change problem ⇑ Corresponding author. E-mail address: [email protected] (S. Xu). https://doi.org/10.1016/j.ijheatmasstransfer.2018.07.017 0017-9310/Ó 2018 Published by Elsevier Ltd.
with heat conduction for the first time. Subsequently, Chatterjee and Chakraborty [6] successfully applied enthalpy updating procedure and one relaxation factor for convergence of iteration. Huber et al. [7] investigated convection-dominated solid-liquid phase change in a cavity by setting number of iterations to be one for calculation efficiency. In this regard, Eshraghi and Felicelli [8] introduced an implicit LBM to handle with the source term, thus avoiding iteration steps. Recently, Huang [9] proposed a new thermal LBM model for solid-liquid phase change by adding latent heat source into transient term, which is characterized by avoiding the iteration process and linear equations. Since 2015, some research focus on the heat and mass transfer process of solidification for droplets based on enthalpy-based LBM [10–14]. Among them, some [10,11] are applied for single component multiphase system while others [12–14] for multicomponent situations. To the best of knowledge, droplet solidification on cold surface with different wettability using multicomponent multiphase LBM has never been studied. In this paper, the novel multicomponent multiphase enthalpy-based LBM model based on Zhao and Cheng [14] is developed for simulating solidification of a liquid droplet surrounded by air and contacting a cold solid surface.
2. Description of enthalpy-based multicomponent multiphase LBM model In this section, we briefly present the new enthalpy-based multicomponent multiphase LBM model for liquid-solid phase change.
137
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
The main characteristics for new LBM model are as follows: (1) Triple distribution functions were used, two for velocity field and one for temperature field; (2) Iteration steps and a group of linear equations are avoided by modifying the equilibrium distribution function of the temperature field; (3) The velocity field and temperature field are linked each other through equilibrium function and buoyancy source term; (4) The bounce-back boundary condition for liquid-solid interface is in form of an additional collision term in LBM equation for velocity field. 2.1. Velocity field The evolution for velocity f i with immersed moving boundary scheme and force form of exact difference method (EDM) is described as [14]: k
k
k
f i ðx þ ei Dt; t þ DtÞ ¼ f i ðx; tÞ kðsÞ
þ Bk Xi
1B
sk
k
kðeqÞ
½f i ðx; tÞ f i
ðx; tÞ
k
þ ð1 Bk ÞDf i ðx; tÞ
ð1Þ
where f i is the density distribution function at node x and time t, ei is the discrete velocity in the i direction, k is the component indicator (1 for liquid phase and 2 for gas phase), the D2Q9 discretization scheme for 2D geometry is used in our model, sk is the relaxation time for component k, calculated by:
sk
tk
1 ¼ 2 þ c s Dt 2
pffiffiffi equaling to c= 3. Bk is a weighting function that relates to the dimensionless relaxation time and liquid volume fraction, represented by:
Bk ¼
ð1 f l Þðsk 0:5Þ f l þ sk 0:5
ð3Þ
where the liquid fraction f l is calculated by:
fl ¼
8 > <0 > :
H 6 Hs
HHs Hl Hs
ð4Þ
H P Hl
with Hs and Hl corresponding to the enthalpies of solidus and liquidus temperature, respectively. kðeqÞ fi
denotes the equilibrium distribution function for the velocity and given by:
"
kðeqÞ
fi
¼ xi qk 1 þ
2
ei ub ðei ub Þ þ c2s 2c4s
u2b 2c2s
#
ð5Þ kðsÞ
where wi is the weighting coefficient, Xi is an additional collision term bouncing back the non-equilibrium part of the distribution function, represented as:
XikðsÞ ¼ f ki ðx; tÞ f ki ðx; tÞ þ f ikðeqÞ ðqk ; us Þ f kðeqÞ ðqk ; ub Þ i
ð6Þ
where us is the velocity of solid phase and i is the opposite direction of i. k Df i
is the force term of EDM proposed by Kupershtokh and Medvedev: k Df i ðx; tÞ
¼
kðeqÞ fi ð
qk ; ub þ Duk Þ
qk ; ub Þ
kðeqÞ fi ð
ð7Þ
where Duk is the velocity difference due to total force and can be given as:
Duk ¼
F k Dt
qk
F k2 ¼ wk ðqk ðxÞÞ
s XX Gkk ðx; x0 Þwk ðqk ðx0 ÞÞðx0 xÞ
ð8Þ
ð9Þ
k
x0
where wk ðqk ðxÞÞ is the effective mass, which is a function of the local density. The density value is used for effective mass in this paper for represent two different fluid components. simplicity. k and k 0 Gkk ðx; x Þ is a coupling constant depending on a Green’s function. The fluid-solid interaction force is given by:
F k2 ¼ wk ðqk ðxÞÞ
X Wðx; x0 Þsðx0 Þðx0 xÞ
ð10Þ
x0
where sðx0 Þ is a none-zero constant at the fluid-solid interface and zero otherwise, Wðx; x0 Þ represents fluid-solid interaction strength and different levels of wettability can be obtained by adjusting the parameter Wðx; x0 Þ. The macroscopic density, bulk velocity and real physical velocity can be obtained by:
qk ¼
X k X X X Dt f i ; ub ¼ qk uk = qk ; ur ¼ qk uk þ F k 2 i k k k
ð11Þ
2.2. Temperature field The passive scalar method [15] is the extensively applied approach for incorporating thermal effects into LBM, especially for multicomponent multiphase problem. In this work, we adopted the passive scalar method by using single LBM temperature equation. The evolution equation for temperature field is given by:
g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ
; Hs < H < Hl
1
and the gravitational force F k3 . The pseudo-potential was applied to represent the particle interaction force:
ð2Þ
tk is the kinematic viscosity, cs is the lattice sound speed
where
where F k is the total force, composed of three parts, including the particle interaction force F k1 , the fluid-solid interaction force F k2
1
sT
½g i ðx; tÞ g eq i ðx; tÞ
ð12Þ
where sT is the relaxation time for temperature field which is given below:
sT ¼
a
c2s Dt
þ
1 2
ð13Þ
where a is thermal diffusivity and defined as a ¼ k=ðC p qÞ (k and C p are thermal conductivity and specific heat capacity at constant pressure, respectively). g eq i ðx; tÞ is the temperature equilibrium distribution function, given as:
g eq i
8 u2 i¼0 < H C p T þ xi C p Tð1 2c 2Þ s h i ¼ 2 2 e u ðe uÞ : xi q 1 þ i þ i u i–0 2c4 2c2 c2 s
s
ð14Þ
s
P where H is the enthalpy by H ¼ ig . The temperature can be obtained as follows:
8 H > T HCs p;s H 6 Hs > > s < Hl H HHs T þ T ; T ¼ Hl Hs s Hl Hs l Hs < H < Hl > > > : T l þ HHl H P Hl C
ð15Þ
p;l
2.3. Computational domain and boundary condition The 2D computational domain with the size of N x N x ¼ 150 150 in lattice unit is used in this work. One droplet with
138
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
the designated radius of 25 lattice units is contacting the cold flat surface, immersing in air. For the velocity field, half-way bounce back conditions are applied on the top and bottom walls with periodic boundary conditions imposed for left and right sides of computation domain. As for the temperature field, the non-equilibrium extrapolation scheme is utilized for four sides. It is noted that dimensionless temperature of top and bottom walls are kept at constant (T ¼ 0) while adiabatic thermal boundary is set for left and right sides. The temperature of droplet is higher than that of cold flat surface and melting point so that the droplet can keep in the form of liquid shape before the interfere of temperature effects. Four cases of different contact angles for droplet interacting with cold solid surface is investigated: h = 160°, 135°, 60°, 20°. The temperature field is activated after velocity field reaches equilibrium state. 3. Results and discussion In our work, the problem for droplet solidification above cold flat surface can be characterized by some dimensionless parameters: Ohnesorge number, Fourier number, Prandtl number, dimensionless temperature number, Stefan number, which are defined by
Oh ¼
ld
ðqd rDd Þ0:5
; Fo ¼
at D2d
; Pr ¼
v
a
; T ¼
T Tm CpT c ; Ste ¼ Td Tm L
ð16Þ
where ud is the liquid dynamic viscosity, qd is the droplet density, r is surface tension, Dd is the droplet radius, a is the thermal diffusivity of ice, m is the kinematic viscosity, T m is the melting temperature, T d is the initial droplet temperature, C p is the specific heat at constant pressure, T c is the characteristic temperature and set T c ¼ T d T m , L is the latent heat. As mentioned before, the equilibrium state of droplet contacting solid surface is preferred followed by intervention of thermal effects. Therefore, the parameter affecting the wettability of solid
surface must be confirmed in advance. The contact angle between the fluid and solid surface can be calculated by Young’s equation. Three parameters need to be known including interfacial tension between two fluid components r12 and between each fluid component and the solid surface rs1 , rs2 , as shown in Eq. (17):
cos h ¼ ðrs2 rs1 Þ=r12 ¼ ðg2w g 1w Þ=ðgðq1 q2 Þ=2Þ
ð17Þ
According to Eq. (17), when the value of g 2w between gas phase and solid wall is greater than the value of g 1w between liquid phase and solid wall, the contact angle h is less than 90 , showing the solid surface is hydrophilic and vice versa. The fluid-fluid interaction coefficient is kept at 0.1 with g ¼ g 12 ¼ g 21 while interaction strength between each fluid and the solid surface varies from 0.05 to 0.05. In addition, fluid-solid interaction strength value is kept as g 1w ¼ g 2w . It is noted that all the units are in lattice units unless otherwise stated specifically. The contact angle can be evaluated by final steady state values of the droplet with height a, droplet length b and droplet radius R: tan h ¼ b=ðaðR aÞÞ. The contact angle between wetting phase and solid surface varies from 7° to 173° when g 1w changes from 0.05 to 0.05. Different level of wetting (non-wetting) conditions can be obtained by adjusting g 1w in the following simulation. Fig. 1 presents 2D simulation results of solidification process for droplet above cold solid surface with different contact angle (h = 160°, 135°, 60°, 20°). The size of whole domain is 150 150 lattice units with initial radius of droplet 25 lattice units. In the simulation, the Ohnesorge number is set to be Oh = 0.13 and Stefan number is Oh = 0.25. Initially, the temperature of droplet and gas is chosen to be T d = 2.0 (which is higher than that of melting temperature point T = 0.5) and T g = 0.0 respectively. The top and bottom wall of calculation domain is set to be T t = 0.0 and T b = 0.0. Other parameters are as follows: C p = 1.0, L = 0.25, viscosity ratio for gas and droplet v g : v l = 1:1, Prandtl number Pr = 0.71 and the thermal diffusivity k = 0.235. The red region denotes the ice region
Fig. 1. Shape evolution of droplet on surface at different wettability with time.
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
and the blue region inside the droplet denotes the liquid state. When the dimensionless temperature of droplet contacting the solid surface drops below T = 0.5, the part of droplet contacting solid surface starts to form ice firstly. The top area of droplet also ices because temperature of gas close to the droplet is low. Note that the droplet area that ices lastly is not the top point of droplet but high area inside droplet. When contact angle is h ¼ 160 , the solidification time is F o = 5.6, which is bigger than that of case with contact angle h ¼ 20 . With the reducing of contact angle (more hydrophilic), the solidification time reduces quickly.
139
To investigate the dimensionless temperature variation along the vertical direction cross the droplet, the results extracted from line ‘‘CD” as indicated in Fig. 2(a) are obtained in Fig. 2. Fig. 2(a) shows the dimensionless temperature contour at F o = 0.4 when droplet forms a contact angle h = 160° with solid surface. The dashed line is the shape of the liquid droplet. The heat transferred from droplet surface to ambient gas leads to temperature gradient and the dimensionless temperature difference between droplet and cold surface changes more quickly because of larger heat flux near the cold surface. To validate this expression, the local heat flux near cold surface denoted by ‘‘EF” line in Fig. 2(a) is calculated [11]. From Fig. 2(c), we can see that heat flux value between droplet and cold surface is almost 8 times larger than that between gas and cold surface. Fig. 2(b) shows the dimensionless temperature along the vertical line ‘‘CD”. At F o = 0.4, the dimensionless temperature of point ‘‘D” is 0.35 larger than that in ‘‘C” point. It is noted that 0.35 difference for dimensionless temperature can be seen during the whole solidification process. In addition, the dimensionless temperature near cold surface decreases quicker than far area due to almost constant low temperature boundary condition seen from two arrows in Fig. 2(b). After Fourier number reaches to F o = 5.6, the temperature profile across the droplet behaves linearly.
4. Conclusions A novel multicomponent multiphase enthalpy-based LBM model for droplet solidification contacting cold solid surface was proposed in this paper. The solidification process of droplet contacting wall with different wettability (h = 160°, 135°, 60°, 20°) is investigated through developed multicomponent multiphase enthalpy-based LBM model. The following conclusions can be made from this paper: (1) The dimensionless temperature in droplet drops quickly when droplet is placed above cold flat surface. Whereas the temperature drops much slower when droplet is during solidification process because of phase change. (2) The dimensionless temperature in droplet near the contact line between cold solid surface and droplet drops more quickly than that near droplet surface. The largest temperature difference around 0.4 exist in droplet during solidification when contact angel is h = 160° while 0.34 in droplet between the zenith and the lowest point. (3) For a droplet with same initial radius, the solidification time for hydrophobic cold surface is almost 5 times larger than that for hydrophilic surface. The freezing time for contact angle (h = 160°, 135°, 60°, 20°) is 5.6, 4.9, 3.0 and 1.1, respectively.
Conflict of interest We have no conflict of interest. We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.
Acknowledgements Fig. 2. The dimensionless temperature contour (a), the dimensionless temperature profile along the vertical diameter of droplet (b) and heat flux distribution (c) along the cold plate at F o = 0.4 with h = 160°.
This work was supported by National Key R&D Plan of China through Grant Nos. 2017YFB0102802.
140
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
References [1] X. Wu, W. Dai, W. Xu, L. Tang, Mesoscale investigation of frost formation on a cold surface, Exp. Therm. Fluid. 31 (2007) 1043–1048. [2] J.P. Hindmarsh, A.B. Russell, X.D. Chen, Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet, Int. J. Heat Mass Transf. 46 (2003) 1199–1213. [3] R. Huang, H.W., An immersed boundary-thermal lattice Boltzmann method for solid-liquid phase change, J. Comput. Phys. 277 (2014) 305–319. [4] W. Miller, S. Succi, D. Mansutti, Lattice Boltzmann model for anisotropic liquidsolid phase transition, Phys. Rev. E 86 (2001) 3578–3581. [5] W. Jiaung, J. Ho, C. Kuo, Lattice Boltzmann method for the heat conduction problem with phase change, Numer. Heat Transf. Part B: Fundam. 39 (2) (2001) 167–187. [6] D. Chatterjee, S. Chakraborty, An enthalpy-based lattice Boltzmann model for diffusion dominated solid–liquid phase transformation, Phys. Lett. A 341 (2005) 320–330. [7] C. Huber, A. Parmigiani, B. Chopard, M. Manga, Lattice Boltzmann model for melting with natural convection, Int. J. Heat Fluid Flow 29 (5) (2009) 1469– 1480.
[8] M. Eshraghi, S.D. Felicelli, An implicit lattice Boltzmann model for heat conduction with phase change, Int. J. Heat Mass Transf. 55 (2012) 2420–2428. [9] R. Huang, H. Wu, P. Cheng, A new lattice Boltzmann model for solid-liquid phase change, Int. J. Heat Mass Transf. 59 (2013) 295–301. [10] X. Li, J. Zhao, P. Cheng, A lattice Boltzmann model for condensation and freezing of dry saturated vapor about a cryogenic spot on an inclined hydrophobic surface, Int. J. Heat Mass Transf. 114 (2017) 628–639. [11] J. Zhao, X. Li, P. Cheng, Lattice Boltzmann simulation of a droplet impact and freezing on cold surfaces, Int. Commun. Heat Mass Transf. 87 (2017) 175–182. [12] X. Zhao, B. Dong, W. Li, Numerical investigation on freezing process of a single freefalling droplet based on the lattice Boltzmann method, Int. J. Heat Mass Transf. 109 (2017) 807–811. [13] W. Xiong, P. Cheng, Numerical investigation of air entrapment in a molten droplet impacting and solidifying on a cold smooth substrate by 3D lattice Boltzmann method, Int. J. Heat Mass Transf. 124 (2018) 1262–1274. [14] J. Zhao, P. Cheng, A lattice Boltzmann method for simulating laser cutting of thin metal plates, Int. J. Heat Mass Transf. 110 (2017) 94–103. [15] M.K. Ikeda, P.R. Rao, L.A. Schaefer, A thermal multicomponent lattice Boltzmann model, Comput. Fluids. 101 (2014) 250–262.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Technical Note
A multicomponent multiphase enthalpy-based lattice Boltzmann method for droplet solidification on cold surface with different wettability Peng Xu, Sichuan Xu ⇑, Yuan Gao, Pengcheng Liu Tongji University, School of Automotive Studies, Caoan Road 4800, Shanghai 201804, PR China
a r t i c l e
i n f o
Article history: Received 8 May 2018 Received in revised form 2 July 2018 Accepted 4 July 2018
Keywords: Lattice Boltzmann method Droplet Solidification Multicomponent multiphase
a b s t r a c t In this paper, a liquid droplet solidification contacting cold solid surface is modeled by a novel multicomponent multiphase enthalpy-based lattice Boltzmann method (LBM). The freezing process of liquid droplet coupled with thermal and mass transfer on cold flat surface with different wettability (h = 160°, 135°, 60°, 20°) is investigated in detail. The solid volume fraction of droplet solidification with Fourier number are shown under different contact angles. The simulation results show that the solidification time for droplet contacting on a hydrophilic surface is much smaller than that of hydrophobic surface when liquid droplet shares the same initial radius. Ó 2018 Published by Elsevier Ltd.
1. Introduction The solidification process is a common phenomenon in many physical fields, such as air conditioning, phase change energy storage, defrosting, and etc. It is a complicated multiphase flow problem coupled with heat and mass transfer and phase change, especially for multicomponent multiphase system. Over the past several decades, many scientists have studied the liquid-solid phase change process for water droplet either by experimental observation or by macroscopic numerical analysis [1,2]. The LBM has recently gained much attention as a powerful tool to simulate complex physical problem due to advantages of capacity for investigating complicated geometries, simple implementation, high computation efficiency, easy implementation of parallel-processing. It is difficult to determine the phase change interface by conventional computational fluid dynamics (CFD) because of complicated heat and mass transfer mechanism with nonlinear solid-liquid interface characteristics. Three types of liquid-solid phase change LBM models have been developed: immersed boundary method [3], phase-field method [4] and enthalpy-based method [5–13]. Among them, the enthalpy-based LBM model is employed widely due to its simplicity and effectiveness for solid-liquid phase change problem with assumption of no solid-phase movement. Jiaung et al. [5] proposed the enthalpy-based LBM model to simulate phase change problem ⇑ Corresponding author. E-mail address: [email protected] (S. Xu). https://doi.org/10.1016/j.ijheatmasstransfer.2018.07.017 0017-9310/Ó 2018 Published by Elsevier Ltd.
with heat conduction for the first time. Subsequently, Chatterjee and Chakraborty [6] successfully applied enthalpy updating procedure and one relaxation factor for convergence of iteration. Huber et al. [7] investigated convection-dominated solid-liquid phase change in a cavity by setting number of iterations to be one for calculation efficiency. In this regard, Eshraghi and Felicelli [8] introduced an implicit LBM to handle with the source term, thus avoiding iteration steps. Recently, Huang [9] proposed a new thermal LBM model for solid-liquid phase change by adding latent heat source into transient term, which is characterized by avoiding the iteration process and linear equations. Since 2015, some research focus on the heat and mass transfer process of solidification for droplets based on enthalpy-based LBM [10–14]. Among them, some [10,11] are applied for single component multiphase system while others [12–14] for multicomponent situations. To the best of knowledge, droplet solidification on cold surface with different wettability using multicomponent multiphase LBM has never been studied. In this paper, the novel multicomponent multiphase enthalpy-based LBM model based on Zhao and Cheng [14] is developed for simulating solidification of a liquid droplet surrounded by air and contacting a cold solid surface.
2. Description of enthalpy-based multicomponent multiphase LBM model In this section, we briefly present the new enthalpy-based multicomponent multiphase LBM model for liquid-solid phase change.
137
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
The main characteristics for new LBM model are as follows: (1) Triple distribution functions were used, two for velocity field and one for temperature field; (2) Iteration steps and a group of linear equations are avoided by modifying the equilibrium distribution function of the temperature field; (3) The velocity field and temperature field are linked each other through equilibrium function and buoyancy source term; (4) The bounce-back boundary condition for liquid-solid interface is in form of an additional collision term in LBM equation for velocity field. 2.1. Velocity field The evolution for velocity f i with immersed moving boundary scheme and force form of exact difference method (EDM) is described as [14]: k
k
k
f i ðx þ ei Dt; t þ DtÞ ¼ f i ðx; tÞ kðsÞ
þ Bk Xi
1B
sk
k
kðeqÞ
½f i ðx; tÞ f i
ðx; tÞ
k
þ ð1 Bk ÞDf i ðx; tÞ
ð1Þ
where f i is the density distribution function at node x and time t, ei is the discrete velocity in the i direction, k is the component indicator (1 for liquid phase and 2 for gas phase), the D2Q9 discretization scheme for 2D geometry is used in our model, sk is the relaxation time for component k, calculated by:
sk
tk
1 ¼ 2 þ c s Dt 2
pffiffiffi equaling to c= 3. Bk is a weighting function that relates to the dimensionless relaxation time and liquid volume fraction, represented by:
Bk ¼
ð1 f l Þðsk 0:5Þ f l þ sk 0:5
ð3Þ
where the liquid fraction f l is calculated by:
fl ¼
8 > <0 > :
H 6 Hs
HHs Hl Hs
ð4Þ
H P Hl
with Hs and Hl corresponding to the enthalpies of solidus and liquidus temperature, respectively. kðeqÞ fi
denotes the equilibrium distribution function for the velocity and given by:
"
kðeqÞ
fi
¼ xi qk 1 þ
2
ei ub ðei ub Þ þ c2s 2c4s
u2b 2c2s
#
ð5Þ kðsÞ
where wi is the weighting coefficient, Xi is an additional collision term bouncing back the non-equilibrium part of the distribution function, represented as:
XikðsÞ ¼ f ki ðx; tÞ f ki ðx; tÞ þ f ikðeqÞ ðqk ; us Þ f kðeqÞ ðqk ; ub Þ i
ð6Þ
where us is the velocity of solid phase and i is the opposite direction of i. k Df i
is the force term of EDM proposed by Kupershtokh and Medvedev: k Df i ðx; tÞ
¼
kðeqÞ fi ð
qk ; ub þ Duk Þ
qk ; ub Þ
kðeqÞ fi ð
ð7Þ
where Duk is the velocity difference due to total force and can be given as:
Duk ¼
F k Dt
qk
F k2 ¼ wk ðqk ðxÞÞ
s XX Gkk ðx; x0 Þwk ðqk ðx0 ÞÞðx0 xÞ
ð8Þ
ð9Þ
k
x0
where wk ðqk ðxÞÞ is the effective mass, which is a function of the local density. The density value is used for effective mass in this paper for represent two different fluid components. simplicity. k and k 0 Gkk ðx; x Þ is a coupling constant depending on a Green’s function. The fluid-solid interaction force is given by:
F k2 ¼ wk ðqk ðxÞÞ
X Wðx; x0 Þsðx0 Þðx0 xÞ
ð10Þ
x0
where sðx0 Þ is a none-zero constant at the fluid-solid interface and zero otherwise, Wðx; x0 Þ represents fluid-solid interaction strength and different levels of wettability can be obtained by adjusting the parameter Wðx; x0 Þ. The macroscopic density, bulk velocity and real physical velocity can be obtained by:
qk ¼
X k X X X Dt f i ; ub ¼ qk uk = qk ; ur ¼ qk uk þ F k 2 i k k k
ð11Þ
2.2. Temperature field The passive scalar method [15] is the extensively applied approach for incorporating thermal effects into LBM, especially for multicomponent multiphase problem. In this work, we adopted the passive scalar method by using single LBM temperature equation. The evolution equation for temperature field is given by:
g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ
; Hs < H < Hl
1
and the gravitational force F k3 . The pseudo-potential was applied to represent the particle interaction force:
ð2Þ
tk is the kinematic viscosity, cs is the lattice sound speed
where
where F k is the total force, composed of three parts, including the particle interaction force F k1 , the fluid-solid interaction force F k2
1
sT
½g i ðx; tÞ g eq i ðx; tÞ
ð12Þ
where sT is the relaxation time for temperature field which is given below:
sT ¼
a
c2s Dt
þ
1 2
ð13Þ
where a is thermal diffusivity and defined as a ¼ k=ðC p qÞ (k and C p are thermal conductivity and specific heat capacity at constant pressure, respectively). g eq i ðx; tÞ is the temperature equilibrium distribution function, given as:
g eq i
8 u2 i¼0 < H C p T þ xi C p Tð1 2c 2Þ s h i ¼ 2 2 e u ðe uÞ : xi q 1 þ i þ i u i–0 2c4 2c2 c2 s
s
ð14Þ
s
P where H is the enthalpy by H ¼ ig . The temperature can be obtained as follows:
8 H > T HCs p;s H 6 Hs > > s < Hl H HHs T þ T ; T ¼ Hl Hs s Hl Hs l Hs < H < Hl > > > : T l þ HHl H P Hl C
ð15Þ
p;l
2.3. Computational domain and boundary condition The 2D computational domain with the size of N x N x ¼ 150 150 in lattice unit is used in this work. One droplet with
138
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
the designated radius of 25 lattice units is contacting the cold flat surface, immersing in air. For the velocity field, half-way bounce back conditions are applied on the top and bottom walls with periodic boundary conditions imposed for left and right sides of computation domain. As for the temperature field, the non-equilibrium extrapolation scheme is utilized for four sides. It is noted that dimensionless temperature of top and bottom walls are kept at constant (T ¼ 0) while adiabatic thermal boundary is set for left and right sides. The temperature of droplet is higher than that of cold flat surface and melting point so that the droplet can keep in the form of liquid shape before the interfere of temperature effects. Four cases of different contact angles for droplet interacting with cold solid surface is investigated: h = 160°, 135°, 60°, 20°. The temperature field is activated after velocity field reaches equilibrium state. 3. Results and discussion In our work, the problem for droplet solidification above cold flat surface can be characterized by some dimensionless parameters: Ohnesorge number, Fourier number, Prandtl number, dimensionless temperature number, Stefan number, which are defined by
Oh ¼
ld
ðqd rDd Þ0:5
; Fo ¼
at D2d
; Pr ¼
v
a
; T ¼
T Tm CpT c ; Ste ¼ Td Tm L
ð16Þ
where ud is the liquid dynamic viscosity, qd is the droplet density, r is surface tension, Dd is the droplet radius, a is the thermal diffusivity of ice, m is the kinematic viscosity, T m is the melting temperature, T d is the initial droplet temperature, C p is the specific heat at constant pressure, T c is the characteristic temperature and set T c ¼ T d T m , L is the latent heat. As mentioned before, the equilibrium state of droplet contacting solid surface is preferred followed by intervention of thermal effects. Therefore, the parameter affecting the wettability of solid
surface must be confirmed in advance. The contact angle between the fluid and solid surface can be calculated by Young’s equation. Three parameters need to be known including interfacial tension between two fluid components r12 and between each fluid component and the solid surface rs1 , rs2 , as shown in Eq. (17):
cos h ¼ ðrs2 rs1 Þ=r12 ¼ ðg2w g 1w Þ=ðgðq1 q2 Þ=2Þ
ð17Þ
According to Eq. (17), when the value of g 2w between gas phase and solid wall is greater than the value of g 1w between liquid phase and solid wall, the contact angle h is less than 90 , showing the solid surface is hydrophilic and vice versa. The fluid-fluid interaction coefficient is kept at 0.1 with g ¼ g 12 ¼ g 21 while interaction strength between each fluid and the solid surface varies from 0.05 to 0.05. In addition, fluid-solid interaction strength value is kept as g 1w ¼ g 2w . It is noted that all the units are in lattice units unless otherwise stated specifically. The contact angle can be evaluated by final steady state values of the droplet with height a, droplet length b and droplet radius R: tan h ¼ b=ðaðR aÞÞ. The contact angle between wetting phase and solid surface varies from 7° to 173° when g 1w changes from 0.05 to 0.05. Different level of wetting (non-wetting) conditions can be obtained by adjusting g 1w in the following simulation. Fig. 1 presents 2D simulation results of solidification process for droplet above cold solid surface with different contact angle (h = 160°, 135°, 60°, 20°). The size of whole domain is 150 150 lattice units with initial radius of droplet 25 lattice units. In the simulation, the Ohnesorge number is set to be Oh = 0.13 and Stefan number is Oh = 0.25. Initially, the temperature of droplet and gas is chosen to be T d = 2.0 (which is higher than that of melting temperature point T = 0.5) and T g = 0.0 respectively. The top and bottom wall of calculation domain is set to be T t = 0.0 and T b = 0.0. Other parameters are as follows: C p = 1.0, L = 0.25, viscosity ratio for gas and droplet v g : v l = 1:1, Prandtl number Pr = 0.71 and the thermal diffusivity k = 0.235. The red region denotes the ice region
Fig. 1. Shape evolution of droplet on surface at different wettability with time.
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
and the blue region inside the droplet denotes the liquid state. When the dimensionless temperature of droplet contacting the solid surface drops below T = 0.5, the part of droplet contacting solid surface starts to form ice firstly. The top area of droplet also ices because temperature of gas close to the droplet is low. Note that the droplet area that ices lastly is not the top point of droplet but high area inside droplet. When contact angle is h ¼ 160 , the solidification time is F o = 5.6, which is bigger than that of case with contact angle h ¼ 20 . With the reducing of contact angle (more hydrophilic), the solidification time reduces quickly.
139
To investigate the dimensionless temperature variation along the vertical direction cross the droplet, the results extracted from line ‘‘CD” as indicated in Fig. 2(a) are obtained in Fig. 2. Fig. 2(a) shows the dimensionless temperature contour at F o = 0.4 when droplet forms a contact angle h = 160° with solid surface. The dashed line is the shape of the liquid droplet. The heat transferred from droplet surface to ambient gas leads to temperature gradient and the dimensionless temperature difference between droplet and cold surface changes more quickly because of larger heat flux near the cold surface. To validate this expression, the local heat flux near cold surface denoted by ‘‘EF” line in Fig. 2(a) is calculated [11]. From Fig. 2(c), we can see that heat flux value between droplet and cold surface is almost 8 times larger than that between gas and cold surface. Fig. 2(b) shows the dimensionless temperature along the vertical line ‘‘CD”. At F o = 0.4, the dimensionless temperature of point ‘‘D” is 0.35 larger than that in ‘‘C” point. It is noted that 0.35 difference for dimensionless temperature can be seen during the whole solidification process. In addition, the dimensionless temperature near cold surface decreases quicker than far area due to almost constant low temperature boundary condition seen from two arrows in Fig. 2(b). After Fourier number reaches to F o = 5.6, the temperature profile across the droplet behaves linearly.
4. Conclusions A novel multicomponent multiphase enthalpy-based LBM model for droplet solidification contacting cold solid surface was proposed in this paper. The solidification process of droplet contacting wall with different wettability (h = 160°, 135°, 60°, 20°) is investigated through developed multicomponent multiphase enthalpy-based LBM model. The following conclusions can be made from this paper: (1) The dimensionless temperature in droplet drops quickly when droplet is placed above cold flat surface. Whereas the temperature drops much slower when droplet is during solidification process because of phase change. (2) The dimensionless temperature in droplet near the contact line between cold solid surface and droplet drops more quickly than that near droplet surface. The largest temperature difference around 0.4 exist in droplet during solidification when contact angel is h = 160° while 0.34 in droplet between the zenith and the lowest point. (3) For a droplet with same initial radius, the solidification time for hydrophobic cold surface is almost 5 times larger than that for hydrophilic surface. The freezing time for contact angle (h = 160°, 135°, 60°, 20°) is 5.6, 4.9, 3.0 and 1.1, respectively.
Conflict of interest We have no conflict of interest. We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.
Acknowledgements Fig. 2. The dimensionless temperature contour (a), the dimensionless temperature profile along the vertical diameter of droplet (b) and heat flux distribution (c) along the cold plate at F o = 0.4 with h = 160°.
This work was supported by National Key R&D Plan of China through Grant Nos. 2017YFB0102802.
140
P. Xu et al. / International Journal of Heat and Mass Transfer 127 (2018) 136–140
References [1] X. Wu, W. Dai, W. Xu, L. Tang, Mesoscale investigation of frost formation on a cold surface, Exp. Therm. Fluid. 31 (2007) 1043–1048. [2] J.P. Hindmarsh, A.B. Russell, X.D. Chen, Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet, Int. J. Heat Mass Transf. 46 (2003) 1199–1213. [3] R. Huang, H.W., An immersed boundary-thermal lattice Boltzmann method for solid-liquid phase change, J. Comput. Phys. 277 (2014) 305–319. [4] W. Miller, S. Succi, D. Mansutti, Lattice Boltzmann model for anisotropic liquidsolid phase transition, Phys. Rev. E 86 (2001) 3578–3581. [5] W. Jiaung, J. Ho, C. Kuo, Lattice Boltzmann method for the heat conduction problem with phase change, Numer. Heat Transf. Part B: Fundam. 39 (2) (2001) 167–187. [6] D. Chatterjee, S. Chakraborty, An enthalpy-based lattice Boltzmann model for diffusion dominated solid–liquid phase transformation, Phys. Lett. A 341 (2005) 320–330. [7] C. Huber, A. Parmigiani, B. Chopard, M. Manga, Lattice Boltzmann model for melting with natural convection, Int. J. Heat Fluid Flow 29 (5) (2009) 1469– 1480.
[8] M. Eshraghi, S.D. Felicelli, An implicit lattice Boltzmann model for heat conduction with phase change, Int. J. Heat Mass Transf. 55 (2012) 2420–2428. [9] R. Huang, H. Wu, P. Cheng, A new lattice Boltzmann model for solid-liquid phase change, Int. J. Heat Mass Transf. 59 (2013) 295–301. [10] X. Li, J. Zhao, P. Cheng, A lattice Boltzmann model for condensation and freezing of dry saturated vapor about a cryogenic spot on an inclined hydrophobic surface, Int. J. Heat Mass Transf. 114 (2017) 628–639. [11] J. Zhao, X. Li, P. Cheng, Lattice Boltzmann simulation of a droplet impact and freezing on cold surfaces, Int. Commun. Heat Mass Transf. 87 (2017) 175–182. [12] X. Zhao, B. Dong, W. Li, Numerical investigation on freezing process of a single freefalling droplet based on the lattice Boltzmann method, Int. J. Heat Mass Transf. 109 (2017) 807–811. [13] W. Xiong, P. Cheng, Numerical investigation of air entrapment in a molten droplet impacting and solidifying on a cold smooth substrate by 3D lattice Boltzmann method, Int. J. Heat Mass Transf. 124 (2018) 1262–1274. [14] J. Zhao, P. Cheng, A lattice Boltzmann method for simulating laser cutting of thin metal plates, Int. J. Heat Mass Transf. 110 (2017) 94–103. [15] M.K. Ikeda, P.R. Rao, L.A. Schaefer, A thermal multicomponent lattice Boltzmann model, Comput. Fluids. 101 (2014) 250–262.