Numerical study of the effect of forced convective flow on dropwise condensation by thermal LBM simulation

The enhancement mechanism of forced convective flow on dropwise condensation over a cold spot is numerically investigated by two-dimensional hybrid thermal lattice Boltzmann (LB) model based on the Shan-Chen pseudopotential LB model. After validating the present LB model, dropwise condensation over a cold spot as the nucleation region is simulated. The well-known power law for the growth of a single condensing droplet is demonstrated. Finally, the simulation of dropwise condensation considering the convection flow or not is carried out in the constant contact radius (CCR) mode. Using the CCR model, the effect of contact angle can be also investigated. The result of streamline field indicates that the forced convectional flow complicates the internal flow of droplet and main flow. The dragging force from main flow changes the size of two symmetric vortices inside the droplet. And the channel flow is also strongly influenced by the suction effect caused by condensation at the three phase contact line. By comparison, the heat transfer enhancement of the superimposed flow is not worth mentioning. The present study illustrates the mechanisms of dropwise condensation under forced convectional flow.


Introduction
Because of the complicated dynamic behavior of droplets, dropwise condensation becomes to be a typical case of multiphase flow involving vapor-liquid phase change, free surface flow, coupled heat and mass transfer, multi-scale feature, and interfacial behaviors. Therefore, even though various experimental and theoretical investigations have been carried out, dropwise condensation has been rarely studied numerically [1,2].
Over the past two decades, compared to the conventional CFD based multiphase models, the LB model has been developed into a powerful simulation method in modeling multiphase flow undergoing phase change. Theories and applications of LB model for multiphase flow and phase-change heat transfer have been deeply reviewed and summarized by Li and his coworkers [3]. For multiphase flow, a popular model used in LB model is the pseudopotential LB method proposed by Shan and Chen [4,5] because of its computational efficiency and straightforward algorithm. In this model, the phase separation is realized spontaneously by the introduced short-range intermolecular interaction based on a potential function. To alleviate the weaknesses of the original Shan-Chen model, several improvements [6][7][8][9][10][11] have been carried out. For phase-change heat transfer, the thermal pseudopotential LB model is developed by devising the temperature distribution function [12][13][14][15][16]. The temperature distribution function is used to solve the macroscopic energy equation. The hybrid thermal LB model, which solves the temperature equation by the finite-difference scheme, is devised for avoiding the error term in energy equation recovered. Using the hybrid thermal LB model, Li and his coauthors have made several successful works [17][18].
Using the hybrid thermal LB model, the present paper is focused on the coupling effect between the fully developed laminar flow and dropwise condensation. After validating the present LB model, dropwise condensation considering the convective flow or not are simulated using the same properties. The main objective is to understand the mechanism of dropwise condensation in the forced convective laminar flow.

Model description
The multiphase single-component pseudopotential LB model proposed by Shan and Chen [4,5] is employed to simulate the multiphase flow for the density and velocity fields. The energy equation is discretized by the finite difference scheme for the temperature field.

The multiphase single-component pseudopotential LB model
Based on the Bhatnagar-Gross-Krook (BGK) collision operator [19], the discretized LB equation is given below for the flow fields: (1) where fα(x,t) denotes the particle distribution function, x is the particle position, t is the time, eα represents the discrete velocity along the α-th direction, τf is called as the non-dimensional relaxation time, the equilibrium density distribution function fα eq (x,t) can be written as below: (2) where wα and cs are the weight coefficients and lattice sound speed, respectively. In this work, two-dimensional 9-velocity scheme (D2Q9) is chosen. The lattice quantities (density ρ and velocity u) are given as: ff        ue (3) In Eq. (1), ∆fα(x,t) is the body force term. In LB model, the hydrodynamic equations are correctly recovered by incorporating the body force. According to the former works [8,20], the exact difference method (EDM) developed by Kupershtokh et al. [21] shows better stability and accuracy. Therefore, the EDM forcing scheme is used in this work: (4) where ∆u=Fδt/ρ is the effect of body force in velocity during each time step. For EDM forcing scheme, the actual fluid velocity v is defined by: 2 t    F vu (5) In pseudopotential LB model, the short-range intermolecular interaction force is significant and responsible for the phase transition. The mixed scheme for the interparticle interaction force is used in this work: x + e x + e e (6) where β is the weighting factor depending on the nonideal equation of state. The Green function G(x+eα) denotes the interaction strength, which can be found in [20]. The pseudopotential ψ(x) is determined by the local density and equation of state and given by: The pseudopotential ψ(x) is linked to the temperature fields through the non-ideal equation of state. The Peng-Robinson (P-R) equation of state is adopted in this work: where a = 0.45724R 2 Tc 2 /pc and b = 0.0778RTc/pc, φ(T) can be found in [6]. In this work, we set a = 2/49, b = 2/21 and R=1. In present work, the saturation temperature is set as Tsat = 0.86Tc, which has the coexistence liquid density ρl = 6.5 and the coexistence vapor density ρv = 0.38. After numerical test, the weighting factor in Eq. (6) is used as β = 1.18 which can give stable and accurate results in a wide temperature range. It should be mentioned that the lattice unit (l.u.) is used for all quantities here. And the gravity is free.

