# Computational Study of Liquid Film Evaporation along a Wavy Wall of a Vertical Channel

Authors: Monssif Najim; M’barek Feddaoui; Abderrahman Nait Alla; Adil Charef
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A numerical study of mixed convection heat and mass transfer along a vertical channel with a wavy wall is performed. The wavy wall is heated by a constant flux, while the other is adiabatic. The discretisation of equations in both liquid and gas phases is realised using an implicit finite difference scheme. Results of simulation compare the effect of multiple parameters, especially amplitude and characteristic length of the curve, on the liquid film evaporation process. The results indicate that heat and mass transfer is enhanced by increasing the amplitude and number of wall waves. Moreover, a very small value of waves amplitude of the wall may reduce the sensible heat and mass transfer.


---

## Body

## 1. Introduction

Corrugated channels are one of several devices proposed by manufacturers in the HVAC industry aiming at the enhancement of heat transfer in many heat exchanger applications. A part of those applications include flowing liquid films. Thus, simultaneous heat and mass transfer between the thin liquid film and gas/vapour stream is widely encountered. At this stage, evaporation process has a great impact on performance and energy efficiency. Those processes involve conjugate heat transfer between liquid and gas with the change of phase due to evaporation or condensation, what makes their physics quite complex.So far, many studies have been realised. Early, Parken et al. [1] have conducted an experimental study of heat transfer in a falling liquid film outside horizontal tube. The authors investigated the influence of multiple inlet parameters like feedwater flow rate and temperature in addition to some configuration parameters like feedwater distribution system and tube diameter on the average heat transfer coefficient. Yan et al. [2] have experimentally studied the role of interfacial heat and mass transfer in the evaporative cooling process of a liquid film in a vertical channel. The results showed the independency of the evaporative cooling on the inlet liquid film temperature and flow rate. The authors identified that effective cooling can be reached by smaller liquid film flow rates and high temperatures. A second part of the same study is carried out numerically by Yan and Lin [3]. The numerical results were very similar to their experiment and showed the dominance of latent heat transfer in the process. Afterwards, many numerical investigations have been done by Yan [4–6]. He et al. [7] realised a numerical study of counter-courant liquid film cooling of a uniformly heated wall. The study revealed two modes of cooling: direct film cooling and evaporative film cooling. The authors defined the convenient conditions to establish each cooling mode. To understand the combined buoyancy effects of thermal and mass diffusion on the turbulent mixed convection flows inside a tube, Feddaoui et al. [8] have conducted a detailed numerical study reporting the best conditions to have a better heat transfer, precisely, a higher gas flow Reynolds number, a higher heat flux, or a lower inlet water flow. Feddaoui et al. [8, 9] extended many of the previous studies of Yan and his coworkers. They reported numerical results for tube and channel’s geometries. The authors reported the convenient conditions for the effective cooling of the liquid film and for heat removal from the wall. They also mentioned the existence of two cooling mechanisms: convective cooling in the case of low liquid film temperature and phase change cooling. Airstream humidification by an evaporating liquid film was discussed by Orfi et al. [10]. The study reported the effect of liquid film characteristics on the heat and mass transfer and observed that the interfacial heat flux may take higher values than the imposed heat flux. A description of convective heat and mass transfer of falling liquid film evaporation in a cavity was reported by Ben Jabrallah et al. [11]. It is shown that intensification of heat and mass transfer occurs with lower feed mass flow rate. Experimental investigation of a wavy liquid film falling down a flat plate was conducted by Zhou et al. [12] in order to measure the thickness of the liquid film. The presented results revealed the influence of Reynolds number and liquid feed on the thickness of the liquid film. Cherif et al. [13] examined experimentally the mixed convective heat and mass transfer along a vertical rectangular channel during falling liquid film evaporation. In addition, they performed a numerical study and simulation results showed a good agreement with their experiment. Aziz et al. [14] analysed the impact of thermocapillarity and radiative heat transfer on flow and heat transfer in a thin liquid film on an unsteady stretching sheet. They found that the thermocapillarity reduces the temperature distribution and the surface shear stress, while it improves the velocity and dimensionless heat flux. In contrast, it was found that radiation parameter helps to enhance the temperature, while reducing the dimensionless heat flux. The authors did not consider the mass transfer process. Recently, Nait Alla et al. [15, 16] have deeply analysed the evaporation of multiple liquid films under different conditions and geometries. The results were presented for glycols and ethanol liquid films evaporation for insulated, heated, and partially heated walls. Heat and mass transfer in thin liquid film over stretching surface is deeply examined by Shkarah et al. [17]. The analytical solutions permitted us to calculate the distribution of the heat flux and the evaporating liquid film thickness in a microchannel. The authors claimed the possibility of neglecting the capillary pressure. The obtained results were compared to numerical solutions. Najim et al. [18] were interested in the study of salinity impact on the evaporation of saline water film inside a vertical tube. The computational study identified the convenient conditions to enhance the heat and mass transfer and showed that increasing the concentration of salinity enhances the evaporation.The problem of flows in corrugated channels was early investigated using analytic methods. Goldstein Jr. and Sparrow [19] were the first to investigate experimentally the effect of triangular waves of a corrugated channel wall. They observed an enhancement of the average heat transfer coefficient comparing to flat channel walls for turbulent regime, but with an increasing demand on the pumping power. Many studies have been done afterward, for various waves shapes, lengths, amplitudes, and flow regimes [20–22]. A very limited number of works investigated mixed convection over wavy walls with conjugated heat and mass transfer. For example, Jang et al. [23, 24] analysed numerically the mixed and natural convection along a vertical wavy wall and showed that higher amplitude-wavelength ratio increases the fluctuation of velocity, temperature, and concentration fields.Through this literature review, a large number of studies examined the liquid film evaporation on smooth surfaces, most of them on a flat plate. However, many solutions on the industry consist of the use of irregular surfaces: textured, wired, or wavy surfaces crossed by a liquid film. To the best knowledge of this work authors, there is no numerical study dedicated to analysing the heat and mass transfer of a flowing liquid film over a wavy wall. The present paper presents a numerical investigation of the heat and mass transfer during liquid film evaporation over a wavy wall of a vertical channel. The aim of this study is to reveal the effects of the amplitude and characteristic length of the wavy surface on the evaporation.


