# Volume of Fluid (VOF) Modeling of Liquid Film Evaporation in Mixed Convection Flow through a Vertical Channel

Authors: Hayat El Baamrani; Lahcen Bammou; Ahmed Aharoune; Abdallah Boukhris
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, the volume of fluid (VOF) method in the OpenFOAM open-source computational fluid dynamics (CFD) package is used to investigate the coupled heat and mass transfer by mixed convection during the evaporation of water-thin film. The liquid film is falling down on the left wall of a vertical channel and is subjected to a uniform heat flux density, whereas the right wall is assumed to be insulated and dry. The gas mixture consists of air and water vapor. The governing equations in the liquid and in the gas areas with the boundary conditions are solved by using the finite volume method. The results which include temperature, velocity, and vapor mass fraction are presented. The effect of heat flux density, liquid inlet temperature, and mass flow rate on the heat and mass transfer are also analyzed. Better liquid film evaporation is noted for the system with a higher heat flux density and inlet liquid temperature or a lower mass flow rate. Therefore, the VOF method describes well the thermal and dynamic behavior during the evaporation of the liquid film.


---

## Body

## 1. Introduction

Falling liquid film evaporation has been widely studied in recent decades due to its interest in improving coupled heat and mass transfer for various engineering applications, such as seawater desalination, distillation, drying of agricultural products, cooling of electronic components, and air conditioning. Many studies linked to the heat and mass transfer characteristics are available in the literature. Classical works including evaporation of pure liquid film have been largely studied [1–7].A mixed convection heat and mass transfer during the liquid film evaporation has been the subject of a large number of numerical or experimental studies. For numerical study, Van and Lin [8] studied the effects of natural convection all along a finite liquid film evaporation in a vertical channel. They showed that the assumption of an extremely thin film is only useful for a system with a small liquid mass flow rate. Yan et al. [9, 10] presented a numerical study of the effect of the wall heat flux and inlet liquid mass flow rate during finite liquid film evaporation of laminar mixed convection in a vertical channel. The results indicate that the heat and mass transfer is dominated by latent mode. Later, Feddaoui et al. [11] presented a numerical study for liquid film evaporation inside an insulated vertical channel. They reported that the evaporation is higher for low liquid flow rate or high gas flow Reynolds number. They also treated the case of liquid film evaporation with turbulent flows [12]. Mezaache and Daguenet [13] presented a numerical study of the effect of inlet conditions on the evaporation process along the gas-liquid interface over an inclined plate. Their results showed that the effect of inclination is highly dependent on the liquid mass flow rate and gas velocity and an increase in the inclination improves the evaporation by increasing the vapor mass flow rate. Laaroussi et al. [14] investigated the mixed convection in a vertical channel of thin film evaporation on wetted walls. They compared the Boussinesq approximation with the density change for air-water vapor and air-hexane vapor mixtures. Mass diffusion induces an important flow acceleration in the boundary layers and the recirculation of flow at the center of the channel for large mass evaporation rates.Cherif et al. [15] presented an experimental study of the evaporation by natural and forced convection in a vertical channel, and a thin liquid film falls on two parallel plates. They showed that the thermal yield of evaporation increases with the down air velocity and the heat flux density in order to improve heat and mass transfer during liquid film evaporation.Cherif et al. [16] investigated experimentally the influence of film evaporation on mixed convective heat and mass transfer in vertical parallel channel wetted by a water film. The results showed that the evaporative cooling increases with small heat flux and large air velocities. Alla et al. [17] studied the evaporation of glycol liquid film in mixed convection. Their results showed that the propylene glycol evaporates more than glycol, and the heat transfer by the latent mode is more important than the heat transfer by sensible mode. Nasr and Alghamdi [18] investigated numerically the evaporation and condensation under mixed convection of binary liquid film. They showed that the water condensation increases by reducing inlet temperature and increasing the amount of liquid film, and an increase of inlet liquid concentration of ethylene glycol enhances the ethylene glycol evaporation and water vapor condensation. Zhang et al. [19] studied the evaporation by laminar mixed convection in a vertical three-dimensional rectangular channel with nitrogen and water vapor. The results indicate that the wall temperature and the inlet temperature improve the heat and mass transfer. Alla et al. [20] have investigated numerically two configurations. In the first case, a uniform heat flow is applied to the wet wall while in the second case, and the same heat flux was applied to preheat the liquid at the inlet. The authors indicate that the evaporation rate increases with the second case. Fang et al. [21] studied the heat transfer coefficient correlation of the multiphase transformation inside and outside the tube. They showed that the heat transfer coefficient increases with the feed flow rate and ethanol content in the evaporation tube. The heat transfer coefficient decreases when the heat transfer temperature difference between the inside and outside of the tube increases.We also find in the literature some other studies that have focus on heat and mass transfer with the evaporation of a liquid film flowing in a porous layer. Leu et al. [22] numerically studied the heat and mass improvement of liquid film evaporation with a porous layer. The results showed that the heat transfer performance of the liquid film was improved by using a porous layer, and the average Nusselt and Sherwood numbers are increased inversely with porosity, porous layer thickness, and relative humidity. Nasr and Al-Ghamdi [23] analyzed the influence of porous layer thickness and porosity on the coupled heat and mass transfer performance during the evaporation of flowing liquid film. The results indicate that the presence of the porous layer improves the heat and mass transfer performance at the liquid-gas interface. Terzi et al. [24] numerically studied the performance of liquid film evaporation by introducing a liquid saturated porous layer along an inclined channel. The results indicated that the use of the porous layer improves the heat and mass transfer.The VOF method has been used in the literature in many problems involving two-phase flows. Piscaglia et al. [25] presented a development and validation of a two-phase dynamic VOF solver for the simulation of internal nozzle flows. The phase change at the interface was modeled using the Schnerr and Saver model. Large Eddy simulations (LES) are used to model turbulence effect. The approach can be considered to be reliable for investigations of the internal nozzle flows during transient operations. Pham and Choi [26] performed a numerical simulation for the water hammering using interFoam solver. The pressure-implicit method is used for continuity, momentum, and energy equations of a volume of fluid model discretization. Implicit Euler and central difference schemes were used for temporal and spatial discretization. The solver is validated with a one-dimensional Stefan problem’s benchmark case. Amidu et al. [27] presented a hybrid VOF- Eulerian multifluid model for solving a hot low and high void fraction subcooled flow boiling. They showed that the best estimation of the radial velocity profile is achieved by using the two-phase rough wall function instead of the classical single-phase wall function.In our knowledge, all the previous studies deal with the heat and mass transfer during the liquid film evaporation with the finite volume method by considering boundary conditions at gas-liquid interface. The main objective of the present work is the use of volume of fluid (VOF) method to investigate numerically the mixed convection heat and mass transfer with liquid film evaporation. To this end, an OpenFoam solver was developed and successfully validated with previous numerical results taken from the literature. The effect of parameters such as the inlet liquid temperature mass flow rate and wall heat flux density on the heat and mass transfer during water film evaporation in a vertical channel was investigated.