Energy equation
After neglecting the viscous heat dissipation, the governing equation for the temperature field can be derived from the local balance law of entropy [17]: (9) with the thermal conductivity λ and the specific heat at constant volume cv. The second-order Runge-Kutta method is employed to implement the time discretization: where K(T) is the right-side term of Eq. (9). Based on the second-order isotropic difference schemes, the spatial gradient and Laplacian of and arbitrary quantity ϕ can be calculated as: x e e (12) x e x (13)

Model validations
To validate the capability of present model for simulating vapor-liquid phase change, several tests are carried out.

Phase transition
The 2D phase transition process at a constant temperature is simulated. A square domain (150 × 150 l.u.) with all periodic boundary conditions is initialized with the local density ρ = ρc+∆ρ, where ∆ρ is a small random perturbation parameter. During simulation, the phase transition is achieved automatically and the droplet arises after reaching the equilibrium state. Then the corresponding vapor and liquid densities can be determined. Over a wide temperature range, the liquid and vapor densities from the simulation show a good agreement with the theoretical Maxwell construction, as shown in Fig. 1.

Laplace test
For Laplace test, a droplet with radius R is initialized at the center of the computational domain having all periodic boundary conditions. Once the simulation reaches to the equilibrium state, the pressure difference ∆p across the droplet interface can be quantified.

Wettabilit test
The equilibrium contact angle (wetting condition) on solid wall is implemented through the geometric formulation proposed by Ding and Spelt [22]. Using the second-order central difference scheme, the discretized form for numerical application is written as: with the prescribed static contact angle θw a , ρi,0 denotes the density of the ghost layer (i,0) beneath the boundary lattices, where the first and second indexes represent the coordinate tangential and normal to the boundary. For θw a = 60 º and 120 º , the numerical simulation obtain the equilibrium contact angle θw = 60.6º and 121.5º respectively. Fig. 3 shows the equilibrium shape of liquid droplet over the solid wall.

2D droplet evaporation
Evaporation of a stationary droplet is considered to verify the accuracy of the present model in simulating liquid-vapor phase change. For droplet evaporation, the classic D2 law states that the square of the droplet radius reduces linearly over time. The simulation is examined on a 200 l.u. × 200 l.u. computational domain with a droplet (radius R = 30 l.u.) placed at the center. At the initial time step, the droplet temperature is set at the saturation temperature 0.86Tc and a uniform temperature 1.0Tc is given to the vapor surrounding the droplet. All boundaries are employed at the periodic boundary condition and constant temperature condition (1.0Tc).
The specific heat at a constant volume is cv = 5.0. According to the D2 law, the thermal conductivity λ should be constant, and the evaporation rate depends linearly on λ. Two cases are carried out: Case A with λ = 2/3, Case B with λ = 3/3. Fig. 4 plots the square of the droplet radius normalized by the initial radius during the evaporation process for these two cases. As expected, the linear relationship between (R/R0) 2 and t is clearly observed for two cases, and the increasing λ accelerates