## 2. Analysis

Let us consider a water liquid film flowing downward a wavy wall of a vertical channel as illustrated in Figure1. The liquid film enters the channel with TL0 at feed rate Γ0, while the humid gas flow enters the channel at Re0, T0, and P0 in contact with the liquid-gas interface and the adiabatic wall. The channel length and initial width are, respectively, L=1.5m and d=0.02m. As the liquid flows down the wavy wall, subjected to uniform heating flux, a quantity of the flow evaporates and the liquid film thickness decreases.Figure 1
Schematic diagram of the physical system.


### 2.1. Governing Equations

To formulate the governing equations, the following assumptions are done:(i)
The flow in the two phases is considered laminar, incompressible, and of Newtonian fluids.(ii)
The problem is two-dimensional and steady.(iii)
The thermodynamic and thermophysical properties of air and air-water vapour mixture are a function of both local temperature and concentration.(iv)
The film surface is in thermodynamic equilibrium and smooth.(v)
Heat transfer by radiation, the effect of the surface tension, Dufour, Soret, and other secondary effects are negligible.(vi)
The axial diffusion of heat, momentum, and mass were assumed to be negligible.(vii)
The vapour-gas mixture is an ideal gas mixture and the saturation conditions are evaluated at the liquid-mixture interface.The pure component data are approximated by polynomials in function of temperature and mass fraction. For further details, the thermophysical properties are available in [25, 26].Under such assumptions, the governing equations in the liquid phase could be written as(1)∂∂xρLuL+∂∂yρLvL=0,(2)uL∂∂xρLuL+vL∂∂yρLuL=-∂P∂x+∂∂yμL∂uL∂y+ρLg,(3)uL∂ρLCPLTL∂x+vL∂ρLCPLTL∂y=∂∂yλL∂TL∂y.Similarly, the conservation equations governing the gaseous phase are(4)∂∂xρGuG+∂∂yρGvG=0,(5)uG∂ρGuG∂x+vG∂ρGuG∂y=-∂P∂x+∂∂yμG∂uG∂y+ρG-ρ0g,(6)uG∂ρGCPGTG∂x+vG∂ρGCPGTG∂y=∂∂yλG∂TG∂y+∂∂yDCPV-CPa∂w∂yTG,(7)ρGuG∂w∂x+vG∂w∂y=∂∂yD∂w∂y.The curvature of the wavy wall is described by the following sinusoidal function:(8)ϖx=asin2πxl,where a is the amplitude of the wave and l its characteristic length. Consequently, the number of waves is defined as n=L/l.Accordingly, the effective width of the channel at eachx-location is(9)b=d-ϖx.


### 2.2. Boundary Conditions

The governing equations are subject to the following boundary conditions:(i)
At the inlet (x=0):(10)TG=T0;uG=u0;w=w0;Hr=Hr0,P=P0;TL=TL0;Γ=Γ0.(ii)
At the heated wall (y=b):(11)uL=0;λL∂TL∂y=QW.(iii)
At the insulated wall (y=0):(12)uG=0;∂TG∂y=0;∂w∂y=0.The matching conditions at the interfacey=b-δ are as follows:(13)uIx=uG,I=uL,I,(14)TIx=TG=TL,I,(15)τI=μ∂u∂yL,I=μ∂u∂yG,I.Energy balance at the gas-liquid interface is(16)λ∂T∂yL,I=λ∂T∂yG,I+m˙hfg,where m˙=-ρGvI=ρGD/1-wI∂w/∂y and hfg is the latent heat of vapourization.In addition, we should add mass conservation equations to compute the unknown pressure gradientdP/dx:(17)Γ0=∫b-δbρudyL-∫0xρGvIdx,(18)b-δ0ρ0u0=∫0b-δρGuGdy+∫0xρGvIdx.After the possession of velocities, temperatures, and water vapour mass fraction fields (u, v, T, and w), the evaluation of heat and mass transfer is proceeded through nondimensional numbers defined below.(i)
The local sensible Nusselt number to evaluate sensible heat transfer at the interface:(19)Nus,x=Qs,IDhλGTI-Tb=DhTI-Tb∂T∂yG,I.(ii)
The local latent Nusselt number to evaluate latent heat transfer at the interface:(20)NuL,x=QL,IDhλGTI-Tb=mI˙hfgDhλGTI-Tb.(iii)
The local Sherwood number to evaluate mass transfer at the interface:(21)Shx=hTDhD=m˙I1-wIDhρGDwI-wb.(iv)
The accumulated evaporation rate Mr:(22)Mr=1Γ0∫0xm˙Idx.The bulk quantities are evaluated by the following calculation:(23)Tb=∫0b-δxρGCpGuGTGdy∫0b-δxρGCpGuGdy,wb=∫0b-δxρGuGwdy∫0b-δxρGuGdy.


