Three-Dimensional Mathematical Model of the Liquid Film Downflow on a Vertical Surface

Abstract

Film downflow 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 film downflow allows for the determination of the hydrodynamic parameters of the liquid film flow 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 fluid density. Additionally, equations for determining the averaged film thickness, the averaged velocity vectors over the film thickness, the longitudinal and vertical velocity components, and the initial angle of streamline deviation from the vertical axis were analytically obtained.

1. 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 the fluid film movement, heat and mass transfer in column equipment, the separation of gas–liquid mixtures, and fluid flow in heat exchangers. Scientists for the above purpose have investigated the heat and mass transfer processes occurring in microchannels [3,4,5,6] and free and induced liquid film flow over different surfaces [7,8,9,10,11,12,13,14,15,16,17,18] and inside of different channels [19,20,21,22,23,24,25]. This flow occurs due to the influence of various forces, which determine the film formation or destruction, and the inclination angle of the surface. At the same time, the developed mathematical models and correlation coefficients proposed for them should consider each of the above process features, and the possibility of changing the physical state of the flow components forms a liquid film (condensation, evaporation).

The liquid film motion inside of microchannels was considered in [3,4,5,6] by the application of numerical and experimental methods as well as developing mathematical models. The research [3] describes the development of the criterion model for the prediction of the critical heat flux premature caused by periodic liquid film movement in microchannels. The research implements the simplifications and assumptions for the Navier–Stokes equations, such as neglecting the gravity and inertia forces, considering the film movement only along the channel, and introducing the lubrication approximation. The authors of [4] studied the falling water film movement inside of microchannel with the “water–carbon dioxide” mixing section. The obtained results consider the mass exchange between two phases by numerical simulation for the case of laminar flow and multiphase VOF (volume of fluid) model with the indirect comparison between the numerical simulation and experimental results. The Lattice Boltzmann method (LBM) is commonly used for liquid flow motion in microchannels [5,6], especially in cases of studying hydrodynamic phenomena during the heat-mass transfer processes [6], which involve formation and breakup process of the liquid film in this channel type.

Studies on falling film flow on vertical and angled surfaces have been carried out in [7,8,9,10,11,12,13] by experimental and numerical methods. These studies use the VOF model. However, the comparison of the predicted and experimental values showed that the VOF model is unsuitable for predicting the hydrodynamic entrance length of Newtonian liquid film over a cylindrical surface. A new empirical model and an improved minimal surface model were developed based on experimental data and used to predict the hydrodynamic entrance length. However, usage of the VOF model for the laminar wavy motions of liquid film flow over flat surfaces in combination with the CSF (Continuum Surface Force) model gives good agreement with the experimental data [8,9].

For better accuracy and further improvement of the results obtained by CFD (Computational Fluid Dynamics) methods, the authors of the study [10] implement the MAC (Marker and Cell) and the moving grid. The application of the VOF model is enough for the study of a non-Newtonian liquid film [11]. The Lattice Boltzmann method can be used not only for liquid film flow in microchannels but also for fluid film flow with different Reynolds numbers on the vertical surface [12]. It should be stressed that, in the case of the Shan–Chen model, usage takes place at the generality limit of the Lattice Boltzmann method. The authors of the paper [13] developed the mathematical model of liquid film flow on flat elements of an aircraft. In this case, film thickness changes along the surface. The dependences of the falling film hydrodynamic parameters on fluid viscosity, surface tension, and primary flow rate were found during the above research. Additionally, the influence of primary flow and phase interaction parameters on the local and general characteristics of the film thickness, its distribution over the surface, the wave formation, and their height were considered. Doubtless, authors describe the influence of surface topology, angle of inclination, and its material on the film flow behavior.

The authors of [14,15,16,17,18] have developed mathematical models for the liquid film flow on vertical surfaces by adopting various simplifications and assumptions into the Navier–Stokes equations. The study [14] complements the abovementioned system by the energy equation. As a result, the authors evaluated the dimensionless temperature field and the evaporation rate. Liquid film evaporation and the irregularity of its heating lead to a temperature gradient at the gas–liquid phase interface, and, consequently, Marangoni convection appears. During the research, the Marangoni effect leads to the film breakdown; on the other hand, the gravity force and surface tension stabilize 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 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.

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