Dropwise condensation on plate
After validating the present model, the growth process of a single droplet is simulated in a rectangle domain (Lx × Ly = 800 × 200 l.u.). The thermal conductivity is taken as λ = ρcvχ with cv = 5.0 and χ = 0.04. The kinematic viscosity υ is set to 0.05. The whole domain is initialized as the vapor having the saturation temperature Tsat. A cold spot of 30 l.u. as a nucleation region is placed at the center of bottom wall and is fixed at Tw = 0.7Tc. The rest of the bottom wall is at Tw = 0.83Tc. The periodic boundary condition is applied in the x-direction. For the bottom surface, the non-slip boundary condition including the forcing term is used. And the convective boundary condition proposed by Lou et al. [23] is specified at the top boundary. The prescribed static contact angle is θw a = 90º and the contact angle hysteresis is considered with a hysteresis window (0º, 90º) [24]. For showing the heat transfer during condensation, the variation of the averaged heat flux over the cold spot with time is plotted in Fig. 5. At the start of condensation, the heat flux increases dramatically because of nucleation. Once the first droplet forms, the increasing conduction resistance of droplet decreases the heat flux. Finally, the heat flux goes to a stable state. It can be concluded that the heat transfer of dropwise condensation is mainly attributed by the small size droplet. Additionally, Fig. 6 presents the time evaluation of droplet radius during condensation. As known, a single droplet growth in condensation can be described by the power law (R ~ t a ). Based on a theoretical analysis, the power law exponent a = 0.5 for a 2D droplet growth has been demonstrated by Ashrafi and Moosavi [25]. As seen in Fig. 6, the present simulated droplet radius obey a power law with a = 0.407. The difference with the theoretical value falls into the normal difference range presented in reference [25]. Following the parabolic velocity profile along ydirection, the velocity inlet boundary condition is implemented. As before, the convective outlet boundary condition is employed at the right boundary. These two cases have the same thermal boundary condition, which set the saturation temperature Tsat for all boundaries in spite of a cold spot 80 l.u. located in the center of bottom wall. As the nucleation region, the temperature of cold spot is set as Tw, and the sub-cooling ∆T = Tsat -Tw is chosen as 0.16Tc. For Case B, the parabolic velocity profile is initialized on the whole computational domain, and the first 20 000 steps are carried out without condensation so that the velocity distribution can reach its equilibrium state. At the beginning of condensation, the wettability effect is not added so the liquid can cover the cold spot freely. When the whole cold spot is covered by the liquid film, the wettability effect is taken into consideration with a hysteresis window of (0º, 180º) for contact angle. Therefore, the droplet grows in the constant contact radius (CCR) mode, which implies that the contact angle increases gradually whereas the contact line is pinned on the solid wall. Using this way, the effect of contact angle can be also researched and the droplet basis can always keep at the constant temperature Tw. In this work, the averaged inlet velocity is U = 0.025. Thus, the Reynolds number is Re = ULy/υ = 100 which corresponds to the laminar flow. Fig. 7 gives the evolution of droplet profile during condensation for two cases. In Case A, the contact angle is increased owing to the constant contact radius. Due to the shearing effect of external flow, the deformation of droplet can be observed and the deformation is enhanced with the growth of droplet.
In order to deeply understand the mechanism of dropwise condensation, it is necessary to study the fluid dynamics. Fig. 8 shows the streamline field inside and outside of the condensing droplet for two cases. As shown in the left side, three different contact angles (CA) for Case A are presented, which is corresponded by Case B in right side at the same condensation time. For Case A, the x-coordinate of showing region span from 320 to 480. In order to display the flow field in the channel for Case B, the range from 310 to 560 is used. That is why the droplet in right side looks smaller than the left side although the droplet size from two cases keeps close at the same condensing time. From the left side, two symmetric vortices (namely, Marangoni flow) inside the droplet can be clearly found in different CA because of the temperature gradient. Also it can be seen that most of the vapor flow goes to the three phase contact line and condensate. However, two symmetric vortices are replaced by two different size vortices in Case B. The left vortex is narrowed by the flow dragging force having the opposite direction, while the same direction of the flow dragging force expand the right vortex. An interesting phenomenon is the suction effect caused by the condensation at the three phase contact line. The pressure difference induced by the suction effect forces the main flow going to the bottom wall, and the main flow is strongly compressed under a big size droplet. At the small contact angle CA = 69.4º, the strong suction effect brings out a strong compression for the main flow. Additionally, the lee-side vortex structure is also replaced by the suction flow in the back region of droplet, and the region of suction flow is enlarged with the increasing of droplet height. About the effect of convection flow in heat transfer, Fig. 9 and Fig. 10 present the temperature distribution and the local heat flux over the cold spot. From Fig. 9, it looks that the convection flow has no influence for the heat transfer inside the droplet. Also it can be found that the maximum temperature gradient locates at the three phase contact line, which can also explain that three phase contact line is the main place for condensation. From the heat flux shown in Fig. 10, the heat flux is decreased with increasing the contact angle, which is attributed to the variation of conduction resistance. Comparing the local heat flux for two cases, the heat transfer is slightly enhanced by the forced convectional flow. Therefore, it can conclude that the effect of forced convectional flow for heat transfer can be neglected compared with the latent heat released by condensation.

Conclusions
In present work, the 2D hybrid thermal LB model is used to simulate dropwise condensation. The fluid flow and temperature field is solved by the pseudopotential LB model and finite difference scheme, respectively. In order to avoid the drawbacks of the original pseudopotential LB model, some recent improvements from literature are utilized in the present work. For validation, some tests are also carried out.
Based on the CCR model, the effect of forced convective flow for dropwise condensation in different contact angles is investigated. For two symmetric vortices inside the droplet in dropwise condensation without convection flow, superimposed flow narrows the left vortex while enlarges the right one. The maximum temperature gradient can explain that most of the vapor flow goes to the three phase contact line and condensates. Condensation at three phase contact line brings out a suction effect. The suction effect is responsible for the flowing down of the main flow. Compared with the latent heat of condensation, the contribution of forced convective flow for heat transfer is no worth of mention. Therefore, the conclusion is that the superimposed convective flow complicates the fluid dynamics for the main flow and internal flow of droplet whereas the enhancement effect for heat transfer is trivial.