Three-Dimensional Mathematical Model of the Liquid Film Downﬂow on a Vertical Surface

: Film downﬂow from captured liquid without wave formation and its destruction is one of the most important aspects in the development of separation equipment. Consequently, it is necessary to create well-organized liquid draining in areas of captured liquid. Thus, the proposed 3D mathematical model of ﬁlm downﬂow allows for the determination of the hydrodynamic parameters of the liquid ﬁlm ﬂow and the interfacial surface. As a result, it was discovered that the interfacial surface depends on the proposed dimensionless criterion, which includes internal friction stress, channel length, and ﬂuid density. Additionally, equations for determining the averaged ﬁlm thickness, the averaged velocity vectors over the ﬁlm thickness, the longitudinal and vertical velocity components, and the initial angle of streamline deviation from the vertical axis were analytically obtained.


Introduction
Liquid film flow is one of the most common phenomena in the chemical and oil/gas industries during the processes of a heat/mass transfer and a separation. It is worth noting that film downflow allows for increasing the intensity and efficiency of the abovementioned processes. The value of the thermal conductivity and heat transfer between the coolants in heat exchangers, the speed of mass transfer between the phases in mass transfer equipment, and the efficiency in draining liquids in separation equipment [1,2] increase using film flows. Therefore, the high efficiency of these processes should be achieved by studying the fluid downflow and its hydrodynamic and heat exchange parameters, the above parameters' dependence on flow regimes and surface topology, as well as the forces that contribute to films forming and destruction.
The use of the general method, which is not appropriate for obtaining reliable data, cannot be used to determine the hydrodynamic parameters of a liquid film and its thickness. Therefore, the numerical identification allows for the determination of new mathematical models for the processes, including liquid film flow. The authors of [15] described fluid flow motion on a vertical surface with evaporation into gas flow by solving the equation system. This system includes equations of continuity, mass transfer, mass balance, and volume fraction of liquid evaporation into the gas phase and Dalton's law. According to the dependencies described in [15], the film thickness and its flow rate decrease along with the surface height as a result of the water evaporation, and the film thickness and speed decrease during cross-interaction are comparable with the backflow and even slightly exceed it. The study [16] includes a considerable amount of new experimental data on a falling film with fully developed waves. The relationship between the wave amplitude and the full range of fluid properties (the Kapitza number (Ka) in a range of 2-3890) and flow rates (the Weber ration (We) in a range of 0.08-58.60) is determined. The obtained data at the constant Kapitza number indicate that the amplitude of waves gradually increases with an increasing flow velocity. Additionally, the authors declare two modes: low-Kapitza fluids (Ka < 200) and high-Kapitza fluids (Ka > 200). These two liquid types have different wave amplitudes, specifically for the abovementioned parameter of liquid film flow for the first type of saturate around three times the flat film thickness and the second around 10 times. It is important to note that the article also discusses the physical meaning for the existence of two Kapitza number ranges for wave amplitudes. The Kapitza number can be expressed as a ratio of up-flow forming capillary forces to viscous resistance forces. The critical value of the ratio is needed to have locally negative flows that can support large-amplitude waves. The data also indicate the dependence of film roughness (normalized standard deviation) on the inverse of the Weber number; herewith, the dependence of film roughness on a Kapitza number is weaker. The authors of [16] have developed a mathematical model for the two-dimensional case, which gives a good agreement with the experimental data. The study [17] compares using three approaches to the description of the wavy film flow on the vertical surfaces, namely the Navier-Stokes equations in their full statement and two integral methods: Shkadov's and the regularized fundamental model. The usage range of each approach was determined depending on the relationship between the dimensional neutral disturbance wavelength and the wavelength, the Reynolds number, and the Kapitza number. The study [18] investigates the falling film flow with periodic waves on a vertical surface. The Reynolds number of the investigated flow is in a range of 5-10. As a result, authors calculate the film waveforms, the time of the regular wave mode forming, and amplitudes of periodic disturbances.
As mentioned above, the liquid film flow is considered not only on vertical surfaces but also inside channels of various configurations [19][20][21][22][23][24][25]. The paper [19] describes the influence of bubble approach velocities on the liquid film drainage process under constant conditions for plane parallel immobile surfaces by the Stefan-Reynolds model. The relations between the bubble approach velocity and the critical thickness of the liquid film were found (the film wetting stability decreases at a high speed). For example, these effects are essential for oil recovery from tar sands. The authors of study [20] carried out the experimental studies of falling liquid films in internally structured tubes with intermediate and steep spiral grooves. Water and propane were used for forming the liquid film. The structured internal surface of the vertical tube increases the stability of the film flow along its length compared to a smooth interior surface of the vertical pipe. The study [21] describes the usage of CFD methods for determining the thickness of a liquid film in a roller crystallizer. The authors have carried out the two-dimensional transient numerical simulation with implementation for describing the two-phase flow (liquid-glycerin) of the VOF model with geometric reconstruction Youngs. Numerical simulations determined the parameters of the liquid film for different speeds of the drum rotation. It was established that the increase of the speed causes an increase in the film thickness at the film's rising side and the minimum film thickness point (top of the rotating drum). In [22], the experimental studies of the gas-liquid flow in the vertical 180 • return bend were carried out to calculate the film thickness by modifying the Weber and Froude criteria. They were chosen because they consider the main active inertia forces and surface tension. The resulting expression for calculating the film thickness is used to determine the liquid film thickness in gas-liquid annular systems with different diameters of pipe in oil and gas production systems. The behavior of the film flow and the deposition characteristics of particles Energies 2020, 13,1938 4 of 15 on the refractory furnace wall were studied in [23] by using a commercial CFD code coupled with required user-defined subroutines. The Eulerian-Lagrangian approach described the multiphase flow, wherein the particle-wall interaction mechanisms were considered by particle capturing entrainment and wall burning sub-models. The available data were selected as input data from the commercial operation of the furnace. The numerical simulations were performed as a transient case. It was defined as a result of a comparison of the initial data with real operational parameters that use the technique described in [23], which can be implemented for predicting the liquid film formation on the furnace walls. The above methodology allows obtaining the film formation in the liquid/liquid systems [24,25].
It is clear from the listed literature that the liquid film flow over surfaces of various configurations can be described only by using the specific approaches with considering the accompanying processes, which occur in the process equipment. One of the most critical aspects in the development of separation equipment is providing the captured liquid film flow mode without wave formation and the creation of its organized removal [26][27][28]. The avoidance of waves on the falling film is necessary for terms of decreasing the possibility of droplet breaking from its top points and further their removal from the separation element. The mathematical models of liquid film flow and particle movement on a curvilinear surface were developed in [26,29] for the case of the strong influence of the interfacial friction forces. The result of this interfacial interaction is a modification of the liquid film velocity vector direction (angle appearance to the gravity force).
It should be emphasized that another method for the decreasing probability of secondary droplet splashing is the usage of the modular dynamic separation devices. Nevertheless, the parameters of the liquid film are critical for the development of the above type of separation devices [30,31].
It is worth mentioning that the above paper does not consider the change in film thickness. Therefore, the main aim of the research is the development of a three-dimensional mathematical model of a liquid film flow considering the change of thickness in both vertical and horizontal directions. Such a model allows for developing draining channels in the separation device and locating them in places with critical film thickness. The determination of the liquid film thickness on length may be useful for predicting the location of filtering and drainage channels in places with critical film thickness. The determination of the liquid film thickness on height can be useful for obtaining the height of the channel to prevent the exceeding of the flow rate of the liquid phase in the separation channels.

