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

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.

1. Introduction

Falling film exchangers are commonly employed in ammonia (NH₃)- 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:

Droplet to droplet-jet mode:

$$R{e}_f=0.0417G{a}^{0.3278}.$$

(1)

Droplet-jet to jet mode:

$$R{e}_f=0.0683G{a}^{0.3204}.$$

(2)

Jet to jet-sheet mode:

$$R{e}_f=0.8553G{a}^{0.2483}.$$

(3)

Jet-sheet to sheet mode:

$$R{e}_f=1.068G{a}^{0.2563}.$$

(4)

Re_f and Ga are defined as:

$$R{e}_f=\frac{4{\Gamma }_{1/2}}{\mu }$$

(5)

$$Ga=\frac{\rho {\sigma }^3}{{\mu }^{4}g},$$

(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,20,21]:

$$\delta =\sqrt[3]{\frac{3{\mu }_l{\Gamma }_{1/2}}{{\rho }_l\left({\rho }_l-{\rho }_g\right)gsin\theta }}.$$

(7)

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 Re_f = 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, 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.

2. 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. Thermophysical properties of LiBr solution and water vapors.

Fluid		Density (kg/m3)	Viscosity (Pa-s)	Surface Tension (N/m)
LiBr solution	C = 0.45	1453	0.002274	0.08
	C = 0.65	1814	0.01102	0.085
Water Vapors		0.0304	0.0000101

.

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:

$$\frac{\partial \rho }{\partial t}+ \nabla ·\left(\rho \stackrel{\rightharpoonup }{V}\right)=0$$

(8)

$$\frac{\partial \left(\rho \stackrel{\rightharpoonup }{V}\right)}{\partial t}+ \nabla ·\left(\rho \stackrel{\rightharpoonup }{V}\stackrel{\rightharpoonup }{V}\right)=-\nabla P+\nabla ·\left(\stackrel{=}{\tau }\right)+\rho \stackrel{\rightharpoonup }{g}+{\stackrel{\rightharpoonup }{F}}_{vol},$$

(9)

where V is the velocity, t is time, P is the pressure, τ is the tensor, and F_vol 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,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.

$$\frac{\partial \left({\alpha }_i{\rho }_i\right)}{\partial t}+ \nabla ·\left({\alpha }_i{\rho }_i\stackrel{\rightharpoonup }{V_i}\right)=0 \left(i=g or l\right)$$

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

$$\rho ={\alpha }_l{\rho }_l+(1-{\alpha }_l){\rho }_g$$

(11)

$$\mu ={\alpha }_l{\mu }_l+(1-{\alpha }_l){\mu }_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:

$${\alpha }_g=1-{\alpha }_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 F_vol. The CSF model improves the characterization of the liquid‒gas interface by accommodating surface tension σ effects:

$${\stackrel{\rightharpoonup }{F}}_{vol}=\sigma K\nabla {\alpha }_l[ \frac{\rho }{\frac{1}{2}\left({\rho }_l+{\rho }_g\right)}].$$

(14)

The surface curvature K is computed by gradient of surface unit normal vector $\widehat{n}$:

$$K=-\nabla ·\widehat{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 Re_f and flow modes for each case are listed in

Table 2. Films Reynolds number Re_f and flow modes for all cases.

C	Γ1/2 (kg/(m·s))	Ref	Ga	Transitional Ref [16]		Mode
				Droplet to Droplet-Jet	Droplet-Jet to Jet
0.45	0.01	17.6	2.84 × 109	52.3	73	Droplet
	0.03	52.8				Droplet-Jet
	0.05	87.9				Jet
0.65	0.01	3.6	7.7 × 106	7.5	11	Droplet
	0.03	10.9				Droplet-Jet
	0.05	18.2				Jet

. The gravity is acting in the negative y-direction with a value of 9.81 m/s².

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.

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

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

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

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

$$\delta \alpha \sqrt[3]{{\Gamma }_{1/2}}.$$

(16)

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

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

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 R_tot is the sum of all thermal resistances in the falling film network, which can be mathematically described as:

$$R_{tot}=R_{int}+R_w+R_{cond}+R_{evap},$$

(17)

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:

$$R_{cond}= \frac{\delta }{k}.$$

(18)

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⁻³ m²·K/W, whereas the thermal resistance is 0.621 × 10⁻³ m²·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.

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