Thermal Lattice Boltzmann Simulation of Evaporating Thin Liquid Film for Vapor Generation

Thin film evaporation (TFE) plays an important role in many industrial applications, such as power generation, cooling, and thermal management. Effective evaporation takes place in the thin liquid film region with relatively low film thickness and low intermolecular forces. In this paper, a numerical approach based on the thermal lattice Boltzmann method (TLBM) is employed to investigate the heat and mass transfer phenomena in TFE. The TLBM approach is validated by simulating some benchmark problems, and is then used to study a vapor generation problem where TFE is involved. Specifically, vapor is generated from evaporating pores, the solid walls of which are hydrophilic. Factors that affect the overall vapor generation efficiency are investigated via the numerical approach. Methods that can improve the overall efficiency are further proposed. Simulations reveal that distributed scenarios (using distributed small pores instead of a big one) and hydrophobic pore ends render more efficient vapor generation.


Introduction
The fast development of microelectronic devices requires advanced thermal management.For devices with large heat flux such as integrated circuits and concentrator solar cells, highly efficient heat transfer approaches are desired.One of the promising solutions that has gained much attention is thin film evaporation (TFE), which can accommodate large heat fluxes owing to liquid-to-vapor phase change.At the contact line between a wetting liquid and solid, there is a thin film region.When the film is ultra thin (usually on the order of 10 nm), the thermal resistance is very small.However, the intermolecular force between liquid and solid is very strong, thus almost no evaporation takes place.This region is called the adhesion film region.On the other hand, the film thickness is relatively large in the intrinsic meniscus region, which results in large thermal resistance and weak evaporation.Hence, most of the evaporation takes place in the region with relatively low thickness and low intermolecular forces.The schematic of evaporating thin film is demonstrated in Figure 1.
Many researchers so far have conducted relevant studies through both theoretical analysis [1][2][3][4] and experimental approaches [5][6][7].From the theoretical point of view, it is troublesome to solve a sophisticated model that involves complicated partial differential equations.However, some assumptions should be made if a simplified problem is considered instead, which weakens the physical meaning.Since TFE is usually considered in micro-scale and even nano-scale devices, measuring the physical quantities, such as the thickness of liquid film and the local temperature, via experimental methods is also challenging.Besides, it is extremely difficult to go deep into the interfacial phenomena where evaporation takes place.Hence, it is imperative to investigate TFE using simulation approaches.Note that TFE is essentially a multiphase flow problem with heat transfer.There are many approaches that can be used in fluid flow simulations, such as the finite volume method [8], smoothed particle hydrodynamics [9], molecular dynamics (MD) [10], lattice Boltzmann method (LBM) [11], etc. Comparisons of various approaches can be found in [12,13].In our study of TFE, LBM is mainly considered, which is a mesoscopic computational fluid dynamics approach.It is known that LBM is good at tackling fluid-solid interaction and complex geometries.
Besides, LBM employs a simple algorithm (collision and streaming), which makes it intrinsically scalable and parallelizable for computations.LBM has been used in many engineering and scientific problems, e.g., [14][15][16][17][18][19][20].In this paper, we will show that LBM, and more specifically thermal LBM (TLBM), can reproduce several important phenomena in TFE, even without very fine mesh grids.Hence, it facilitates our understanding of the physics behind TFE and further device-level design.In the last decade, TLBM has been widely investigated [21].There are several TLBM approaches, which can be classified as the multispeed method [22], the double distribution method [23], and the hybrid method [24].A more detailed discussion of various approaches can be found in the excellent review paper [25].It is noticed that TLBM has been applied in the study of many thermal flow problems, e.g., natural convection [26,27], and boiling [28][29][30].More recently, the authors of [31] reported results on TFE at a hydrophilic microstructured surface using a double distribution TLBM approach.The surface consists of micro pillars.The influence of wettability and heat flux on liquid film distribution and liquid-vapor interface morphology were investigated.In addition, the dependence of dry-out heat flux and wettability was discussed.In this paper, we also propose a double distribution TLBM approach that is based on a pseudo-potential function.The approach is considered in the study of TFE.We mainly focus on the interfacial phenomena, which may be hardly observed through experimental methods.By introducing some flow and thermal boundary conditions, the proposed method is further applied in the investigation of a practical problem, which is direct vapor generation.Specifically, vapor is generated in evaporating pores where heat is concentrated.TFE is involved in the problem due to the fact that the solid walls of pores are hydrophilic.Factors that affect the overall efficiency of vapor generation are investigated, including the pore distribution and wettability of the solid walls.With the help of TLBM simulations, a scenario that can achieve higher vapor generation efficiency is proposed.As a result, the TLBM approach may facilitate our understanding of the phenomena that relate to TFE, and our design of devices where TFE is involved.Before presenting the main results, the TLBM approach together with some validations are discussed first in the following section.

