CFD Analysis of Falling Film Hydrodynamics for a Lithium Bromide (LiBr) Solution over a Horizontal Tube

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.


Introduction
Falling film exchangers are commonly employed in ammonia (NH 3 )-and lithium bromide (LiBr)-based absorption refrigeration systems [1][2][3][4][5][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][10][11][12][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 Re f and modified Galileo number Ga for flow mode recognition and transition [14][15][16][17]. The correlations developed by Roques et al. [16] are as follows:

Problem Statement
In order to characterize the film distribution, one 25.4-mm tube was selected with an impingement height h i 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.
Energies 2020, 13, x FOR PEER REVIEW 3 of 15 minimum thickness to lie in the valley area, 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 this study, a 2D CFD model for an aqueous LiBr solution falling over a horizontal tube was modeled, which is used in the refrigeration and multi-effect desalination industries [28][29][30][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.

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.

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].

Conservation Equations
The mass and momentum balances can be found by continuity and Navier-Stokes equations, as shown below: where V is the velocity, t is time, P is the pressure, τ is the tensor, and F vol is the volumetric force.

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][37][38][39][40]. They modified the continuity equation as in Equation (8) by introducing void fraction α to solve the mass balance for both phase separately.
The subscripts g and l represent gaseous and liquid phases, respectively. The density and viscosity are evaluated on the basis of the void fraction.
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:

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 F vol . The CSF model improves the characterization of the liquid-gas interface by accommodating surface tension σ effects: Energies 2020, 13, 307

of 15
The surface curvature K is computed by gradient of surface unit normal vectorn:

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 Re f 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/s 2 .

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.
Energies 2020, 13, x FOR PEER REVIEW 5 of 15 The surface curvature K is computed by gradient of surface unit normal vector ̂:

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/s 2 .

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 georeconstruct 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.

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, a mesh with 29,400 elements was selected for the numerical calculations. In addition, a time step of 10 µs was selected after the time dependency test.

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, a mesh with 29,400 elements was selected for the numerical calculations. In addition, a time step of 10 μs was selected after the time dependency test.

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°.

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 Re f = 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 • .

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, a mesh with 29,400 elements was selected for the numerical calculations. In addition, a time step of 10 μs was selected after the time dependency test.

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°.   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.  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 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.  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 timeaveraged film thickness profile is represented by the solid line in Figure 6a,b. The jet mode is usually considered steady with weak perturbations.  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.

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 after θ = 150° due to the accumulation of liquid. As the liquid load is increased from 0.01 kg/(m·s) to 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: With the same liquid load value, the film thickness 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 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.

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 after θ = 150° due to the accumulation of liquid. As the liquid load is increased from 0.01 kg/(m·s) to 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: With the same liquid load value, the film thickness 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

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 after θ = 150 • due to the accumulation of liquid. As the liquid load is increased from 0.01 kg/(m·s) to 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: With the same liquid load value, the film thickness 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 respectively. The film thickness in the impingement and detachment zones, as shown in Figure 6a, is much greater than the in-between zone, i.e., θ = 20 • -160 • . Therefore, the average film thickness δ avg is evaluated for θ = 20-160 • and is presented in Figure 9b, with liquid load and concentration variation. 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.
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, respectively. The film thickness in the impingement and detachment zones, as shown in Figure 6a, is much greater than the in-between zone, i.e., θ = 20°-160°. Therefore, the average film thickness δavg is evaluated for θ = 20-160° and is presented in Figure 9b, with liquid load and concentration variation. 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.   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 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, respectively. The film thickness in the impingement and detachment zones, as shown in Figure 6a, is much greater than the in-between zone, i.e., θ = 20°-160°. Therefore, the average film thickness δavg is evaluated for θ = 20-160° and is presented in Figure 9b, with liquid load and concentration variation. 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.   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 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 V n and tangential velocity V t 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 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).

Velocity Distribution
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 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).

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: 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.

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 R tot is the sum of all thermal resistances in the falling film network, which can be mathematically described as: where R int is the thermal resistance of internal flow, R w is the tube wall thermal resistance, R cond is the conduction thermal resistance of liquid film, and R evap is the thermal resistance by evaporation. From a hydrodynamics study, R cond can be computed by the following equation: 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 m 2 ·K/W, whereas the thermal resistance is 0.621 × 10 −3 m 2 ·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.

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 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

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 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. 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.