### 2.3. Solution Method

To resolve the parabolic equations ((1) to (8)), a fully implicit scheme is adopted. The numerical solution is achieved using finite difference method. Each finite difference equation’s system forms a tridiagonal matrix, which can be solved using the TDMA method (Patankar [27]). It is still necessary to satisfy the global mass flow constraint. This is done by pressure gradient and axial velocity profile correction at each axial step, according to Raithby and Schneider [28] method. Moreover, to ensure the accuracy of the numerical computation, a nonuniform grid in both axial and radial directions is employed. Accordingly, the grid is refined at the interface and the tube inlet.


#### 2.3.1. Marching Procedure

After the initialization of liquid and gas inlet values, the numerical computation is solved systematically. The following steps are involved in the developed numerical model.ForI=1, the liquid film longitudinal velocity is obtained by solving the momentum equation neglecting the inertia and convection terms. The initial liquid film thickness is then calculated by (24)δ0=3μLΓ0ρL2g1/3.For I=2 to the last x-direction node, we have the following:(1)
For the first iteration, an arbitrary pressure gradientdp/dx is assumed.(2)
An arbitrary value of the liquid film thicknessδxit is assumed for the I+1 section.(3)
Calculate the longitudinal velocities in the liquid and vapour phases in theI+1 section using (2) and (3). The interface velocity is obtained from (15).(4)
Integrate continuity equation to obtain the transversal velocities in both phases.(5)
The temperatures in the liquid film and gas flow along with mass fraction in the gas phase are then calculated from (3), (6), and (7).(6)
Calculate the relative errorELit of the liquid film mass balance using (17).(7)
A better approximation to the liquid film thickness is then obtained using the secant method. Thus,(25)δxit+1=δxit-δxit-δxit-1ELit-ELit-1ELit.The convergence criteria used areELit=10-6. Usually five to seven iterations suffice to obtain converged solution.(8)
The mass flow in vapour phase is then calculated and for a given inlet vapour velocity the error in the mass flowEGit is determined by (18).(9)
A better approximation to the pressure gradient is then obtained using the secant method. Thus,(26)dpdxit+1=dpdxit-dp/dxit-dp/dxit-1EGit-EGit-1EGit.Five to six iterations are usually sufficient to satisfy the convergence criteria for the error in the mass balance(10)
Check the convergence of all calculated variables: longitudinal velocities, species mass fractions, and temperatures. If the convergence criteria (Err≤10-6) are satisfied, then the new thermophysical properties are calculated; if not, the steps (1) to (9) are repeated for a new iteration (it+1).The calculations are then advanced in the longitudinal directionx.


#### 2.3.2. Velocity-Pressure Coupling

In order to satisfy the global mass flow constraint, the correction of the pressure gradient and axial velocity profile at each axial position is achieved using the method proposed by Raithby and Schneider [28], described by Tannehill et al. [29].To illustrate, we will letS=-dp/dx. We make an initial guess for the pressure gradient noted -dP/dx∗ and calculate the corresponding guessed gas mass flow rate M˙i∗:(27)M˙i∗=∫0b-δxρui,j∗dy.Due to the linearity of the momentum equation with frozen coefficients, the correct velocity at each point from an application of Newton’s method would be(28)ui,j=ui,j∗-∂ui,j∂S∂P∂x-∂P∂x∗=ui,j∗+UΔS,where ΔS=-∂P/∂x-∂P/∂x∗ and U=∂ui,j/∂S.Δ S is the change in the pressure gradient required to satisfy the global mass flow constraint. The differential equations are actually discretised with respect to the pressure gradient (S) to obtain difference equations for U. The Thomas algorithm is used to solve the system of algebraic equations for U. The coefficients for the unknowns in these equations will be the same as for the original implicit discretised equations for ui,j (momentum equation).The boundary conditions onU must be consistent with the velocity boundary conditions. On boundaries where the velocity is specified, U=0, while ∂U/∂n=0 (n normal boundary) where a velocity gradient is specified. The solution for ui,j is then used to compute ΔS by noting that UΔS is the correction in velocity at each point required to satisfy the global mass flow constraint. Thus, we can write(29)ΔS=M˙i-M˙i∗∫0b-δxρui,jdy,where the integral is evaluated by numerical means. The M˙i in the previous equation is the known value specified by the initial conditions. The required value of ΔS is determined from (29). The correct values of velocity uj can then be determined from (28). The continuity equation is then integrated to determinevjn+1.


#### 2.3.3. Grid Dependency Test