Double Distribution TLBM
The model we consider in our study is the so-called double distribution thermal lattice Boltzmann model (TLBM).It was firstly proposed in [23], where density and energy are considered in different distribution functions.A TLBM approach was further developed in [32] with quantitative analysis, which also inspires the modeling employed in this work.

Density Distribution Function
The D2Q9 lattice Boltzmann (LB) model has been widely used in simulations of many complex flows [33].It consists of two parts, namely streaming and collision.A D2Q9 model with Bhatnagar-Gross-Krook approximation of the collision term is given by where f i (x, t) is the particle distribution function in the i-th discretized velocity direction at position x and time t, ∆x, ∆t and τ f are the lattice spacing, time step and relaxation time, respectively.The equilibrium distribution function is obtained by and one has with w 0 = 4/9, w 1,2,3,4 = 1/9, w 5,6,7,8 = 1/36, c s = √ 3c the sound speed, c = ∆x/∆t, u = ∑ 8 i=0 f i e i , and the velocity set e i is given by Letting F be an external force, then the macroscopic properties such as density and velocity can be obtained by Furthermore, the kinematic viscosity ν relates to the relaxation time by ν = (τ − 0.5)c 2 s ∆t.The exact difference method shown as above has good performance in simulating multiphase flows due to its simplicity and accuracy [34].Moreover, the unphysical influence of relaxation time can be avoided.In the TFE problem, the external forces F include cohesive force (between fluid nodes), adhesion force (between fluid and solid nodes), and body force (e.g.gravity).They will be discussed in the following.
A single component multiphase LB model is considered to investigate TFE.Numerous results on the multiphase LB model have been reported.Readers may refer to [21,25,35] and the references therein for more detailed information.A pseudopotential approach inspired by [32,36] is considered in this work.The cohesive force is described as the gradient of a psudopotential function, which has the following form where P is the pressure, G and β are two tuning parameters.As studied in [32], the parameter β controls the resolution of liquid-vapor interface.Note that high resolution, however, suffers from numerical instability.The parameter G affects the density ratio of liquid and vapor.The model reduces to the one discussed in [32] as G equals 1.It is worth mentioning that the pseudopotential function can accommodate different equation of states (EOS).In our study, the normalized Van der Waals EOS is considered, that is In the above equation, the dimensionless physical variables are defined as where ρ, P, and T denote the density, pressure, and temperature in physical units, respectively, while ρc , Pc , and Tc correspond to the critical point.Throughout this paper, the normalized variables which are intrinsically dimensionless will be used in the lattice Boltzmann model.The pseudopotential function φ(ρ, T) is now obtained by substituting the pressure expression (7) into Equation (6).The cohesive force between fluid nodes is then expressed as Note that in the computation of F int , a discretized realization of gradient is required.The two-dimensional gradient ∇φ can be approximated by where λ 1,2,3,4 = 4/12 and λ 5,6,7,8 = 1/12 [32].On the other hand, we may define ψ as 2φ.Then ∇φ can also be approximated as Based on Equations ( 10) and (11), an improved method is developed [34], In the above approximation of gradient, A is a weighting factor.For Van der Waals EOS, it is suggested that A = 0.55 gives good performance [34].Note that when the neighbor of a fluid node is solid, the adhesion force should be considered.A convenient way is to assign the density and temperature for solid nodes, denoted as ρ w and T w .Now the cohesive force and the adhesion force can be jointly given by If a body force F g exists, e.g., gravity, then the total force is It is worth mentioning that solid nodes are assigned with ρ w only when we calculate the adhesion force.That is to say, ρ w does not affect the density distribution function in the streaming and collision processes.