Material and Methods
One of the most widely used separation devices is vane-type mist extractors because they offer the right balance among efficiency, operating range, pressure drop requirement, and installed cost. This type of impingement mist extractor consists of a series of baffles, vanes, or plates commonly installed by the parallel bet; thus, the cross-section changes along with the gas-liquid flow movement. The gas-liquid mixture passes through the vane-type mist extractor; herein, particles from gas flow do not change their trajectories due to sufficient momentum to break through the gas streamlines [26]. The fluid film downflow is considered on a vertical surface of the separation channel. It falls under the action of gravity force considering the interaction with a wall and with the main gas-liquid flow. The governing equations for describing incompressible liquid film flow are the system of the Navier-Stokes equations [32] and the continuity equation [33] with three velocities' component V = (u, v, w), m/s: Energies 2020, 13, 1938 5 of 15 where ρ-liquid density, kg/m 3 ; µ-dynamic viscosity of the liquid, Pa·s; p-dynamic pressure of a fluid, Pa; g-gravitational acceleration, m/s 2 . The solution of the above system of equations is carried out using the following simplifications and assumptions: - The liquid is incompressible, and the process of its flow is isothermal (ρ = const, and therefore divV = 0); - The liquid film flow is considered on the vertical surface of the separation channel considering frictional interaction between phases; accordingly, the film does not flow vertically downward, but at some angle to the gravity direction. Therefore, film flow is presented as a two-dimensional case in xand z-directions. These axes are directed downward and along the wall, respectively. The design scheme is shown in Figure 1; where ρ-liquid density, kg/m 3 ; µ-dynamic viscosity of the liquid, Pa·s; p-dynamic pressure of a fluid, Pa; g-gravitational acceleration, m/s 2 .
The solution of the above system of equations is carried out using the following simplifications and assumptions: -The liquid is incompressible, and the process of its flow is isothermal (ρ = const, and therefore divV = 0); -The liquid film flow is considered on the vertical surface of the separation channel considering frictional interaction between phases; accordingly, the film does not flow vertically downward, but at some angle to the gravity direction. Therefore, film flow is presented as a twodimensional case in x-and z-directions. These axes are directed downward and along the wall, respectively. The design scheme is shown in Figure 1;  The next step after the implementation of all the introduced simplifications and assumptions into the system (1) is the evaluation of the components included in the obtained system [34]. The following values U H 2 , U δ 2 , U L 2 are obtained for comparison, U δ 2 has the highest order of them. That is why only ∂ 2 u ∂y 2 stays on the right side of the equation.
The second system equation is considered by analogy to the components order evaluation of the first system equation. The variables are replaced by their values on the liquid film surface and, as a result, the following components 10 −2 , respectively. This data shows that it is necessary to consider only the following equation components since they have a lesser order of smallness: ∂p ∂z and µ ∂ 2 w ∂y 2 . As a result of the implementation all the above simplifications and assumptions for the components, the following system of differential equations are obtained: Considering the first equation allows for finding the equation for the velocity component u. Considering the boundary condition of zero velocity on the vertical wall ( u| y = 0 = 0), the velocity component u(x, y, z) in the general case can be approximately expressed as a third-degree polynomial [35] with respect to the coordinate y: The following boundary conditions allow determining unknown coefficients in expression (3): (1) Due to the gas in the vane-type, separation equipment flows along the z-axis and the vertical velocity component allows for writing the following physical boundary condition on the interfacial surface [36]: τ yx = µ ∂u ∂y y = δ(x,z) = 0; (2) A no-slip boundary condition is applied on the wall surface. Therefore, the velocity components u and v in the first system Equation (2) are equal to zero.
The expression of the velocity component u with obtained coefficients and terms rearrangement has the form: where ν-kinematic viscosity, m 2 /s. The second system Equation (2) allows for finding the velocity component w in the direction of the z-axis. In the case of laminar flow [37] with the absence of wave formation, when the motion occurs under the inertial and gravity forces, the following equation can be written: where ρ mix -gas-liquid mixture density, kg/m 3 . The substitution of the abovementioned expression to the second Equation (2) with the consideration of the dimensionless density ratio k = 1 − ρ mix ρ , consequent integration and transformation allow obtaining the following expression for the velocity component: Energies 2020, 13, 1938 7 of 15 It should also be noted that the (6) expression considers the no-slip boundary condition w| y = 0 = 0 and the interfacial friction condition on the film surface τ yz = µ ∂w The substitution of the obtained expressions for the lengthwise (6) and downward (4) velocities components into the last system Equation (2) (the continuity equation) allows for finding the film thickness equation. The result of the transformations are as follows: The unknown coefficient c included in (4) and (7) is determined by considering the liquid film flow thickness depends on the x, z coordinates ( Figure 1). The auxiliary functions U(x) and W(z) allow obtaining the dependence of δ on x and z: The result of the two functions multiplication (8) is introduced instead of the film thickness in expression (7): The algebraic sum in Equation (9) is equal to zero, so it allows dividing it into two separate equations. The first of it has the following form: According to g νU dU dx = −k g νW dW dz = const, the above equation can be written as a result of the two constants multiplication −g ν c 3 ; one of them (c 3 ) is unknown. After dividing into two differential equations, the thickness is as follows: where δ 0 -liquid film thickness at the origin. The second equation obtained after dividing into two parts of the expression (9) with its negative terms have the form: The next step is (12) integration, after variables' separation and U x U = c 3 accounting: where c 0 (z)-the value of the constant at the trapped liquid film start, included in the vertical velocity component expression. The selection of an elementary area with dimensions dx×dz and momentum equation with respect to the x-axis [38] allows evaluating unknown constant c 3 from formulas (11) and (13): Energies 2020, 13, 1938 8 of 15 Equation (14) is integrated with respect to the y-coordinate and differentiated by the x-coordinate considering the abovementioned expression for the vertical velocity component u and the value of its partial derivative with respect to the y coordinate on the wall ∂u ∂y y = 0 = gδ ν − 3cδ 2 : The above equation allows for the getting of the quadratic equation relative to the parameter c 0 by the substituting of the film thickness δ with the constant c from the expression of the vertical velocity component with the subsequent averaging of the expression over the wall height H and accounting for the 2nd remarkable limit [39]: This equation has two roots. Only one case (3ν << 427 60 g ν δ 4 0 c 3 e 4c 3 z k ) is considered for further model development because it leads to the significant simplification of this equation's analytical expressions roots. Additionally, the equation roots (16) take real values due to the physical content of the problem. Thus, it is necessary to use the condition z ≥ k c 3 ln 0.92 c 3 δ 0 . As a result, the condition c 3 δ 0 ≥ 0.92 must be satisfied with the positive z-coordinate in the flow direction (starting from the initial section). One of the corresponding to the above requirements for roots of the Equation (16) takes the indicated below form after its linearization with respect to the z coordinate, for the case c 3 δ 0 >> 0.92: The final equation for the vertical velocity component u has the following form: The constant c 3 is unknown, as can be seen from the expressions for the vertical velocity component and film thickness of the flowing liquid. The balances' equation for the flow rates [40] allows writing the following equality: where Q u0 = w(x, y)dydx-the liquid flow rate, which is entering the film horizontally at the wall end along which fluid flows (Figure 1).
The above expressions for liquid flow rate can be integrated and substituted into (19). As a result, the transcendental equation makes the following form: where R-the error function-∆Q u = Q u1 − Q u0 is the difference of the liquid flow rates in a vertical direction, and ∆Q w = Q w1 − Q w0 is the difference of the liquid flow rates in the horizontal direction. In this case, the differences of the liquid flow rates are equal in module value according to Equation (20): Energies 2020, 13, 1938 9 of 15 The exclusion of the terms of higher-order smallness from the balances' equation for the flow rates (18) allows for obtaining the following expressions: In this case, the (19) equation has the following solution: This equation emphasizes that the liquid film shape depends on the constant c 3 . Notably, its sign indicates directions of increasing or decreasing film thickness. This fact makes it necessary to introduce a dimensionless film shape criterion Cr = τ ρgL , the limit value of which is 0.23. Thus, the following individual cases for film motion can be highlighted: The final expressions for the film thickness and the vertical and lengthwise velocity components are listed below in respect that the latest obtained equations: u(x, y, z) = g ϑ y w(x, y, z) = where the liquid volume concentration c p determines the initial thickness δ 0 = 1 2 c p B of the liquid film. To prove the reliability of the proposed mathematical model, a numerical example has been calculated for the following parameters: air blower FPZ K04 MS; nominal flow rate Q n = 137.0 m 3 /h; density ρ = 1 × 10 3 kg/m 3 ; liquid volume concentration c p = 0.012; dimensionless density ratio k = 0.9; kinematic viscosity ν = 1 × 10 −4 m 2 /s; gravitational acceleration g = 9.81 m/s; channel width B = 0.165 m; channel height H = 0.050 m; channel length L = 0.200 m.
As a result of numerical calculations, the charts for film thickness at its boundaries are shown in Figure 2 for the above initial data and Equation (23). Figure 2a shows the shape of the gas-liquid interfacial surface at the start and end of the film. Figure 2b shows the shape of the interfacial surface at the top and bottom of the film. Overall, Figure 2c presents the 3D graph of the film thickness function.
at the top and bottom of the film. Overall, Figure 2c presents the 3D graph of the film thickness function.
Additionally, the values of the velocity components on the interface have a scientific and practical interest. Excluding the variable y from the expression for the film thickness δ(x, z) allows for evaluating the following velocity components on the interfacial surface: It should be additionally noted that the liquid film downflows at a certain angle to the vertical sedimentation surface. This angle is determined considering that the streamlines on the phase interface matched with the phase velocities at the interface surface (u S = dx dt , w S = dz dt ) [41]. The expression for tgα = dz dx = k 1 + 2 k τ ρgδ is the result of considering (26) expression, and considering its component x decreases exponential function. Thus, the angle is equal to α = arctg[2L/(δ 0 ·Cr)] for the case Cr >> 0.5kδ 0 /L. This angle is about 28 • at the inlet and about 46 • at the outlet of the channel. The velocity components average values [42] are determined by averaging of the expressions (24), (25) over the film thickness: Additionally, the values of the velocity components on the interface have a scientific and practical interest. Excluding the variable y from the expression for the film thickness δ(x, z) allows for evaluating the following velocity components on the interfacial surface: Energies 2020, 13, 1938 11 of 15 Additionally, for the inlet cross-section, Figure 3a presents the components of the average film flow velocity, and Figure 3b shows the elements of the film flow velocity on the interfacial surface.
Additionally, for the inlet cross-section, Figure 3a presents the components of the average film flow velocity, and Figure3b shows the elements of the film flow velocity on the interfacial surface. Additionally, the proposed methodology allows for developing the draining channels in the separation device. Notably, the calculated Reynolds number for the liquid film = 2 / is about 31, which is more than critical value Recr = 24. This effect allows for avoiding the use of drainage channels into the design of the separation device. Notably, for the abovementioned case study, the total cross-sectional area of the draining channel = ( − ) is equal to 6.3·10 -6 m 2 , which particularly corresponds to 8 draining holes of the diameter 1 mm.