A grid dependency test is done to avoid convergence problems due to thin grids use and to find the optimum solution between computational time and results precision. Several grid sizes have been tested (Table1). It was found that in all grid sizes the difference in local Nusselt and Sherwood numbers is always less than 5%. The grid with 501×(201+101) is chosen because it gives results close enough to those of the thin grid and sufficiently accurate to describe the heat and mass transfer and to generate a good wall curvature.Table 1
Test of grid dependency by comparing the local interfacial Nusselt and Sherwood numbers of various grid arrangements (conditionsQw=5000 W·m−2; Γ0=0.02 kg·m−1·s−1; T0=25°C; TL0=60°C; Hr=0%; Re=1000, a=0.0004m, n=20).
Grid size NI × (NG + NL) x ( m ) 0.25 0.5 0.75 1 1.5 201 × ( 101 + 51 ) Nus,x 1.146 1.756 2.227 2.560 3.156 NuL,x 16.132 29.640 45.943 62.807 108.028 Shx 1.339 2.186 2.918 3.446 4.316 201 × ( 201 + 101 ) Nus,x 1.196 1.816 2.298 2.635 3.236 NuL,x 16.638 30.170 46.459 63.190 107.952 Shx 1.392 2.253 3.000 3.539 4.428 501 × ( 101 + 51 ) Nus,x 1.121 1.728 2.187 2.518 3.100 NuL,x 15.832 29.293 45.201 62.113 107.360 Shx 1.318 2.165 2.881 3.409 4.258 501 × ( 201 + 101 ) Nus,x 1.176 1.796 2.265 2.601 3.190 NuL,x 16.394 29.885 45.768 62.509 107.115 Shx 1.377 2.240 2.972 3.511 4.382 701 × ( 101 + 51 ) Nus,x 1.116 1.726 2.178 2.507 3.085 NuL,x 15.776 29.332 45.095 61.973 107.295 Shx 1.313 2.166 2.873 3.398 4.239 701 × ( 201 + 101 ) Nus,x 1.173 1.796 2.258 2.593 3.178 NuL,x 16.353 29.939 45.665 62.349 106.923 Shx 1.374 2.243 2.966 3.503 4.367


## 2.1. Governing Equations

To formulate the governing equations, the following assumptions are done:(i)
The flow in the two phases is considered laminar, incompressible, and of Newtonian fluids.(ii)
The problem is two-dimensional and steady.(iii)
The thermodynamic and thermophysical properties of air and air-water vapour mixture are a function of both local temperature and concentration.(iv)
The film surface is in thermodynamic equilibrium and smooth.(v)
Heat transfer by radiation, the effect of the surface tension, Dufour, Soret, and other secondary effects are negligible.(vi)
The axial diffusion of heat, momentum, and mass were assumed to be negligible.(vii)
The vapour-gas mixture is an ideal gas mixture and the saturation conditions are evaluated at the liquid-mixture interface.The pure component data are approximated by polynomials in function of temperature and mass fraction. For further details, the thermophysical properties are available in [25, 26].Under such assumptions, the governing equations in the liquid phase could be written as(1)∂∂xρLuL+∂∂yρLvL=0,(2)uL∂∂xρLuL+vL∂∂yρLuL=-∂P∂x+∂∂yμL∂uL∂y+ρLg,(3)uL∂ρLCPLTL∂x+vL∂ρLCPLTL∂y=∂∂yλL∂TL∂y.Similarly, the conservation equations governing the gaseous phase are(4)∂∂xρGuG+∂∂yρGvG=0,(5)uG∂ρGuG∂x+vG∂ρGuG∂y=-∂P∂x+∂∂yμG∂uG∂y+ρG-ρ0g,(6)uG∂ρGCPGTG∂x+vG∂ρGCPGTG∂y=∂∂yλG∂TG∂y+∂∂yDCPV-CPa∂w∂yTG,(7)ρGuG∂w∂x+vG∂w∂y=∂∂yD∂w∂y.The curvature of the wavy wall is described by the following sinusoidal function:(8)ϖx=asin2πxl,where a is the amplitude of the wave and l its characteristic length. Consequently, the number of waves is defined as n=L/l.Accordingly, the effective width of the channel at eachx-location is(9)b=d-ϖx.


## 2.2. Boundary Conditions

The governing equations are subject to the following boundary conditions:(i)
At the inlet (x=0):(10)TG=T0;uG=u0;w=w0;Hr=Hr0,P=P0;TL=TL0;Γ=Γ0.(ii)
At the heated wall (y=b):(11)uL=0;λL∂TL∂y=QW.(iii)
At the insulated wall (y=0):(12)uG=0;∂TG∂y=0;∂w∂y=0.The matching conditions at the interfacey=b-δ are as follows:(13)uIx=uG,I=uL,I,(14)TIx=TG=TL,I,(15)τI=μ∂u∂yL,I=μ∂u∂yG,I.Energy balance at the gas-liquid interface is(16)λ∂T∂yL,I=λ∂T∂yG,I+m˙hfg,where m˙=-ρGvI=ρGD/1-wI∂w/∂y and hfg is the latent heat of vapourization.In addition, we should add mass conservation equations to compute the unknown pressure gradientdP/dx:(17)Γ0=∫b-δbρudyL-∫0xρGvIdx,(18)b-δ0ρ0u0=∫0b-δρGuGdy+∫0xρGvIdx.After the possession of velocities, temperatures, and water vapour mass fraction fields (u, v, T, and w), the evaluation of heat and mass transfer is proceeded through nondimensional numbers defined below.(i)
The local sensible Nusselt number to evaluate sensible heat transfer at the interface:(19)Nus,x=Qs,IDhλGTI-Tb=DhTI-Tb∂T∂yG,I.(ii)
The local latent Nusselt number to evaluate latent heat transfer at the interface:(20)NuL,x=QL,IDhλGTI-Tb=mI˙hfgDhλGTI-Tb.(iii)
The local Sherwood number to evaluate mass transfer at the interface:(21)Shx=hTDhD=m˙I1-wIDhρGDwI-wb.(iv)
The accumulated evaporation rate Mr:(22)Mr=1Γ0∫0xm˙Idx.The bulk quantities are evaluated by the following calculation:(23)Tb=∫0b-δxρGCpGuGTGdy∫0b-δxρGCpGuGdy,wb=∫0b-δxρGuGwdy∫0b-δxρGuGdy.


## 2.3. Solution Method

