Celerity-Amplitude Relations for Solitary Waves on Vertically Falling

ARTICLE pubs.acs.org/IECR

Celerity-Amplitude Relations for Solitary Waves on Vertically Falling Films Cesar E. Meza and Vemuri Balakotaiah* Department of Chemical and Biomolecular Engineering, University of Houston, Houston, Texas 77204-4004, United States ABSTRACT: Experimental data relating the amplitude and celerity of solitary waves on vertically falling and naturally excited films are presented over a wide of range of fluid properties and flow rates (2 < Ka < 3883, 0.08 < We < 35). It is found that for a fixed set of fluid properties (Ka) and flow rate (1/We), the celerity-amplitude relationship for all finite amplitude waves that exist on the film is nearly linear with a slope that decreases with the flow rate (1/We). For viscous fluids (2 < Ka < 129), the observed normalized maximum wave amplitude increased monotonically with the flow rate (1/We) and saturated at values around 3 (hmax ≈ 3) and found to be nearly independent of the fluid viscosity (Ka). For less viscous fluids (206 < Ka < 3883) the waves were found to be continuously evolving with maximum normalized amplitudes of about 10. A recently developed two equation h
q model is used to analyze the spatio-temporal dynamics of waves on the film. The model predictions corroborate the hypothesis that the existence of very large amplitude waves for less viscous fluids is due to the phenomenon of flow reversal (or existence of up-flows) near the wall. The model predicts a critical Ka value below which no flow reversal occurs. It is also shown that once the wave amplitude exceeds a critical value, the amplitude-celerity relationship is nearly linear for all solitary waves that exist on the film and is independent of the inlet forcing frequency or amplitude or initial conditions. Local bifurcation analysis in the traveling wave coordinate, computational results on extended domains with periodic inlet forcing, and computational results on finite domains with periodic boundary conditions are used to explain the experimentally observed wave celerity-ampliude relations on naturally excited films.

I. INTRODUCTION Wave formation on vertically falling films has been studied extensively since the early seminal work of Kapitza and Kapitza1 and Benjamin2 and Yih.3 Wavy films appear in many applications such as process equipment used in heat and mass transfer operations and chemical reactors. Understanding of the hydrodynamic characteristics of the wavy film is essential for predicting and developing accurate correlations for the heat and mass transfer rates in systems in which film flows occur. In addition to these practical aspects, the wavy film is also of theoretical interest as it is one of the very few known examples of fluid flow in which instability occurs at zero Reynolds number. The vertically falling wavy film is often cited as an example that displays complex spatio-temporal behavior at small Reynolds numbers. For these reasons, there is an extensive literature describing the experimental as well as the theoretical results on the behavior of vertically falling films. We give below a brief review of some experimental and theoretical results and refer to the monographs by Alekseenko et al.4 and Chang and Demekhin5 for a more detailed review of the relevant literature and results. Since the early work of Kapitza and Kapitza,1 experimental data on the wavy film behavior was presented by Krantz and Goren,6 Alekseenko et al.,7 Liu and Gollub,8 Liu et al.,9 Nosoko and Miyara,10 and many others. In addition to verifying the linear growth characteristics of the waves near the inception region, these studies have provided a wealth of data on the spatial and temporal evolution of finite amplitude periodic and solitary waves at low to moderate Reynolds numbers. More recently, Tihon et al.11,12 and Dietze et al.13 verified the existence of upflows in the front regions near the wall of large amplitude solitary waves on a vertically falling film. Alekseenko et al.14,15 presented r 2011 American Chemical Society

data on the characteristics of 3-D stationary solitary waves (also called horseshoe waves) evolving from a localized disturbance. On the theoretical side, the linear stability analysis of Benjamin2 and Yih3 established that the vertically falling film is unstable at all Reynolds numbers (i.e., the critical Reynolds number is zero). Since this early work, many simplified models were derived from the Navier
Stokes equations using boundary layer and long wave (lubrication) scaling. The purpose of these models is to obtain a simplified low-dimensional description of the film in terms of the streamwise coordinate (x) and time (t). These models may be classified as single evolution equations16,17 for the local film thickness h(x,t) or a pair of evolution equations18
21 for the film thickness and the local flow rate q(x,t). The analysis of single evolution equations has shown that their range of validity is limited to very small wave amplitudes (typically less than 20% of the flat film thickness). In recent work, Ruyer-Quil and Manneville18 and Mudunuri and Balakotaiah19 have reasoned that at least two independent local degrees of freedom are needed to describe the evolution of finite amplitude waves for a realistic range of parameters. These authors presented several two equation h
q models by extending the earlier approach of Shakdov20 and others. The accuracy of these models was also verified by comparing the predicted wavy film characteristics with the available numerical solutions of the two-dimensional Navier
Stokes equations.19 Special Issue: Ananth Issue Received: March 22, 2011 Accepted: June 9, 2011 Revised: June 6, 2011 Published: June 09, 2011 13258

dx.doi.org/10.1021/ie2005803 | Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research On the computational side, the two-dimensional Navier
Stokes equations describing vertically falling film were solved by many authors including Ramaswamy et al.,22 Malamataris and Vlachogiannis,23 Gao et al.,24 and Nosoko and Miyara.10 More recently, Malamataris and Balakotaiah25 and Dietze et al.13 presented computations describing the flow field underneath the large amplitude waves, including the regions in which upflows occur. The occurrence of stationary as well as oscillatory cellular patterns with up-flows in a vertically falling film is a surprising new phenomenon discovered only in the past decade, even though it was conjectured by Massot et al.26 nearly four decades ago. It should be pointed out that the numerical simulations in the literature were limited to the two-dimensional case and small domain sizes in the streamwise direction, typically less than one thousand flat film thicknesses. In spite of a large number of published theoretical and experimental studies reviewed above, at present we are still limited in our ability to predict the hydrodynamic characteristics of large amplitude waves that exist on vertically falling films for realistic values of flow rates and fluid physical properties. For example, most prior experimental and theoretical/computational studies of waves were restricted to cases in which the wave amplitudes did not exceed about two times the flat (Nusselt) film thickness. However, it is known that for less viscous fluids such as water, the wave amplitudes can be as large as ten times the flat film thickness.27 In addition, a few of the many questions that are relevant but not answered satisfactorily in the literature are as follows: (i) For a given set of fluid properties and flow rates, is there an upper bound on the wave amplitude (which is a measure of the nonlinearity)? If so, how is the amplitude of this largest wave related to the physical parameters characterizing the film? Under what transient or steady-state conditions is this largest wave realized? (ii) Is there a map or phase diagram in the parameter space that can describe the film complexity in different regimes? What is the boundary in the parameter space beyond which up-flows (or flow reversal) can occur? (iii) For periodic inlet forcing, how is the wave amplitude related to the forcing frequency and amplitude, especially for very low forcing frequencies? At what downstream distance does the wavy film reach an asymptotic state in which the noise insensitive film statistics such as the maximum amplitude, probability density function, and standard deviation of the film profile do not vary? (iv) For naturally excited waves (or noise driven structures) how are the wave amplitudes and film roughness related to the physical parameters? (v) What is the relationship between two-dimensional and three-dimensional solitary wave amplitudes and celerities? (vi) What are the boundaries in the parameter space crossing at which the shear wave instabilities start interacting with the free surface instabilities? (vii) What is the boundary in the parameter space beyond which there is no smooth solution to the film flow equations (i.e., the free surface curvature goes to infinity leading to wave breaking or drop formation)? The present work is an attempt to answer some of these questions partially. In the first part of this work, we present experimental data relating the amplitude and velocity of solitary waves on vertically falling and naturally excited films over a wide of range of fluid properties and flow rates. We show that for a given fluid and flow rate, there is a nearly linear relationship between the wave amplitude and celerity for all the finite amplitude waves that exist on the film. We also present data on the maximum observed wave amplitude as a function of flow rate (or reciprocal of the Weber number) and show that there is a clear distinction

ARTICLE