$$\{\begin{array}{c}\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right)=\rho g-\frac{\partial p}{\partial x}+\mu {\nabla }^{2}u+\frac{1}{3}\mu \frac{\partial }{\partial x}divV;\\ \rho \left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}\right)=-\frac{\partial p}{\partial y}+\mu {\nabla }^{2}v+\frac{1}{3}\mu \frac{\partial }{\partial y}divV;\\ \rho \left(\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\right)=-\frac{\partial p}{\partial z}+\mu {\nabla }^{2}w+\frac{1}{3}\mu \frac{\partial }{\partial z}divV;\\ \frac{\partial \rho }{\partial t}+\frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial \rho w}{\partial z}=0,\end{array}$$

(1)

where ρ—liquid density, kg/m³; µ—dynamic viscosity of the liquid, Pa·s; p—dynamic pressure of a fluid, Pa; g—gravitational acceleration, m/s².

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 x- and z-directions. These axes are directed downward and along the wall, respectively. The design scheme is shown in Figure 1;

-

liquid flow is stationary $\left(\frac{\partial u}{\partial t}=0;\frac{\partial v}{\partial t}=0;\frac{\partial w}{\partial t}=0;\frac{\partial \rho }{\partial t}=0\right)$, and the film flows without wave formation;

-

the liquid film thickness δ varies along the height H and length L of the channel; herewith, it is quite small (δ << L, and δ << H). Due to the abovementioned, the liquid film velocity component in the direction normal to the surface is small in comparison with the components in the main flow directions;

-

the pressure change is slow along with the height of the separation channel due to the specific shape of the vane type of the mist extractors (separation devices). It means that $\frac{\partial p}{\partial x}=0$.

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]. For the first equation of the system:

-

