Numerical Simulation of Falling Liquid Film Flow on a Vertical Plane by Two-Phase Lattice BoltzmannMethod

Falling liquid �lm �ow is widely used inmany processes. Supplementary to experimental studies, Navier-Stokes-basedmodels have been employed for describing �lm �ow phenomena. ese models are o�en 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 �ow by lattice Boltzmann methods. erefore, the well-known Shan-Chen model (Shan and Chen, 1993) has been employed to an isothermal falling liquid �lm. e 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. e computational results show that lattice Boltzmann methods are capable to model falling liquid �lm �ow and that the �owmorphology is in qualitatively good agreement with other numerical and experimental works.


Introduction
Falling liquid �lm �ow is a critical operation in many engineering processes, such as densi�cation 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 �lm becomes wavy from certain Reynolds numbers Re = ( on (see [1][2][3][4][5]) and that waviness generally leads to a heat transfer augmentation.
e theoretical and numerical description of falling liquid �lms is still a challenge.Existing models are either limited in their generality or need enormous computational resources.ese are evolution equations of the �lm 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 �lm 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 �lm �owing down an inclined plane and investigated the stability of the �lm for certain conditions.In the last decades, additional physical effects have been implemented, such as phase change, thermocapillarity, non-Newtonian �uids, electromagnetic interactions, in�uence of ultrasound, and fully temperature-dependent properties (for an overview, see [10]).Shkadov, 1967 [11] was the �rst 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 �lm thickness, the volumetric �ow rate, and the wall shear stress.is model is employable for moderate Reynolds numbers (Re ≲ 100) and large Kapitza numbers (Ka = ( 4 ( 3  < 8 × 10 −3 ) 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 �lm 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 �lm being subjected to nonhomogeneous heating.Miyara, 1999 [17]; Miyara, 2000 [18]; Miyara, 2001 [19] investigated sinusoidally disturbed falling liquid �lms 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 �ow.e governing equations were discretised by employing the �nite 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.e 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 �lm on a grooved surface with the Volume-of-Fluid method for interface tracking.e author considered thermocapillarity and also disruption of the �lm.
Lattice Boltzmann methods are a relatively new approach for simulating �ow phenomena.Based on molecular gas dynamics, �rst 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][29][30][31][32]), have been mostly employed for describing drops or �ow within porous media.To the best knowledge of the present authors, lattice Boltzmann methods have not yet been employed for simulating falling liquid �lm �ows.
e aim of this study is to provide an approach for falling liquid �lm �ow modelling by employing lattice Boltzmann methods.For this purpose, we assume isothermal two-phase �ow 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 �ow morphology from wave-less to laminar-wavy liquid �lm �ow correctly.e 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 �lm �ow are shown.

Numerical Investigation
For solving the hydrodynamics of a falling liquid �lm, the Shan-Chen model [29] has been employed, which is an isothermal lattice Boltzmann model for multiphase and multicomponent (with  components) �ows.

Governing Equations. e 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 Δ.e 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].us, with   =    2 s for the relaxation time,  and  s for the kinematic viscosity and speed of sound, respectively, and  eq  for the equilibrium particle density distribution function (the exponent  will be dropped in the following, since only single-component �ows with  = 1 are being investigated here), which can be derived from the Maxwell distribution by Taylor series and reads Herein,   is weighting factors of the speci�c velocity directions in a D2Q9 lattice (D2Q9 implies two space dimensions and nine velocity directions) (see Figure 1) with the values e vector  eq represents the macroscopic equilibrium velocity and the vectors   are the microscopic velocity vectors in the respective directions, that read Herein,  = 1 lu ts −1 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   =  P / LB , velocity   =  P / LB , and time   =   /  ).e variable  is the mass density at a site  and can be expressed by the zeroth moment of the distribution function with the molecular mass .e equation for the macroscopic equilibrium velocity  eq consists of two parts.e �rst 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 �eld.us, e vector of all forces  is the superposition of three forces, that is, wherein  INT is the interphase force,  FW is the �uid-wall force and  GR is the gravitational force.e interphase force can be calculated with where  is the cohesion parameter and  is the interaction potential that is de�ned by e 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 �rst term represents the ideal gas law and the second term represents a nonlinear part being shaped like the van-der-Waals equation.e �uid-wall forces can be calculated by wherein  ads is the adsorption parameter and the variable e interrelation of the cohesion parameter  and the adsorption parameter  ads is represented by whereby the contact angle  at the three-phase point can be correlated with the density  [34] as e 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].