Conclusions
The mathematical model for a liquid film flow in the three-dimensional formulation has been developed with consideration of the change in its thickness in both vertical and horizontal directions by the implementation of the simplifications, assumptions, and order evaluation of the components included in the system of the Navier-Stokes equations and continuity equation. The reliability of this model has been proved by the numerical calculation of the flow rate. Notably, the average value of the axial velocity component on the interfacial surface corresponds to the nominal value for the air blower with the allowable relative error.
The proposed model can be useful for the development of draining channels inside of the separation equipment, the design improvement of the heat and mass transfer equipment (e.g., plate absorbers, and heat exchangers). The draining channels' development can be realized in future research by using the above model for dynamic and inertia-filtrating separation devices. The main parameters for the abovementioned purpose are the expressions for the film thickness in the dependence of the wall length on its top and bottom, for the velocity components and the averaged values. According to the obtained numerical results, the average value of the axial velocity component on the interfacial surface w av = 1 H H 0 w(x, 0)dx is equal to 0.017 m/s. Consequently, the average air velocity w a = w s + 0.5Bµ(dw/dy) y = δ /µ a is equal to 4.4 m/s, and the flow rate of the air blower Q 0 = w a BH is equal to 132.7 m 3 /h. On the other hand, the average velocity for the nominal flow rate of the air blower w n av = Q n /(BH) is equal to 4.6 m/s. Consequently, the relative error of determining the velocity ε w = w av −w n av w n av is equal to 4.3%. Additionally, the proposed methodology allows for developing the draining channels in the separation device. Notably, the calculated Reynolds number for the liquid film Re w = 2δw av /ν is about 31, which is more than critical value Re cr = 24. This effect allows for avoiding the use of drainage channels into the design of the separation device. Notably, for the abovementioned case study, the total cross-sectional area of the draining channel S = νH 2u (Re w − Re cr ) is equal to 6.3 × 10 −6 m 2 , which particularly corresponds to 8 draining holes of the diameter 1 mm.