comparison of the elements $\nu \frac{\partial u}{\partial y}$ and $w\frac{\partial u}{\partial z}$ from the left side of the equation shows that the component $w\frac{\partial u}{\partial z}$ can be neglected. The order of magnitude of the following components are the following: film thickness δ ($\partial y$) is about 10⁻³ m, the channel length L($\partial z$) is about 10⁻¹ m. The velocity component directed normal to the wall—10⁻² m/s. The above allows for the conclusion that the order of the $\frac{w·\raisebox{1ex}{$\partial u$}\!\left/ \!\raisebox{-1ex}{$\partial z$}\right.}{v·\raisebox{1ex}{$\partial u$}\!\left/ \!\raisebox{-1ex}{$\partial y$}\right.}\cong \frac{W}{V}\frac{\raisebox{1ex}{$U$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}{\raisebox{1ex}{$U$}\!\left/ \!\raisebox{-1ex}{$\delta $}\right.}\cong \frac{W}{V}\frac{\delta }{L}\cong \frac{{10}^{-1}}{{10}^{-2}}\frac{{10}^{-4}}{{10}^{-1}}\cong {10}^{-2}$. Since $\frac{V}{W}\cong \frac{{10}^{-2}}{{10}^{-1}}\cong {10}^{-1}$, the velocity component v can be neglected too.

-

the components from the right side of the equation are supplemented by replacing the variables with their values on the liquid film surface u → U, v → V, w → W, x → H (the wall height along which the liquid flows), y → δ, z → L. The following values $\frac{U}{H^2}$, $\frac{U}{{\delta }^2}$, $\frac{U}{L^2}$ are obtained for comparison, $\frac{U}{{\delta }^2}$ has the highest order of them. That is why only $\frac{{\partial }^{2}u}{\partial 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 $U\frac{W}{H}$, $V\frac{W}{\delta }$, $W\frac{W}{L}$, $\rho \frac{U^2}{H}$, $\mu \frac{W}{H^2}$, $\mu \frac{W}{{\delta }^2}$, $\mu \frac{W}{L^2}$ are obtained. The orders are equal to about ${10}^{-1}\frac{{10}^{-1}}{{10}^{-1}}\cong {10}^{-1}$, ${10}^{-2}\frac{{10}^{-1}}{{10}^{-1}}\cong {10}^{-2}$, ${10}^{-1}\frac{{10}^{-1}}{{10}^{-1}}\cong {10}^{-1}$, ${10}^3\frac{{10}^{-2}}{{10}^{-1}}\cong {10}^2$, ${10}^{-3}\frac{{10}^{-1}}{{10}^{-2}}\cong {10}^{-2}$, ${10}^{-3}\frac{{10}^{-1}}{{10}^{-6}}\cong {10}^2$, and ${10}^{-3}\frac{{10}^{-1}}{{10}^{-2}}\cong {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$: \frac{\partial p}{\partial z}$ and $\mu \frac{{\partial }^{2}w}{\partial 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:

$$\{\begin{array}{c}\rho \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=\rho g+\mu \frac{{\partial }^{2}u}{\partial y^2};\\ 0=-\frac{\partial p}{\partial z}+\mu \frac{{\partial }^{2}w}{\partial y^2};\\ \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0.\end{array}$$

(2)

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:

$$u\left(x,y,z\right)=a\left(x,z\right)y+b\left(x,z\right)y^2+c\left(x,z\right)y^3.$$

(3)

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]: ${{\tau }_{yx}=\mu \frac{\partial u}{\partial y}|}_{y=\delta \left(x,z\right)}=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:

$$u=\frac{g}{\nu }y\left(\delta -\frac{y}{2}\right)+cy\left(y^2-3{\delta }^2\right),$$

(4)

where ν—kinematic viscosity, m²/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:

$$\frac{\partial p}{\partial z}=-\left(1-\frac{{\rho }_{mix}}{\rho }\right)g,$$

(5)

where ρ_mix—gas–liquid mixture density, kg/m³.

The substitution of the abovementioned expression to the second Equation (2) with the consideration of the dimensionless density ratio k = $\left(1-\frac{{\rho }_{mix}}{\rho }\right)$, consequent integration and transformation allow obtaining the following expression for the velocity component:

$$w=\left[\frac{\tau }{\mu }+\frac{kg}{\nu }\left(\delta -\frac{y}{2}\right)\right]y.$$

(6)

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 ${\tau }_{yz}=\mu {\frac{\partial w}{\partial y}|}_{y=\delta }=\tau \left(x,z\right)$.

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:

$$\frac{g}{\nu }\frac{\partial \delta }{\partial x}-3\frac{\partial }{\partial x}\left(c{\delta }^2\right)+\frac{kg}{\nu }\frac{\partial \delta }{\partial z}=0.$$

(7)

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:

$$\delta \left(x,z\right)=U\left(x\right)W\left(z\right).$$

(8)

The result of the two functions multiplication (8) is introduced instead of the film thickness in expression (7):

$$\frac{g}{\nu U}\frac{dU}{dx}-3UW\frac{\partial c}{\partial x}-6c\frac{dU}{dx}W+k\frac{g}{\nu W}\frac{dW}{dz}=0.$$

(9)

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:

$$\frac{g}{\vartheta }\frac{U_x^{\prime }}{U}=-k\frac{g}{\vartheta }\frac{W_z^{\prime }}{W}.$$

(10)

According to $\frac{g}{\nu U}\frac{dU}{dx}=-k\frac{g}{\nu W}\frac{dW}{dz}=const$, the above equation can be written as a result of the two constants multiplication $\frac{-g}{\nu }{c}_3$; one of them (c₃) is unknown. After dividing into two differential equations, the thickness is as follows:

$$\delta ={\delta }_0{e}^{c_3\left(\frac{z}{k}-x\right)}.$$

(11)

where ${\delta }_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:

$$3UW\frac{\partial c}{\partial x}=-6c\frac{dU}{dx}W.$$

(12)

The next step is (12) integration, after variables’ separation and $\frac{U_x{}^{\prime }}{U}=c_3$ accounting:

$$c\left(x,z\right)=c_0\left(z\right)e^{-2c_{3}x},$$

(13)

where c₀(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₃ from formulas (11) and (13):

$$\frac{\partial }{\partial x}{{\displaystyle \int }}_0^{\delta }{u}^{2}dy-g\delta =-\nu {\frac{\partial u}{\partial y}|}_{y=0}.$$

(14)

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 ${\frac{\partial u}{\partial y}|}_{y=0}=\frac{g\delta }{\nu }-3c{\delta }^2$:

$$\left(\frac{2}{3}\frac{g^2}{\nu }{\delta }^2-\frac{61}{10}c\frac{g{\delta }^3}{\nu }+\frac{68}{5}{c}^2{\delta }^4\right)\frac{\partial \delta }{\partial x}+\left(\frac{136}{35}c{\delta }^5-\frac{61}{60}\frac{g{\delta }^4}{\nu }\right)\frac{\partial c}{\partial x}=3\nu c.$$

(15)

The above equation allows for the getting of the quadratic equation relative to the parameter c₀ 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]:

$$c_0^2\frac{612}{35}{\delta }_0^5{e}^{5c_3\frac{z}{k}}{c}_3+c_0\left(3\nu -\frac{427}{60}\frac{g}{\upsilon }{\delta }_0^4{c}_3{e}^{4c_3\frac{z}{k}}\right)+\frac{2}{3}\frac{g^2}{{\nu }^2}{\delta }_0^2{e}^{2c_3\frac{z}{k}}=0.$$

(16)

This equation has two roots. Only one case ($3\nu <<\frac{427}{60}\frac{g}{\nu }{\delta }_0^4{c}_3{e}^{4c_3\frac{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\ge \frac{k}{c_3}ln\frac{0.92}{c_3{\delta }_0}$. As a result, the condition $c_3{\delta }_0\ge 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{\delta }_0>>0.92$:

$$c_0\left(z\right)\approx 0.093\frac{g}{\nu c_3{\delta }_0^2}\left(1-\frac{2c_3}{k}z\right).$$

(17)

The final equation for the vertical velocity component u has the following form:

$$u\left(x,y,z\right)=\frac{g}{\nu }y\left[\delta -\frac{y}{2}+\frac{0.093}{c_3{\delta }_0^2}\left(y^2-3{\delta }^2\right)\left(1-\frac{2c_3}{k}z\right)e^{-2c_{3}x}\right].$$

(18)

The constant c₃ 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:

$$Q_{u0}+Q_{w0}=Q_{u1}+Q_{w1},$$

(19)

where $Q_{u0}={{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{\delta \left(o,z\right)}u\left(0,y,z\right)dydz$—the liquid flow rate, which is entering the film vertically from the top; $Q_{w0}={{\displaystyle \int }}_0^H{{\displaystyle \int }}_0^{\delta \left(x,0\right)}w\left(x,y\right)dydx$—the liquid flow rate, which is entering the film horizontally at the origin; $Q_{u1}={{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{\delta \left(H,z\right)}u\left(H,y,z\right)dydz$—the liquid flow rate, which is entering the film vertically from the bottom; $Q_{w1}={{\displaystyle \int }}_0^H{{\displaystyle \int }}_0^{\delta \left(x,0\right)}w\left(x,y\right)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:

$$R\left(c_3\right)=\Delta Q_u+\Delta Q_w=0,$$

(20)

where R—the error function—$\Delta Q_u=Q_{u1}-Q_{u0}$ is the difference of the liquid flow rates in a vertical direction, and $\Delta 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): $\left|\Delta Q_u\right|=\left|\Delta Q_w\right|$.

The exclusion of the terms of higher-order smallness from the balances’ equation for the flow rates (18) allows for obtaining the following expressions:

$$\Delta Q_u\approx -0.058\frac{gL{\delta }_0^2}{\upsilon c_3}{e}^{\frac{4c_{3}L}{k}};\Delta Q_w\approx \frac{\tau {\delta }_0^2}{4\mu c_3}{e}^{\frac{2c_{3}L}{k}}.$$

(21)

In this case, the (19) equation has the following solution:

$$c_3=\frac{k}{2L}ln\frac{4.31\tau }{\rho gL}.$$

(22)

This equation emphasizes that the liquid film shape depends on the constant c₃. Notably, its sign indicates directions of increasing or decreasing film thickness. This fact makes it necessary to introduce a dimensionless film shape criterion $Cr=\frac{\tau }{\rho gL}$, the limit value of which is 0.23.

Thus, the following individual cases for film motion can be highlighted:

-

if $\frac{\tau }{\rho gL}>0.23$, the liquid film thickness decreases in the direction of the positive x-coordinate (vertically) and increases in the direction of the positive z-coordinate (horizontally);

-

if $\frac{\tau }{\rho gL}=0.23$, the liquid film thickness does not change over the entire surface of its flow δ = δ₀;

-

if $\frac{\tau }{\rho gL}<0.23$, the fluid film thickness increases in the direction of the positive x coordinate (vertically) and decreases in the direction of the positive z-coordinate (horizontally).

The final expressions for the film thickness and the vertical and lengthwise velocity components are listed below in respect that the latest obtained equations:

$$\delta ={\delta }_0{\left(\frac{4.31\tau }{\rho gL}\right)}^{\frac{z-kx}{2L}},$$

(23)

$$\begin{array}{c}u\left(x,y,z\right)=\frac{g}{\vartheta }y\{{\delta }_0{\left(\frac{4.31\tau }{\rho gL}\right)}^{\frac{z-kx}{2L}}-\frac{y}{2}+\frac{0.186L}{k{\delta }_0^{2}ln\left(\frac{4.31\tau }{\rho gL}\right)}[y^2-\\ -3{\delta }_0^2{\left(\frac{4.31\tau }{\rho gL}\right)}^{\frac{z-kx}{L}}]\left[1-\frac{2z}{k}{\left(\frac{4.31\tau }{\rho gL}\right)}^{\frac{z-kx}{2L}}\right]{\left(\frac{4.31\tau }{\rho gL}\right)}^{-\frac{kx}{L}}\},\end{array}$$

(24)

$$w\left(x,y,z\right)=\left\{\frac{\tau }{\mu }+\frac{kg}{\vartheta }\left[{\delta }_0{\left(\frac{4.31\tau }{\rho gL}\right)}^{\frac{z-kx}{2L}}-\frac{y}{2}\right]\right\}y,$$

(25)

where the liquid volume concentration c_p determines the initial thickness ${\delta }_0=\frac{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³/h; density ρ = 1 × 10³ kg/m³; liquid volume concentration c_p = 0.012; dimensionless density ratio k = 0.9; kinematic viscosity ν = 1 × 10⁻⁴ m²/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.

The obtained results show that the film thickness varies at its start δ(x, 0) along with the sedimentation surface height in the range of 1.00–1.05 mm (increases by 4.8%). Additionally, the film thickness at the end δ(x, L) varies in a range of 0.81–0.85 mm (it increases by 4.8% ). At the top of the liquid film δ(0, z) along the sedimentation surface length, the thickness is in a range of 1.00–0.81 mm (it decreases by 18.8%). Finally, at the bottom, the thickness δ(H, z) is in a range of 1.05–0.85 mm (it decreases by 18.8%).

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=\frac{dx}{dt},$ $w_S=\frac{dz}{dt}$) [41]. The expression for tgα = $\frac{dz}{dx}=k\left(1+\frac{2}{k}\frac{\tau }{\rho g\delta }\right)$ is the result of considering (26) expression, and considering its component $\frac{0.372\delta }{C_3{\delta }_0^2}{e}^{-2C_{3}x}$ decreases exponential function. Thus, the angle is equal to α = arctg[2L/(δ₀·Cr)] for the case Cr >> 0.5kδ₀/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:

$$\overline{u}\left(x,z\right)=\frac{1}{\delta \left(x,z\right)}{{\displaystyle \int }}_0^{\delta \left(x,z\right)}u\left(x,y,z\right)dy=\frac{g{\delta }^2}{3\upsilon }\left(1-\frac{0.349}{C_3\delta }{e}^{-2C_{3}x}\right),$$

(26)

$$\overline{w}\left(x,z\right)=\frac{1}{\delta \left(x,z\right)}{{\displaystyle \int }}_0^{\delta \left(x,z\right)}w\left(x,y,z\right)dy=k\frac{g{\delta }^2}{3\upsilon }\left(1+\frac{3}{2k}\frac{\tau }{\rho g\delta }\right).$$

(27)

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:

$$\begin{array}{c}{u}_{\mathrm{s}}\left(x,z\right)=u\left(x,\delta \left(x,y\right),z\right)=\frac{g{\delta }^2}{2\nu }\cdot \left(1-\frac{0.372\delta }{C_3{\delta }_0^2}{e}^{-2C_{3}x}\right),\\ w_{\mathrm{s}}\left(x,z\right)=w\left(x,\delta \left(x,y\right),z\right)=k\frac{g{\delta }^2}{2\nu }\left(1+\frac{2}{k}\frac{\tau }{\rho g\delta }\right).\end{array}$$

(28)

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.

According to the obtained numerical results, the average value of the axial velocity component on the interfacial surface ${\overline{w}}_{av}=\frac{1}{H}{{\displaystyle \int }}_0^H\overline{w}\left(x,0\right)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₀ = w_aBH is equal to 132.7 m³/h. On the other hand, the average velocity for the nominal flow rate of the air blower ${\overline{w}}_{av}^n$ = Q_n/(BH) is equal to 4.6 m/s. Consequently, the relative error of determining the velocity ${\epsilon }_w=\left|\frac{{\overline{w}}_{av}-{\overline{w}}_{av}^n}{{\overline{w}}_{av}^n}\right|$ 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 $R{e}_w=2\delta {\overline{w}}_{av}/\nu$ 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=\frac{\nu H}{2\overline{u}}\left(R{e}_w-R{e}_{cr}\right)$ is equal to 6.3 × 10⁻⁶ m², which particularly corresponds to 8 draining holes of the diameter 1 mm.

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