between the “viscous” and “less viscous” fluids characterized by the Kapitza number. In the former case, the maximum observed wave amplitudes are limited to about three times the flat (Nusselt) film thickness while in the latter case the maximum observed wave amplitudes are about ten times the flat film thickness. In the second part of this work, we present a corroborative theory using a simplified low-dimensional two partial differential equation (h
q) model along with local bifurcation analysis and numerical simulations. Even though the model examined is valid only in the visco-capillary regime, or when inertial effects in the film are not dominant, and covers only part of the experimental data, we show that it captures the main experimental trends, has acceptable accuracy within its region of validity, and provides insight on the behavior of the wavy film. This paper is organized as follows. In the next section, we provide a brief description of our experimental system and present the data on maximum observed wave amplitudes as a function of fluid properties and flow rate. We also present the data on wave velocity-wave amplitude for a selected set of fluids at different flow rates. In the following section, we present a corroborative theory, analysis and numerical simulations using a two equation h
q model and explain our experimental results. In the last section, we summarize the main contributions of this work along with the limitations of our experimental as well as theoretical/computational results and discuss some possible improvements or extensions.

II. EXPERIMENTAL RESULTS ON WAVE AMPLITUDES The main goal of the experiments was to measure the amplitude of the fully developed and saturated waves over a wide range of fluid viscosities and flow rates. Glycerin and water solutions of different concentrations were used to vary the viscosity in the range 1 to 250 cP. The densities and surface tension varied in much smaller ranges of 1.00 to 1.23 g/cm3 and 72.0 to 63.7 dyn/cm, respectively. Details of the experimental apparatus, procedures, accuracy of the conductance probe method used in measuring the film thickness, and data collection methods may be found in an earlier publication27and also in the thesis work of Meza.28 The wave amplitudes were measured in a tube of ID 38.1 mm, the length of the tubing from feed to the first set of probes (conductance probes 1 and 2) was 5.23 m, and the probes were separated by a distance of 2.54 cms (1 in.). The second set of conductance probes (probes 3 and 4) were located 1.74 m below the first set. This total length of about 7 m is larger than that used in most earlier studies of vertically falling films. The reason for using two sets of measuring points was for comparison of wave amplitudes and to verify if they have reached saturation values. In the case of high viscosity fluids the wave amplitudes were found to be saturated within experimental error [For the conductance probe technique, this error depends on the physical thickness of the film, and decreases with increasing film thickness. The typical error is about 5 to 10% of the measured value28]. In some of the experiments with low viscosity (high Ka) fluids, the amplitudes were found to be not saturated even at 7 m from inlet. In all of our experiments, the waves were found to be three-dimensional (or there is no azimuthal symmetry) but the wavelength in the azimuthal direction was higher than that in the streamwise direction, and both of these length scales are much higher than the wave amplitudes. Thus, the data reported here on amplitudes and wave velocities (celerities) are for three-dimensional naturally excited solitary waves with long wavelength in both the streamwise and the azimuthal directions. 13259

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Summary of Overall Experimental Fluid Properties and Flow Ranges Re (range)

We (range)

