energies Article CFD Analysis of Falling Film Hydrodynamics for a Lithium Bromide (LiBr) Solution over a Horizontal Tube Furqan Tahir 1, * , Abdelnasser Mabrouk 1,2 and Muammer Koç 1 1 2 * Division of Sustainable Development (DSD), College of Science & Engineering (CSE), Hamad Bin Khalifa University (HBKU), Education City, Doha 34110, Qatar; aaboukhlewa@hbku.edu.qa (A.M.); mkoc@hbku.edu.qa (M.K.) Qatar Environment and Energy Research Institute (QEERI), Hamad Bin Khalifa University (HBKU), Education City, Doha 34110, Qatar Correspondence: ftahir@mail.hbku.edu.qa; Tel.: +974-30074879 Received: 30 November 2019; Accepted: 6 January 2020; Published: 8 January 2020



 Abstract: Falling film evaporators are used in applications where high heat transfer coefficients are required for low liquid load and temperature difference. One such application is the lithium bromide (LiBr)-based absorber and generator. The concentration of the aqueous LiBr solution changes within the absorber and generator because of evaporation and vapor absorption. This causes the thermophysical properties to differ and affects the film distribution, heat, and mass transfer mechanisms. For thermal performance improvement of LiBr-based falling film evaporators, in-depth analysis at the micro level is required for film distribution and hydrodynamics. In this work, a 2D numerical model was constructed using the commercial CFD software Ansys Fluent v18.0. The influence of the liquid load corresponding to droplet and jet mode, and the concentration, on film hydrodynamics was examined. It was found that the jet mode was more stable at a higher concentration of 0.65 with ±0.5% variation compared to lower concentrations. The recirculation was stronger at a low concentration of 0.45 and existed until the angular position (θ) = 10◦ , whereas at 0.65 concentration it diminished after θ = 5◦ . The improved heat transfer is expected at lower concentrations due to lower film thickness and thermal resistance, more recirculation, and a higher velocity field. Keywords: falling film; film thickness; CFD; horizontal tube; lithium bromide; LiBr; transient analysis; VOF 1. Introduction Falling film exchangers are commonly employed in ammonia (NH3 )- and lithium bromide (LiBr)-based absorption refrigeration systems [1–6]. In falling film exchangers, a thin layer of liquid film is formed around the tube and heat is transferred to/from the fluid flowing inside the tube. The heat and mass transfer coefficients strongly depend on the film hydrodynamics. Extensive experimental, analytical, and numerical works were performed to describe the falling film flow on bare and enhanced tubes [7,8]. The falling film characterization comprises of (a) film thickness, (b) wave formation, (c) wettability, (d) flow development regions, such as stagnation, (e) impingement, (f) developing, developed detachment, and necking zones, (g) flow modes, such as droplet, jet, and sheet, and (h) flow regimes, which include laminar, wavy laminar, and turbulent [7,9–13]. The falling film modes were experimentally distinguished: droplet, droplet-jet, jet, jet-sheet, and sheet modes. Researchers have developed empirical correlations based on films’ Reynolds number Ref and modified Galileo number Ga for flow mode recognition and transition [14–17]. The correlations developed by Roques et al. [16] are as follows: Energies 2020, 13, 307; doi:10.3390/en13020307 www.mdpi.com/journal/energies Energies 2020, 13, 307 2 of 15 Droplet to droplet-jet mode: Re f = 0.0417Ga0.3278 . (1) Re f = 0.0683Ga0.3204 . (2) Re f = 0.8553Ga0.2483 . (3) Re f = 1.068Ga0.2563 . (4) Droplet-jet to jet mode: Jet to jet-sheet mode: Jet-sheet to sheet mode: Ref and Ga are defined as: Re f = Ga = 4Γ1/2 µ ρσ3 µ4 g , (5) (6) where Γ1/2 represents the liquid load on one side of the tube, µ is the viscosity, g is the acceleration due to gravity, ρ is the density, and σ is the surface tension. In 1916, Nusselt established the film thickness δ correlation as a function of angular position θ, liquid load Γ1/2 , and thermophysical properties [18]. The Nusselt solution, as shown in Equation (7), is widely used for thickness estimation and comparison purposes [9,19–21]: s 3µl Γ1/2  δ= 3  . (7) ρl ρl − ρ g gsinθ Rogers and Goindi [22] measured the film thickness on an aluminum tube of 132 mm outer tube diameter D at three different angular positions, θ = 45◦ , 90◦ , and 135◦ . They found a 30% deviation on averaged when compared with Equation (7). Gstoehl et al. [23] applied the laser method to analyze the film thickness around the tube. The film thicknesses θ = 22◦ –62◦ and 112◦ –152◦ on a 19.05-mm copper tube were recorded for different intertube spacing s. They concluded that Equation (7) predicts the film thickness well for the upper region (θ = 0◦ –90◦ ), but overestimates the thickness for the lower region (θ > 90◦ ). Hou et al. [24] conducted experiments to measure film thickness using a displacement micrometer from θ = 15◦ –165◦ for D = 20–32 mm, s = 10–40 mm, and Ref = 150–800. They reached a similar conclusion, that Equation (7) overestimates the film thickness in the lower portion, and proposed a modified version of Equation (7) by incorporating intertube spacing in the tube diameter ratio (s/D). In addition, it was observed that the minimum film thickness lies in the range of 90◦ –110◦ instead of at a fixed location, i.e., θ = 90◦ as per Equation (7). Chen et al. [25] used the laser technique to capture film thickness in the longitudinal and circumferential directions. In their findings, the maximum film thickness was found to be at the crest location, where it was twice the minimum film thickness in the axial direction. Furthermore, higher film thickness for seawater was recorded as compared to water. Since there are limitations in experimental procedures covering the entire parameters and capturing other details such as operating conditions, velocity, and temperature profiles within the film, computational fluid dynamics (CFD) has become an important tool for in-depth analysis to design optimal conditions for various applications. Tahir et al. [11] analyzed the liquid load distribution for each column and concluded that the triangular pitch configuration exhibits more maldistribution than a square pitch arrangement. Ji et al. [19] performed CFD simulations by implementing a volume of fluid (VOF) model to track the liquid/gas interface, with LiBr solution as the working fluid. They found the minimum thickness near 120◦ , so adjusted Equation (7) by replacing sin(θ) with sin(3/4 θ), such that it predicts the minimum thickness at 120◦ . Lin et al. [26] performed experiments and numerical calculations for a three-dimensional geometry, to examine the film distribution. They validated their VOF model with the experimental data and observed the minimum thickness to lie in the valley area, Energies 2020, 13, 307 3 of 15 in which the heat transfer performance was found to poor. Wang et al. [27] also studied film thickness numerically, with a three-dimensional geometry of two 25.4-mm diameter tubes. They found that the minimum film thickness location depends on both the angular position and the axial displacement. Based on their CFD results, they proposed a film thickness correlation similar to Hou et al. [24], but with different parameters for three sections, namely impingement, transition, and departure regions. Hosseinnia et al. [20] carried out a numerical analysis for mass transfer mechanism in the LiBr-based absorber. The VOF model and h-refinement (a mesh adaptation technique) were implemented in order to solve conservation of mass, momentum, and energy equations. Their analysis showed that the mass transfer mechanism is greatly affected by the flow mode, as the absorption rate was decreased more than 10-fold when the flow mode was changed from droplet to jet. For this reason, most of the absorber operates under droplet, droplet-jet, and modes or low liquid load [9,20]. The film distribution around the tube can be divided into stagnation, impingement, developing, developed, detachment, and necking regions. The majority of the heat and mass transfer takes place in the first three regions [9], thus the film distribution affects the heat and mass transfer mechanisms. In – modeled, this study, a 2D CFD model for an aqueous LiBr solution falling over a horizontal tube was which is used in the refrigeration and multi-effect desalination industries [28–31].The LiBr concentration changes in an absorber and generator due to vapor absorption and evaporation. As the concentration varies, the thermophysical properties differ and alter the magnitude of forces and hence the film hydrodynamics. The influence of LiBr concentration on the film hydrodynamics and conduction thermal resistance (which is often neglected) is lacking in the literature. Therefore, the effects of LiBr concentration, with corresponding liquid loads for droplet and jet mode, on film distribution were analyzed. First, the CFD model was validated using experimental data from the literature. Then, the transient behavior of falling film, film thickness distribution, average film thickness, and residence time were quantified and discussed. Afterwards, the tangential and normal velocities distribution for impingement and developing regions for low and high LiBr concentrations C were examined. 2. Problem Statement In order to characterize the film distribution, one 25.4-mm tube was selected with an impingement height hi of 5 mm, as shown in Figure 1. The domain was equipped with a 2-mm orifice at the top of the tube. The width of the domain was 33 mm; since the nature of the problem is symmetric, only half of the domain is selected for CFD calculations. The aqueous LiBr solution is in the liquid phase and water vapors are considered as the gaseous phase. Two LiBr concentrations are chosen as per the common working range, i.e., low concentration C = 0.45 and high concentration C = 0.65. The thermophysical properties of the liquid and gaseous phases were evaluated as per [32] at 30 ◦ C (the lower working temperature of the absorber [33]), as shown in Table 1. Mass inlet Velocity inlet (1 mm orifice) (0 m/s) Symmetry D = 25.4 mm tube wall y 35.4 mm 16.5 mm x Pressure outlet Figure 1. 2D geometry with boundary conditions. Energies 2020, 13, 307 4 of 15 Table 1. Thermophysical properties of LiBr solution and water vapors. Fluid Density (kg/m3 ) Viscosity (Pa-s) Surface Tension (N/m) 1453 1814 0.0304 0.002274 0.01102 0.0000101 0.08 0.085 C = 0.45 C = 0.65 Water Vapors LiBr solution 3. CFD Model For numerical model simplification, it is assumed that the flow is laminar for all cases. There is no fluid spreading in the z-direction. The thermophysical properties are constant throughout the distribution. There is no heat and mass transfer (adiabatic condition). The flow is symmetric along the y-axis. There is no vapor flow, and the wall adhesion contact angle θw is 0◦ (purely hydrophilic) to ensure the complete wetting of the tube [34]. 3.1. Conservation Equations The mass and momentum balances can be found by continuity and Navier-Stokes equations, as shown below:  ⇀ ∂ρ + ∇· ρV = 0 (8) ∂t  ⇀  ⇀⇀ ∂ ρV = ⇀ ⇀ + ∇· ρV V = −∇P + ∇· τ + ρ g + F vol , (9) ∂t where V is the velocity, t is time, P is the pressure, τ is the tensor, and Fvol is the volumetric force. 3.2. Volume of Fluid Model (VOF) Model In falling film evaporators, the liquid and gaseous phases are separated by distinct boundaries. Hirt and Nichols developed a multiphase model named the volume of fluid (VOF) model to accurately capture the liquid phase interface [35]. The VOF model has been widely implemented in many applications such as bubble dynamics, reactors, falling films, capillary flows, flashing, and splashing in tanks to numerically distinguish the liquid/gas interface [36–40]. They modified the continuity equation as in Equation (8) by introducing void fraction α to solve the mass balance for both phase separately.  ⇀ ∂(αi ρi ) + ∇· αi ρi Vi = 0 (i = g or l) ∂t (10) The subscripts g and l represent gaseous and liquid phases, respectively. The density and viscosity are evaluated on the basis of the void fraction. ρ = αl ρl + (1 − αl )ρ g (11) µ = αl µl + (1 − αl )µ g , (12) where α is the void/volume fraction representing how much liquid is present in a liquid-gas mixture. The liquid void fractions αg and αl are defined as follows: α g = 1 − αl . (13) 3.3. Continuum Surface Force (CSF) Model Brackbill et al. [41] developed a continuum surface force (CSF) model as shown in Equation (14) to calculate the volumetric force Fvol . The CSF model improves the characterization of the liquid-gas interface by accommodating surface tension σ effects: ⇀ F vol ρ  ]. = σK∇αl [  1 ρ + ρ g l 2 (14) Energies 2020, 13, 307 5 of 15 𝑛̂ The surface curvature K is computed by 𝐾 gradient = −𝛻. 𝑛̂of surface unit normal vector n̂: K = −∇·n̂. (15) 3.4. Operating and Boundary Conditions For the boundary conditions, (a) a mass inlet for liquid entrance is selected at the top, (b) a velocity inlet with 0 m/s is provided for water vapors, (c) a wall is selected for a tube with a zero contact angle, (d) a pressure outlet is for the outlet at the bottom, and (e) a symmetric boundary condition is chosen for the rest, as shown in Figure 1. The mass flow inlet varied from 0.01 to 0.5 kg/(m·s), with 0.02 kg/(m·s) intervals for both C = 0.45 and 0.65 concentrations, corresponding to droplet, droplet-jet, and jet mode. The films Reynolds number Ref and flow modes for each case are listed in Table 2. The gravity is acting in the negative y-direction with a value of 9.81 m/s2 . Table 2. Films Reynolds number Ref and flow modes for all cases. Γ C Γ 1/2 (kg/(m·s)) Ref 0.01 0.03 0.05 0.01 0.03 0.05 17.6 52.8 87.9 3.6 10.9 18.2 0.45 0.65 Ga Transitional Ref [16] Droplet to Droplet-Jet Droplet-Jet to Jet 2.84 × 109 52.3 73 7.7 × 106 7.5 11 Mode Droplet Droplet-Jet Jet Droplet Droplet-Jet Jet 3.5. Methodology The domain is meshed quadrilateral dominant elements, and fine boundary layers are provided along the y-axis and tube wall in order to accurately capture the film distribution, as shown in Figure 2. The minimum thickness of the boundary layer is 0.02 mm adjacent to the wall. The conservation of mass and momentum are discretized across each cell in a commercial CFD code, i.e., Ansys Fluent v18.0 (Ansys Inc., Canonsburg, [Pa], United States). A second-order upwind scheme is used for Equation (9), semi-implicit method for pressure linked equations (SIMPLE) algorithm is applied for pressure-velocity coupling and pressure staggering option (PRESTO) is employed for pressure approximation. For void fraction gradient, a piecewise linear interpolation method known as geo-reconstruct scheme is ‒ implemented. The gradients are evaluated at the center of cell using a Green-Gauss cell-based approach. The deduced equations are solved for each cell in a first-order implicit routine for an unsteady solution. 0.02 mm Figure 2. Mesh with fine boundary layers for CFD analysis. Energies 2020, 13, 307 6 of 15 3.6. Mesh Dependency Test To obtain better results with the least computational effort, a mesh dependency test was performed. Three different meshes with 21,800, 29,400, and 36,100 elements were chosen and Γ Γ solved for Γ1/2 = 0.01 kg/(m·s) and C = 0.65. Figure 3 shows the film thickness distribution when the film completely covers the tube surface, i.e., t = 2.79 s. As is evident from the film distribution, when the mesh size increased from 21,800 to 29,400 elements, the distribution changed. However, a further increase in the mesh size led to minimal deviation. Hence, to save computational time,ofa mesh 10 μs with 29,400 elements was selected for the numerical calculations. In addition, a time step of of 10 10 μs µs was selected after the time dependency test. Film Thickness, δ (mm) δ 0.36 21,800 29,400 0.32 36,100 0.28 0.24 0.2 20 40 60 80 100 120 140 θ Angular Position, θ ( ) 160 Γ Figure 3. Film thickness distribution for ΓΓ1/2 = 0.01 kg/(m·s) and C = 0.65 at t = 2.79 s. 3.7. Validation The CFD model was compared with the experimental data by Hou et al. [24] for a 25.4-mm tube with 10 mm intertube spacing under Ref = 574, as shown in Figure 4a. The CFD model predicts well the θ –◦ film thickness distribution, with a maximum error of 16.9% in the range of θ = θ 30◦ –160 – and 25.6% for θ Γ θ =θ20◦ near the impingement. In Figure 4b, the time-averaged film thickness for ΓΓ1/2 = 0.01 kg/(m·s) and C = 0.65 is equated with the Nusselt film thickness and found to be in good agreement, with a θ maximum error of 13.3% at θ θ = 20◦ . 0.75 0.36 0.65 δ δ (mm) Film Thickness, δ δ (mm) Film Thickness, CFD Hou et al. 0.55 0.45 0.35 0.25 0.28 0.24 0.2 20 (a) Nusslet CFD 0.32 40 60 80 100 120 Angular Position, θ ( ) θ 140 160 20 (b) 40 60 80 100 120 140 160 Angular Position, θ ( ) θ Figure 4. CFD model comparison of (a) experimental data by Hou et al. [24] and (b) Nusselt solution as in Equation (7). 4. Results and Discussion 4.1. Transient Analysis Figure 5 shows the volume fraction contours of the LiBr solution (blue region) with a liquid load of Γ1/2 = 0.01 kg/(m·s) and concentration of C = 0.65. At t = 0.4 s, the liquid falls freely with gravity as Energies 2020, 13, 307 7 of 15 Γ the driving force and is opposed by the surface tension, which attempts to form a droplet shape, as shown at t = 0.8 s. At t = 1.2 s, the liquid hits the top of the tube and starts spreading around the tube wall. The impact region is known as the stagnation point; as the liquid spreads it enters stagnation and then the developing region. At t = 1.6 s, the liquid continues to spread, forming a thin layer of LiBr solution. At this stage, the liquid film is in the developing region, as the velocity profiles within the film are not fully developed. From t = 2 s to 2.4 s, the tube surface is mostly covered by liquid film and the flow is fully developed. At t = 2.8 s, the tube is entirely covered with a thin layer of LiBr solution and the liquid starts gathering underneath the tube. The adhesion and surface tension forces hold the liquid at the bottom and prevent it from falling. The liquid keeps on gathering until the inertial force dominates and the liquid separates from the tube surface and falls freely. t = 0.4 s t = 0.8 s t = 1.2 s t = 1.6 s t = 2.0 s t = 2.4 s t = 2.8 s t = 3.2 s Figure 5. Volume fraction contours for Γ1/2 = 0.01 kg/(m·s) and C = 0.65 w.r.t. time. Γ Figure 6a,b shows the film thicknesses δins at an angular position ranging from 5◦ to 175◦ for different time instants under the operating conditions of Γ1/2 = 0.01 kg/(m·s) and C = 0.45/0.65. It is clear δ that the film thickness pattern is not fixed, but changes withΓrespect to time. Therefore, for comparison and analysis, the time-averaged film thickness δ is considered in this work. The time-averaged film thickness profile is represented by the solid line in Figure 6a,b. The jet mode is usually considered δ steady with weak perturbations. Figure 7 displays the transient behavior of film thickness at θ = 90◦ for a liquid load corresponding to the jet mode with C = 0.45 and 0.65. For C = 0.45, the instantaneous film θ thickness varies over time after the tube is completely wet and becomes stable with an average thickness of 0.245 mm. However, there remain some fluctuations, with a maximum deviation of 14.1%. At C = 0.65, the density and viscosity are 24.8% and 384% higher than at C = 0.45. The dense, viscous fluid rapidly becomes stable, with very low fluctuations of +0.5%. The perturbations positively affect the heat and mass transfer, and are present for a low-density and viscous fluid such as C = 0.45. Energies 2020, 13, 307 δ tavg t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 0.8 0.6 0.4 0.2 0 0 (a) 30 60 90 120θ 150 1.8 δ 1 Ins. Film Thickness, δ ins (mm) Ins. Film Thickness, δ ins (mm) 8 of 15 t2 t3 1.2 t4 t5 0.9 t6 t7 0.6 t8 t9 0.3 t10 0 180 0 (b) Angular Position, θ (°) t1 1.5 30 60 90 120 Angular Position, θ (°) θ tavg 150 180 Figure 6. Film thickness at different times after the tube is completely covered and time-averaged film Γ thickness profiles for Γ1/2 = 0.01 kg/(m·s) and (a) C = 0.45, and (b) C = 0.65. Γ Ins. Film Thickness,δ δ ins (mm) 0.4 0.35 C = 0.45 C = 0.65 0.3 0.25 0.2 0.15 0 0.1 0.2 0.3 0.4 0.5 Time from Tube is Completely Covered, t (s) Figure 7. Instantaneous film thickness atθθ = 90◦ for Γ Γ1/2 = 0.05 kg/(m·s). θ Γ 4.2. Film Thickness and Residence Time The effect of liquid load on film thickness distribution is shown in Figure 8a,b. The film thickness decreases as the fluid moves from top to bottom, until it reaches to a minimum value, which is around θ = 90◦ . After that, the film thickness again begins to rise, sharply near the bottom of the tube, which is θ ◦ due to the accumulation of liquid. As the liquid load is increased from 0.01 kg/(m·s) to after θ = 150 θ θ 0.03 kg/(m·s), the film distribution shifts upwards for both C = 0.45 and C = 0.65. However, when the θ liquid load is raised from 0.03 kg/(m·s) to 0.05 kg/(m·s), the growth in film thickness is smaller than that of 0.01 kg/(m·s) to 0.03 kg/(m·s). This is because it follows the same trend as in Equation (7), which implies that the film thickness is directly proportional to the cubic root of liquid load, keeping the other parameters constant: p (16) α 3√3𝛤Γ1/2 . 𝛿δ 𝛼 1/2 3 With the same liquid load value, the film thickness 𝛿𝛼 √𝛤1/2 of C = 0.65 is higher than that of C = 0.45 because the density and viscosity vary with concentration. Figure 9a represents the time required for the tube surface to be completely covered by the liquid film. The higher film thickness for C = 0.65 with the same mass flux indicates a smaller fluid velocity field, which means a dense, viscous fluid needs more time for tube surface coverage. This indicates that the residence time for a higher concentration is more than that for a lower concentration, which means that more heat transfer is possible due to the higher residence time. On the other hand, the higher thickness results in increased thermal resistance, which is undesirable for heat transfer. At Γ1/2 = 0.01 kg/(m·s), the time needed for C = 0.65 and C = 0.45 is 2.79 s and 1.52 s, respectively. With the increase in liquid load, the fluid moves faster with a higher velocity field. Thus, the coverage time decreases to 0.88 s and 0.5 s for C = 0.65 and C = 0.45, Γ Γ Energies 2020, 13, 307 9 of 15 respectively. The film thickness in the impingement and detachment zones, as shown in Figure 6a, θ – ◦ ◦ is much greater the average film thickness δavg – thickness δ than the in-between θ zone, – i.e., θ = 20 –160 .θTherefore, ◦ is evaluated in Figure 9b, with liquid load and concentration variation. thicknessfor δ θ = 20–160 and isθpresented – The average thickness for C = 0.65 is 43.3% at 0.01 kg/(m·s) and 44% at 0.05 kg/(m·s), higher than that for C = 0.45. The higher average film thickness is due to a change in thermophysical properties with concentration. 0.45 0.65 0.4 Film Thickness, δ (mm) δ Film Thickness, δ (mm) δ 0.01 kg/(m-s) 0.03 kg/(m-s) 0.35 0.05 kg/(m-s) 0.3 0.03 kg/(m-s) 0.5 0.05 kg/(m-s) 0.4 0.35 0.2 0.15 0.3 0.25 0.1 (a) 0.55 0.45 0.25 20 0.01 kg/(m-s) 0.6 40 θ120 Angular Position, θ ( ) 60 80 100 0.2 140 160 20 40 (b) 60 80 100 θ120 140 160 Angular Position, θ ( ) δ δavg (mm) Average Film Thickness, Figure 8. Film thickness distribution w.r.t. liquid load for (a) C = 0.45 and (b) C = 0.65. Time to cover tube, t (s) 3 C=0.45 2.5 C=0.65 2 1.5 1 0.5 0 0.01 0.02 0.03 0.04 0.05 0.45 C=0.45 0.4 C=0.65 0.35 0.3 0.25 0.2 0.15 0.01 0.02 0.03 0.04 0.05 Liquid Load, Γ1/2 (kg/(m∙s)) Liquid Load, Γ1/2 (kg/(m∙s)) (b) Γ ∙ Γ ∙ Figure 9. (a) Time required for the complete wetting of tube and (b) average film thickness for different liquid loads and concentrations. (a) 4.3. Velocity Distribution Figure 10a,b show velocity vectors at 0.05 kg/(m·s) for low and high concentrations. As the fluid falls and hits the tube, recirculation within the film begins. These recirculation formations enhance the fluid mixing and convection effects in the impingement and developing zones. Higher recirculation results in improved thermal performance [42]. The velocity vectors for C = 0.45 exhibit stronger recirculation than is seen for C = 0.65. The higher density and viscosity result in higher inertial force, which opposes recirculation. The film thickness has been normalized to dimensionless film thickness, as shown in Figure 10c, for comparison purposes. Figure 11a,b shows the normal Vn and tangential velocity Vt components for C = 0.45 and angular position ranging from 5◦ to 60◦ . The normal velocity profile at 5◦ displays fluid flow in the positive, normal direction. Moreover, it is clear from the tangential velocity profile that the fluid near the wall is moving downward, while some part is flowing in the opposite direction. This behavior shows recirculation phenomena in the film. The maximum normal and tangential velocity opposite to the mainstream flow at 5◦ were found to be 0.024 m/s and 0.06 m/s, respectively. The recirculation can be seen up to 10◦ , which is evident in the normal and tangential velocity components. From 20◦ , the flow in the film starts developing, the tangential component keeps rising because of the gravitational and inertial force, and the tangential component profile gets fully developed from 40◦ . The maximum tangential velocity at 60◦ is found Energies 2020, 13, 307 10 of 15 – to be 0.189 m/s. The normal velocity magnitude decreases as the flow develops, which is apparent from the normal velocity profiles for 20–60◦ . Figure 11c,d shows the normal and tangential velocity components for C = 0.65. In this case, the velocity profiles reflect the weak recirculation. The maximum tangential velocity opposite to the flow is 0.0021 m/s, which is 96.5% lower than for C = 0.45. Here, the flow enters the developing region from 10◦ and the tangential velocity component increases. The flow is developed at 40◦ , with tangential velocity increasing 30◦ onwards. The maximum tangential velocity at 60◦ is 0.108 m/s, which is 42.8% less than for C = 0.45. Hence, the less dense and viscous flow (C = 0.45) possesses strong recirculation and higher velocity, which may enhance the heat and mass transfer process when matched with a fluid of high density and viscosity (C = 0.65). Recirculation region δ (b) (a) n η = n/δ (c) Figure 10. Velocity vector for Γ1/2 = 0.05 kg/(m·s) and (a) C = 0.45 (b) C = 0.65 and (c) dimensionless film thickness η schematic. Γ thickness η schematic 0.2 0.03 5 deg 10 deg 20 deg 30 deg 40 deg 50 deg 60 deg 0 -0.01 -0.02 5 deg 10 deg 20 deg 30 deg 40 deg 50 deg 60 deg -0.03 -0.04 -0.05 -0.06 0 Normal velocity, Vn (m/s) (a) 0.2 0.15 0.1 0.05 0 -0.05 -0.1 0.4 0.6 0.8 0 1 (b) Dimensionless film thickness, η 0.12 0 0.1 -0.005 -0.01 -0.015 5 deg 10 deg 20 deg 30 deg 40 deg 50 deg 60 deg -0.02 -0.025 -0.035 0 0.2 0.2 0.4 0.6 0.8 1 0.8 1 Dimensionless film thickness, η 0.005 -0.03 (c) Tangential velocity, Vt (m/s) 0.01 Tangential velocity, Vt (m/s) Normal velocity, Vn (m/s) 0.02 5 deg 10 deg 20 deg 30 deg 40 deg 50 deg 60 deg 0.08 0.06 0.04 0.02 0 -0.02 0.4 0.6 Dimensionless film thickness, η 0.8 1 0 (d) 0.2 0.4 0.6 Dimensionless film thickness, η Figure 11. Velocity profile for ΓΓ1/2 = 0.05 kg/(m·s) and C = 0.45 (a) normal component and (b) tangential component; and for Γ1/2Γ= 0.05 kg/(m·s) and C = 0.65 (c) normal component and (d) tangential component. 𝑅𝑡𝑜𝑡 = 𝑅𝑖𝑛𝑡 +𝑅𝑤 +𝑅𝑐𝑜𝑛𝑑 + 𝑅𝑒𝑣𝑎𝑝 Energies 2020, 13, 307 11 of 15 4.4. Thermal Resistance The thermal performance of falling film evaporators depends on the overall heat transfer coefficient or total thermal resistance, which is the reciprocal of the first. The total thermal resistance Rtot is the sum of all thermal resistances in the falling film network, which can be mathematically described as: Rtot = Rint + Rw + Rcond + Revap , (17) where Rint is the thermal resistance of internal flow, Rw is the tube wall thermal resistance, Rcond is the conduction thermal resistance of liquid film, and Revap is the thermal resistance by evaporation. From a hydrodynamics study, Rcond can be computed by the following equation: Rcond = δ . k (18) Conduction Thermal Resistance Rcond, (x10 -3 m2·K/W) Figure 12 shows the average conduction thermal resistance variation with liquid load and concentration for θ = 20◦ –160◦ . The average thermal resistance at 0.01 kg/(m·s) and C = 0.45 is 0.361 × 10−3 m2 ·K/W, whereas the thermal resistance is 0.621 × 10−3 m2 ·K/W for C = 0.65, which is 72% higher. Furthermore, the thermal resistance increases with the liquid load for both concentrations. The main factors are the film thickness, which is higher for C = 0.65, and the thermal conductivity, which decreases with an increasing LiBr concentration. At 0.05 kg/(m·s), the thermal resistance for C = 0.65 is 72.9% higher than that of C = 0.45. It can be deduced that, with increasing liquid load and concentration, the heat transfer rate may decrease because of the higher thermal resistance. 1.25 C=0.45 C=0.65 1 0.75 0.5 0.25 0.01 0.02 0.03 0.04 0.05 Liquid Load, Γ1/2 (kg/(m∙s)) Figure 12. Average conduction thermal resistance variation w.r.t. liquid load. 5. Conclusions The CFD model of an aqueous LiBr solution falling over a horizontal tube was developed in this work. The effects of LiBr concentration and liquid load on falling film distribution were analyzed. It was found that jet mode, corresponding to 0.05 kg/(m·s), was more stable at a higher concentration, i.e., C = 0.65, with small fluctuations of +0.5%. In contrast, low concentration (C = 0.45) flow with 0.05 kg/(m·s) had maximum fluctuations of 14.1%. The film thickness increased with concentration and liquid load, and the residence time increased with higher concentration and lower liquid loads. The average film thickness and residence time at 0.05 kg/(m·s) and C = 0.65 were found to be 44% and 76% higher than for C = 0.45, respectively. Weak recirculation was observed for higher concentration; the recirculation exists until 5◦ for C = 0.65 and 10◦ for C = 0.45. In addition, the maximum tangential velocity when the flow was developed was 0.108 m/s for C = 0.65, vs. 0.189 m/s for C = 0.45. The flow θ Γ Energies 2020, 13, 307 12 of 15 was fully developed at θ = 40◦ for both concentrations. The thermal resistance analysis demonstrated that the thermal resistance rises with the concentration and liquid load. For Γ1/2 = 0.05 kg/(m·s), the thermal resistance for C = 0.65 was 72.9% higher than that of C = 0.45. It may be deduced from this hydrodynamics study that the heat transfer performance is expected to be better at lower concentrations because of the lower film thickness, more recirculation, higher velocity, and lower thermal resistance. In practical applications, such as refrigeration and desalination, this implies that regions in falling film evaporators with low concentration would have better thermal performance as compared to high concentration regions. Author Contributions: Conceptualization, A.M. and M.K.; methodology, F.T.; software, F.T.; validation, F.T.; formal analysis, M.K.; investigation, F.T. and A.M.; resources, A.M. and M.K.; data curation, A.M.; writing—original draft preparation, F.T.; writing—review and editing, A.M. and M.K.; visualization, F.T.; supervision, A.M. and M.K.; project administration, A.M. and M.K.; funding acquisition, F.T. All authors have read and agreed to the published version of the manuscript. Funding: The publication of this article was funded by the Qatar National Library (QNL). Acknowledgments: The authors would like to acknowledge the resources and support provided by the College of Science and Engineering (CSE) and the Qatar Environment and Energy Research Institute (QEERI) of Hamad Bin Khalifa University (HBKU). The publication of this article was funded by the Qatar National Library. Conflicts of Interest: The authors declare no conflict of interest. Nomenclature C CFD CSF D F g Ga hi k K LiBr PRESTO R Ref s SIMPLE t V VOF Greek letters α σ µ Γ1/2 ρ δ θ θw η LiBr mass concentration (%) computational fluid dynamics continuum surface force tube diameter (m) force (N) acceleration due to gravity (m/s2 ) modified Galileo number impingement height (m) thermal conductivity (W·m/K) curvature lithium bromide pressure staggering option thermal resistance (m2 ·K/W) film Reynolds number intertube spacing (m) semi-implicit method for pressure linked equations time (s) velocity (m/s) volume of fluid void fraction surface tension (N/m) viscosity (Pa·s) liquid load, on one side of tube (kg/(m·s)) density (kg/m3 ) film thickness (m) angular position (◦ ) wall adhesion contact angle (◦ ) dimensionless film thickness Energies 2020, 13, 307 Subscripts avg cond evap g int ins l n t tot w vol 13 of 15 time averaged from θ = 20–160◦ conduction evaporation gaseous phase internal instantaneous liquid phase normal tangential total wall volumetric References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. Danilova, G.N.; Burkin, V.G.; Dyundin, V.A. Heat transfer in spray-type refrigerator evaporators. Heat Transf. Res. 1976, 8, 105–113. Lorenz, J.J.; Yung, D. Film Breakdown and Bundle-Depth Effects in Horizontal-Tube, Falling-Film Evaporators. J. Heat Transf. 1982, 104, 569–571. [CrossRef] Gonzalez, G.J.M.; Jabardo, J.M.S.; Stoecker, W.F. Falling Film Ammonia Evaporators; University of Illinois: Champaign, IL, USA, 1992. Zeng, X.; Chyu, M.-C.; Ayub, D.Z.H. Evaporation heat transfer performance of nozzle-sprayed ammonia on a horizontal tube. ASHRAE Trans. 1995, 101, 136–149. Hoffmann, L.; Greiter, I.; Wagner, A.; Weiss, V.; Alefeld, G. Experimental investigation of heat transfer in a horizontal tube falling film absorber with aqueous solutions of LiBr with and without surfactants. Int. J. Refrig. 1996, 19, 331–341. [CrossRef] Kwon, K.; Jeong, S. Effect of vapor flow on the falling-film heat and mass transfer of the ammonia/water absorber. Int. J. Refrig. 2004, 27, 955–964. [CrossRef] Ribatski, G.; Jacobi, A.M. Falling-film evaporation on horizontal tubes—A critical review. Int. J. Refrig. 2005, 28, 635–653. [CrossRef] Fernández-Seara, J.; Pardiñas, Á.Á. Refrigerant falling film evaporation review: Description, fluid dynamics and heat transfer. Appl. Therm. Eng. 2014, 64, 155–171. [CrossRef] Tahir, F.; Mabrouk, A.; Koc, M. Review on CFD analysis of horizontal falling film evaporators in multi effect desalination plants. Desalin. Water Treat. 2019, 166, 296–320. [CrossRef] Tahir, F.; Mabrouk, A.; Koc, M. CFD Analysis of Falling Film Wettability in MED Desalination plants. In Proceedings of the Qatar Foundation Annual Research Conference Proceedings, Ar-Rayyan, Qatar, 19–20 March 2018; Volume 2018, p. EEPD650. Tahir, F.; Mabrouk, A.A.; Koc, M. CFD analysis of spray nozzle arrangements for multi effect desalination evaporator. In Proceedings of the 3rd Thermal Fluids Engineering Conference, Ft. Lauderdale, FL, USA, 4–7 March 2018; pp. 935–941. Mabrouk, A.A.; Bourouni, K.; Abdulrahim, H.K.; Darwish, M.; Sharif, A.O. Impacts of tube bundle arrangement and feed flow pattern on the scale formation in large capacity MED desalination plants. Desalination 2015, 357, 275–285. [CrossRef] Mitrovic, J. Flow structures of a liquid film falling on horizontal tubes. Chem. Eng. Technol. 2005, 28, 684–694. [CrossRef] Hu, X.; Jacobi, A.M. The Intertube Falling Film: Part 1—Flow Characteristics, Mode Transitions, and Hysteresis. J. Heat Transf. 1996, 118, 616–625. [CrossRef] Hu, X.; Jacobi, A.M. Flow characteristics of liquid droplets and jets falling between horizontal circular tubes. In Experimental Heat Transfer, FLuid Mechanics and Thermodynamics; PISA, Edizioni ETS: Paris, France, 1997; pp. 1295–1302. Roques, J.F.; Dupont, V.; Thome, J.R. Falling Film Transitions on Plain and Enhanced Tubes. J. Heat Transf. 2002, 124, 491–499. [CrossRef] Energies 2020, 13, 307 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 14 of 15 Roques, J.-F.; Thome, J.R. Falling Film Transitions Between Droplet, Column, and Sheet Flow Modes on a Vertical Array of Horizontal 19 FPI and 40 FPI Low-Finned Tubes. Heat Transf. Eng. 2003, 24, 40–45. [CrossRef] Nusselt, W. Die oberflachenkondensation des wasserdamphes. VDI-Zs 1916, 60, 541–546. Ji, G.; Wu, J.; Chen, Y.; Ji, G. Asymmetric distribution of falling film solution flowing on hydrophilic horizontal round tube. Int. J. Refrig. 2017, 78, 83–92. [CrossRef] Hosseinnia, S.M.; Naghashzadegan, M.; Kouhikamali, R. CFD simulation of water vapor absorption in laminar falling film solution of water-LiBr—Drop and jet modes. Appl. Therm. Eng. 2017, 115, 860–873. [CrossRef] Ding, H.; Xie, P.; Ingham, D.; Ma, L.; Pourkashanian, M. Flow behaviour of drop and jet modes of a laminar falling film on horizontal tubes. Int. J. Heat Mass Transf. 2018, 124, 929–942. [CrossRef] Rogers, J.T.; Goindi, S.S. Experimental laminar falling film heat transfer coefficients on a large diameter horizontal tube. Can. J. Chem. Eng. 1989, 67, 560–568. [CrossRef] Gstoehl, D.; Roques, J.F.; Crisnel, P.; Thome, J.R. Measurement of Falling Film Thickness Around a Horizontal Tube Using a Laser Measurement Technique. Heat Transf. Eng. 2004, 25, 28–34. [CrossRef] Hou, H.; Bi, Q.; Ma, H.; Wu, G. Distribution characteristics of falling film thickness around a horizontal tube. Desalination 2012, 285, 393–398. [CrossRef] Chen, X.; Shen, S.; Wang, Y.; Chen, J.; Zhang, J. Measurement on falling film thickness distribution around horizontal tube with laser-induced fluorescence technology. Int. J. Heat Mass Transf. 2015, 89, 707–713. [CrossRef] Lin, S.; Liu, X.; Li, X. The spatial distribution of liquid film thickness outside the horizontal falling film tube. Int. J. Heat Mass Transf. 2019, 143, 118577. [CrossRef] Wang, J.; Chen, X.; Lu, T.; Chen, X.; Shen, S.; Liu, B. Three-dimensional film thickness distribution of horizontal tube falling film with column flow. Appl. Therm. Eng. 2019, 154, 140–149. [CrossRef] Tahir, F.; Baloch, A.A.B.; Ali, H. Resilience of Desalination Plants for Sustainable Water Supply in Middle East. In Sustainability Perspectives: Science, Policy and Practice, Strategies for Sustainability; Khaiter, P.A., Erechtchoukova, M.G., Eds.; Springer Nature Switzerland AG: Basel, Switzerland, 2020; pp. 303–329. ISBN 978-3-030-19550-2. Tahir, F.; Mabrouk, A.; Koc, M. Energy and exergy analysis of thermal vapor compression (TVC) and absorption vapor compression (AVC) for multi effect desalination (MED) plants. In Proceedings of the 8th Global Conference on Global Warming (GCGW), Doha, Qatar, 22–25 April 2019; p. 25. Mabrouk, A.; Abotaleb, A.; Abdelrehim, H.; Tahir, F.; Koc, M.; Abdelrashid, A.; Nasralla, A. High Performance MED Desalination Plants, Part II: Novel Integration MED-Absorption Vapor Compression. In Proceedings of the IDA 2017 World Congress on Water Reuse and Desalination, Sao Paulo, Brazil, 15–20 October 2017. Narváez-Romo, B.; Chhay, M.; Zavaleta-Aguilar, E.W.; Simões-Moreira, J.R. A critical review of heat and mass transfer correlations for LiBr-H2 O and NH3 -H2 O absorption refrigeration machines using falling liquid film technology. Appl. Therm. Eng. 2017, 123, 1079–1095. [CrossRef] Pátek, J.; Klomfar, J. A computationally effective formulation of the thermodynamic properties of LiBr–H2 O solutions from 273 to 500K over full composition range. Int. J. Refrig. 2006, 29, 566–578. Kaushik, S.C.; Arora, A. Energy and exergy analysis of single effect and series flow double effect water–lithium bromide absorption refrigeration systems. Int. J. Refrig. 2009, 32, 1247–1258. [CrossRef] Khan, S.A.; Tahir, F.; Baloch, A.A.B.; Koc, M. Review of Micro–Nanoscale Surface Coatings Application for Sustaining Dropwise Condensation. Coatings 2019, 9, 117. [CrossRef] Hirt, C.; Nichols, B. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–225. [CrossRef] Baloch, A.; Ali, H.; Tahir, F. Transient Analysis of Air Bubble Rise in Stagnant Water Column Using CFD. In Proceedings of the ICTEA: International Conference on Thermal Engineering, Doha, Qatar, 25–28 February 2018; Volume 2018. Nigim, T.H.; Eaton, J.A. CFD prediction of the flashing processes in a MSF desalination chamber. Desalination 2017, 420, 258–272. [CrossRef] Ozmen-Cagatay, H.; Kocaman, S. Dam-Break Flow in the Presence of Obstacle: Experiment and CFD Simulation. Eng. Appl. Comput. Fluid Mech. 2011, 5, 541–552. [CrossRef] Energies 2020, 13, 307 39. 40. 41. 42. 15 of 15 Lopes, R.J.G.; Quinta-Ferreira, R.M. Assessment of CFD−VOF Method for Trickle-Bed Reactor Modeling in the Catalytic Wet Oxidation of Phenolic Wastewaters. Ind. Eng. Chem. Res. 2010, 49, 2638–2648. [CrossRef] Gupta, R.; Fletcher, D.F.; Haynes, B.S. On the CFD modelling of Taylor flow in microchannels. Chem. Eng. Sci. 2009, 64, 2941–2950. [CrossRef] Brackbill, J.; Kothe, D.; Zemach, C. A continuum method for modeling surface tension. J. Comput. Phys. 1992, 100, 335–354. [CrossRef] Yang, L.; Liu, Y.; Yang, Y.; Shen, S. Microscopic mechanisms of heat transfer in horizontal-tube falling film evaporation. Desalination 2016, 394, 64–71. [CrossRef] © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).