## 2. Problem Formulation

### 2.1. Definition of Geometry

The present work deals with a numerical study of liquid film evaporation by mixed convection in a vertical channel, as illustrated in Figure1. The geometry under consideration is two vertical and parallel plates with height H and width b. The first plate is maintained at the uniform heat flux qe and wetted by a thin liquid film, with inlet mass flow B0 and temperature TL0, while the second one is dry and insulated. The gas enters in the channel with inlet temperature T0, a water-vapor concentration W0, and velocity V0.Figure 1
Schematic diagram of geometry.


### 2.2. Assumptions

The following simplifying assumptions are taken into account:Flows and transfers in the two phases are laminar and incompressibleThe effect of the superficial tension is negligibleThe liquid and gas flows are two dimensional and transientThe channel walls are supposed very thin


### 2.3. Governing Equations

In the present work, the volume of the fluid VOF method [28] is used to track the liquid-gas interface. In this method, a single momentum equation is shared by the fluids, and the volume fraction of each of the fluids in each computational cell is tracked throughout the domain [29, 30]. The volume fraction αL indicating the ratio of the liquid volume in a cell to the volume of the cell is equal to one if the cell is completely filled with liquid and is equal to zero if the cell is totally filled with gas:(1)αLx,t=VLiquidV=0,inside gas phase,0<αL<1,inside transitional region,1,inside liquid phase.The fluid density and viscosity of a single mixture phase in a computational cell are defined as follows:(2)ρ=αLρL+1−αLρv,μ=αLμL+1−αLμv.The continuity equation is(3)∇⋅ρU⟶=0.Taking into account the effect of gravity, the momentum equation is defined as(4)∂∂tρU⟶+∇⋅ρU⟶U⟶=−∇p+∇⋅μ∇U⟶+∇U⟶T+ρg⟶.The energy equation is(5)∂∂tρE+∇⋅U⟶ρE+p=∇⋅k∇T+Sh.Lee’s model [31] is adopted in the VOF method, and it is widely used to simulate mass transfer during evaporation and condensation processes. In this model, the mass transfer depends mainly on the saturate and shared temperatures. The mass transfer can be described as follows:(6)IfTL>Tsatm˙LV=coeff∗αLρLTL−TsatTsat,IfTV<Tsatm˙VL=coeff∗αVρVTsat−TVTsat.