Energy Distribution Function
In the double distribution LB model, an energy distribution function is considered independent of the density distribution in order to mimic the energy transport.The following energy distribution function is employed in our study where g i (x, t) is the energy distribution in the i-th discretized direction, and τ g is the relaxation time for energy.The equilibrium distribution g eq i is defined as In the above equations, c v is the heat capacity at constant volume (assumed to be 1 in our study).This model differs from the one proposed in [23], where both viscous dissipation and the compression work done by the pressure were taken into account.The influence of viscous dissipation can be evaluated via the Brinkman number Br, which is defined as Br = νu 2 /(αT 0 − αT), where T 0 and T represent temperature of wall and fluid, respectively.The importance of compression work is reflected by the Mach number, that is Ma = |U| /c s .In the problem of TFE, both Br and Ma are usually very small.Hence, viscous dissipation and compression work are negligible in our study.The resulting energy distribution function is significantly simplified.Furthermore, the second order truncation error can be absorbed in the relaxation time τ f and τ g .It is noticed that for the multiphase flow, the density difference between liquid and vapor phases is large, which deteriorates the numerical stability in the thermal flow simulations.Hence, another variable g is introduced, that is g = g/ρ.Now g is considered in the steaming and collision steps.Furthermore, the equilibrium distribution of g i can easily be extracted from Equations ( 16)- (18).The temperature T is then computed via The thermal diffusivity α is expressed as α = 2c 2 s (τ g − 0.5)∆t.

Boundary Conditions
Boundary conditions have great influence on the numerical simulations.Abundant achievements have been reported for the flow field.The classic bounce-back scheme is a simple treatment for the non-slip boundary condition.Bounce back combined with specular refection can generate slip boundary condition [11].Other popular boundary conditions include Zou and He's non-equilibrium bounce-back method [37] and Guo's extrapolation method [38], etc.The results relating to thermal boundary conditions are relatively fewer.In the following part, we will first discuss the thermal boundary conditions.After that, a replenishment boundary condition is introduced, which will be employed in our study of the vapor generation problem.