To resolve the parabolic equations ((1) to (8)), a fully implicit scheme is adopted. The numerical solution is achieved using finite difference method. Each finite difference equation’s system forms a tridiagonal matrix, which can be solved using the TDMA method (Patankar [27]). It is still necessary to satisfy the global mass flow constraint. This is done by pressure gradient and axial velocity profile correction at each axial step, according to Raithby and Schneider [28] method. Moreover, to ensure the accuracy of the numerical computation, a nonuniform grid in both axial and radial directions is employed. Accordingly, the grid is refined at the interface and the tube inlet.


### 2.3.1. Marching Procedure

After the initialization of liquid and gas inlet values, the numerical computation is solved systematically. The following steps are involved in the developed numerical model.ForI=1, the liquid film longitudinal velocity is obtained by solving the momentum equation neglecting the inertia and convection terms. The initial liquid film thickness is then calculated by (24)δ0=3μLΓ0ρL2g1/3.For I=2 to the last x-direction node, we have the following:(1)
For the first iteration, an arbitrary pressure gradientdp/dx is assumed.(2)
An arbitrary value of the liquid film thicknessδxit is assumed for the I+1 section.(3)
Calculate the longitudinal velocities in the liquid and vapour phases in theI+1 section using (2) and (3). The interface velocity is obtained from (15).(4)
Integrate continuity equation to obtain the transversal velocities in both phases.(5)
The temperatures in the liquid film and gas flow along with mass fraction in the gas phase are then calculated from (3), (6), and (7).(6)
Calculate the relative errorELit of the liquid film mass balance using (17).(7)
A better approximation to the liquid film thickness is then obtained using the secant method. Thus,(25)δxit+1=δxit-δxit-δxit-1ELit-ELit-1ELit.The convergence criteria used areELit=10-6. Usually five to seven iterations suffice to obtain converged solution.(8)
The mass flow in vapour phase is then calculated and for a given inlet vapour velocity the error in the mass flowEGit is determined by (18).(9)
A better approximation to the pressure gradient is then obtained using the secant method. Thus,(26)dpdxit+1=dpdxit-dp/dxit-dp/dxit-1EGit-EGit-1EGit.Five to six iterations are usually sufficient to satisfy the convergence criteria for the error in the mass balance(10)
Check the convergence of all calculated variables: longitudinal velocities, species mass fractions, and temperatures. If the convergence criteria (Err≤10-6) are satisfied, then the new thermophysical properties are calculated; if not, the steps (1) to (9) are repeated for a new iteration (it+1).The calculations are then advanced in the longitudinal directionx.


### 2.3.2. Velocity-Pressure Coupling

In order to satisfy the global mass flow constraint, the correction of the pressure gradient and axial velocity profile at each axial position is achieved using the method proposed by Raithby and Schneider [28], described by Tannehill et al. [29].To illustrate, we will letS=-dp/dx. We make an initial guess for the pressure gradient noted -dP/dx∗ and calculate the corresponding guessed gas mass flow rate M˙i∗:(27)M˙i∗=∫0b-δxρui,j∗dy.Due to the linearity of the momentum equation with frozen coefficients, the correct velocity at each point from an application of Newton’s method would be(28)ui,j=ui,j∗-∂ui,j∂S∂P∂x-∂P∂x∗=ui,j∗+UΔS,where ΔS=-∂P/∂x-∂P/∂x∗ and U=∂ui,j/∂S.Δ S is the change in the pressure gradient required to satisfy the global mass flow constraint. The differential equations are actually discretised with respect to the pressure gradient (S) to obtain difference equations for U. The Thomas algorithm is used to solve the system of algebraic equations for U. The coefficients for the unknowns in these equations will be the same as for the original implicit discretised equations for ui,j (momentum equation).The boundary conditions onU must be consistent with the velocity boundary conditions. On boundaries where the velocity is specified, U=0, while ∂U/∂n=0 (n normal boundary) where a velocity gradient is specified. The solution for ui,j is then used to compute ΔS by noting that UΔS is the correction in velocity at each point required to satisfy the global mass flow constraint. Thus, we can write(29)ΔS=M˙i-M˙i∗∫0b-δxρui,jdy,where the integral is evaluated by numerical means. The M˙i in the previous equation is the known value specified by the initial conditions. The required value of ΔS is determined from (29). The correct values of velocity uj can then be determined from (28). The continuity equation is then integrated to determinevjn+1.


### 2.3.3. Grid Dependency Test

A grid dependency test is done to avoid convergence problems due to thin grids use and to find the optimum solution between computational time and results precision. Several grid sizes have been tested (Table1). It was found that in all grid sizes the difference in local Nusselt and Sherwood numbers is always less than 5%. The grid with 501×(201+101) is chosen because it gives results close enough to those of the thin grid and sufficiently accurate to describe the heat and mass transfer and to generate a good wall curvature.Table 1
Test of grid dependency by comparing the local interfacial Nusselt and Sherwood numbers of various grid arrangements (conditionsQw=5000 W·m−2; Γ0=0.02 kg·m−1·s−1; T0=25°C; TL0=60°C; Hr=0%; Re=1000, a=0.0004m, n=20).
Grid size NI × (NG + NL) x ( m ) 0.25 0.5 0.75 1 1.5 201 × ( 101 + 51 ) Nus,x 1.146 1.756 2.227 2.560 3.156 NuL,x 16.132 29.640 45.943 62.807 108.028 Shx 1.339 2.186 2.918 3.446 4.316 201 × ( 201 + 101 ) Nus,x 1.196 1.816 2.298 2.635 3.236 NuL,x 16.638 30.170 46.459 63.190 107.952 Shx 1.392 2.253 3.000 3.539 4.428 501 × ( 101 + 51 ) Nus,x 1.121 1.728 2.187 2.518 3.100 NuL,x 15.832 29.293 45.201 62.113 107.360 Shx 1.318 2.165 2.881 3.409 4.258 501 × ( 201 + 101 ) Nus,x 1.176 1.796 2.265 2.601 3.190 NuL,x 16.394 29.885 45.768 62.509 107.115 Shx 1.377 2.240 2.972 3.511 4.382 701 × ( 101 + 51 ) Nus,x 1.116 1.726 2.178 2.507 3.085 NuL,x 15.776 29.332 45.095 61.973 107.295 Shx 1.313 2.166 2.873 3.398 4.239 701 × ( 201 + 101 ) Nus,x 1.173 1.796 2.258 2.593 3.178 NuL,x 16.353 29.939 45.665 62.349 106.923 Shx 1.374 2.243 2.966 3.503 4.367