### 2.4. Boundary Conditions

The inlet conditions (x = 0) are as follows.In the liquid phase,(7)TL=TL0,B=B0.In gas phase,(8)TG=T0,VG=V0,WG=W0.At the wet wall (y = 0),(9)VL=0,−λ∂TL∂y=qe.At the dry wall (y = b),(10)VG=0,∂TG∂y=0,∂WG∂y=0.


## 2.1. Definition of Geometry

The present work deals with a numerical study of liquid film evaporation by mixed convection in a vertical channel, as illustrated in Figure1. The geometry under consideration is two vertical and parallel plates with height H and width b. The first plate is maintained at the uniform heat flux qe and wetted by a thin liquid film, with inlet mass flow B0 and temperature TL0, while the second one is dry and insulated. The gas enters in the channel with inlet temperature T0, a water-vapor concentration W0, and velocity V0.Figure 1
Schematic diagram of geometry.


## 2.2. Assumptions

The following simplifying assumptions are taken into account:Flows and transfers in the two phases are laminar and incompressibleThe effect of the superficial tension is negligibleThe liquid and gas flows are two dimensional and transientThe channel walls are supposed very thin


## 2.3. Governing Equations

In the present work, the volume of the fluid VOF method [28] is used to track the liquid-gas interface. In this method, a single momentum equation is shared by the fluids, and the volume fraction of each of the fluids in each computational cell is tracked throughout the domain [29, 30]. The volume fraction αL indicating the ratio of the liquid volume in a cell to the volume of the cell is equal to one if the cell is completely filled with liquid and is equal to zero if the cell is totally filled with gas:(1)αLx,t=VLiquidV=0,inside gas phase,0<αL<1,inside transitional region,1,inside liquid phase.The fluid density and viscosity of a single mixture phase in a computational cell are defined as follows:(2)ρ=αLρL+1−αLρv,μ=αLμL+1−αLμv.The continuity equation is(3)∇⋅ρU⟶=0.Taking into account the effect of gravity, the momentum equation is defined as(4)∂∂tρU⟶+∇⋅ρU⟶U⟶=−∇p+∇⋅μ∇U⟶+∇U⟶T+ρg⟶.The energy equation is(5)∂∂tρE+∇⋅U⟶ρE+p=∇⋅k∇T+Sh.Lee’s model [31] is adopted in the VOF method, and it is widely used to simulate mass transfer during evaporation and condensation processes. In this model, the mass transfer depends mainly on the saturate and shared temperatures. The mass transfer can be described as follows:(6)IfTL>Tsatm˙LV=coeff∗αLρLTL−TsatTsat,IfTV<Tsatm˙VL=coeff∗αVρVTsat−TVTsat.