Geometry and Mesh.
A falling liquid �lm is investigated, �owing down an inclined wall with an angle  against horizon, whereupon it is assumed that the �lm only develops along the wall.is implies that the domain can be reduced to two dimensions as shown in Figure 2. In this �gure, also the boundary conditions are marked.More details about the boundary conditions are given in Section 2.3.
e applied discretisation scheme (D2Q9) can be used only for uniform lattices, wherein the spacing is limited by two dependencies as follows.(i) �ecessary resolution of �lm thickness (�uzmin and Derksen, 2011 [36] showed that the �lm thickness should be at least twice as thick as the interface thickness.Here, we have an interface thickness of �ve grid points.) (ii) For neglecting compressibility effects, By virtue of Obviously, the demand for a small Mach number undergoes the demand for su�cient resolution of the �lm� thus, the critical space resolution is e extent of the computational domain depends upon the length of �lm under consideration and upon the prospective wave magnitude.e dimensions in both coordinate directions are given as a multiple of the inlet �lm thickness ℎ 0 , wherein the length is  = 1000ℎ 0 and the width is  = ℎ 0 .is results in a mesh of 30000 × 10 nodes, since the �lm thickness is resolved by 30 nodes.e grid independence test has been performed and the result for Re = 20 is shown in Figure 3.

Boundary and Initial
Conditions. e boundary conditions are shown in Figure 2. ere are three solid wall boundaries (A, C, D) modelled by bounce-back, a velocity inlet (B), and a zero gradient outlet.e velocity given in B is constant in space but may be varied in time with where  0 is the initial liquid velocity,  and  are the relative amplitude and frequency of the disturbance, respectively, whereas the velocity component in  direction remains   = 0. e missing incoming distribution functions are calculated by using the Zou and He boundary condition Zou and He, 1997 [37].e 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 �uid.
For initialising the domain, the �lm was supplied with a constant thickness ℎ 0 and a uniform velocity  0 = (0,  0 )  .e gas was assumed to be at rest.

Implementation and Validation of the Numerical Model.
e 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 �oiseuille �ow and �ow around a cylinder which showed an excellent performance of the code for single-phase �ows.e 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 threephase points [38].
e square-to-circle computations have been performed by de�ning squares of different sizes being transformed to circles, whose radii have been measured.e pressure difference between inside and outside of the drop has been also determined.In Figure 4, the pressure difference Δ is plotted versus the curvature 1/.It is shown that the results can be �tted 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      mu lu ts − .

Test Series and
where  F,P is the kinematic viscosity of the �lm and  P is the magnitude of the gravitational acceleration with an angle  against horizon.e Kapitza number, a dimensionless number, only depending on properties of the �uid and on gravitational acceleration, cannot be prede�ned in this model, since both �uid 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 �ow at Re ∈ {, 8, } (free transition) (ii) sinusoidally disturbed �ow at Re  .e parameters for the computations are shown in Table 1.

Results and Discussion
3.1.Nondisturbed Film.e hydrodynamics of falling liquid �lms without arti�cial disturbance at the inlet have been studied for three different Reynolds numbers representing speci�c changes in the �ow 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 �elds and the characteristic wave velocities.