## 2.3.1. Marching Procedure

After the initialization of liquid and gas inlet values, the numerical computation is solved systematically. The following steps are involved in the developed numerical model.ForI=1, the liquid film longitudinal velocity is obtained by solving the momentum equation neglecting the inertia and convection terms. The initial liquid film thickness is then calculated by (24)δ0=3μLΓ0ρL2g1/3.For I=2 to the last x-direction node, we have the following:(1)
For the first iteration, an arbitrary pressure gradientdp/dx is assumed.(2)
An arbitrary value of the liquid film thicknessδxit is assumed for the I+1 section.(3)
Calculate the longitudinal velocities in the liquid and vapour phases in theI+1 section using (2) and (3). The interface velocity is obtained from (15).(4)
Integrate continuity equation to obtain the transversal velocities in both phases.(5)
The temperatures in the liquid film and gas flow along with mass fraction in the gas phase are then calculated from (3), (6), and (7).(6)
Calculate the relative errorELit of the liquid film mass balance using (17).(7)
A better approximation to the liquid film thickness is then obtained using the secant method. Thus,(25)δxit+1=δxit-δxit-δxit-1ELit-ELit-1ELit.The convergence criteria used areELit=10-6. Usually five to seven iterations suffice to obtain converged solution.(8)
The mass flow in vapour phase is then calculated and for a given inlet vapour velocity the error in the mass flowEGit is determined by (18).(9)
A better approximation to the pressure gradient is then obtained using the secant method. Thus,(26)dpdxit+1=dpdxit-dp/dxit-dp/dxit-1EGit-EGit-1EGit.Five to six iterations are usually sufficient to satisfy the convergence criteria for the error in the mass balance(10)
Check the convergence of all calculated variables: longitudinal velocities, species mass fractions, and temperatures. If the convergence criteria (Err≤10-6) are satisfied, then the new thermophysical properties are calculated; if not, the steps (1) to (9) are repeated for a new iteration (it+1).The calculations are then advanced in the longitudinal directionx.


## 2.3.2. Velocity-Pressure Coupling

In order to satisfy the global mass flow constraint, the correction of the pressure gradient and axial velocity profile at each axial position is achieved using the method proposed by Raithby and Schneider [28], described by Tannehill et al. [29].To illustrate, we will letS=-dp/dx. We make an initial guess for the pressure gradient noted -dP/dx∗ and calculate the corresponding guessed gas mass flow rate M˙i∗:(27)M˙i∗=∫0b-δxρui,j∗dy.Due to the linearity of the momentum equation with frozen coefficients, the correct velocity at each point from an application of Newton’s method would be(28)ui,j=ui,j∗-∂ui,j∂S∂P∂x-∂P∂x∗=ui,j∗+UΔS,where ΔS=-∂P/∂x-∂P/∂x∗ and U=∂ui,j/∂S.Δ S is the change in the pressure gradient required to satisfy the global mass flow constraint. The differential equations are actually discretised with respect to the pressure gradient (S) to obtain difference equations for U. The Thomas algorithm is used to solve the system of algebraic equations for U. The coefficients for the unknowns in these equations will be the same as for the original implicit discretised equations for ui,j (momentum equation).The boundary conditions onU must be consistent with the velocity boundary conditions. On boundaries where the velocity is specified, U=0, while ∂U/∂n=0 (n normal boundary) where a velocity gradient is specified. The solution for ui,j is then used to compute ΔS by noting that UΔS is the correction in velocity at each point required to satisfy the global mass flow constraint. Thus, we can write(29)ΔS=M˙i-M˙i∗∫0b-δxρui,jdy,where the integral is evaluated by numerical means. The M˙i in the previous equation is the known value specified by the initial conditions. The required value of ΔS is determined from (29). The correct values of velocity uj can then be determined from (28). The continuity equation is then integrated to determinevjn+1.


## 2.3.3. Grid Dependency Test