## 2.4. Boundary Conditions

The inlet conditions (x = 0) are as follows.In the liquid phase,(7)TL=TL0,B=B0.In gas phase,(8)TG=T0,VG=V0,WG=W0.At the wet wall (y = 0),(9)VL=0,−λ∂TL∂y=qe.At the dry wall (y = b),(10)VG=0,∂TG∂y=0,∂WG∂y=0.


## 3. Solution Method

The system of equations with boundary conditions are solved numerically using the open-source library OpenFOAM [32], which is based on the finite volume method. The icoReactingMultiphaseInterFoam solver is used in this study based on the VOF method. The pressure-velocity coupling is handled by solving the momentum equations and a pressure equation derived from a combination of the continuity and the momentum equations. This coupling has been ensured by the PIMPLE algorithm (pressure-implicit method for pressure-linked equations). Figure 2 shows the PIMPLE flowchart, which is a combination of the pressure implicit with splitting of operators’ (PISO) algorithm and the Semi-implicit method for pressure-linked equations’ (SIMPLE) algorithm. Although this method is formulated for the pseudo-transient simulations with a large courant number, it permits an under-relaxation in the equation in order to assist in stabilizing the numerical convergence whilst maintaining numerical accuracy. The spatial schemes used to approximate the convective and diffusive fluxes are both of the second-order accuracy. The gradient, the divergence, and the Laplacian terms are approximated by the Gauss Linear, the Gauss linearUpwind, and the Gauss linear orthogonal schemes, respectively. These second-order schemes give very good precision in the results despite their difficulty in achieving convergence. The Euler scheme is employed for the time derivative discretization. Adjustable time step is used to ensure the courant number limit to 1. The convergence criteria are satisfied for all calculated variables when (Err≤10−7).Figure 2
PIMPLE flowchart.The courant number is given by(11)Co=V⋅ΔtΔx.


### 3.1. Mesh Independence Validation

The mesh quality affects the stability and the calculation accuracy. In order to improve the mesh quality, the verification of the grid independence is necessary. All domains are meshed using the structured mesh, and finer grids are used in the regions near the inlet and the wall, as shown in Figure3. The average cell length in the region of the liquid film is 0.1 mm.Figure 3
Mesh structure of the physical model.To examine the independence of mesh size, three different grid sizes are used which are 118959, 263940, and 598410, respectively. Figure4 shows the wall temperature variation for three different grids. The result shows that there is a little difference between the simulation with 263940 and 598410 grids. Hence, the size of grid is 263940 in this study.Figure 4
Variation of wall temperature in three different grids.


### 3.2. Numerical Validation

In the present study, to check the validation and the adequacy of the used numerical schemes, we compared the results obtained by the code with results available in the literature. The numerical solver we developed was validated with the results of Yan [9] study, concerning the evaporation of the pure water film in a vertical channel which is uniformly heated on the left wall. The comparison has been done for T0 = 20°C, W0 = 50%, TL0 = 20°C, Re = 2000, and qe = 1000 W/m2. The comparison illustrates the evolution of the interfacial temperature (Figure 5(a)), the axial velocity profile (Figure 5(b)), and the interfacial vapor mass fraction (Figure 5(c)). A good agreement between our results and those obtained by Yan [9] is observed. Therefore, the present numerical method and the used grid layout are considered to obtain accurate results for practical purpose.Figure 5
Comparison of the present study with that of Yan [9]. Validation of calculated (a) interfacial temperature, (b) axial velocity profile, and (c) interfacial vapor mass fraction.
(a)(b)(c)


## 3.1. Mesh Independence Validation

