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

Authors: Andreas Hantsch; Ulrich Gross
Journal: Journal of Engineering
 (2013)

---

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


---

## Body

## 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 numbersRe=(uh)/ν 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 (Re≲100) and large Kapitza numbers (Ka=(η4g)/(ρσ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 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
(1)fiα(x+eiΔt,t+Δt)-fiα(x,t)=Ωα,
where fiα is the particle density distribution function in the velocity direction i at a point in time t and a time step Δt. 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,
(2)Ωα=-1τα+0.5(fiα(x,t)-fiα,eq(x,t)),
with τα=ναcs2 for the relaxation time, ν and cs for the kinematic viscosity and speed of sound, respectively, and fiα,eq for the equilibrium particle density distribution function (the exponent α will be dropped in the following, since only single-component flows with α=1 are being investigated here), which can be derived from the Maxwell distribution by Taylor series and reads
(3)fieq(x,t)=wiρ[1+3ei·ueqc2+92(ei·ueq)2c4-32(ueq)2c2].
Herein, wi 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
(4)wi={49,fori=0,19,fori=1,…,4,136,fori=5,…,8.Figure 1
Velocity directionsi and velocity vectors ei in D2Q9 lattice.The vectorueq represents the macroscopic equilibrium velocity and the vectors ei are the microscopic velocity vectors in the respective directions, that read
(5)eic=(010-101-1-110010-111-1-1)fori∈{0,1,…,8}.
Herein, c=1luts-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 ζl=LP/LLB, velocity ζu=uP/uLB, and time ζt=ζl/ζu). The variable ρ is the mass density at a site x and can be expressed by the zeroth moment of the distribution function
(6)ρ(x,t)=m∑i=08‍fi,
with the molecular mass m. The equation for the macroscopic equilibrium velocity ueq 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 F that determine the velocity field. Thus,
(7)ueq=mρ∑i=08‍eifi+τFρ.
The vector of all forces F is the superposition of three forces, that is,
(8)F=FINT+FFW+FGR,
wherein FINT is the interphase force, FFW is the fluid-wall force and FGR is the gravitational force. The interphase force can be calculated with
(9)FINT(x,t)=-Gψ(x,t)∑i=18‍wiψ(x+eiΔt,t)ei,
where G is the cohesion parameter and ψ is the interaction potential that is defined by
(10)ψ(x,t)=ψ0exp(-ρ0ρ)withρ=ρ(x,t),(11)ψ0=4.0,ρ0=200.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
(12)p=ρ3+G6ψ2(x,t),
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(13)FFW(x,t)=-Gadsψ(x,t)∑i=18‍wis(x+eiΔt)ei,
wherein Gads is the adsorption parameter and the variable s(14)s(x+eiΔt)={1,ifx+eiΔtisawallnode,0,ifx+eiΔtisnotawallnode.
The interrelation of the cohesion parameter G and the adsorption parameter Gads is represented by
(15)Gads=Gψ(ρ(x)),
whereby the contact angle θ at the three-phase point can be correlated with the density ρ [34] as
(16)ψ(ρ(x))={ψ(ρF)forθ=0∘,0.5[ψ(ρF)+ψ(ρG)]forθ=90∘,ψ(ρG)forθ=180∘.The computation of the gravitational force can be performed with(17)FGR(x,t)=ρ(x,t)g,
where g 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 2.3.Figure 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.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.)
(18)δx≲110h0,h0=1.0lu.(ii)
For neglecting compressibility effects,(19)Ma=ucS≤0.1=3u=3δtδx≈3ηtδx2δx=3ηtδxlu≤0.1.
By virtue of
(20)u=δtδx=ηtδxwithδt∝δx2.Obviously, the demand for a small Mach number undergoes the demand for sufficient resolution of the film; thus, the critical space resolution is
(21)δx≲0.1(3ηtlu)-1withηt=2.0.
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 h0, wherein the length is L=1000h0 and the width is H=5h0. This results in a mesh of 30000×150 nodes, since the film thickness is resolved by 30 nodes. The grid independence test has been performed and the result for Re=20 is shown in Figure 3.Figure 3
Grid independence test forRe=20 without sinusoidal disturbances with the given resolutions at t=0.2s.


### 2.3. Boundary and Initial Conditions

The boundary conditions are shown in Figure2. 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
(22)uy=u0[1+εsin(2πFt)],
where u0 is the initial liquid velocity, ε and F are the relative amplitude and frequency of the disturbance, respectively, whereas the velocity component in x direction remains ux=0. 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 allfi at the boundary node equal to the corresponding value of fi at the neighbouring node in the fluid.For initialising the domain, the film was supplied with a constant thicknessh0 and a uniform velocity u0=(0,u0)T. 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 Figure4, the pressure difference Δp is plotted versus the curvature 1/R. 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 σ=(12.45±0.32)muluts-2.Figure 4
Evaluation of the pressure difference depending on curvature of the drop.


### 2.5. Test Series and Parameters

The characteristic velocityu0,P and the initial film thickness h0,P may be calculated with
(23)u0,P=13Re2νF,PgPsinβ3,h0,P=3ReνF,P2gPsinβ3,
respectively, following Nusselt, 1916, [39]. Herein, Re is the Reynolds number defined as
(24)Re=u0,Ph0,PνF,P,
where νF,P is the kinematic viscosity of the film and gP is the magnitude of the gravitational acceleration with an angle β against horizon. The Kapitza number,
(25)Ka=μF4gρFσ3,
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 G 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 (β=90° as follows):(i)
non-disturbed flow atRe∈{5,8,20} (free transition)(ii)
sinusoidally disturbed flow atRe=20.1.The parameters for the computations are shown in Table1.Table 1
Parameters for computation.
Re = 5
Re = 8
Re = 20
Re = 20.1
h
0
,
P/mm
0.115
0.135
0.183
0.184
u
0
,
P/(mm/s)
43.50
59.51
109.5
109.9
H
P/mm
0.578
0.676
0.917
0.919
L
P/mm
115.5
135.2
183.3
183.8
t
0
,
P/s
2.655
2.272
1.674
1.672
ε
0
0
0
0.05
F/(1/s)
0
0
0
45
g
P/(m/s2)
2.452
2.452
4.660
4.415
Ka/10−3
28.58
3.544
0.050
0.049
β/°
90.0
ν
F
,
P/(m2/s)
1.004809 × 10−6
G
−120.0
G
ads
−187.16
ρ
F/(mu/lu2)
528.8
ρ
G/(mu/lu2)
86.4
σ/(mu lu/ts2)
14.33


## 2.1. Governing Equations

The lattics Boltzmann equation of the particle density distribution function of a componentα reads
(1)fiα(x+eiΔt,t+Δt)-fiα(x,t)=Ωα,
where fiα is the particle density distribution function in the velocity direction i at a point in time t and a time step Δt. 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,
(2)Ωα=-1τα+0.5(fiα(x,t)-fiα,eq(x,t)),
with τα=ναcs2 for the relaxation time, ν and cs for the kinematic viscosity and speed of sound, respectively, and fiα,eq for the equilibrium particle density distribution function (the exponent α will be dropped in the following, since only single-component flows with α=1 are being investigated here), which can be derived from the Maxwell distribution by Taylor series and reads
(3)fieq(x,t)=wiρ[1+3ei·ueqc2+92(ei·ueq)2c4-32(ueq)2c2].
Herein, wi 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
(4)wi={49,fori=0,19,fori=1,…,4,136,fori=5,…,8.Figure 1
Velocity directionsi and velocity vectors ei in D2Q9 lattice.The vectorueq represents the macroscopic equilibrium velocity and the vectors ei are the microscopic velocity vectors in the respective directions, that read
(5)eic=(010-101-1-110010-111-1-1)fori∈{0,1,…,8}.
Herein, c=1luts-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 ζl=LP/LLB, velocity ζu=uP/uLB, and time ζt=ζl/ζu). The variable ρ is the mass density at a site x and can be expressed by the zeroth moment of the distribution function
(6)ρ(x,t)=m∑i=08‍fi,
with the molecular mass m. The equation for the macroscopic equilibrium velocity ueq 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 F that determine the velocity field. Thus,
(7)ueq=mρ∑i=08‍eifi+τFρ.
The vector of all forces F is the superposition of three forces, that is,
(8)F=FINT+FFW+FGR,
wherein FINT is the interphase force, FFW is the fluid-wall force and FGR is the gravitational force. The interphase force can be calculated with
(9)FINT(x,t)=-Gψ(x,t)∑i=18‍wiψ(x+eiΔt,t)ei,
where G is the cohesion parameter and ψ is the interaction potential that is defined by
(10)ψ(x,t)=ψ0exp(-ρ0ρ)withρ=ρ(x,t),(11)ψ0=4.0,ρ0=200.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
(12)p=ρ3+G6ψ2(x,t),
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(13)FFW(x,t)=-Gadsψ(x,t)∑i=18‍wis(x+eiΔt)ei,
wherein Gads is the adsorption parameter and the variable s(14)s(x+eiΔt)={1,ifx+eiΔtisawallnode,0,ifx+eiΔtisnotawallnode.
The interrelation of the cohesion parameter G and the adsorption parameter Gads is represented by
(15)Gads=Gψ(ρ(x)),
whereby the contact angle θ at the three-phase point can be correlated with the density ρ [34] as
(16)ψ(ρ(x))={ψ(ρF)forθ=0∘,0.5[ψ(ρF)+ψ(ρG)]forθ=90∘,ψ(ρG)forθ=180∘.The computation of the gravitational force can be performed with(17)FGR(x,t)=ρ(x,t)g,
where g 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 2.3.Figure 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.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.)
(18)δx≲110h0,h0=1.0lu.(ii)
For neglecting compressibility effects,(19)Ma=ucS≤0.1=3u=3δtδx≈3ηtδx2δx=3ηtδxlu≤0.1.
By virtue of
(20)u=δtδx=ηtδxwithδt∝δx2.Obviously, the demand for a small Mach number undergoes the demand for sufficient resolution of the film; thus, the critical space resolution is
(21)δx≲0.1(3ηtlu)-1withηt=2.0.
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 h0, wherein the length is L=1000h0 and the width is H=5h0. This results in a mesh of 30000×150 nodes, since the film thickness is resolved by 30 nodes. The grid independence test has been performed and the result for Re=20 is shown in Figure 3.Figure 3
Grid independence test forRe=20 without sinusoidal disturbances with the given resolutions at t=0.2s.


## 2.3. Boundary and Initial Conditions

The boundary conditions are shown in Figure2. 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
(22)uy=u0[1+εsin(2πFt)],
where u0 is the initial liquid velocity, ε and F are the relative amplitude and frequency of the disturbance, respectively, whereas the velocity component in x direction remains ux=0. 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 allfi at the boundary node equal to the corresponding value of fi at the neighbouring node in the fluid.For initialising the domain, the film was supplied with a constant thicknessh0 and a uniform velocity u0=(0,u0)T. 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 Figure4, the pressure difference Δp is plotted versus the curvature 1/R. 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 σ=(12.45±0.32)muluts-2.Figure 4
Evaluation of the pressure difference depending on curvature of the drop.


## 2.5. Test Series and Parameters

The characteristic velocityu0,P and the initial film thickness h0,P may be calculated with
(23)u0,P=13Re2νF,PgPsinβ3,h0,P=3ReνF,P2gPsinβ3,
respectively, following Nusselt, 1916, [39]. Herein, Re is the Reynolds number defined as
(24)Re=u0,Ph0,PνF,P,
where νF,P is the kinematic viscosity of the film and gP is the magnitude of the gravitational acceleration with an angle β against horizon. The Kapitza number,
(25)Ka=μF4gρFσ3,
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 G 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 (β=90° as follows):(i)
non-disturbed flow atRe∈{5,8,20} (free transition)(ii)
sinusoidally disturbed flow atRe=20.1.The parameters for the computations are shown in Table1.Table 1
Parameters for computation.
Re = 5
Re = 8
Re = 20
Re = 20.1
h
0
,
P/mm
0.115
0.135
0.183
0.184
u
0
,
P/(mm/s)
43.50
59.51
109.5
109.9
H
P/mm
0.578
0.676
0.917
0.919
L
P/mm
115.5
135.2
183.3
183.8
t
0
,
P/s
2.655
2.272
1.674
1.672
ε
0
0
0
0.05
F/(1/s)
0
0
0
45
g
P/(m/s2)
2.452
2.452
4.660
4.415
Ka/10−3
28.58
3.544
0.050
0.049
β/°
90.0
ν
F
,
P/(m2/s)
1.004809 × 10−6
G
−120.0
G
ads
−187.16
ρ
F/(mu/lu2)
528.8
ρ
G/(mu/lu2)
86.4
σ/(mu lu/ts2)
14.33


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

Figures5 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 k, whereupon the respective real time is the product of k and a characteristic time step. Every last subfigure (k=9), that is, the uppermost one, shows the film interface at time
(26)tt0,P=tu0,PLDh0,P=2withLD=LPh0,P=1000,
which implies that a fluid particle with the velocity u0,P has crossed the complete domain twice. This time corresponds to the time in discrete units of t=900000ts.Figure 5
Evolution of the interface withRe=5 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=5.310s, whereby tk=k·0.590s, with k∈{0,1,2,…,9}.The wave evolution forRe=5 is shown in Figure 5. It can be observed that for k=1 (tP=0.590s), first waves occur from y/h0≲800 on with a small amplitude. In the next following points in time, the position of wave initiation flows downstream, and it stabilises at y/h0≈200 for times greater than tP≳2.950s (k≳5). The waves occurring here are limited in their amplitude to approximately half of the initial film thickness h0. 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 Re≈6 from a laminar wave-less shape to small waves.In Figure6, the evolution of the film surface is shown for Re=8. Contrary to Re=5, 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 Re=5. According to Gross et al., 2009, [5] a water film has so-called “small waves” at Re=8, which also can be observed here.Figure 6
Evolution of the interface withRe=8 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=5.544s, whereby tk=k·0.504s, with k∈{0,1,2,…,9}.Figure 7
Evolution of the interface withRe=20 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=3.348s, whereby tk=k·0.372s, with k∈{0,1,2,…,9}.The wave evolution forRe=20 is shown in Figure 7. Similar to Re=8, a first solitary wave establishes whose magnitude increases until the top of the wave reaches the upper side of the domain (k=1 at y/h0≈340). 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 Δ(y/h0)≈150 in k=1 to Δ(y/h0)≈400 in k≳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 h/h0=2. Gross et al., 2009 [5] observed waves with increased magnitude and the occurrence of some bow waves (see Figure 7, k=9 at y/h0≈410 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 parameterG 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 Figures8 and 9 for Re=20 at k=9 (t=3.348s). 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
Section of the velocity field forRe=20 (Figure 7, k=9, t=3.348s).Figure 9
Section of the velocity field forRe=20 (Figure 7, k=9, t=3.348s).Figure8 shows a section of the transition region from laminar wave-less to laminar wavy film. There are sinusoidal waves with a wavelength of approximately 10h0 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 (x/h0=0). 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 1lu/ts whilst the gas bulk velocity is approximately 0.1lu/ts.In Figure9, 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 Re=20 which can also be seen here. Obviously, both of these waves have a length of approximately 25h0 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 Table2, 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 400mm/s for Re=20. For lower Reynolds numbers no experimental data are available.Table 2
Wave velocity of solitary waves (solitary w.) and beginning wave production (b. w. production).
Re
u
wave in mm/s
Solitary w.
b. w. production
5
—
32.59
8
126.86
42.74
20
387.50
≈0It can be seen in Table2 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.


### 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 Figure10, 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 (We=ρh0u2/σ), 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.Figure 10
Comparison of the present results forRe=20.1, F=45s-1 and ε=0.05 with Gao et al., 2003 [41].


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

Figures5 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 k, whereupon the respective real time is the product of k and a characteristic time step. Every last subfigure (k=9), that is, the uppermost one, shows the film interface at time
(26)tt0,P=tu0,PLDh0,P=2withLD=LPh0,P=1000,
which implies that a fluid particle with the velocity u0,P has crossed the complete domain twice. This time corresponds to the time in discrete units of t=900000ts.Figure 5
Evolution of the interface withRe=5 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=5.310s, whereby tk=k·0.590s, with k∈{0,1,2,…,9}.The wave evolution forRe=5 is shown in Figure 5. It can be observed that for k=1 (tP=0.590s), first waves occur from y/h0≲800 on with a small amplitude. In the next following points in time, the position of wave initiation flows downstream, and it stabilises at y/h0≈200 for times greater than tP≳2.950s (k≳5). The waves occurring here are limited in their amplitude to approximately half of the initial film thickness h0. 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 Re≈6 from a laminar wave-less shape to small waves.In Figure6, the evolution of the film surface is shown for Re=8. Contrary to Re=5, 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 Re=5. According to Gross et al., 2009, [5] a water film has so-called “small waves” at Re=8, which also can be observed here.Figure 6
Evolution of the interface withRe=8 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=5.544s, whereby tk=k·0.504s, with k∈{0,1,2,…,9}.Figure 7
Evolution of the interface withRe=20 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=3.348s, whereby tk=k·0.372s, with k∈{0,1,2,…,9}.The wave evolution forRe=20 is shown in Figure 7. Similar to Re=8, a first solitary wave establishes whose magnitude increases until the top of the wave reaches the upper side of the domain (k=1 at y/h0≈340). 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 Δ(y/h0)≈150 in k=1 to Δ(y/h0)≈400 in k≳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 h/h0=2. Gross et al., 2009 [5] observed waves with increased magnitude and the occurrence of some bow waves (see Figure 7, k=9 at y/h0≈410 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 parameterG 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 Figures8 and 9 for Re=20 at k=9 (t=3.348s). 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
Section of the velocity field forRe=20 (Figure 7, k=9, t=3.348s).Figure 9
Section of the velocity field forRe=20 (Figure 7, k=9, t=3.348s).Figure8 shows a section of the transition region from laminar wave-less to laminar wavy film. There are sinusoidal waves with a wavelength of approximately 10h0 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 (x/h0=0). 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 1lu/ts whilst the gas bulk velocity is approximately 0.1lu/ts.In Figure9, 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 Re=20 which can also be seen here. Obviously, both of these waves have a length of approximately 25h0 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 Table2, 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 400mm/s for Re=20. For lower Reynolds numbers no experimental data are available.Table 2
Wave velocity of solitary waves (solitary w.) and beginning wave production (b. w. production).
Re
u
wave in mm/s
Solitary w.
b. w. production
5
—
32.59
8
126.86
42.74
20
387.50
≈0It can be seen in Table2 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.


## 3.1.1. Evolution of the Interface

Figures5 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 k, whereupon the respective real time is the product of k and a characteristic time step. Every last subfigure (k=9), that is, the uppermost one, shows the film interface at time
(26)tt0,P=tu0,PLDh0,P=2withLD=LPh0,P=1000,
which implies that a fluid particle with the velocity u0,P has crossed the complete domain twice. This time corresponds to the time in discrete units of t=900000ts.Figure 5
Evolution of the interface withRe=5 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=5.310s, whereby tk=k·0.590s, with k∈{0,1,2,…,9}.The wave evolution forRe=5 is shown in Figure 5. It can be observed that for k=1 (tP=0.590s), first waves occur from y/h0≲800 on with a small amplitude. In the next following points in time, the position of wave initiation flows downstream, and it stabilises at y/h0≈200 for times greater than tP≳2.950s (k≳5). The waves occurring here are limited in their amplitude to approximately half of the initial film thickness h0. 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 Re≈6 from a laminar wave-less shape to small waves.In Figure6, the evolution of the film surface is shown for Re=8. Contrary to Re=5, 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 Re=5. According to Gross et al., 2009, [5] a water film has so-called “small waves” at Re=8, which also can be observed here.Figure 6
Evolution of the interface withRe=8 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=5.544s, whereby tk=k·0.504s, with k∈{0,1,2,…,9}.Figure 7
Evolution of the interface withRe=20 at different points in time; the subfigures 0 to 9 show the evolution for increasing time until t9=2t0,P=3.348s, whereby tk=k·0.372s, with k∈{0,1,2,…,9}.The wave evolution forRe=20 is shown in Figure 7. Similar to Re=8, a first solitary wave establishes whose magnitude increases until the top of the wave reaches the upper side of the domain (k=1 at y/h0≈340). 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 Δ(y/h0)≈150 in k=1 to Δ(y/h0)≈400 in k≳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 h/h0=2. Gross et al., 2009 [5] observed waves with increased magnitude and the occurrence of some bow waves (see Figure 7, k=9 at y/h0≈410 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 parameterG 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 Figures8 and 9 for Re=20 at k=9 (t=3.348s). 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
Section of the velocity field forRe=20 (Figure 7, k=9, t=3.348s).Figure 9
Section of the velocity field forRe=20 (Figure 7, k=9, t=3.348s).Figure8 shows a section of the transition region from laminar wave-less to laminar wavy film. There are sinusoidal waves with a wavelength of approximately 10h0 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 (x/h0=0). 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 1lu/ts whilst the gas bulk velocity is approximately 0.1lu/ts.In Figure9, 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 Re=20 which can also be seen here. Obviously, both of these waves have a length of approximately 25h0 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 Table2, 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 400mm/s for Re=20. For lower Reynolds numbers no experimental data are available.Table 2
Wave velocity of solitary waves (solitary w.) and beginning wave production (b. w. production).
Re
u
wave in mm/s
Solitary w.
b. w. production
5
—
32.59
8
126.86
42.74
20
387.50
≈0It can be seen in Table2 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.


## 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 Figure10, 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 (We=ρh0u2/σ), 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.Figure 10
Comparison of the present results forRe=20.1, F=45s-1 and ε=0.05 with Gao et al., 2003 [41].


## 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 ofRe∈{5,8,20} 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.



---

*Source: 484137-2013-01-17.xml*