Conclusions
The mathematical model for a liquid film flow in the three-dimensional formulation has been developed with consideration of the change in its thickness in both vertical and horizontal directions by the implementation of the simplifications, assumptions, and order evaluation of the components included in the system of the Navier-Stokes equations and continuity equation. The reliability of this model has been proved by the numerical calculation of the flow rate. Notably, the average value of the axial velocity component on the interfacial surface corresponds to the nominal value for the air blower with the allowable relative error.
The proposed model can be useful for the development of draining channels inside of the separation equipment, the design improvement of the heat and mass transfer equipment (e.g., plate absorbers, and heat exchangers). The draining channels' development can be realized in future research by using the above model for dynamic and inertia-filtrating separation devices. The main parameters for the abovementioned purpose are the expressions for the film thickness in the dependence of the wall length on its top and bottom, for the velocity components and the averaged values.
One of the most useful results is the dimensionless film shape criterion. It shows that the film shape depends on the stresses along the film length, the wall length, the fluid density, and the gravitational acceleration. The limiting value of this criterion is 0.23. This criterion indicates areas with a maximum fluid concentration on the walls. According to this, designing the local draining channels can be proposed. In particular, if the criterion is greater/lower than its limiting value, the liquid film thickness decreases/increases in the direction of the positive x-coordinate (vertically) and increases/decreases in the direction of the positive z-coordinate (horizontally), respectively. For the particular case study, the film thickness arises by 4.8% at the inlet along with the sedimentation surface height. It decreases by 4.8% and 18.8% at the end and top of the liquid, respectively. The comparison of the calculated results with the existing data for the chosen air blower proves the reliability of the proposed model. In particular, the relative error of determining the average velocity component is about 4.3%.