Abstract
Falling liquid film flow is widely used in many processes. Supplementary to experimental studies, Navier-Stokes-based models have been employed for describing film flow phenomena. These models are often disadvantageous since they are either strongly limited in their generality or need enormous computational resources. In this investigation, a new approach is proposed for modelling flow by lattice Boltzmann methods. Therefore, the well-known Shan-Chen model (Shan and Chen, 1993) has been employed to an isothermal falling liquid film. The validity of the implementation has been checked against some single-phase and two-phase reference cases. Test series have been conducted for three different Reynolds numbers without external disturbances and for one Reynolds number with sinusoidally pulsating inlet velocity. The computational results show that lattice Boltzmann methods are capable to model falling liquid film flow and that the flow morphology is in qualitatively good agreement with other numerical and experimental works.
1. Introduction
Falling liquid film flow is a critical operation in many engineering processes, such as densification of suspensions or emulsions or the separation of mixtures, liquefaction in the condenser of a power plant or heat and mass transfer in a heat pipe. Experiments and previous simulations have shown that the film becomes wavy from certain Reynolds numbers on (see [1–5]) and that waviness generally leads to a heat transfer augmentation.
The theoretical and numerical description of falling liquid films is still a challenge. Existing models are either limited in their generality or need enormous computational resources. These are evolution equations of the film thickness (also known as Benney equation, long-wave equation), (weighted residual) integrated boundary layer equations (IBL, WIBL, WRIBL), and the full Navier-Stokes equations, with increasing computational demand.
Based on the two-dimensional Navier-Stokes equations, evolution equations of the film thickness have been designed, by introducing a stream function and assuming long disturbance wavelengths and Reynolds number of the order of one. First research has been performed by Yih, 1955 [6]; Yih, 1963 [7], Benjamin, 1957 [8], and Benney, 1966 [9], who derived the evolution equation for a film flowing down an inclined plane and investigated the stability of the film for certain conditions. In the last decades, additional physical effects have been implemented, such as phase change, thermocapillarity, non-Newtonian fluids, electromagnetic interactions, influence of ultrasound, and fully temperature-dependent properties (for an overview, see [10]).
Shkadov, 1967 [11] was the first to develop the integrated boundary layer equations from Navier-Stokes equations and qualitatively good agreement with experiments in the nonlinear regime was obtained. It fails, however, to represent the critical behaviour. Ruyer-Quil and Manneville, 1998 [12]; Ruyer-Quil and Manneville 2000 [13]; Ruyer-Quil and Manneville, 2002 [14] extended the model of Shkadov and created a new one, based on weighted residual integrated boundary layer equations (WRIBL), which represent experimental data much better. Nguyen and Balakotaiah, 2000 [15] gave a literature review and designed a model consisting of three partial differential equations for the film thickness, the volumetric flow rate, and the wall shear stress. This model is employable for moderate Reynolds numbers () and large Kapitza numbers () and in this range, it is much better than other models as viscous effects of higher order were considered together with pressure corrections across the film thickness. Scheid, 2003 [16] developed a second-order WRIBL model which was capable to simulate the transition from two-dimensional to three-dimensional waves and the reaction of a film being subjected to nonhomogeneous heating.
Miyara, 1999 [17]; Miyara, 2000 [18]; Miyara, 2001 [19] investigated sinusoidally disturbed falling liquid films by solving the governing equations with the marker and cell algorithm for different Prandtl, Reynolds, Weber, and Froude numbers as well as for condensation. Groß, 2007 [20] and Soemers, 2008 [21] created a model for studying three-dimensional multiphase flow. The governing equations were discretised by employing the finite element method. In addition to the solution of the governing equations, the discretisation method offers an error estimation. Banerjee et al., 2004 [22], Fulgosi et al., 2003 [23], and Lakehal et al., 2003 [24] simulated turbulence and transport of scalar quantities across interfaces with direct numerical simulations for carbon dioxide in oceans, as an example. The authors found the turbulence development close to the interface to be similar to the turbulence development close to rigid walls but less anisotropic. Gambaryan-Roisman, 2007 [25] investigated the behaviour of a liquid film on a grooved surface with the Volume-of-Fluid method for interface tracking. The author considered thermocapillarity and also disruption of the film.
Lattice Boltzmann methods are a relatively new approach for simulating flow phenomena. Based on molecular gas dynamics, first models were presented in the late 1980s [26] (for a background, see [27]). Multiphase and multicomponent models, based on different approaches (e.g., see [28–32]), have been mostly employed for describing drops or flow within porous media. To the best knowledge of the present authors, lattice Boltzmann methods have not yet been employed for simulating falling liquid film flows.
The aim of this study is to provide an approach for falling liquid film flow modelling by employing lattice Boltzmann methods. For this purpose, we assume isothermal two-phase flow that can be described by the model of Shan and Chen, 1993 [29]. Numerical tests are performed at low Reynolds numbers in order to investigate the ability of the model to predict the change in flow morphology from wave-less to laminar-wavy liquid film flow correctly. The paper is organised as follows. In Section 2 the numerical model is shown and details of the mathematical model, the geometry and mesh, the boundary and initial conditions, some preliminary tests and the parameters of the calculations are presented. In Section 3, the results of the calculations for disturbed and nondisturbed film flow are shown.
2. Numerical Investigation
For solving the hydrodynamics of a falling liquid film, the Shan-Chen model [29] has been employed, which is an isothermal lattice Boltzmann model for multiphase and multicomponent (with components) flows.
2.1. Governing Equations
The lattics Boltzmann equation of the particle density distribution function of a component reads where is the particle density distribution function in the velocity direction at a point in time and a time step . The collision term represents the local interaction of the particles and is modelled with the well-known BGK approximation that leads to a linearisation of the collision term [33]. Thus, with for the relaxation time, and for the kinematic viscosity and speed of sound, respectively, and for the equilibrium particle density distribution function (the exponent will be dropped in the following, since only single-component flows with are being investigated here), which can be derived from the Maxwell distribution by Taylor series and reads Herein, is weighting factors of the specific velocity directions in a D2Q9 lattice (D2Q9 implies two space dimensions and nine velocity directions) (see Figure 1) with the values
The vector represents the macroscopic equilibrium velocity and the vectors are the microscopic velocity vectors in the respective directions, that read Herein, is the reference velocity in lattice units (lattice units can be converted to physical units either by nondimensional numbers or by conversion factors. E., g., for the length , velocity , and time ). The variable is the mass density at a site and can be expressed by the zeroth moment of the distribution function with the molecular mass . The equation for the macroscopic equilibrium velocity consists of two parts. The first part depends on the particle density distribution function at the current time and the second part depends on the effect of external forces that determine the velocity field. Thus, The vector of all forces is the superposition of three forces, that is, wherein is the interphase force, is the fluid-wall force and is the gravitational force. The interphase force can be calculated with where is the cohesion parameter and is the interaction potential that is defined by
The interphase force leads to the effects of phase separation and surface tension. From (9) and (10), the equation of state can be calculated with wherein the first term represents the ideal gas law and the second term represents a nonlinear part being shaped like the van-der-Waals equation.
The fluid-wall forces can be calculated by wherein is the adsorption parameter and the variable The interrelation of the cohesion parameter and the adsorption parameter is represented by whereby the contact angle at the three-phase point can be correlated with the density [34] as
The computation of the gravitational force can be performed with where is the vector of the gravitational acceleration.
By using a Chapman-Enskog expansion (see, e.g., [35]), Shan and Chen, 1993 [29] have derived the macroscopic balancing equations at the limits of low disturbance frequency and long wave lengths. Some thermodynamic properties of the model have been additionally discussed by He and Doolen, 2002 [30] and Chen et al., 2007 [26].
2.2. Geometry and Mesh
A falling liquid film is investigated, flowing down an inclined wall with an angle against horizon, whereupon it is assumed that the film only develops along the wall. This implies that the domain can be reduced to two dimensions as shown in Figure 2. In this figure, also the boundary conditions are marked. More details about the boundary conditions are given in Section .
The applied discretisation scheme (D2Q9) can be used only for uniform lattices, wherein the spacing is limited by two dependencies as follows.(i)Necessary resolution of film thickness (Kuzmin and Derksen, 2011 [36] showed that the film thickness should be at least twice as thick as the interface thickness. Here, we have an interface thickness of five grid points.) (ii)For neglecting compressibility effects, By virtue of Obviously, the demand for a small Mach number undergoes the demand for sufficient resolution of the film; thus, the critical space resolution is The extent of the computational domain depends upon the length of film under consideration and upon the prospective wave magnitude. The dimensions in both coordinate directions are given as a multiple of the inlet film thickness , wherein the length is and the width is . This results in a mesh of nodes, since the film thickness is resolved by nodes. The grid independence test has been performed and the result for is shown in Figure 3.
2.3. Boundary and Initial Conditions
The boundary conditions are shown in Figure 2. There are three solid wall boundaries (A, C, D) modelled by bounce-back, a velocity inlet (B), and a zero gradient outlet. The velocity given in B is constant in space but may be varied in time with where is the initial liquid velocity, and are the relative amplitude and frequency of the disturbance, respectively, whereas the velocity component in direction remains . The missing incoming distribution functions are calculated by using the Zou and He boundary condition Zou and He, 1997 [37].
The outlet (E) was assumed to be free of gradients, which is realised by setting all at the boundary node equal to the corresponding value of at the neighbouring node in the fluid.
For initialising the domain, the film was supplied with a constant thickness and a uniform velocity . The gas was assumed to be at rest.
2.4. Implementation and Validation of the Numerical Model
The previously described model has been implemented in a self-written serial C++-program, which has been validated with single- and two-phase benchmark tests. Computations have been conducted for Poiseuille flow and flow around a cylinder which showed an excellent performance of the code for single-phase flows. The two-phase benchmarks were the transition from a liquid square to a circle and the following determination of the surface tension as well as the investigation of the contact angle formation at the three-phase points [38].
The square-to-circle computations have been performed by defining squares of different sizes being transformed to circles, whose radii have been measured. The pressure difference between inside and outside of the drop has been also determined. In Figure 4, the pressure difference is plotted versus the curvature . It is shown that the results can be fitted by a linear regression line and that the deviation is relatively high for small curvatures, which occurs due to low pressure differences. From the regression line, the surface tension can be calculated to be .
2.5. Test Series and Parameters
The characteristic velocity and the initial film thickness may be calculated with respectively, following Nusselt, 1916, [39]. Herein, is the Reynolds number defined as where is the kinematic viscosity of the film and is the magnitude of the gravitational acceleration with an angle against horizon. The Kapitza number, a dimensionless number, only depending on properties of the fluid and on gravitational acceleration, cannot be predefined in this model, since both fluid densities and the surface tension are coupled by the cohesion parameter and the viscosity is linked to the Reynolds number.
By employing the present implementation, the following test cases have been investigated for a vertical wall ( as follows):(i)non-disturbed flow at (free transition)(ii)sinusoidally disturbed flow at .
The parameters for the computations are shown in Table 1.
3. Results and Discussion
3.1. Nondisturbed Film
The hydrodynamics of falling liquid films without artificial disturbance at the inlet have been studied for three different Reynolds numbers representing specific changes in the flow morphology according to Gross et al., 2009. [5].
In the following sections, the wave shapes depending on time and space are shown as well as the velocity fields and the characteristic wave velocities.
3.1.1. Evolution of the Interface
Figures 5 to 7 illustrate the evolution of the interface with time for different Reynolds numbers. Every figure is composed of ten subfigures in which the film thickness is plotted versus the strongly compressed axial position in the film. The subfigures are shown in a subsequent order and have a specific value for a variable , whereupon the respective real time is the product of and a characteristic time step. Every last subfigure (), that is, the uppermost one, shows the film interface at time which implies that a fluid particle with the velocity has crossed the complete domain twice. This time corresponds to the time in discrete units of .
The wave evolution for is shown in Figure 5. It can be observed that for (), first waves occur from on with a small amplitude. In the next following points in time, the position of wave initiation flows downstream, and it stabilises at for times greater than (). The waves occurring here are limited in their amplitude to approximately half of the initial film thickness . The observations of these calculations correspond qualitatively good to experimental data obtained by Gross et al., 2009, [5], who found that a water film in a vertical tube has a first transition at from a laminar wave-less shape to small waves.
In Figure 6, the evolution of the film surface is shown for . Contrary to , a first solitary wave occurs flowing downstream faster than the other waves. This wave has its origin in the starting of the calculation. During the whole computation process no new solitary wave was observed for this Reynolds number. It can be seen that small waves appear only behind the solitary wave and that the beginning of the wave production is also moving downstream. The small wave amplitudes are higher than the ones occurring at . According to Gross et al., 2009, [5] a water film has so-called “small waves” at , which also can be observed here.
The wave evolution for is shown in Figure 7. Similar to , a first solitary wave establishes whose magnitude increases until the top of the wave reaches the upper side of the domain ( at ). This wave and also the waves further downstream appear due to the starting of the computation. Upstream the solitary wave a small wave producing region establishes whose position remains the same during the whole computation period, but with increasing length from in to in . From the wave producing region, waves are released whose magnitude is about twice the magnitude of the small waves, and some of them grow to more than . Gross et al., 2009 [5] observed waves with increased magnitude and the occurrence of some bow waves (see Figure 7, at and Figure 9).
The figures shown here demonstrate that it is possible to model the hydrodynamics of falling liquid films by lattice Boltzmann methods, but it should be mentioned that the direct comparison of these results with experimental data is difficult, due to the limitations of the model; the fluid densities are coupled with the surface tension by the cohesion parameter which reduces the degrees of freedom for realistic cases.
3.1.2. Velocity Field
For better insight, two sections of the flow field are shown in Figures 8 and 9 for at (). In both figures, the film (grey) is indicated and the velocity field is represented by arrows. Similar to the previously shown figures, the coordinate in flow direction is strongly compressed.
Figure 8 shows a section of the transition region from laminar wave-less to laminar wavy film. There are sinusoidal waves with a wavelength of approximately but with increasing amplitude, which corresponds to the experimental findings, where initially small waves occur whose amplitude increases with time. The velocity field shows a cocurrent motion of film and gas with the gas velocity being higher. Additionally, two effects can be seen. First, there are large velocity vectors directing in the gas phase and second, the film seems to flow towards the wall (). Both effects appear to be due to the interphase and fluid-wall forces modelled by interaction potentials. These spurious currents have been previously described by Chen et al., 2007 [26], and Sukop and Thorne jr. 2007 [34] who found out that the currents at the interface do not represent a mass transfer through the interface. The magnitude of the largest parasitic currents at the interface is approximately whilst the gas bulk velocity is approximately .
In Figure 9, a wave with a bow wave in front of it is shown. Gross et al., 2009 [5] have discovered that water films exhibit bow waves for which can also be seen here. Obviously, both of these waves have a length of approximately and a much larger mightiness than the waves shown in Figure 8.
3.1.3. Wave Velocity
Both the wave velocity of a solitary wave and the velocity of a wave producing region can be evaluated by measuring wave positions at two consecutive points in time. In Table 2, the wave velocities are presented for the simulations being performed here. Good qualitatively agreement is obtained by comparing the results with experimental data shown by Philipp et al., 2006 [40] who reported an approximate wave velocity of for . For lower Reynolds numbers no experimental data are available.
It can be seen in Table 2 that with increasing Reynolds number, the solitary wave velocity also increases but the velocity of the wave producing region almost remains constant, and for it is approximately zero.
3.2. Sinusoidally Disturbed Film
Apart from falling liquid films with a natural transition to wavy interface, computations have been performed with a sinusoidally disturbed film, and the results are compared to those obtained by Gao et al., 2003 [41] with a volume of fluid approach and Navier-Stokes equations.
In Figure 10, the results of the present computations are presented and also those from literature. It can be seen that both approaches show large waves with smaller ones in-between. However, there is only a qualitative agreement between both. The reason for this can be found in the fact that the Shan-Chen model allows only limited variation of the Weber number (), and hence, the Weber number in our case is whilst Gao et al., 2003 [41] have a value of .
4. Conclusion
In the present investigation, falling liquid films on a vertical plane have been modelled with the Shan-Chen model, which is an isothermal multiphase and multicomponent lattice Boltzmann model based on interaction potentials. The implementation of the model has been validated against various single- and two-phase flow problems and proved its validity.
The simulation of falling liquid films has been conducted for cases with Reynolds numbers of which showed that the employed Shan-Chen model is capable to simulate falling liquid film flows and that there is a good qualitative agreement with experimental data. Characteristic flow morphologies have been calculated for their respective Reynolds numbers. Direct comparison of the wave profiles is difficult, since limitations of the Shan-Chen model prevent that the Kapitza number can be freely chosen.
To conclude, it can be stated that the lattice Boltzmann methods are capable to simulate falling liquid films but the constraints given by the Shan-Chen model limit their generality. The investigation of liquid film flow with heat transfer, larger Reynolds numbers, and other geometrical parameters will be conducted in the future.
Acknowledgment
The authors greatly acknowledge the valuable comments and suggestions given by the unknown reviewers.