A grid dependency test is done to avoid convergence problems due to thin grids use and to find the optimum solution between computational time and results precision. Several grid sizes have been tested (Table1). It was found that in all grid sizes the difference in local Nusselt and Sherwood numbers is always less than 5%. The grid with 501×(201+101) is chosen because it gives results close enough to those of the thin grid and sufficiently accurate to describe the heat and mass transfer and to generate a good wall curvature.Table 1
Test of grid dependency by comparing the local interfacial Nusselt and Sherwood numbers of various grid arrangements (conditionsQw=5000 W·m−2; Γ0=0.02 kg·m−1·s−1; T0=25°C; TL0=60°C; Hr=0%; Re=1000, a=0.0004m, n=20).
Grid size NI × (NG + NL) x ( m ) 0.25 0.5 0.75 1 1.5 201 × ( 101 + 51 ) Nus,x 1.146 1.756 2.227 2.560 3.156 NuL,x 16.132 29.640 45.943 62.807 108.028 Shx 1.339 2.186 2.918 3.446 4.316 201 × ( 201 + 101 ) Nus,x 1.196 1.816 2.298 2.635 3.236 NuL,x 16.638 30.170 46.459 63.190 107.952 Shx 1.392 2.253 3.000 3.539 4.428 501 × ( 101 + 51 ) Nus,x 1.121 1.728 2.187 2.518 3.100 NuL,x 15.832 29.293 45.201 62.113 107.360 Shx 1.318 2.165 2.881 3.409 4.258 501 × ( 201 + 101 ) Nus,x 1.176 1.796 2.265 2.601 3.190 NuL,x 16.394 29.885 45.768 62.509 107.115 Shx 1.377 2.240 2.972 3.511 4.382 701 × ( 101 + 51 ) Nus,x 1.116 1.726 2.178 2.507 3.085 NuL,x 15.776 29.332 45.095 61.973 107.295 Shx 1.313 2.166 2.873 3.398 4.239 701 × ( 201 + 101 ) Nus,x 1.173 1.796 2.258 2.593 3.178 NuL,x 16.353 29.939 45.665 62.349 106.923 Shx 1.374 2.243 2.966 3.503 4.367


## 3. Results and Discussion

Before tackling the numerical results of the simulations, it is paramount to test the accuracy and validity of the computations. Our code results are compared with the case of the evaporation of pure water and ethanol liquid films in a vertical channel with respect to the same conditions given in the experimental and numerical studies. We first compared our result with the numerical study of Yan [5] as shown in Figure 2(a). The comparison shows the variation of interface temperature along the channel for different inlet mass flows. A very satisfying agreement is found between the current computational study and the simulation of Yan (maximum relative error of 1.9%). Moreover, the comparison is done with the evaporative cooling experiment of Yan et al. [2] in an adiabatic vertical channel plotted in Figure 2(b). A good agreement between our simulations and the experiment curves is found with a maximum relative error of 4.7%. Accordingly, the present numerical algorithm is validated and can be considered accurate for the practical purpose. The difference found for Γ0=0.01 kg·m−1·s−1 in Figure 2(b) is due to the difficulty in maintaining null ethanol vapour fraction in the gas phase during the experiment as stated by the authors [3].Figure 2
(a) Local temperature distribution along the liquid film interface, and (b) local temperature distribution along the heated channel wall, comparisons of simulations with previous numerical studies.
(a)
Pure water (b)
EthanolIn this study, calculations are specifically performed for pure water film evaporation over a wavy vertical channel. The following sets of conditions are selected in the computation: the dry air enters atT0=25°C; Re = 2000 and 1 atm. The liquid film flows over the wavy wall at TL0=60°C and Γ0=0.02 kg·m−1·s−1. The wall is subjected to a uniform heat flux Qw=4000 W·m−2.A first overview of the results shows that all plotted parameters take a wavy form when the wall is wavy. However, the wavy form of those parameters does not give clear information of the studied parameters variation in most cases. To overcome this issue, we treat the curved parameters as a signal. Thus, we determine for each curve the upper and lower envelope and finally the mean of both envelopes. The obtained value (with subscriptxc) represents the equilibrium line of the curves evolution at each location.To compare the effect of the wall wavelength and amplitude on the heat and mass transfer, the variations of liquid film thickness and local nondimensional parameters are plotted. It is noted that small amplitude values are chosen to ensure smooth interface condition and the adherence of the liquid film to the wall. In fact, very large amplitudes amplify turbulent effects and may create a discontinuity of the liquid film. From an economic point of view, large amplitudes increase the weight of the evaporator and consequently the process cost. Figure3 shows the effect of the waves parameters on the liquid film thickness along the channel. For all cases, the liquid film thickness decreases along the channel due to the evaporation. It is also observed that the liquid film thickness has a periodical variation and takes also crests and troughs along the channel, but their amplitudes decrease from the inlet to x=0.25m and then increase along the channel. Moreover, the liquid film is less thick for the largest number and the greatest amplitude of the wall waves. This indicates that the evaporation of water could be enhanced by incrementing the waves number or magnitude. The effect of wall wave and amplitude is clearly shown in Figure 3(b) where the equilibrium lines of liquid film thickness are compared to the flat wall case. The liquid film thickness is smaller for n=40 and a=0.0004m. In fact, doubling the waves number increases the surface of the wall more than doubling the wave amplitude, although having a larger surface enhances the heat exchange between the liquid film and the wall, and thus the liquid film evaporation.Figure 3
Evolution of local liquid film thickness along the channel.To improve our understanding of the heat and mass transfer occurring at the interface of liquid-vapour phases, the variations of local sensible and latent Nusselt numbers in addition to the local Sherwood number at the interface are presented in Figures4, 5, and 6. The local sensible Nusselt number oscillates along the channel with a decreasing wave amplitude downstream the channel. The latent Nusselt number has also an oscillating trend but unlike Nus,x, the amplitude of the wave increases going to the outlet. Moreover, the Nusselt number decreases continuously along the channel for Nus,x, while it decreases and then increases slightly for NuL,x (Figures 4(b) and 5(b)). The upper and lower peaks of the Nusselt and Sherwood numbers are observed in the upper and lower peaks of the wall, respectively. The larger Nus,x and NuL,x are observed for large amplitude and waves numbers. It is worth noticing that the heat exchange by evaporation is more effective at the interface since latent Nusselt number is much important than Nus,x. In fact, the liquid film is heated gradually along the channel and the evaporation increases gradually from the inlet to the outlet. Consequently, the heat transfer by latent heat is incremented on the one hand, while the sensible heat flux decreases on the other hand. Concerning local Sherwood number (Figure 6), it can be seen that its variation is similar to the variation of local sensible Nusselt number along the channel. To explain the higher and lower values of Nus,x, NuL,x, and Shx located at the upper and lower peaks of the wavy wall, respectively, one must refer to the velocity of the flowing fluids. In fact, the section of the channel takes the smallest value at the upper peaks, which means that, for a given gas flow Reynolds number, the average velocity and velocity gradient increase at those locations leading to the enhancement of heat transfer. In contrast, at the lower peaks locations, the section of the channel has the smaller value and hence the velocity and velocity gradient decrease. Those effects become more important by increasing the wave amplitude and waves number. However, the highest local values of Nus,x and Shx correspond to the first wave. Afterward, the flows begin to be periodically fully developed. As a result, Nus,x and Shx start to take uniform periodic values at the rest of the channel. Oppositely, the highest local latent Nusselt number takes place at the last wave. This can be explained by the fact that as the liquid film flows downstream the channel, its temperature approaches the saturation temperature and gives rise to the evaporation at the interface.Figure 4
Evolution of local sensible Nusselt number at the interface for multiple waves amplitude and number along the channel.Figure 5
Evolution of local latent Nusselt number at the interface for multiple waves amplitude and number along the channel.Figure 6
Evolution of local Sherwood number at the interface for multiple waves amplitude and number along the channel.Figure7 shows the variation of accumulated evaporation rate along the channel for various wavelengths and amplitudes. The accumulated evaporation rate increases continuously along the channel due to the continuous increasing evaporation. It is observed that the higher evaporation rate occurs for the large wave amplitude and bigger waves number. In fact, this trend of Mr confirms the findings about the liquid film thickness. When evaporation increases the Mr increases and thus liquid films thickness decreases.Figure 7
Evolution of accumulated evaporation rate at the interface for multiple waves amplitude and number along the channel.After describing the effect of wavelength and amplitude on the local thermal and mass characteristics, the next step is to investigate the effect of both wavelength and amplitude, on the average Nusselt and Sherwood numbers, and compare with the case of flat wall presented by dashed line (NuL,av∗=85.73). It is clearly shown in Figure 8 that for any wavelength and amplitude the average latent Nusselt number is always higher than NuL,av∗. However, it can be admitted that for the channel with small wave amplitude, the NuL,av is slightly larger than the flat channel and the waves number does not contribute significantly on increasing it. Therefore, with the increase of the waves amplitude, the effect of waves number becomes significant.Figure 8
Average latent Nusselt number at the interface in function of waves number for multiple waves amplitudes.The behaviour of average sensible Nusselt number toward wavelength and amplitude is completely different fromNuL,av as shown in Figure 9. The Nus,av for low waves number and small amplitudes is under the average sensible Nusselt number of flat channel Nus,av∗=3.49. It can be seen that, for a=0.0001m, the average sensible Nusselt number increases slightly by increasing the waves number but does not reach the Nus,av∗ of the flat channel, while for channels with a=0.0002m and a=0.0004m, their Nus,av exceed Nus,av∗ for n>26 and n>12, respectively. Concerning the average sensible number presented in Figure 10, it is always high than the flat channel’s average Sherwood number Shav∗=4.52 for the largest waves amplitude. For the small value of wave amplitude, there is no significant difference and a very slight increase of Shav by the increase of n compared to flat channel. For a=0.0002m, the Sherwood number is increased only for n>14.Figure 9
Average sensible Nusselt number at the interface in function of waves number for multiple wave amplitudes.Figure 10
Average Sherwood number at the interface in function of waves number for multiple wave amplitudes.