The mesh quality affects the stability and the calculation accuracy. In order to improve the mesh quality, the verification of the grid independence is necessary. All domains are meshed using the structured mesh, and finer grids are used in the regions near the inlet and the wall, as shown in Figure3. The average cell length in the region of the liquid film is 0.1 mm.Figure 3
Mesh structure of the physical model.To examine the independence of mesh size, three different grid sizes are used which are 118959, 263940, and 598410, respectively. Figure4 shows the wall temperature variation for three different grids. The result shows that there is a little difference between the simulation with 263940 and 598410 grids. Hence, the size of grid is 263940 in this study.Figure 4
Variation of wall temperature in three different grids.


## 3.2. Numerical Validation

In the present study, to check the validation and the adequacy of the used numerical schemes, we compared the results obtained by the code with results available in the literature. The numerical solver we developed was validated with the results of Yan [9] study, concerning the evaporation of the pure water film in a vertical channel which is uniformly heated on the left wall. The comparison has been done for T0 = 20°C, W0 = 50%, TL0 = 20°C, Re = 2000, and qe = 1000 W/m2. The comparison illustrates the evolution of the interfacial temperature (Figure 5(a)), the axial velocity profile (Figure 5(b)), and the interfacial vapor mass fraction (Figure 5(c)). A good agreement between our results and those obtained by Yan [9] is observed. Therefore, the present numerical method and the used grid layout are considered to obtain accurate results for practical purpose.Figure 5
Comparison of the present study with that of Yan [9]. Validation of calculated (a) interfacial temperature, (b) axial velocity profile, and (c) interfacial vapor mass fraction.
(a)(b)(c)


## 4. Results and Discussion