(4% (unitless)

(7% (unitless)

3883

86.6
809

0.80
33.2

3699

92.0
3127

0.08
28.8

1.23

2506

65.1
2461

0.08
34.6

density

surface tns.

viscosity

glycerin (vol %)

(.4% (g/cm3)

(.8% (dyn/cm)

(mPa s) (.8% (cP)

0.0

0.99

72.0

0.90

0.0

1.00

72.0

0.93

10.0

1.03

70.3

Ka (.1.6% (unitless)

20.0

1.06

69.1

1.69

1635

54.3
1849

0.09
30.4

30.0

1.09

68.2

2.63

901

57.2
1261

0.09
15.4

40.0

1.11

67.6

3.79

564

32.5
881

0.10
24.4

50.0 55.0

1.13 1.13

67.2 67.0

5.50 7.94

337 206

33.9
639 20.1
430

0.10
13.7 0.12
20.1

60.0

1.15

66.8

11.3

65.0

1.16

66.5

16.2

129 79.6

23.5
294

0.15
9.72

12.0
186

0.19
18.3

70.0

1.17

66.2

22.8

50.6

12.0
125

0.23
11.7

75.0

1.18

65.8

28.4

37.6

6.05
90.8

0.29
27.5

78.0

1.19

65.6

33.8

29.8

10.5
39.1

0.94
8.45

80.0

1.19

65.3

48.2

18.5

4.65
42.0

0.53
20.4

83.0 85.0

1.20 1.21

65.0 64.8

72.8 93.0

10.6 7.65

4.09
16.6 2.07
51.3

1.42
14.5 0.16
34.0

87.0

1.21

64.5

104

6.58

2.73
12.5

1.42
17.8

90.0

1.22

64.0

172

3.35

2.23
10.9

0.89
12.6

92.0

1.23

63.7

248

2.05

1.17
8.75

0.81
22.7

Before presenting the experimental results, we note that there are three independent dimensionless groups that characterize the vertically falling film on the inside wall of a cylindrical pipe (under the assumption that the gas phase is stationary). These are the Reynolds number Re ¼

4Q 4Γ 4u0 h0 ¼ ¼ N N πDν ν ν

ð1Þ

σ Ka ¼
1=3 μ νg

ð2Þ

the Kapitza number

and the ratio of flat (Nusselt) film thickness to the pipe radius ξ¼

h0N R

ð3Þ

Another important dimensionless group is the Weber number, We, which is not independent but can be used in place of the Reynolds number. For the case of the one-dimensional flat film in the interior of a tube, it may be shown that
3 4g h0N Re ¼ f ðξÞ ð4Þ 3ν2    3 1 2 4 4 f ðξÞ ¼ 1
4ð1
ξÞ + 3ð1
ξÞ + 4ð1
ξÞ ln 1
ξ 16ξ3 ð5Þ We ¼ ¼

σ Fðu0N Þ2 h0N

¼

31=3 45=3 Ka Re
5=3

14:53Ka Re
5=3

f ðξÞ1=3

ð6Þ f ðξÞ1=3 Here, Γ is the flow rate per unit width, μ(ν) is the viscosity (kinematic viscosity), σ is the surface tension, and g is the

acceleration due to gravity. The flat (Nusselt) film thickness and average velocity of the film are denoted by hN0 and uN0 , respectively. The function f(ξ) accounting for the pipe wall curvature varies from unity (at ξ = 0) to 3/16 (at ξ = 1). For ξ , 1, f(ξ) ≈ 1
ξ, which implies that for very small ξ, say ξ < 0.01, we can ignore the effect of pipe wall curvature and reduce the number of dimensionless groups to two, namely, We (or Re) and Ka. This limiting case corresponds to that of a falling film on a vertical plate. In all of our experiments, ξ < 0.1 and hence the error in neglecting the curvature effect is of the same order of magnitude as the error in the measurement of the film thickness. In a few cases, where the maximum wave amplitudes are about 20% of the pipe radius, the curvature effect is not negligible, and the data reported is not corrected for this effect . Tables 1 and 2 give a summary of experimental conditions for which the results are reported here. Table 1 identifies the range of fluid properties and flow rates of all the experiments, and Table 2 summarizes the experimental conditions for four different fluids that are of the main focus of this work. The fifth column of Table 1 gives the Kapitza number of the fluid and last two columns give the range of Reynolds and Weber numbers for which the wave amplitudes were measured for that fluid. It should be noted that the Reynolds number range for each series of experiments with a fixed viscosity (or Kapitza number) is different. The lowest Reynolds number in the experiments was determined by the minimum volumetric flow when fluid can wet all the inner tube surface and form a film. Therefore, for low viscosity fluids, accurate and meaningful experimental data could be obtained (with the conductance probe method) only for higher Reynolds numbers where the film was at least 0.2 mm thick. For higher viscosity fluids, the upper Reynolds number limit for our experiments was determined by the pump capacity. The higher the viscosity, the lower the volumetric flow rate that could be reached on the particular pump and hence the lower the maximum Reynolds number achievable in our experiments. 13260

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

Table 2. Dimensionless Groups, Nusselt Film Thickness, and Velocity for Selected Experiments calculated

calculated

v (0.6%

Γ (4%

Ka (1.6%

Re (4%

We (7%

exp #

glycerin (vol %)

(mm2/s) (cSt)

(cm2/s)

(unitless)

(unitless)

(unitless)

1

87.0

85.7

0.59

6.58

2.73

17.8

1.16

2

87.0

85.7

0.72

6.58

3.37

12.7

1.23

5.82

3

87.0

85.7

0.89

6.58

4.16

8.88

1.33

6.72

4

87.0

85.7

1.13

6.58

5.26

5.96

1.44

7.87

5

87.0

85.7

1.43

6.58

6.66

4.04

1.55

6

87.0

85.7

2.10

6.58

9.86

2.13

1.76

11.9

7 8

87.0 65.0

85.7 14.0

2.67 0.42

6.58 79.6

12.5 12.0

1.42 18.3

1.91 0.56

14.0 7.45 14.5

hN0 (1.5% (mm)

uN0 (3.1% (cm/s) 5.09

9.20

9

65.0

14.0

1.14

79.6

32.8

3.47

0.79

10

65.0

14.0

2.04

79.6

58.2

1.32

0.96

21.3

11

65.0

14.0

4.54

79.6

131

0.35

1.24

36.4

12

65.0

14.0

6.55

79.6

186

0.19

1.41

46.4

13

40.0

3.41

0.28

554

14

40.0

3.41

0.96

554

114

3.04

0.46

20.8

15 16

40.0 40.0

3.41 3.41

1.89 4.36

554 554

223 513

0.99 0.25

0.58 0.77

32.5 56.8

17

40.0

3.41

7.56

554

881

18

0.0

0.93

0.21

3699

19

0.0

0.93

0.46

3699

196

20

0.0

0.93

0.77

3699

333

21

0.0

0.93

1.48

3699

22

0.0

0.93

2.56

23 24

0.0 0.0

0.93 0.93

4.19 7.30

32.5

24.4

9.02

0.93

81.7

0.18

11.7

8.09

0.24

19.4

3.36

0.28

27.6

635

1.15

0.35

42.5

3699

977

0.56

0.40

56.6

3699 3699

1806 3127

0.20 0.08

0.49 0.59

85.2 123.1

92.0

0.10

0.31

28.8

A. Experimental Maximum Wave Amplitude-Flow Rate Relations. The Nusselt’s flat film solution is used to scale and

present the experimental data in dimensionless form. We neglect curvature effect and analyze the data by assuming that the film is on a flat wall and is characterized by the two independent dimensionless groups, namely, the Kapitza and the Weber numbers. We take the Nusselt film average velocity uN0 and thickness hN0 as the characteristic velocity and length scales. [Remark: The thickness and mean velocity of the Nusselt flat film given in Table 1 are calculated from the flow rate and fluid properties using hN0 = (3νΓ/g)1/3and uN0 = (gΓ2/3ν)1/3, respectively]. Figure 1 shows the maximum measured normalized wave amplitudes (Hmax = (h0 max)/(hN0 )
1) and the normalized standard deviation (hSD) or film roughness, where hSD ¼

h0SD

h0SD h0N

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N
u 1 X 2 ¼t h0i
hh0 i ; ðN
1Þ i ¼ 1

ð7Þ

hh 0 i ¼

N 1X h0 ð8Þ N i¼1 i

as a function of the reciprocal of the Weber number for the experiments with low Ka fluids. The maximum wave amplitude h0 max was determined by taking the average of five largest wave amplitudes that are within 10% of each other in the recorded wave profile over 200 s. [If the peak height of a wave that exists on the film is denoted by h0 peak, then h0 max is taken as the average of

Figure 1. Plots of normalized maximum average amplitude (Hmax) and standard deviation (hSD), as a function of 1/We for low-Ka fluids.

the five highest values of h0 peak that are within 10% in the recorded trace]. The following observations may be made from 13261

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

Figure 2. Plots of normalized maximum average amplitude (Hmax) and standard deviation (hSD), as a function of 1/We for high-Ka fluids.

the data shown in Figure 1: (a) the maximum normalized dimensionless wave amplitude Hmax increases monotonically with 1/We and saturates at values between 1 and 2 (or hmax values between 2 and 3) (b) For small wave amplitudes or large Weber numbers (1/We < 0.1), Hmax increases as R/We, where the constant R is about 7, (c) the film roughness or the normalized standard deviation (hSD) also increases monotonically with 1/We with the saturation occurring at around 0.27, and (d) The dependence of Hmax and hSD on the Kapitza numbers appears to be weak, that is, the maximum wave amplitude and film roughness appear to depend only on the Weber number and not on the Kapitza number, at least in the range of the experiments (2 < Ka < 129). [Remark: Observation (b) implies that for a given fluid (fixed Ka), Hmax increases with the 5/3rd power of the flow rate, Γ]. The scatter in the vertical direction in the data shown in Figure 1 may be due to several reasons: (i) the variation of the Kapitza number, which varies from a low of 2 to a high of about 129, and (ii) the combined experimental error in the conductance probe technique, curvature correction, and in the measurement of other variables. For the sake of clarity of the experimental points, we have not shown the error bars in Figures 1 and 2. As stated earlier, the error in the wave amplitudes is estimated to be about 10%, except in the region of very small wave amplitudes (when the physical film thicknesses 1.3) for any fluid at a fixed flow rate and that (b) the slope of the Ce
hpeak line decreases with increasing flow rate (or decreasing We) and increasing Ka (decreasing viscosity). To summarize, our experimental data presented in Figures 1 to 6 lead to several important observations. The first of these is that

the behavior of the viscous wavy film is different from that of less viscous fluids. Equivalently, there appears to be a cutoff value for the Kapitza number below which wave amplitudes are limited to about three times the flat film thickness and the amplitude is insensitive to the fluid properties (and depending mainly on the flow rate or the Weber number). In contrast, for less viscous fluids such as a water, the wave amplitudes can be as high as ten times the flat film thickness and sensitive to fluid properties. On the basis of prior work that examined the flow patterns underneath the large amplitude waves,25 it may be conjectured that the existence of large amplitude waves in less viscous (or high surface tension) fluids is due to the occurrence of up-flows near the wall region, that is, the clockwise rotating vortices in the front capillary regions support the large amplitude waves. In the next section, we provide a corroborative theory for this hypothesis. The second experimental observation of a nearly linear relation between wave velocity and amplitude is already established in the literature but only for solitary waves generated using periodic inlet forcing and small amplitude waves.4,7,8 In this work, we have shown that this relation is also valid for naturally excited waves and over a wide range of flow rates and fluid properties, that is, the celerity-amplitude relation is independent of the mechanism of generation of these finite amplitude solitary waves. This result was already corroborated in our prior theoretical and computational studies on extended domains.27 In the next section, we discuss a new method for obtaining the celerity-amplitude relation using the h
q model with periodic boundary conditions (BCs) and show that the linear relationship between celerity and amplitude is only valid for finite amplitude solitary waves. 13263

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

The third main experimental observation is the relationship between the normalized deviation wave amplitude (Hmax) and (1/We). Specifically, our experiments with high viscosity fluids (Figure 1) show that for small amplitude waves (Hmax < 1), we have Hmax = R/We, where the constant R is about 7, and this relation appears to be independent of the Kaptiza number. A similar relation was presented by Alekseenko et al.4 for the case of high Ka fluids. On the basis of these observations, it may be conjectured that for the case of very thin films (corresponding to large values of We), the wave amplitude depends only on the Weber number and is independent of the Kapitza number. In the next section, we present local bifurcation results corroborating this conjecture.

III. ANALYSIS OF VERTICALLY FALLING FILMS USING A LOW-DIMENSIONAL MODEL The solution of the full N
S equations gives information on the wavy film behavior, but the presence of unbounded nonlinearities in the BCs and fast changing interface shape makes the computational task difficult. While it is certainly feasible with the present day computers, few numerical results are available in the literature for the full 3-D film flow equations for realistic values of the physical parameters (e.g., when hpeak exceeds 2). As reviewed in the introduction, the case of 2-D film flow has been solved by several authors.10,22
25 However, even in this 2-D case, the domain length (in the streamwise direction) over which the full N
S equations are solved is typically less than 1000 film thicknesses. This should be compared with the length of the test section in our experiments (7 m), which is between 4000 and 35,000 film thicknesses, depending on the fluid properties and flow rate. As demonstrated in our earlier work,27 the use of short domain lengths in the numerical solution may not capture the dynamics of slow growing long waves (low forcing frequencies) which have the largest amplitudes at saturation. In fact, to our knowledge, at present, there are no simulations of either 2-D or 3-D waves of large amplitudes (e.g., hpeak > 3). In contrast, as shown in Figures 2, 5 and 6, experimental wave amplitudes can be as high as ten times the Nusselt film thickness. Simplified low-dimensional models are useful in reducing the effort of integrating in two (or three) spatial dimensions by eliminating the transverse spatial coordinate. As reviewed in the introduction, earlier theoretical work has shown that integral boundary layer models having two independent degrees of freedom, namely, the local film thickness, h(x,t), and the local flow rate, q(x,t), can predict the qualitative features of the wavy film in the viscocapillary regime (We > 1).19,21 Here, we use one of the h
q models developed in our earlier work to explain our experimental results. Defining the dimensionless local flow rate q(x,t) as Z hðx, tÞ uðx, y, tÞ dy ð9Þ qðx, tÞ ¼ 0

the model in dimensionless form (for vertically falling film) is given by19 h t + qx ¼ 0 ð10Þ 2 q 18 q2 hx 34 qqx Re qt ¼ 4h
4 2 + Re 2
Re 5 h 35 35 h h 2 339 hx qx 9 2 373 qhx 41 2
hhx +
h hxx
70 h 7 70 h2 84 2837 qhxx 449 1 qxx + ReWehhxxx +
420 h 70 3

ð11Þ

[Here, the streamwise coordinate and time are scaled using hN0 and hN0 /uN0 and the subscript t (or x) denotes partial derivative. With this scaling, h = q = 1 corresponds to the Nusselt film]. The accuracy of this model is not known except for the fact that it becomes more accurate for low values of Ka and high values of We. Comparison of the numerical results of this h
q model for a few sets of parameter values corresponding to those used by Malamataris et al.25 and Nosoko and Miyara10 in the solution of the full 2-D film flow equations showed excellent agreement, including the prediction of the regions in which up-flows occur. In the h
q model, the local flow rate q(x,t) becomes negative in the up-flow regions, and this usually occurs at the minimum of the capillary ripples preceding a big wave. It should also be pointed out that the phenomenon of flow reversal in the front regions of large amplitude waves was predicted by using the h
q models and validated recently computationally as well as experimentally.11
13,18,25,27 Thus, the above h
q model is expected to be accurate in the visco-capillary regime (We > 1), which is the focus of the numerical and local bifurcation studies presented below. A. Flow Reversal and Cellular Patterns in Vertically Falling Films. As can be expected, the behavior of the wavy film is least complex for highly viscous fluids (small values of Ka) and low flow rates (small values of 1/We) leading to thin films in which inertial effects may be significant but not dominant). As discussed in the next section, in this limit, the h
q model may be analyzed using local bifurcation theory to obtain quantitative relationships between wave amplitude, wave celerity, and the two dimensionless numbers (Ka,We) characterizing the film. As stated above, Hmax is limited to about 2 for low Ka fluids while it can be as high as 8 or 9 for high Ka fluids. The physical reason for the existence of a critical Kapitza number for flow reversal (and cellular pattern formation in the film leading to large amplitude waves) may be explained by noting that Ka may be interpreted as a ratio of pattern forming capillary forces (σ/hN0 ) to pattern resisting viscous forces (μ(νg)1/3/hN0 ). When Ka is small, the viscous forces dominate over the capillary forces and there can be no flow reversal, while the converse is true for high Ka fluids. Thus, as in other systems in which cellular patterns are formed (e.g., Rayleigh-Benard or Marangoni convection), a critical value of the ratio (or Ka value) is needed to form cellular patterns (flow reversal) leading to large amplitude waves. The above intuitive explanation is supported by experimental data as well as computational results using the h
q model (and numerical calculations using the Navier
Stokes equations). For example, there are no capillary ripples in front of the waves when Ka is sufficiently small. However, the number of these capillary ripples (and hence the number of regions of flow reversal) increases rapidly with increasing Ka and decreasing We (for high Ka). Figure 7 shows three stationary wave profiles computed using the h
q model [Remark: A stationary wave is one that has a constant shape and speed]. Figure 8 shows the corresponding flow rate (q(x,t)) profiles. The first case corresponds to a viscous fluid in one of our experimental conditions with Ka = 6.5 and We = 8.9. Here, the wave amplitude is limited to about hpeak = 2.5 (and has a celerity of 5.74), the substrate thickness (h0) is 0.91 and there are no capillary ripples in the front region of the wave (except for a very small dip with hmin = 0.89). The maximum flow rate at the wave peak is about 9.7 while the minimum flow rate is 0.64 compared to a substrate flow rate (q0 = h30) of 0.75. The second case with Ka = 210 and We = 21.38 corresponds to recent experimental results of Alekseenko et al.15 In this case, the maximum wave amplitude is hpeak = 2.69, and the wave has 13264

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

Figure 7. Plots showing variation of the wave profile with the Kapitza and Weber numbers.

about four capillary ripples. The flow rate at the first minimum is negative, while it is almost zero at the second, indicating the presence of counter rotating cells. [Remark: For a detailed description of the streamlines/flow patterns in the up-flow regions of the waves, we refer to recent articles by Malamataris and Balakotaiah25 and Dietze et al.13]. The third case with Ka = 3449.3 and We = 33.38 corresponds to the famous experiments of Kapitza and Kapitza with water.1 In this case, the peak wave amplitude is about 3.06 (celerity of 5.1) and there are about 17 capillary ripples in the front with six up-flow regions. The flow rate at the first minimum (qmin =
0.85) is higher in magnitude than that in the substrate (q0 = 0.47). It should also be pointed out that the wavelength in this case is about 300 film thicknesses, which is comparable to the domain size used in most numerical simulations of the two-dimensional N
S equations. In addition, the time required for the wave to reach a steady profile is of the order of 700 dimensionless time units for an initial perturbation of magnitude 30% in the flow rate. Again, this time is higher than that used in most direct simulations of the N
S equations and for shorter times, unsteady and oscillating wave profiles are obtained.25 For less viscous fluids such as water, we have computed as many as 20 regions of up-flow in the front region of a large wave having a wavelength of about 500 and an amplitude of 4.5 for We values of 30 or higher. [Remark: All the numerical calculations using the h
q model were done using Mathematica 8 software. Eighth order finite difference method with fixed but variable mesh size (depending on Ka value) was used with infinite precision].

ARTICLE

Figure 8. Plots of local flow rates for the wave profiles shown in Figure 7.

We have used the h
q model to determine the approximate boundary in the (Ka, We) plane at which flow reversal occurs. As stated earier, this boundary is important as it demarcates the region in which large amplitude waves can be observed. To determine this boundary, for each Ka value, the Weber number is decreased from a large value until a large amplitude wave with a capillary wave in the front, in which the local flow rate at the minimum is close to zero, is obtained. The sensitivity of this minimum flow rate with the domain length (and hence the wave amplitude) is used to determine the largest possible wave (whose substrate thickness is very close to that of Nusselt film), and this value is taken as the critical Weber number. Figure 9 shows the approximate boundary determined using this procedure. The following observations can be made from the results shown in Figure 9: (i) For Ka values below 15, there is no critical Weber number in the region We > 1 at which flow reversal occurs. Thus, there appears to be a critical Kapitza number (around 15) below which there is no flow reversal. (ii) For high Ka values (Ka > 100), the critical We at which up-flows appear, increases approximately as (Ka)1/2. While this increase can be expected intuitively, we do not have an explanation for this (approximate) square root dependence. (iii) For low viscosity fluids such as water with a Ka value of 3700, the critical Weber number at which up-flows occur is around 450, implying that flow reversal and hence large amplitude waves can be obtained for We.1 [Remark: Our definition of Re = 4Γ/ν includes a factor 4 for consistency with circular pipes, in which most experiments on falling films are carried out. However, most literature definitions of Re exclude this factor]. 13265

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

This relation suggests that the local flow rate in a wave moving with a steady speed is related to the local film thickness linearly with the slope equal to the celerity (Ce) of the wave. This relation is useful in determining the celerity of the different waves that exist on a film. [Also, if the relation between q(x,t) and h(x,t) is not linear, it implies that the wave is not stationary and still evolving]. The constant of integration determines the different classes of waves. Imposing that the flat film h(z) = h0, q(z) = q0 = h30 be a solution gives c0 ¼ h30
Ceh0

ð13Þ

In the literature, the flat film (or substrate) is taken as the Nusselt film or q0 = h0 = 1 and hence the constant c0 = 1
Ce. However, as we show in the next section, this assumption restricts the class of solutions that can exist. From a physical point of view, the substrate thickness on which waves exist could be less than that of the Nusselt film, as the waves carry most of the fluid mass and momentum. In the steady traveling wave coordinate, the averaged momentum balance, eq 11 leads to the following third order nonlinear ordinary differential equation describing the wave profile in the z-coordinate (with q = Ceh + h30
Ceh0): hzzz ¼ Figure 9. Approximate boundary of the region in which upflows can occur in a vertically falling film.

The results in Figure 9 may used to provide a qualitative explanation of the maximum wave amplitude data shown in Figures 1 and 2. For low Ka (viscous) fluids, the wave amplitude is mainly determined by the balance between inertial/gravity and viscous forces, while the capillary forces play a minor role. However, for high Ka (less viscous) fluids, capillary effects can be dominant and create regions of up-flows. As stated earlier, these up-flows (with counter-rotating vortices) can support large amplitude waves and were confirmed in recent experiments11
13 but only when the wave amplitudes were relatively small (about 2). The h
q models as well as recent calculations with the 2-D Navier
Stokes equations25 predict that for high Ka fluids such as water, there can be many such up-flow regions (>20) whose width can oscillate in time, leading to cellular patterns and breathing (expanding and contracting) vortices in the film. Again, these predictions corroborate the experimental observation that for low viscosity fluids, the solitary wave profiles (and their celerities) are not stationary even at very large distances from the inlet. B. Amplitude-Celerity Relations for Solitary Wave Families Using Local Bifurcation Analysis. In this section, we explain the experimentally observed wave amplitude data at large Weber numbers (low flow rates) using the results of local bifurcation theory. We extend the earlier local bifurcation results on the amplitude and celerity of fully developed waves for We.1. In the traveling wave coordinate z = x
Cet, where Ce is the (real celerity of the wave), the evolution equations of the h
q model can be reduced to a set of ordinary differential equations (ODEs) if we assume that a steady waveform is obtained (i.e., Ce is constant and the wave shape depends only on z, h(x,t) = h(z)). These nonlinear ODEs can be analyzed using dynamical systems theory. Integration of the continuity eq 10 in the steady traveling coordinate leads to q
Ceh ¼ c0

ð12Þ

3 4q 2 ½
4h + 2
ReCeqz hReWe h 5 18Re q2 hz 34Re qqz 339 hz qz 9 2 + + hhz
+ 35 h2 35 h 70 h 7 373 qh2z 41 2 2837 qhzz 449 qzz 

+ h hzz + 70 h2 84 420 h 70

ð14Þ

This third order nonlinear differential equation can be written as three first order equations, and the established local and global bifurcation techniques can be applied to analyze the various bifurcations in the parameter space (h0,Ce,Ka,We). Since the details of the local bifurcation analysis are already presented elsewhere,19,27,29
31 we report here only the new results. Also, we do not consider here the global bifurcations (which can be very complex for high Ka and/or low We values) and refer to these earlier works. To explain our local bifurcation and numerical results, we review here very briefly the properties of small amplitude long waves as well as the different bifurcations that can occur from the flat film. First, we note that the leading order (long wave) approximation of eq 11 leads to the relation q = h3, which upon substitution into the continuity equation leads to the hyperbolic equation ∂h ∂h + 3h2 ¼ 0 ∂t ∂x

ð15Þ

Equation 15 leads to the following observations: (i) the celerity of small amplitude long waves is Ce = 3h2, (ii) infinitesimal amplitude long waves on the Nusselt film (h = 1) have a (linear) celerity of three, (iii) depression (trough) waves move slower than elevation (hump) waves, (iv) when the substrate thickness is less than that of the Nusselt film, the linear celerity can be smaller than three. The fixed points (steady-states) of eq 14 are given by hss1 ¼ h0 hss2 ¼


h0 + 2

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 Ce
h20 4

ð16Þ

Figure 10 shows a schematic diagram of these two branches which meet at point D corresponding to the linear celerity, 13266

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

extended domains, and these largest waves (having Nusselt film as the substrate) were denoted as the Γ2 wave by Meza and Balakotaiah.27 [Remark: These waves were referred to as the “tsunami waves” by Meza and Balakotaiah as these are the largest waves that can exist on the film for any fixed Ka and We. However, we do not use that nomenclature here and simply refer to them as the Γ2 waves]. Weakly nonlinear analysis of the h
q model near the double zero point ((1/We) f 0 and Ce f 3h20) results in the following expressions for the peak height and celerity of the γ2 family: hpeak
h0 ¼ Figure 10. Schematic diagram showing the different wave families and bifurcations in the traveling wave coordinate.

Ce = 3h20. While these steady-state branches are independent of Ka and We, the various bifurcations that occur on these branches depend on Ka and We. Figure 10 shows the bifurcations that occur for a low Ka fluid and We.1. A similar diagram for a high Ka fluid or We < 1 can be extremely complex and contain several period doublings (corresponding to multipeak waves) as well as chaotic windows. However, in the limit of large Weber numbers ((1/We) f 0), the various bifurcations corresponding to different wave families can be analyzed.19,27,29
31 The slow moving γ1 wave family is characterized by deeper troughs than crests and celerities bounded by the Hopf bifurcation point A and the homoclinic bifurcation point B (Figure 10). The homoclinic (or infinite period bifurcation) point correpsonds to a solitary wave of the γ1 wave family. Local bifurcation analysis leads to the following expressions relating the amplitude, and celerity of the γ1 solitary wave (near the double zero point corresponding to (1/We) f 0 and Ce f 3h20): hpeak
h0 ¼ and

2:52h60 ; We

h0
hmin ¼2 hpeak
h0


 Ce ¼ 3h20
6h0 hpeak
h0

ð17Þ

ð18Þ

As stated above, these relations are generalizations of those presented in earlier studies which assumed h0 = 1. By allowing h0 to be less than unity, the above relation can describe solitary waves on the “substrate” whose thickness is less than the Nusselt film. We note that the above celerity-amplitude relation is linear for each fixed h0. Even though there is only one solitary wave of the γ1 family for each fixed h0, the celerity-amplitude relation is also applicable for all periodic waves whose celerities vary between the Hopf bifurcation point A (corresponding to small amplitude periodic waves) and the homoclinic point B (corresponding to solitary waves).29,30 When h0 is varied, the above relations describe the solitary waves of the γ1 family having different substrate thicknesses. When eqs 17 and 18 are plotted by taking h0 as a parameter, the relation between peak amplitude and celerity is very close to linear. An important point to note about the γ1 family of waves is that the celerity decreases as the amplitude increases. The second class of waves are those moving with higher celerity and referred to as the γ2 family, characterized by higher crests and shallow troughs. When h0 = 1 these waves are the largest waves that can be observed in time simulations over

and

7:56h60 ; We

hmin ¼ h0


 Ce ¼ 3h20 + 2h0 hpeak
h0

ð19Þ

ð20Þ

Again, for any fixed h0, the celerity-peak amplitude relation is linear and also applicable to all periodic waves of the γ2 family whose celerities vary between the Hopf bifurcation point C and the homoclinic point E (Figure 10). An important point to note about the γ2 family of waves is that they have high celerities and high amplitudes, and the celerity increases as the amplitude increases. Again, it should be pointed out that the largest of these waves, the Γ2 wave has the highest celerity, amplitude, and wavelength, and represents the largest wave that could exist for those specific Ka and We values (because of nonlinear saturation effects). The third family of waves (not shown in Figure 10) emanating from the flat film are the long waves corresponding to wavenumber R = 0 branch of the neutral stability curve. These are difficult to study experimentally as they correspond to very low forcing frequencies and require a very long spatial domain or very long times to develop. The bifurcations that occur from the branch R = 0 and Ce = 3 can be studied in the traveling wave coordinate.19,27 These solitary waves (obtained in the limit of low forcing frequencies) have a substrate that is slightly lower than the Nusselt film thickness and travel with speed that is close to but exceeds 3. We shall refer to these waves as the γ*2 family as they break-up to become the γ2 (and Γ2) waves.27 Weakly nonlinear analysis of the h
q model near the double zero point (1/We) f 0 and Ce f 3 results in the following expressions for hpeak, the peak height, and Ce, the celerity of the solitary γ*2 wave: hpeak
1 ¼

5:12 ; We

hpeak
1 ¼2 1
hmin

ð21Þ

and Ce
3 ¼

2 3
hpeak
1 4

ð22Þ

An interesting observation is that for all three wave families, the wave peak height increases linearly with the inverse of the Weber number (for (1/We) f 0) and independent of the Kapitza number (confirming the experimental results). However, for the γ*2 family, the relation between celerity and peak height is not linear or close to linear but is quadratic. Thus, these long waves have celerities very close to the linear value of 3 even when their amplitudes are not small. For example, a wave with a hpeak of 1.5 has a celerity of only 3.1875. We believe that this is one reason for the nonexistence of correlations (and scatter of experimental data) between wave amplitude and celerity for small amplitude waves. 13267

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research It is important to note again that the above local bifurcation results have only a limited region of validity, namely, for wave celerities near the linear limit of 3h20 (for γ1 and γ2 families) or 3 (for γ*2 family) and large Weber numbers or small peak amplitudes. These local bifurcation results can be extended by numerical analysis of the global bifurcations to determine the celerity-amplitude relations for finite values of We and Ka using the full h
q model in the traveling wave coordinate. We do not pursue this here. Instead, we determine the celerity-amplitude relations for finite amplitude waves using time simulations of the h
q model. C. Wave Evolution on Extended Domains. In this section, we review briefly the connection between the results of transient simulation of the wavy film on extended domains and the study of individual waves with periodic BCs. In our earlier work,19,27 we have presented some numerical simulations using the h
q model on extended domains with monochromatic pulsing experiments by taking the inlet conditions:
 ð23Þ qð0, t Þ ¼ 1 + ε sin 2πft , hð0, t Þ ¼ qð0, t Þ1=3 and initial conditions, h(x,0) = q(x,0) = 1, corresponding to the Nusselt flat film. The other BCs used were, qx(0,t) = 0, hx(0,t) = 0, and hx(L,t) = 0, where L is the length of the computational domain. Here, f is the dimensionless inlet forcing frequency and is related to the wavenumber (R) by f = RCe/2π. The main observations from these simulations over extended domains (some as large as 50,000 film thicknesses) may be summarized as follows: (i) For high inlet forcing frequencies (slightly below the upper neutral/cutoff value), the γ1 family of waves emanate, grow, and then transform to the γ2 family and then eventually to Γ2 waves. (ii) For very low inlet forcing frequencies (near the lower cutoff value of zero), the long waves of the γ*2 family emerge, they grow (steepen), and give rise to packets of γ2 waves, which disperse and grow eventually to Γ2 waves. (iii) For finite/ short domain lengths, the waves may exit the domain and only either the saturated or unsaturated (still growing) γ1, γ2, and γ*2 waves (depending on the forcing frequency, amplitude, and domain length) may be observed. (iv) For long domain lengths, waves of different families coexist on the film. The relation between wave peak amplitude and wave celerity for all the finite amplitude waves that exist on the film is very close to linear, and this relation is independent of the inlet forcing frequency, and (v) the above observations are also valid for multifrequency inlet forcing of the form qð0, tÞ ¼ 1 +

n X


 εi sin 2πfi t ,

hð0, tÞ ¼ qð0, t Þ1=3 ð24Þ

i¼1

The above observations imply that the celerity-amplitude relations determined using a single inlet forcing frequency are also applicable to naturally excited waves. The main advantage of the simulations over extended domain is that they show the spatial evolution of waves as a function of position and inlet frequency, which can be compared with monochromatic forcing experiments.4,8 The main disadvantage is the long computational time needed (even with the h
q models). Also, if the wave celerity-peak amplitude relation (which is independent of inlet forcing) is the main interest, this can be determined by using periodic BCs. This is the main computational focus of this work and is discussed in the next section. D. Determination of Celerity-Amplitude Relations Using Periodic BCs. Periodic BCs are commonly used to study the

ARTICLE

bifurcation behavior of waves. While the wave evolution on a vertically falling film is not periodic, the use of periodic BCs in a model may be justified if the focus is on the waves of a given wavelength and how they evolve in time or downstream distance. Further, using periodic BCs, we can investigate systematically the relationship between wave amplitude and celerity, and this relationship is expected to be valid in the asymptotic region in which the statistics of the waves on the film reach a pseudo steady-state. The ultimate justification of this assumption is the comparison of the predicted celerity-amplitude relations with the experimental data and simulations using extended domains with inlet forcing. The h
q model was solved on a domain of size L using the BCs, that is, ∂h ∂h ð0, tÞ ¼ ðL, tÞ; ∂x ∂x 2 2 ∂h ∂h ð0, tÞ ¼ 2 ðL, tÞ ∂x2 ∂x

hð0, tÞ ¼ hðL, tÞ;

qð0, tÞ ¼ qðL, tÞ;

L¼

2π R

where R is a wavenumber between 0 and the neutral value Rc. The initial flow rate in the film and the surface film profile was taken to be of the form  
1=3 2πnx ð25Þ qðx, 0Þ ¼ q0 + ε sin ; hðx, 0Þ ¼ qðx, 0Þ L where ε is the perturbation amplitude, q0 is the flow rate in the substrate film, and n is the number of waves in the initial perturbation. In the numerical simulations, q0 = 1 and ε is varied between 0.01 and 0.95 while n is taken as unity in most cases. From the linear stability analysis of the h
q model with h0 = q0 = 1, the neutral wavenumber and celerity may be found using the relations: Ce
3 ¼ R2c We ¼

R2c ð507
449CeÞ 280

6 ð7Ce2
17Ce + 9Þ 35

ð26Þ

For example, for We = 10, the above relations give Rc = 0.440 (or critical wavelength as Lc = 14.3 and neutral celerity Ce = 2.56). For any fixed domain size L > Lc, the size of the initial perturbation (ε) was varied to obtain a finite amplitude wave, and its celerity was determined by the slope of the q
h plot using the relation q
Ceh = constant. The domain length and the perturbation amplitude were varied to obtain waves belonging to different families. The results are presented below for a few cases with fixed Kapitza and Weber numbers. Case 1: Ka = 6.5 and We = 10. These parameter values are selected to compare the celerityamplitude relations obtained using periodic BCs with the same data obtained by numerical simulations on extended domains and global bifurcation results of the h
q model presented in earlier work.26 Figure 11 shows the wave celerity-amplitude data obtained by simulation of single waves with periodic BCs with ε = 0.1. The open circles indicate the different waves obtained as the domain size is increased from the critical value of 14.3 to 200. For small domain sizes, the waves obtained are of γ1 type, and as 13268

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

Figure 11. Computed celerity-amplitude points for different wave families.

the domain size is increased, we get waves of the γ2 family. The largest of these is obtained for L = 200 and has an amplitude of 2.54 and celerity of 6.05. Increasing the domain size beyond 200 led to two or more smaller γ2 waves. For example, L = 400 gives two smaller waves of amplitude 2.53 while L = 500 gave three smaller waves. The Γ2 wave with a wavelength of 35, amplitude of 3.11, and celerity of 7.76 was easily obtained in the time simulations over extended domains (with Nusselt substrate as the initial condition) and is more difficult to obtain in the periodic simulations with initial conditions given by eq 25. To obtain waves of larger amplitude, we have to change the initial condition, eq 25, so that q0 exceeds unity. Now, the Γ2 wave is obtained even with L = 100 by taking q0 = 1.7. The physical reason for requiring higher base flow rate is that these large amplitude waves carry most of the mass (and momentum), and hence the local flow rate in these waves can be much higher than that corresponding to that in the Nusselt film. Thus, if we start with an initial condition corresponding to the Nusselt film, we end up with a smaller wave having a substrate much thinner than the Nusselt film as the growing wave sucks up most of the fluid in the film. For very small initial perturbations (ε e 0.05) of the Nusselt film and longer domain sizes we obtained the waves of the γ*2 family whose celerity-amplitude relation is shown by the open squares in Figure 11. In contrast, smaller domain sizes and large initial perturbations (ε g 0.9) lead to waves whose celerities were below those of the γ2 family (shown by open triangles in Figure 11). Thus, by varying the size of the domain (or wavelength) and magnitude of the initial perturbation, different celerity-amplitude combinations in the triangular region denoted by AHB in Figure 11 could be obtained. As stated above, the determination of the celerity-amplitude relation for the entire γ2 family (up to the Γ2 wave denoted by point G in Figure 11) requires varying also the base flow rate. On the basis of our local bifurcation analysis and detailed calculations for this specific case with periodic BCs as well as numerical simulations of wave evolution on extended domains, the following observations can be made: (i) a multitude of celerity-amplitude combinations are possible for small amplitude waves, whose amplitude is smaller than that defined by the homoclinic point H of the γ*2 family; (ii) for large

Figure 12. Comparison of measured and computed celerity-amplitude relationship for Ka = 6.58, We = 8.88 (top), and We = 4.04 (bottom).

amplitude γ2 waves, defined in this case as those with amplitudes between the homoclinic point H and the homoclinic point of the G of the γ2 family, the celerity-amplitude relation appears to be unique and linear (or very close to linear) and the computed celerities fall very close to or on the line joining the homoclinic points H and G; (iii) the same celerity-amplitude relation is obtained from time simulations over extended domains or using periodic BCs; and (iv) the large amplitude waves (whose amplitudes are closer to the Γ2 wave) are difficult to observe either in numerical simulations or experiments as they require special initial or inlet conditions. This last observation is important in the explanation of our (as well as other) experimental results. Case 2: Ka = 6.58 and We = 8.88, 4.04 This case corresponds to a fluid used in one of our experiments. The parameter values corresponding to the higher flow rate (We = 8.88) are very close to case 1 discussed above. Figure 12 compares the computed celerities with experimental values for the two flow rates selected. The experimental and computed data include both the γ1 and the γ2 waves. The following further observations may be made from this comparison: (i) the agreement is good at the higher flow rate and fair at the lower flow rate, and (ii) the computations show that the wave amplitudes can be much higher than observed in the experiments. For example, for We = 4.04, the Γ2 wave has an amplitude of 3.7, but the maximum amplitude observed experimentally is about 2.8. The line shown through the experimental points is the best fit linear regression and has a smaller slope than that predicted by the simulations. As explained above, two possible reasons for not observing higher amplitude waves in the experiments are as follows: (i) they require special inlet forcing or initial conditions and can not be seen on naturally excited films, and (ii) they exist but only during start-up and have already passed the measuring section before our data was taken. 13269

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

Figure 13. Comparison of measured and computed celerity-amplitude relationship for Ka = 79.6 and We = 18.3.

Figure 14. Comparison of measured and computed celerity-amplitude relationship for Ka = 554 and We = 24.2.

Case 3: Ka = 79.6 and We = 18.3. Figure 13 compares the computed celerities with experimental data, and the agreement here is good. Again, the experimental data covers only γ2 waves with amplitudes of 2.2 or smaller while the simulations predict higher amplitude waves exist at the given flow rate. Case 4: Ka = 554 and We = 24.2. Figure 14 compares the measured and computed celerityamplitude data for this case. Once again, we see that the agreement is fair, but the experimental data at this flow rate does not include waves of amplitude exceeding 3 while the numerical simulations predict that the amplitude can be as high as 4. In addition, as in case 2, the slope of the computed line is higher than that given by the best fit regression line of the data. Case 5: Ka = 3699 and We = 28.8. This case corresponds to our experimental data for water at the smallest flow rate. Figure 15 shows a comparison of the predicted and measured celerities. Again, the agreement is fair, and the computed slope is higher than that of the regression line of the experimental data. On the basis of the above comparisons, it may be concluded that the h
q model predicts the qualitative features of the wavy films and has acceptable quantitative accuracy.

IV. SUMMARY AND DISCUSSION One main contribution of this work is the presentation of experimental celerity-amplitude relationships for finite amplitude waves on vertically falling films over a wide range of fluid properties (2 < Ka < 3883) and flow rates (0.08 < We < 35).

ARTICLE

Figure 15. Comparison of measured and computed celerity-amplitude relationship for Ka = 3699 and We = 28.8.

The data presented show some clear trends: (a) For a given fluid (fixed Ka), the wave amplitude increases monotonically with the flow rate (as can be expected intuitively) and reaches an asymptotic value. (b) The behavior of the low Kapitza fluids (Ka < 129) is found to be different from that of the high Kapitza fluids (Ka > 200). In the former case, the wave amplitudes saturated around 3 times the Nusselt film thickness, and the film profile reaches a pseudo steady-state within the experimental domain (short distances/time). In the case of high Kapitza fluids, the wave amplitudes saturate at much higher values of around 10 times the Nusselt film thickness, and waves are dynamic, that is, continuously evolving by merging and splitting though the wave statistics such as peak height appear to reach an asymptotic limit. (c) For a fixed fluid and flow rate (fixed Ka and We), the celerityamplitude relationship is nearly linear for all finite-amplitude waves that exist on the film. The slope of this relationship increases with Ka (decrease in viscosity) and We (decrease in flow rate). (d) The data also shows that the film roughness (standard deviation normalized by the Nusselt film thickness) is also correlated to the inverse of the Weber number and has a weaker dependence on the Kapitza number. A second contribution of this work is the presentation of a corroborative theory on the wavy film dynamics using a lowdimensional h
q model. One main result of the computations is the determination of the approximate boundary in the (Ka,We) plane that demarcates the region of up-flows and cellular patterns in the film. It is reasoned that large wave amplitudes can occur only in these up-flow regions that exist for high Ka fluids and/or at higher flow rates. The model predictions corroborate the experimental observations on the existence of a critical Kapitza number below which there is no flow reversal. However, the predicted value of 15 is much lower than the experimentally observed value of about 120. A second result of our computations is the demonstration of the relationship between wave amplitude and celerity for all the different (evolving and saturated) solitary waves observed on the film. We have also shown that this relationship becomes linear in the limit of large wave amplitudes. We have also used the results from local bifurcation theory to show that the celerity-amplitude relation is close to linear for waves with finite wavelength but becomes quadratic for waves with very long wavelength. The results of the wave amplitudes computed using the h
q model are in good qualitative and fairly good quantitative agreement with the experimental results presented here. We now discuss some limitations and possible extensions of the theoretical/computational results presented here. As pointed 13270

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research out earlier, one limitation of the theory presented here is the restriction to the case of 2-D films while the experimental results presented are for 3-D films. However, the good agreement between predicted and measured wave amplitudes in the viscocapillary regime (We > 1) indicates that the spanwise wavelength of the various wave families is about the same as that in the streamwise direction. This is indeed confirmed by the experimental results of Alekseenko et al.4,13,14 and theoretical analysis of Scheid et al.20 Thus, the results of the 2-D theory such as the maximum wave amplitude being dependent only Ka and We (and not on the domain length or inlet forcing frequency) and the linear amplitude-celerity relationships may be applicable to the case of 3-D waves on the film, but this needs to be validated. A second limitation of the h
q model studied here is that the peak wave amplitude is limited to about 5, while experimental data shows higher peak amplitudes. We believe that such large amplitude waves can only be described by using more local degrees of freedom in the reduced order models. A third limitation of the current theoretical/modeling studies is that they are mostly limited to describing stationary (or steady traveling) wave families. Our calculations and experiments show that for high Ka fluids, the wave behavior is nonstationary even at moderately high Weber numbers. The characterization of the film behavior, both experimentally and theoretically, in the inertial regime for high Ka fluids is still a challenging task.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT The senior author (V.B) would like to dedicate this manuscript to Professor M. S. Ananth, his teacher and first research mentor. This work was partially supported by a grant (NNX10AL37G) from NASA Glenn Research Center. ’ NOMENCLATURE Ce = dimensionless wave velocity (celerity), uw/uN0 g = accleration due to gravity h0 = dimensional film thickness h = dimensionless (normalized) film thickness, h0 /hN0 hpeak = normalized peak height of a wave, h0 peak/hN0 hmax = maximum (normalized) value of hpeak in the observed trace hN0 = Nusselt flat film thickness Hmax = maximum normalized wave amplitude deviation (hmax
1) Ka = Kapitza number, (σ)/(μ(υg)1/3) q(x,t) = local flow rate (dimensionless) Re = Reynolds number, 4Γ/ν uN0 = average Nusselt velocity, g(hN0 )2/3ν uw = dimensional wave velocity We = Weber number, (σ)/(F(uN0 )2hN0 ) Greek Letters

R = dimensionless wavenumber or frequency, (2πhN0 )/(λ) ε = perturbation amplitude(dimensionless) Γ = volumetric flow rate per unit width λ = wavelength ν = kinematic viscosity μ = viscosity

ARTICLE

F = density σ = surface tension

’ REFERENCES (1) Kapitza, P. L.; Kapitza, S. P. Wave flow of thin layers of a viscous fluid. Zhurn. Eksper. Teor. Fiz 1949, 19 (2), 105.English translation in Collected Papers of P. L. Kapitza; D. Ter Haar, Ed.; Pergamon Press: New York, 1965; Vol.II, p 690. (2) Benjamin, T. B. Wave formation in laminar flow down an inclined plane. J. Fluid. Mech. 1957, 2, 554. (3) Yih, C.-S. Stability of liquid flow down an inclined plane. Phys. Fluids 1963, 6 (3), 321. (4) Alekseenko, S. V., Nakoryakov, V. Y. Pokusaev, B. G. Wave flow of liquid films; Begell House: New York, 1994. (5) Chang, H. C.; Demekhin, E. A. Complex wave dynamics on thin films; Elsevier: New York, 2002. (6) Krantz, W. B.; Goren, S. L. Stability of thin liquid films flowing down a plane. Ind. Eng. Chem., Fundam. 1971, 10, 91. (7) Alekseenko, S. V.; Nakoryakov, V. Y.; Pokusaev, B. G. Wave formation on vertical falling liquid films. AIChE J. 1985, 31, 1446. (8) Liu, J.; Gollub, J. P. Solitary wave dynamics of film flows. Phys. Fluids 1994, 6, 1702. (9) Liu, J.; Schneider, J. B.; Gollub, J. P. Three-dimensional instabilities of film flows Phys. Fluids 1995, 7, 55. (10) Nosoko, T.; Miyara, A. The evolution and the subsequent dynamics of waves on a vertically falling liquid film. Phys. Fluids. 2004, 16, 1118. (11) Tihon, J.; Tovchigrechko, V.; Sobolik, V; Wein, O. Electrodiffusion detection of the near wall flow reversal in liquid films at the regions of solitary waves. J. Appl. Electrochem. 2003, 33, 577. (12) Tihon, J.; Serifi, K.; Argyriadi, K.; Bontozoglou, V. Solitary waves on inclined films: their characteristics and the effects on wall shear stress. Exp. Fluids 2006, 41, 79. (13) Dietze, G. F.; Leefken, A.; Kneer, R. Investigation of back flow phenomenon in falling liquid films. J. Fluid Mech. 2008, 595, 435. (14) Alekseenko, S. V.; Antipin, V. A.; Guzanov, V. V.; Kharlamov, S. M.; Markovich, D. M. Three-dimensional solitary waves on falling liquid film at low Reynolds numbers. Phys. Fluids 2005, 17, 12704. (15) Alekseenko, S. V.; Guzanov, V. V.; Markovich, D. M.; Karlamov, S. M. Characteristics of solitary three-dimensional waves on vertically falling films. Tech. Phys. Lett. 2010, 36 (11), 1024. (16) Benney, D. J. Long waves on liquid films. J. Math. Phys. 1966, 45, 150. (17) Chang, H.-C. Wave evolution on a falling film. Ann. Rev. Fluid Mech. 1994, 26, 103. (18) Ruyer-Quil, C.; Manneville, P. Improved modeling of flows down inclined planes. Eur. Phys. J. B 2000, 15, 357. (19) Mudunuri, R. R.; Balakotaiah, V. Solitary waves on thin falling films in the very low forcing frequency limit. AIChE J. 2006, 12, 3995. (20) Shkadov, V. Wave flow regimes of a thin layer of viscous fluid subject to gravity. Izv. Ak. Nauk SSSR Mekh. Zhi. Gaz 1967, 2, 43–51.(English transl. Fluid Dynamics 2; Faraday Press: New York, 1970; p 29). (21) Scheid, B.; Ruyer-Quil, C.; Manneville, P. Wave patterns in film flows: modelling and three-dimensional waves. J. Fluid Mech. 2006, 562, 183. (22) Ramaswamy, B.; Chippada, S.; Joo, S. W. A full-scale numerical study of interfacial instabilities in thin-films. J. Fluid Mech. 1996, 325, 163. (23) Malamataris, N. A.; Vlachogiannis, M.; Bontozoglou, V. Solitary waves on inclined films: Flow structure and binary interactions. Phys. Fluids 2002, 14, 1082. (24) Gao, D.; Morley, N. B.; Dhir, V. Numerical simulation of wavy falling film flow using VOF method. J. Comput. Phys. 2003, 192, 624. (25) Malamataris, N. A.; Balakotaiah, V. Flow Structure under the large amplitude waves of a vertically falling film. AIChE J. 2008, 54, 1725. (26) Massot, C.; Irani, F.; Lightfoot, E. N. Modified description of wave motion in a falling film. AIChE J. 1966, 12, 445. 13271

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272

Industrial & Engineering Chemistry Research

ARTICLE

(27) Meza, C; Balakotaiah, V. Modeling and Experimental Studies of large amplitude waves on vertically falling films. Chem. Eng. Sci. 2007, 63, 4704. (28) Meza, C., Modeling and experimental studies of large waves on free-falling wavy thin films. Doctoral dissertation, University of Houston, Houston, TX, 2007. (29) Nguyen, L. T.; Balakotaiah, V. Modeling and Experimental Studies of wave evolution on free falling viscous films. Phys. Fluids 2000, 12, 2236. (30) Panga, M. K. R.; Mudunuri, R. R.; Balakotaiah, V. Longwavelength equation for vertically falling films. Phys. Rev. E 2005, 71, 036310. (31) Chang, H.-C.; Demekhin, E. A.; Kopelevich, D. I. Nonlinear evolution of waves on a vertically falling film. J. Fuid Mech. 1993, 250, 433.

13272

dx.doi.org/10.1021/ie2005803 |Ind. Eng. Chem. Res. 2011, 50, 13258–13272