## 4. Conclusion

The study investigated the problem of convective heat and mass transfer inside a vertical channel during liquid film evaporation on a vertical wavy channel. The effects of the wall waves amplitudes and number on the interfacial Nusselt and Sherwood numbers are numerically analysed. The results indicated that all studied parameters vary periodically and similarly to the wall shape. Consequently, maximum values are observed at the crests of the wall and minimum values at the troughs. It was found that the local sensible Nusselt and Sherwood numbers tend to take periodically uniform values near the outlet since the flow becomes to be periodically fully developed except for the local Nusselt number. The variation of average values of Nusselt and Sherwood number with waves number and amplitude was also investigated. The results highlighted that, for very small waves amplitudes of the wavy wall, the average sensible Nusselt and Sherwood numbers may be smaller than the flat wall values. In this case, the increase of waves numbers has no significant contribution to the sensible heat and mass transfer.For the lowest wave amplitudes (a=0.0001m and a=0.0002), the average latent Nusselt number and average Sherwood number were enhanced by a maximum of 1.1% compared to the flat wall. For the high wave amplitude, the average latent Nusselt number enhancement reached a maximum of 10%, while the average Sherwood number enhancement reached a maximum of 21.5%. However, the average sensible Nusselt number was exclusively enhanced when (n > 26) and (n > 12) for a=0.0002 m and 0.0004 m, respectively. For a=0.0001m, the average sensible Nusselt number of the wavy wall is slightly reduced in comparison to the flat wall. Accordingly, the averages of the studied dimensionless numbers increase with the increase of the number of wall waves.



---

*Source: 4208059-2018-04-02.xml*