In this study, we present the effects of flow conditions on heat and mass transfer during the evaporation of the pure water film in a vertical channel whose left wall is uniformly heated. The ranges of physical parameters are listed in Table1. The relative humidity of the air is assigned to 50% for the air/water vapor mixture along a vertical channel.Table 1
The ranges of physical parameters.
TL0 (°C)20T0 (°C)20B0 (kg/m.s)0.01, 0.02, 0.04V0 (m/s)2.1qe (W/m2)1000, 1500, 2000H (m)0.75b (m)0.015Re2000Figure6 illustrates the temperature distributions along the y-axis for various x positions. Figure 6(a) reveals that the temperature decreases monotonically with Y at specific axial locations, which means that the direction of heat transfer is from the liquid film through to the gas stream. We also observe that, for Y between 0 and 0.5, the temperature rises as the flow progresses, while for Y > 0.5, the temperature profile is not affected throughout the flow. Figure 6(b) shows that an increase of the heat flux density applied to the humid wall qe induces an increase of the temperature in the liquid film and the region in the gas flow, where the temperature profile is unaffected and increases by starting from Y = 0.4. Figure 6(c) shows the impact of the rise in mass flow rate B0 on the temperature profile. It is clearly observed that an increase of the mass flow rate B0 decreases the temperature at a given X position. In this case, we observe also that the region in the gas flow, where the temperature profile is not affected, increases and starts at Y = 0.3.Figure 6
Evolution of axial temperature profiles for (a)B0 = 0.02 kg/m.s, qe=1000W/m2, TL0 = 20°C, and Re = 2000, (b) qe=1500W/m2, and (c) B0 = 0.04 kg/m.
(a)(b)(c)Figure7 shows the effects of the applied heat flux on the temperature profile development along the gas-liquid interface. It is observed that the temperature increases from the inlet to the outlet which is due to the heat transferred from the wall to the liquid film. As expected, the interfacial temperature is raised for higher applied heat flux qe.Figure 7
Effects of wall heat flux on the interfacial temperature along the channel forB0 = 0.02 kg/m.s,TL0 = 20°C, and Re = 2000.In order to investigate the importance of the inlet liquid temperature on the heat and mass transfer, we present in Figure8 the evolution of the temperature along the gas-liquid interface for different inlet liquid temperature. The results disclose that the temperature decreases for a higher inlet liquid temperature. This is due to the evaporation process because the evaporation of the liquid film is greater as the liquid temperature rises. Consequently, a large amount of energy is taken from the liquid film to ensure latent heat of evaporation at the interface.Figure 8
Effects of inlet liquid temperature on the interfacial temperature along the channel forB0 = 0.02 kg/m.s,qe=1000W/m2, and Re = 2000.Figure9 shows the distribution of the liquid volume fraction within the channel. From this distribution, we can see that the liquid-gas interface is almost smooth and has no waves all along the channel. It is also observed from this distribution that the liquid film thickness in the lower part of the channel is thinner than that in the upper part, which is simply due to the evaporation of the liquid film.Figure 9
Liquid volume fraction distribution along the channel forB0 = 0.02 kg/m.s,qe=1000W/m2, and Re = 2000.Figure10 presents the evolution of the axial velocity profiles under various conditions. The velocity profile is almost uniform near the inlet. As the downstream flow, the velocity profile becomes parabolic. The maximum velocity value takes place in the outlet and displaces to the humid wall. The change of the velocity peak can be explained by the film evaporation which is in accord with the results provided by Yan [9]. Comparing cases in Figures 10(a)–10(d), a very small variation of velocity is observed by the change in the inlet liquid temperature, mass flow rate, and heat flux density. We also observe from these figures that the interfacial velocity always remains constant along the channel.Figure 10
Evolution of axial velocity profiles for (a)B0 = 0.02 kg/m.s qe=1000W/m2, TL0 = 20°C, and Re = 2000, (b) TL0 = 40°C, (c) B0 = 0.01 kg/m.s, and (d) qe=1500W/m2.
(a)(b)(c)(d)Figure11 presents the interfacial vapor mass fraction profiles along the gas-liquid interface. It is easy to observe that the interfacial vapor mass fraction follows the same trend as the interfacial temperature due to the thermodynamic equilibrium at the liquid-gas interface. A careful inspection in Figure 11 discloses that, as flow moves downstream, the vapor mass fraction increases from the inlet to the outlet of the channel due to the evaporation of the liquid film as a result of the increase in liquid film temperature. It is also observed that the increase in the mass flow rate of the liquid film causes a decrease in the interfacial vapor mass fraction.Figure 11
Development the interfacial vapor mass fraction along the channel forqe=1000W/m2, TL0 = 20°C, and Re = 2000.Figures12(a) and 12(b) show, respectively, the influence of the heat flux density and the inlet liquid film temperature on the interfacial vapor mass fraction along the channel. It can be noted that a higher heat flux leads to a higher latent heat transfer. This results in more water vaporization. On the contrary, the interfacial vapor mass fraction increases as the inlet temperature of the liquid film increases.Figure 12
Effects of heat flux density and inlet liquid temperature on the interfacial vapor mass fraction profiles.
(a)(b)


## 5. Conclusion

In this work, the heat and mass transfer during the evaporation of thin water film in a mixed convection flow along a vertical channel are investigated using the volume of fluid (VOF) method within OpenFOAM CFD package. The effects of the inlet liquid mass flow rate, the heat flux density applied to the wet wall, and the inlet liquid temperature on the heat and mass transfers are analyzed and discussed in detail. The obtained results for water film evaporation show good agreement compared to the results previously published in the literature which indicates that the interface tracking by the VOF method respects well the thermodynamic equilibrium at the liquid-gas interface.The summary of the major results are as follows:(1)
The interfacial temperature is higher for systems with greater heat flux densityqe and also with a lower mass flow rate B0(2)
The velocity profile through the channel is slightly affected by the change in the inlet liquid temperature, mass flow rate, and heat flux density(3)
The increase of imposed heat flux densityqe and inlet liquid temperature TL0 leads to an increase of the evaporation phenomenon(4)
A decrease of mass flow rateB0 improves the water evaporationAll these results show that the VOF method describes well the thermal and dynamic behavior during the evaporation of the liquid film without imposing equilibrium conditions at the liquid-gas interface. This insists us to develop this method in this configuration to benefit from its interface tracking advantages.



---

*Source: 9934593-2021-05-24.xml*