Thermal Boundary Conditions
We mainly focus on two types of thermal boundary conditions, namely constant temperature and constant heat flux.The constant temperature boundary condition is considered first.The south wall shown in Figure 2 is taken as an example, the temperature of which is constant, denoted as T 0 .
The energy distributions g 2 , g 5 and g 6 are unknown after the streaming process.It is assumed that g 2 , g 5 and g 6 follow equilibrium distributions.The idea comes from the "counter-slip" approach proposed in [39].Similar consideration can be found in [26,[40][41][42].It is mentioned in [26] that this treatment provides significant numerical stability and accuracy.For the three unknown distributions it holds where T is a virtual temperature and u w is the wall velocity.On the other hand we have 8 which yields Let γ 2 , γ 5 and γ 6 be the weights in Equations ( 20)-( 22) such that Then it holds In a special case that the wall is stationary, we have γ 2 = 1/6 and γ 5 = γ 6 = 1/12.Now we consider the constant heat flux boundary condition at the normal direction of the south wall.The heat flux can be represented as [41] q(x, t) = Assume that the heat flux is q 0 .Then one has where u wn is normal velocity (in the y direction) of the wall.Thus It is also assumed that the unknown ones follow an equilibrium distribution.Similarly, g 2 , g 5 and g 6 can be obtained from Note that the boundary conditions developed above are not restricted to solid nodes.They are also applicable to fluid nodes.Apart from the approaches presented above, some other thermal boundary conditions can also be found in the literature, for example, the non-equilibrium extrapolation method [43].The detailed flow and thermal boundary conditions that are employed in our study will be mentioned in specific problems.

Replenishment Boundary Condition
Without replenishment, the liquid level of an evaporating pore will decrease when the solid walls of the pore are constantly heated.The liquid will dry out eventually.As shown in Figure 3, the mass flux at the vapor outlet also decreases during this process.However, in a practical application that involves TFE, the liquid is usually replenished from liquid bulk via the capillary force.We propose to use a ghost layer outside the liquid inlet to replenish the liquid.In other words, the ghost layer serves as the reservoir to supply the liquid.Denote the total mass in the pore as m p .At each time step in the simulation, the loss mass caused by evaporation is estimated by where m t and m l are the total mass and loss mass, respectively.Then the ghost layer is assigned with the loss mass.Specifically, for each node of the ghost layer, denoted as x gh , the unknown distributions (assume to be f 2 , f 5 and f 6 ) are given by The effectiveness of the proposed strategy is revealed in Figure 4.It is observed that at the steady state, both total mass and evaporating mass flux are constant.At this point, we have derived a TLBM approach that is able to investigate TFE in an evaporating pore with liquid replenishment.

Contact Angle and Density Ratio
In this subsection, we will show that both the contact angle and density ratio can be modified in the model.Note that contact angle and density ratio in multiphase LBM approaches have been intensively studied [21].More recent work can be found as [44], where high density ratio (liquid to vapor) is achieved using the multiphase model of Lee [45].The main difference of our work is that temperature field is included which brings extra difficulties, for example, thermal fluid simulations are more susceptible to the numerical instability.This section shows that our TLBM approach is able to modify the contact angle and the density ratio even when the flow and temperature field are coupled.The contact angle is investigated first, which relates to the surface wettability.By changing the density of solid nodes, a different contact angle is obtained.Specifically, hydrophilic surface and hydrophobic surface are achieved by choosing ρ w that is close to the liquid density and vapor density, respectively.Figure 5 demonstrates how contact angle changes with ρ w , where the temperature T is set as 0.7.It is observed that the contact angle can be modified within a large range.Different from contact angle, density ratio of liquid and vapor can be influenced by several factors, for example, G, α and T. The dependencies of density ratio on G and α are illustrated in Figure 6.It is observed that the density ratio is mainly determined by the parameter G.With the increase of G, density ratio increases significantly.Compared with G, the parameter α contributes slightly to the density ratio.However, it plays a critical role in the sharpness of liquid-vapor interface.A larger α renders sharper interface, but also makes it vulnerable to numerical instability.Thus a moderate α is suggested in the simulations.Another factor that affects the density ratio is the temperature T. When T increases, both liquid and vapor will have smaller density.Note that the vapor is much more sensitive in the density change with respect to the temperature compared with liquid.As a result, the density ratio is also dependent on temperature.As a conclusion, the thermal lattice Boltzmann model used in our study is able to mimic different wettabilities.The solid surface can be hydrophilic or hydrophobic.Moreover, the density ratio of liquid and vapor ranges from less than 10 to several hundreds, depending on the parameters of the model.

Validations of the TLBM Approach
There are many applications in which TFE occurs in microscale pores.In this section we mainly investigate TFE in a pore using the TLBM approach.The structure illustrated in Figure 7 represents a pore in our simulations.The solid walls are hydrophilic, which ensures that there is a thin-film evaporating region.TFE occurs when the solid walls are heated.We will investigate the thin-film evaporating pore in several aspects, including validations from both reported experiments and theories.

Wall Temperature Distribution
It is known that heat can be removed through evaporation.Since stronger evaporation occurs in the thin-film evaporating region, a local cooling effect can be expected.Based on this idea, the authors of [5] investigated the wall-temperature distribution of an evaporating micro channel.The wall was heated at constant power in their experiments.The structure shown in Figure 7 is employed to simplify the problem.The domain size in the computation is 67 × 200 (x × y) in lattice unit.At the inlet of the liquid, we assume that the net mass flow is zero.To this end, we simply use f 2 = f 4 , f 5 = f 7 , f 6 = f 8 .At the vapor outlet, the Zou-He boundary condition is employed for the flow field.In the flow field, for the solids nodes, i.e., nodes representing the solid walls, bounce-back rule is adopted, which also applies to all the simulations presented in this paper.For the temperature field, the non-equilibrium extrapolation method is used in both inlet and outlet.Besides, the solid walls are heated at the heat flux 8 × 10 −4 , and thus the constant heat flux boundary condition is considered.Other parameters in the model are G = 0.92, β = 0.06, and ρ s = 1.9 (the resulting intrinsic contact angle is about 15 • ).The comparison of TLBM simulations and the experimental data extracted from [5] is given in Figure 8, where X denotes the normalized distance ranging from intrinsic meniscus region (X = 0) to adsorption region (X = 1).Thus the thin-film evaporating region locates in the middle.The temperature is normalized via (T − T min )/(T max ), where T min and T max are the minimum and maximum values in the temperature field, respectively.The results shows that the TLBM simulations agree well with the reported experimental data.It is noticed that there is a sudden drop of temperature in the wall temperature distribution, which implies that the intensive evaporation occurs in this region.In our simulations, it is validated that for the same computational domain, mesh resolution larger than 50 × 150 is fine enough to obtain grid-independent numerical results.

Evaporating Mass Flux
Evaporating mass flux indicates the capability of evaporation heat transfer.In this part, the mass flux in a thin-film evaporating pore is investigated using TLBM, and is then compared with the theoretical predictions.Schrage proposed that the interfacial mass flux can be estimated by the following equation [46,47], where σ is an accommodation ratio, M is the molecular weight, R is the universal gas constant, T lv is the interface temperature, P v eq is the equilibrium vapor pressure associating with the temperature T lv , P v and T v are the pressure and temperature of vapor, respectively.The above equation was derived from the kinetic theory.Denote J as the following expression Then mass flux is proportional to J [31], i.e., In the TLBM simulations investigating the evaporating mass flux, the schematic and boundary conditions introduced in the previous part are still employed.The heat flux applied to the side walls changes from 1 × 10 −4 to 8 × 10 −4 .Hence, different evaporation mass fluxes are achieved.The temperature and pressure in the centre line of the channel are evaluated in order to obtain J as well as the theoretical prediction of m .On the other hand, the mass flux of vapor exiting the pore can also be calculated via where W is the channel width, e y,i is the component of e i at the y direction, and x out is the set of nodes locating at the outlet layer.The above equation evaluates the evaporating mass flux directly from the numerical simulations.The comparisons of simulations and theoretical predictions are illustrated in Figure 9, which yields good agreement.The mass flux discussed above corresponds to an ensemble property of the evaporating pore.It is of interest to look into the details of the liquid-vapor interface, which is an advantage of TLBM.Denote the unit normal vector of the interface as ni .The local mass flux across the interface is estimated by The simulations for interfacial mass flux is shown in Figure 10, where the red solid line represents the contact line throughout adhesion, evaporating and intrinsic meniscus regions.Maximum interfacial mass flux is observed in the evaporating region compared with the adhesion region and the intrinsic meniscus region.It coincides with the classic understanding of TFE.

Simulation of Direct Vapor Generation
In this section we will focus on the characterization and optimal design of a vapor generation device that utilizes solar energy.The schematic is illustrated in Figure 11.A device-level design has been accomplished in Massachusetts Institute of Technology (MIT) [48].The device is a cylinder with a small slot cut in the middle.The absorber is on the top surface of the device, which collects the solar energy under one sun illumination.The bottom of the device is connected to the bulk liquid where water is replenished.As the pore size reduces, solar energy in the pores is more concentrated, which makes it easier for evaporation.Thermal concentration C th which is defined as the ratio of the whole surface area to the pore area can be used in the evaporating pore problem.Specifically, C th = 1/(1 − Γ), where Γ is the ratio of absorber area to the whole surface area (usually Γ ≥ 0.98).
In order to facilitate the evaporation, hydrophilic solid wall is considered for an evaporating pore, which renders TFE.The TLBM approach introduced in this paper is employed to investigate the vapor generation problem.To simplify the problem, it is assumed that the solar power collected by the absorber, denoted as Q a , is described by where q in is the incident power density, S is the total surface area (absorber plus pore).Furthermore, we assume that solar energy absorbed by the absorber is concentrated to the solid walls of pores.In the following, we will first consider the influence of pore size to a single evaporating pore.

Single Pore Scenario
The problem is investigated using the TLBM approach.In the simulations, the depth of pore is fixed as 400∆x, and the width of the pore varies according to the pore size.It is assumed that the total surface is 2 × 10 4 ∆x.Hence, Γ > 0.98 implies that the pore width is less than 400∆x.In our simulations, the heat flux density q in is set to be 1 × 10 −4 .The constant heat flux boundary condition is applied to the solid walls of the pore.By changing the pore width and adopting the corresponding heat flux at the solid walls according to Equation ( 41), scenarios with various pore width are realized.Since the vapor outlet is connected to the ambient, a constant density and temperature corresponding to the ambient are assumed for the boundary layer of the outlet.At the inlet of the pore, the replenishment boundary condition is considered for the flow field, and a constant temperature boundary condition is employed for the temperature field.We mainly focus on the overall efficiency, which is then compared with the experimental data reported in [48].
The overall efficiency η is an important performance index for vapor generation devices.It is defined as the ratio of enthalpy change in the generated vapor to the incident solar power, that is where ṁ is the mass flow rate and h f g is the specific enthalpy change of water (from liquid to vapor).When the dimensionless Van der Waals EOS is considered, the enthalpy h is expressed as [32] And h f g is obtained via its definition, that is The comparison between the experimental data reported in [48] and simulation results is shown in Figure 12.Note that in experiments, the estimation of efficiency takes into account the conduction and radiation losses of the device, which are ignored in simulations.In other words, TLBM simulation studies the ideal case.It is observed that the efficiency decreases as the thermal concentration keeps increasing.According to Equation (42), efficiency is mainly dominated by the mass flow rate ṁ, since the input power Q a and the enthalpy change h f g does not vary much.Hence, we have the following relationship where f e is the mass flux of vapor, and A e is the surface area for evaporation.In our two-dimensional study, A e is the channel width.The above relationship is also plotted in Figure 12 as a comparison with both simulations and experimental data.

Multiple Pores Scenario
In the previous part, the case of a single evaporating pore is mainly studied.It is shown that the pore size plays a crucial role in the overall efficiency.Note that the thin-film evaporating region is more efficient in vapor generation compared with the bulk region.Efficiency can thus be further improved if we increase the portion of the thin-film evaporating region.To this end, scenarios with multiple pores are considered in the following part.Specifically, the total pore size is fixed (fixed thermal concentration), while the pore number varies.In other words, more pores imply a smaller size for each one.We still focus on the two-dimensional model.Let L be the width of the absorber, and C as the thermal concentration.Then the total pore size is L/C.In the case that n identical pores are considered, the width for each pore is L/(Cn).Similar to the single pore scenario case, L is chosen as 2 × 10 4 ∆x.For example, the thermal concentration is 50, i.e., the width is 400 ∆x when n = 1.In order to investigate how the pore number affects the overall efficiency, a consistent heat flux boundary condition is considered in the simulations of different scenarios.The overall efficiency is given by where ṁn is the mass flow rate of the vapor generated in each pore.The results are shown in Figure 13.
When the number of pores increases, the overall efficiency will increase first and then decrease.In other words, maximum efficiency exists which corresponds to a moderate number of pores.For more distributed scenarios, the area of thin-film evaporating region is larger, though the total pore area is the same.As a result, the overall efficiency increases.When the number of pores keeps increasing, the pore size becomes very small.Although the thin-film evaporating region still increases, the water supplement could be slugged due to the suppressed viscosity, which reduces the evaporating mass flow rate.This phenomenon was reported in [49], where it is mentioned that, pores with diameters less that 300 nm can hinder water evaporation.Conditions with various incident power are also shown in Figure 13, specifically, q in = 1 × 10 −4 and q in = 1.2 × 10 −4 are demonstrated.It is revealed that when the incident power is larger, the overall efficiency is also higher.The difference is more obvious for scenarios with more distributed pores.However, it is also noticed that the peak position keeps the same, specifically, n = 11 renders the maximum efficiency in our simulations.Another factor that may affect the overall efficiency of vapor generation is the wettability of solid walls.The pores discussed above have purely hydrophilic solid walls.The most efficient region of TFE is the thin-film evaporating region, the formation of which relies on the hydrophilicity of the wall.Beyond this region, specifically, near the vapor outlet, the liquid film becomes even thinner, the liquid there can hardly transfer to vapor through evaporation.A section of hydrophobic solid wall can be employed at the vapor outlet, the length of which is 50∆x in our simulations.Hence, pore walls have heterogeneous hydrophobicities.The performances of pores with and without hydrophobic sections are shown in Figure 14.The hydrophobic sections contribute to an efficiency increment for distributed scenarios.The advantage of hydrophobic sections is more obvious as the number of pores increases.To look into the reason, the flow filed inside a pore is investigated.Figure 15 illustrates the flow fields of both cases.In case (a), the solid walls are purely hydrophilic.It is observed that some vapor flows back to the pore near the corners of the vapor outlet.However, this phenomenon almost disappears when hydrophobic sections are included, which implies that it is easier for vapor generated in the pore to escape.

Conclusions
A thermal multiphase lattice Boltzmann model with normalized Van der Waals EOS was developed in this paper to simulate heat and mass transfer in TFE.Validations of the approach were carried out by simulating some benchmark problems and comparing with published results.One contribution of this paper is that it shows that TLBM can provide details of the liquid-vapor interface in TFE, for example, the maximal mass flux is achieved in the thin liquid film region.
The detailed information may help people in the phenomenon-level understanding and device-level design.The TLBM approach was further used to study a vapor generation problem where TFE is involved.To this end, a replenishment boundary condition was proposed, which ensured the balanced mass in the evaporating pore.Two main factors that influence the overall vapor generation efficiency were mainly studied.One is the number of pores that are used for vapor generation given a fixed thermal concentration, and the other is the hydrophilicity/hydrophobicity of the solid walls.Based on the TLBM simulations, scenarios leading to higher overall efficiency were suggested, which is another contribution of our work.Specifically, distributed evaporating pores with hydrophobic pore ends are shown to have the highest efficiency for vapor generation.However, when the pores are too dense, the overall efficiency will decrease due to the suppressed viscosity.
Our future work will extend the approach to three-dimensional cases.For the vapor generation problem, the three-dimensional TLBM approach will further enable us to consider the factor of pore shape, and thus will be of more practical value.

Figure 2 .
Figure 2. Illustration of south wall and discretized velocity directions.

Figure 3 .
Figure 3. Evolution of the evaporating mass flux at the vapor outlet: mini plots demonstrate the density distribution at different time instants; the red part is liquid and the blue part is vapor; the computational domain is 101 × 200, and the heat flux imposed to the walls is 8 × 10 −4 .

Figure 4 .
Figure 4. Evolution of the total mass in the pore and the evaporating mass flux at the outlet: the blue solid line is mass flux; the green dotted dash line is the total mass; the implementation of ghost layer boundary condition starts from time 4 × 10 4 .

Figure 5 .
Figure 5.The equilibrium contact angle (CA) changes with ρ w : The red region is liquid; the blue region is vapor; the bottom is solid wall; the whole field is isothermal with T = 0.7.

Figure 6 .
Figure 6.The influence of parameters (G and α) on liquid-vapor density ratio: subfigure (a) gives the relationship between density ratio and the parameter G; subfigure (b) gives the relationship between density ratio and the parameter α; the whole field is isothermal with T = 0.7.

Figure 7 .
Figure 7. Illustration of side walls and wetting liquid: black lines are side walls, blue region represents liquid, above white region is vapor; the two side walls are heated for evaporation.

Figure 8 .
Figure 8.The wall temperature distribution under constant flux: black crosses are experimental data from [5]; red line is obtained by the thermal lattice Boltzmann method (TLBM) simulations; X and T * are normalized location and temperature, respectively; density profile of the half channel is presented in the mini plot.

Figure 9 .
Figure 9. Mass flux under different heating conditions: the black cross markers are the simulations, and the red line is the theoretical prediction by Equation (36), both J and m are normalized.

Figure 10 .
Figure 10.Mass flux around the interface: the red line is the interface of liquid and vapor, the vapor is at the left side and the liquid is at the right side; the black line is the solid wall; the vectors with arrows indicate the local mass fluxes that are normal to the interface, and the vector norms are proportional to the mass fluxes.

Figure 11 .
Figure 11.Microscale pores for vapor generation: the yellow layer is the absorber that absorbs solar energy; under the absorber is the insulator; liquid is shown as the blue region of the figure; vapor is generated from evaporation.

Figure 12 .
Figure12.The normalized overall efficiency vs. thermal concentration: red squares are experimental data reported from[48]; black solid line represents simulation results; efficiency is normalized for comparison, specifically, η * = 1 corresponds to C th = 53, and η * = 0 corresponds to C th = 486.

Figure 13 .
Figure13.Overall efficiency of multiple pores obtained from TLBM simulation: blue squares represent the case of q in = 1.2 × 10 −4 ; red circles represent the case of q in = 1 × 10 −4 ; efficiency is normalized with respect to the case that T wall = 1 × 10 −4 and n = 1.

Figure 14 .
Figure 14.Overall efficiency of multiple pores with hydrophobic sections: blue squares represent the case without any hydrophobic sections; red circles represent the case with hydrophobic sections; the wall temperature is set as 0.72 for both cases; efficiency is normalized with respect to the case that there is only one purely hydrophilic pore.

Figure 15 .
Figure 15.The steady-state density profile (red for liquid and blue for vapor) and streamlines in two cases: (a) the solid wall is purely hydrophilic (contact angle is about 15 • ); (b) the top section of the solid wall is hydrophobic (contact angle is about 145 • , section length is 50∆x) while the remaining part of the wall is still hydrophilic.