Evolution of the
which implies that a �uid particle with the velocity  ,P has crossed the complete domain twice.is time corresponds to the time in discrete units of     ts.e wave evolution for Re   is shown in Figure 5.It can be observed that for    ( P   s), �rst waves occur from ℎ  ≲ 8 on with a small amplitude.In the next following points in time, the position of wave initiation �ows downstream, and it stabilises at ℎ  ≈  for times greater than  P ≳  s ( ≳ ).e waves occurring here are limited in their amplitude to approximately half of the initial �lm thickness ℎ  .e observations of these calculations correspond qualitatively good to experimental data obtained by Gross et al., 2009, [5], who found that a water �lm in a vertical tube has a �rst transition at Re ≈ 6 from a laminar wave-less shape to small waves.
In Figure 6, the evolution of the �lm surface is shown for Re  8. Contrary to Re  , a �rst solitary wave occurs �owing downstream faster than the other waves.is 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.e small wave amplitudes are higher than the ones occurring at Re  .According to Gross et al., 2009, [5] a water �lm has so-called �small waves� at Re  8, which also can be observed here.
e wave evolution for Re   is shown in Figure 7. Similar to Re  8, a �rst solitary wave establishes whose magnitude increases until the top of the wave reaches the upper side of the domain ( = 1 at  0 ≈ 340).is 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 Δ( 0 ) ≈ 150 in  = 1 to Δ( 0 ) ≈ 400 in   4. 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  0 = 2. Gross et al., 2009 [5] observed waves with increased magnitude and the occurrence of some bow waves (see Figure 7,  = 9 at  0 ≈ 410 and Figure 9).e �gures shown here demonstrate that it is possible to model the hydrodynamics of falling li�uid �lms by lattice Boltzmann methods, but it should be mentioned that the direct comparison of these results with experimental data is di�cult, due to the limitations of the model� the �uid densities are coupled with the surface tension by the cohesion parameter  which reduces the degrees of freedom for realistic cases.( = 3.348 s).�n both �gures, the �lm (grey) is indicated and the velocity �eld is represented by arrows.Similar to the previously shown �gures, the coordinate in �ow direction is strongly compressed.Figure 8 shows a section of the transition region from laminar wave�less to laminar wavy �lm.ere are sinusoidal waves with a wavelength of approximately 10 0 but with increasing amplitude, which corresponds to the experimental �ndings, where initially small waves occur whose amplitude increases with time.e velocity �eld shows a cocurrent motion of �lm 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 �lm seems to �ow towards the wall ( 0 = 0).Both effects appear to be due to the interphase and �uid�wall forces modelled by interaction potentials.ese spurious currents have been previously described by Chen et al., 2007 [26], and Sukop and orne jr.2007 [34] who found out that the currents at the interface do not represent a mass transfer through the interface.e magnitude of the largest parasitic currents at the interface is approximately 1 luts whilst the gas bulk velocity is approximately 0.1 luts.In Figure 9, a wave with a bow wave in front of it is shown.Gross et al., 2009 [5] have discovered that water �lms exhibit bow waves for Re = 20 which can also be seen here.Obviously, both of these waves have a length of approximately 25ℎ 0 and a much larger mightiness than the waves shown in Figure 8.

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 400 mm/s for Re = 20.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 Re = 20 it is approximately zero.

Sinusoidally Disturbed Film.
Apart from falling liquid �lms with a natural transition to wavy interface, computations have been performed with a sinusoidally disturbed �lm, and the results are compared to those obtained by Gao et al., 2003 [41] with a volume of �uid 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.e reason for this can be found in the fact that the Shan-Chen model allows only limited variation of the Weber number (We = ℎ 0  2 /), and hence, the Weber number in our case is We = 2.4 whilst Gao et al., 2003 [41] have a value of We = 0.24.

Conclusion
�n the present investigation, falling liquid �lms 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.e implementation of the model has been validated against various single-and two-phase �ow problems and proved its validity.e simulation of falling liquid �lms has been conducted for cases with Reynolds numbers of Re ∈ {5, 8, 20} which showed that the employed Shan-Chen model is capable to simulate falling liquid �lm �ows and that there is a good qualitative agreement with experimental data.Characteristic �ow morphologies have been calculated for their respective Reynolds numbers.Direct comparison of the wave pro�les 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 �lms but the constraints given by the Shan-Chen model limit their generality.e investigation of liquid �lm �ow with heat transfer, larger Reynolds numbers, and other geometrical parameters will be conducted in the future.

7 F 1 :
Velocity directions  and velocity vectors   in D2Q9 lattice.

F 2 :
Geometry of the computational domain (black box); A: wall, B: liquid inlet velocity, C and D: walls, E: gradient-free outlet; : inclination angle against horizon.

F 3 :
Grid independence test for Re = 20 without sinusoidal disturbances with the given resolutions at  = 02 s.

2 F 4 :
R in lu−1  Pressure difference ∆p in mu ts −Evaluation of the pressure difference depending on curvature of the drop.

F 7 :
Evolution of the interface with Re = 20 at different points in time� the sub�gures 0 to 9 show the evolution for increasing time until  9 = 2 0,P = 3.348 s, whereby   =   0.32 s, with   0, , 2,  , 9.
Parameters. e characteristic velocity  ,P and the initial �lm thickness ℎ ,P may be calculated with Interface.Figures 5 to 7 illustrate the evolution of the interface with time for different Reynolds numbers.Every �gure is composed of ten sub�gures in which the �lm thickness is plotted versus the strongly compressed axial position in the �lm.e sub�gures are shown in a subsequent order and have a speci�c value for a variable , whereupon the respective real time is the product of  and a characteristic time step.Every last sub�gure (  ), that is, the uppermost one, shows the �lm interface at time