Nonlinear waves on liquid film surfaces—I. Flooding in a vertical tube

Chemical En&eering Science, Vol. Printed in Great Britain.

41, No.

NONLINEAR

10. pp. 2463-2476,

1986.

WAVES FLOODING

000%2509/86 Pergamon

S3.OOt0.00 Journals Ltd.

ON LIQUID FILM SURFACES-X IN A VERTICAL TUBE

HSUEH-CHIA

CHANG

Department of Chemical Engineering, University of Houston, University Park, Houston, TX 77004, U.S.A.

(Received 31 August 1984) Abstract-The recently available data on flooding in two-phase annular flow are examined via a simplified version of the Navier-Stokes equation. Instantaneous stress exerted on the interface by the gas phase is described by validated empirical correlations. Analysis of the equation then explains why the tit-film thickness is found to be in agreement with time-averaged values under total downflow and total upflow conditions. Wave celerity obtained from the equation is also found to be in excellent quantitative agreement with experimental data. More importantly, the simplicity of the equation allows analysis on the growth of nonlinear, finite-amplitude waves from the surface. A mechanism explaining the observed irregular wave behaviour near flooding is suggested.

INTRODUCTION

“Flooding”

of falling

liquid films in vertical

tubes has

been studied for over 40 years [see the review by Dukler (1976)]. Nevertheless, some important details have only recently been uncovered with the help of more sophisticated equipment. Measurements of instantaneous interfacial stress, local film height and wall stress have only recently become available. In particular, we refer to the report of Smith et al. (1984) and the Ph.D. thesis of Zabaras (1985). (The low liquid Reynolds number experiment, W, = 0.0126 kg/s or Re, = 310, in their studies will be chosen as the standard case to which our theoretical results will be compared.) From these recent experiments, several unique characteristics of annular flow have been revealed. At low or zero upward gas flow rate, the film surface is observed to be smooth near the entrance region. A ripply surface then develops downstream with smallamplitude and high-frequency waves. Further downstream (N 1 m), the structure of the waves changes drastically, with the surface dominated by large “roll” waves with long wavelengths (much larger than the film thickness) separated by relatively flat regions known as “substrates”. Strangely, despite these drastically different surfaces, the local time-averaged film thickness at various column locations is found to be identical within 10 %. For freely falling films, this value agrees with or is slightly less than the Nusselt thickness (Nusselt, 1916), and at low flow rates the local timeaveraged thickness of highly irregular surfaces is in agreement with the corresponding flat-film solution! [In Zabaras’s experiment (1985), a slight negative deviation from the flat-film solution was observed which is attributed to the 7 0/0 difference predicted by Kapitza (1964).] The celerity of individual waves, while increasing downstream, does not change signficantly over the entire column of 2m. However, different waves seem to travel at different celerities with some correlation between the wave amplitude and the celerity (Zabaras, 1985). As mentioned before, the

amplitudes of the waves, in contrast, grow drastically downstream, forming large “roll” waves. There are speculations regarding the formation of mixing vortices within the humps of these roll waves (Dukler, 1976). At extremely high flow rates, total upflow occurs and the liquid film moves upwards from the feed (see Fig. 2). Analogous to the total downflow situation described above, there is also a smooth region immediately above the feed with a fat-film thickness. Again, similar to the downflow case, small and then large waves grow downstream (upwards) but the timeaveraged thickness is still close to the flat-film solution. The other phenomena are also qualitatively similar. Dramatic differences, however, are observed in the transition region between total downflow and total upflow. In this flooding region, the time-averaged thickness at every column location undergoes a discontinuous jump as the gas velocity increases beyond the Aooding point. This is also accompanied by a discontinuous pressure gradient and wall shear. The time-averaged film thickness is still in agreement at different locations of the column but its value does not correspond to any of the various possible flat-film solutions. More interestingly, there is a distinct maximum in the RMS film thickness variation in time at flooding. This large scattering is also reflected by the wide probability density of the film thickness over time at flooding which encompasses waves whose amplitudes are four or five times the mean value. Before and beyond the flooding point, the probability density is narrow with a small variance. Visually, the wave structure at flooding reveals extremely high waves undergoing “chaotic” motion even at the entrance region. (A region of smooth film adjacent to the entrance is no longer present.) It is hence reasonable to speculate that flooding is a highly unstable dynamic phenomenon which cannot be described by flat-film steady-state theory or linear stability theory. The fact that this irregularity in time is only present during flooding indicates that the randomness is not due to 2463

2464

HSUEH-CHIA

external noise but is an inherent property of flooding. The best experimental data on these chaotic patterns of the film surface near flooding are reported by Zabaras (1985), who cross-plotted the instantaneous interfacial stress against the instantaneous film height over a long period of time and obtained a totally uncorrelated pattern. It is interesting to note that the aforementioned unstable behaviour near flooding has inspired a whole class of models. In these models, the flooding point is considered as the condition under which the surface becomes linearly unstable. Individual mechanisms differ in the assumptions employed. Imura ef al. (1977) and Chung (1978) assume potential flow in both phases. Cetinbudaklar and Jameson (1969) assume a viscous laminar velocity profile in the liquid and inviscid gas Bow in the base state. An Orr-Sommerfeld stability analysis is then carried out in all the models to determine the onset of instability which is then assigned the flooding point. One physical argument is that the growth of large waves from the unstable interface will then “bridge” the tube, as first suggested by Duffey et al. (1978) and invalidated by the experiments of Smith et al. (1984) and Zabaras (1985) whose maximum wave amplitude was never close to the necessary value. The other suggested flooding mechanism is that these finite-amplitude waves will carry the liquid upwards. However, the experiments of Zabaras (1985) clearly indicate that the waves are downwardtravelling at and slightly beyond flooding (see Fig. 8). There are other criticisms on the validity of this class of models. Interfacial shear stress on the liquid film, which has to be the underlying physical cause for flooding, is not included. More importantly, onset of linear instability of the flat-film base state at flooding is unreasonable since the flat film is unstable under al1 conditions, as evidenced by the appearance of finiteamplitude waves in freely falling films, total upflow and total downflow. It is hence our contention then that flooding can only be understood through a nonlinear dynamic theory which includes interfacial stress exerted by the gas phase and delineates finite-amplitude waves beyond the linear region. It would be ideal to consider the full Navier-Stokes equation in both phases for this purpose. However, this is not feasible at present for many reasons. Firstly, the highly turbulent gas flow precludes satisfactory computation in this region. (Potential flow cannot be assumed owing to the viscous boundary layer near the interface which exerts the necessary stress for flooding.) Even if freely falling films are considered, which allows the negligence of interfacial stress, solution of the Navier-Stokes equation in the liquid phase is still formidable. A major breakthrough in this problem has recently been achieved by Bach and Villadsen (1984), who employed a Lagrangian finite-elements technique to verify the existence of a sohton solution. Their approach, however, cannot simulate mixing vortices in “roll” waves which have been speculated. Also, non-solitary solutions are not studied.

CHANC

In our present approach, we simplify the NavierrStokes equation by restricting our analysis to long waves. A simpler nonlinear partial differential equation for the film height, as a function of time and spatial coordinates, results from the perturbation analysis. The effects of the gas phase are conveniently studied through empirical correlations for the interfacial stress. In this manner, analysis of the highly turbulent gas phase is circumvented. This parallels the approach of Miles (1957) in his study of windgenerated plane waves. The effects of ripples which ride on the long waves are considered indirectly by using different correlations for downflow and upflow. Qualitatively, ripples roughen the surface of the waves and increase the friction factors of the correlations. Although instantaneous interfacial stress correlations are not available in the literature, we use two timeaveraged correlations which have been shown by Zabaras to agree with his instantaneous data away from the flooding region (for total downflow and upflow). The absence of valid stress correlations for the flooding region does not totally prevent us from studying the flooding phenomenon. Results independent of the stress correlation indicate that finiteamplitude waves reach their maximum amplitudes near the gas velocity corresponding to flooding (see Part II of this series). Away from the flooding region, our theory explains many of the experimental observations of Smith er al. (1984) and Zabaras (1985). The agreement between the time-averaged film thickness at all column locations and the flat-film value is rigorously demonstrated. The computed wave celerity is also in quantitative agreement with their data. The long-wave lubrication theory approach taken here extends the work of Benney (1966), Krishna and Lin (1977), Sivashinsky and Michelson (1980) and Shlang and Sivashinsky ( 1982) for freely falling films to general annular flow with interfacial stress. The major advantage is that the nonlinear character of the Navier-Stokes equation is captured by a single equation. Hence, the growth of nonlinear waves can be moreconveniently studied. We will subsequently show that, upon proper scaling, our equation reduces to the one derived by earlier workers for freely failing film, ~r+4~&+&,+&,,,=O

(1)

where ++J is the reduced wave amplitude and t and < are the scaled time and space variables, respectively. Standard techniques for the well-studied Korteweg de Vries equation (Jeffrey and Kakutani, 1972) are not applicable to eq. (1) owing to the numerically destabilizing fourth-order derivative term. Consequently, there seems to be controversy regarding the dynamic behaviour of solutions to this equation. Atherton (1972) and Tougou (1981) obtained results that indicate that a “saturated” finite-amplitude wave is reached at large time from any initial condition, whereas Sivashinsky and Michelson (1980) obtained chaotic behaviour indefinitely without convergence to a steady wave profile. This controversy will be resolved in Part II of this series.

Nonlinear waves on liquid film surfaces-I

There has been little other theoretical work done on vertical film flow with interfacial stress caused by an upward gas flow. Solevev ef al. (1967) obtained the flatfilm steady-state solution as a function of interracial stress. Since the overall liquid film flow can be either upwards or downwards, three solution values on two branches (corresponding to up and down flows) are obtained at the same interfacial stress value beyond a critical shear stress. Maron and Dukler (1984) added a third branch by assuming a split-flow configuration. Consequently, a maximum of five steady-state solutions are possible at a given shear stress! Smith (1970) found that any downward flow flat-film solution is linearly unstable, as expected, and for large wavelength disturbances to flat surfaces without interfacial stress, the wave speed of the disturbance is c = U*(H*)--

Air

Fig. 1. Coordinate system for the liquid film. The air flow is against gravitational pull.

where

2v

u=

u

0

(5)

V

and U is the component in the x-direction. The boundary conditions for eqs (3) and (4) are, at Y = 0, u=

v=o

(6)

and at the interface, Y = H, kinematic condition and stress balances apply,

v=g+ug CTsl = 0

EQUATIONS

The dimensional governing equations for the film flow problem are the Navier-Stokes and the continuity equations for the liquid phase, the no-slip conditions at the wall, and the kinematic, normal and stress balances on the film surface. We will assume that gas phase is at quasi-steady state, which exerts on the film an interfacial shear stress the value of which is determined instantaneously by the film thickness. We will also neglect the curvature of the tube so that we do not need to work in cylindrical coordinates. (The correction of the flat-film solution due to tube curvature can be which is a shown to be of the order 2H*‘/3D, negligible 0.4% correction for Smith et al.‘s system.) However, the correlations we use for the interfacial stress will be specifically for flow in vertical tubes with explicit dependence on the film thickness. Hence, these correlations restrict our analysis to film flows in tubes. We will also neglect the cross-stream dimension of our system in order to reduce the number of spatial As a result, bifurcation of twodimensions. dimensional waves with azimuthal variation, which has been observed, must await a later study. Using the coordinates and geometry depicted in Fig. 1, the dimensional equations for the liquid phase are f?U = -VP+pV~U+p ; p ~+UmJ (3) > 0

(

v.u=o

Flow

gH*=

where U *(H *) is the interfacial velocity of the flat-film solution. In all these studies, the interfacial (shear) stress is treated as an independent parameter which can be freely varied. While this may be true for flat surfaces as considered by Smith, the interfacial stress in a tube is a function of the film thickness which changes theeffective tube radius [see Wallis (1969)]. Hence, the interfacial stress is not an independent parameter but a system variable that needs to be computed like the film thickness. We will show that this consideration can change the previous results drastically. GOVERNING

2465

(4)

[T,]

+

GHxx

(1+H:)3’2

(8) (9)

=’

where [ ] denotes the difference across Y = H and T,, T,, the interfacial tangential and normal stresses, T, = p(g+g)cos2$

+fl(g-$$)

sin2* (10)

where tan $ = H,

(12)

is the slope of the film surface relative to the vertical wall. For the gas phase, one simply replaces the parameter values in eqs (3), (41, (5), (10) and (11) by the gas-phase values; but since we will use correlations for gas-phase interfacial stress, these replacements will not be carried out explicitly. The stress exerted by the gas flow arises from the complex interaction between the highly turbulent gas flow and the small ripples on the surface of the long waves or flat films. In making this statement, we have invoked the two-wave structure observed in experiments. The long waves are about 5cm long and they

2466

HSUEH-CHIA CHANG

carry the bulk of the liquid flow. The small ripples (_ 0.1 cm) appear on the crests of long waves and the thin troughs (substrates) separating the crests. By making our long-wave approximation, we will, in effect, neglect the ripples and concentrate on the allimportant long waves. However, the ripples must be taken into account in the evaluation of the interfacial stress, as we have mentioned. In essence, this is how the long waves interact with the ripples. At the present time, even with the help of high-power computers, it is still impracticable to compute the highly turbulent gas flow field and its interaction with the ripples in a small boundary layer at the interface to produce the interfacial shear stress. One reason for the computational difficulty is that the gas cannot be considered inviscid in the boundary layer. Hence, we follow Miles (1957) in his studies of windgenerated plane waves by using correlations of interfacial stresses. Unlike his correlations for infinite planes, however, ours contain an explicit dependence on the film thickness because of the tube configuration. Also, since the ripples may undergo bifurcations of their own, the surface “roughness” can drastically increase or decrease. Hence, these correlations may not be a continuous function of the film thickness, H. In general, the normal stress is negligible and the tangential stress can be written as r$ = 7-5 (H/D, us).

(13)

It is important to realize that the literature correlations represented by eq. (13) are usually obtained from time-averaged stress and film-thickness measurements. In our formulation, whichassumes a passive gas phase behaving at quasi-steady state, we will assume that the insranraneous stress for a time-dependent long wave is also described by eq. (13). Experimental evidence for this has been supplied by Zabaras (1985) who finds this assumption to be invalid only near the chaotic flooding point. There are two correlations for upflow. The Wallis version is given by (Wallis, 1969) T$ = To - T1 g

correlation I

(14)

where TIJ= 0.005 ps ug2/2 T1 = 300 x 0.005 p, u,2/2.

(15)

The correlation by Bharathan (1978) has no linear term, T$=

2 -To-T,

$ 0

correlation II

(16)

where Tz = 406 p, u;/2

height (Zabaras, 1985) for total downflow and total upflow configurations. A stress correlation for the flooding region is not available. Consequently, we shall continue to use correlation I or II in the flooding region for lack of any convenient correlation. In Zabaras’s experiments, both correlations are found to underestimate the friction factor (viz. stress) in the flooding region. Insofar as interfacial stress destabilizes the flat-film base state in the flooding region, our subsequent study yields aconservative lower bound on the highly nonlinear wave behaviour near flooding, viz. one should expect even higher amplitude waves and even more chaotic behaviour in this region than what we predict. Outside of the flooding regions, in the total upflow and total downflow regions, these two correlations are quantitatively accurate and so are our results, as we will demonstrate in the next section. Another important point regarding correlations I and 11 is the implicit assumption that the stress is the same on both the leeward and windward portions of the wave with equal local film thickness. Consequently, we are imposing the stipulation that flow separation of the gas phase does not occur behind the wave. Although gas flow separation is a necessary and important factor of another model for flooding (Shearer and Davidson, 1965), we feel that the small curvature of the long waves (amplitude/wavelength ratio < 1) supports this contention. Nevertheless, a certain degree of asymmetry probably exists in actual waves. For downward film flow, eq. (14) or (16) is still applicable to a critical gas flow rate, u,. For vg < uc in the downward case, the stress decreases drastically to a value independent of H (Wallis, 1969) which corresponds to the zeroth-order in eqs (14) and (16), Tf =

and we have rounded off the power of (H/D) in the actual correlation (2.04) to an integer (2) for simplicity. Both correlations have been validated by Zabaras’s measurements of instantaneous stress and local film

To.

(18)

The actual drop-off point seems to be a complex function of liquid flow rate, tube diameter and other fluid and equipment properties [see Fig. 11.2 of WalIis (1969)]. It seems to correspond to the disappearance of ripples on the surface. The constant stress To corresponds to a Moody friction factor of 0.005, which is the value for a smooth pipe. The higher-order terms in eqs (14) and (16) correspond to the increase in stress due to the increase in film thickness and the appearance of ripples. The To value, when compared to the Moody correlation, corresponds to a grain size of 0.001-0.03 D for the “sand roughness” factor. Since the long waves are much longer than these values, this is another indication that they do not contribute to the instantaneous interfacial stress. We will lump eq. (18) and eq. (14) or (16) together: Tf (H, ue) =

(17)

-

- T, - z#‘(up - To - Ti

- u,) for downflow for upflow

(19)

where &’ is the Heaviside function and Ti is equal to T, (H/D) or T2 (H/D)L, depending on the correlation used. Since there is no available correlation for the critical gas flow rate u, at the flooding point, it will be

Nonlinear waves on liquid film surfaces--I treated as a free parameter which fitting the experimental data.

is determined

by

In subsequent analysis, we will use the gas flow rate us (or the equivalent gas mass flow rate f) as our control parameter instead of the interfacial stress Tf used by Solevov et al. (1967), Maron and Dukler (1984) and Smith (1976). Note that from eqs (14)-(19), the relationship between 2-f and 9 is not one-to-one and, more importantly, the interfaclal stress is dependent on H, a system variable. Hence, Tf should not be used as an independent parameter. FLAT-FILM

BASE STATES

We shall solve eqs (3~(9) for the time-independent flat-film solution (x-independent) by using eq. (19). Again, since eq. (19) includes the ripple effects, by “flat film” we mean flat with respect to the long waves. Ripples appear on these surfaces and they modify the stress according to eq. (19). They are, however, invisible in the flat-film solution of our long-wave approximation. The formal perturbation expansion will be carried out in a later section. If the gas phase is driven by a constant pressure gradient in the x-direction, AP/L, then for the flat-film solution the liquid pressure is equal to the gas pressure from eq. (9) and hence P is only a function of x and dP p=-. ax

For the air-water system of Smith et al. (1984), the gravitational term, pg. is much larger than the pressure gradient term, which is of the order psg. Hence, eqs (3)-(9) reduce to d’U* --+pg=o (21) ’ dY2 and u*(o) = 0 (22a)

$!&

(H*)

Note that these equations are insufficient to determine the flat-film thickness H *. To do so, we need to introduce the normalization conditions stipulated by liquid flow rate Q (volumetric flow rate per unit circumference). There are, however, three possible configurations-upflow, downflow and splitflow. Upand downflow correspond to the cases where all the liquid flows up and down, respectively, from the feed. Splitflow is the unusual case near flooding when a portion of the feed moves up while the rest flows down. These different configurations, depicted in Fig. 2, may be possible at the same us. causing steady-state multiplicity. The associated normalization conditions are H*

s0 H*

dY = - Q

U*

s

for downflow

(24a)

for upflow

(24b)

and Jb” U d~+~~~~

for splitfiow (24~)

U dY1 = Q

where Q > 0 is the liquid flow rate per unit crossstream direction (cm*/s), and from eq. (23),

y 0 =*

PC?

(TbtH*)+e,,

cc

P

>

t244

is the point at which U * is identically zero. The three flow configurations are depicted in Fig. 2. Substituting eq. (23) into eqs (24a)-(24c) and introducing the dimensionless variables R = H/H, Fo =

(25)

Tol(Q~c/Hk) = 37.oI~gH,

F,

= 3TJ2pgD

F,

= 3 T2H,/2

(26a) (26b)

pg D*

t2w

where

= Tf(H*)

“3

which yields the flat-film velocity profile u*

U*dY=Q

AP L

2467

(27)

is the Nusselt thickness which corresponds to downward liquid flow with no interfacial stress (F, = FL

=

Down

UP

Split

Fig. 2. Schematic depictions of the three flow configurations. Flow profiles A and B are the limiting cases of splitflow.

2468

HSUEH-CHIA CHANG

= F2 = 0), one obtains the following three algebraic equations for R: 2(1 - F,)R3

- FoR2

2(1 -Fi)R3-FFoR2

= 2

downflow

(28a)

=

upflow

(28b)

-2

and [8(3-2FJ3+54Fi-541 + [27 - 8 Ffj -

R3

24(3 - 2 Fi)2] F,R2 54 = 0

+ 24(3 - 2 Fi) F;R

splitflow

R

(28~)

where F,(R)

= F,

(29a)

if Wallis’s correlation is used [correlation I of eq. (14)] and F,(R)

= F,R

(29b)

if Bharathan’s is used [correlation II in eq. (16)]. For a given liquid-gas system, the liquid Aow rate only affects H,. The other parameters are fixed at Smith et d’s experimental values: D = 5.1 cm, ps = IO- 3 g/cm”. p = 0.01 g/cm s and p = 1 g/cm3. For liquid flow rate per unit circumference, Q, of 0.7864cm2/s (corresponding to an overall liquid flow rate of W, = 0.0126 kg/s in Smith et al.‘s and Zabaras’s experiments), H, is found to be 0.0289 cm. We use this value as our reference value (standard case). Using eqs (15), (17) and (26), one obtains

2

-0 Air

Flow

f Rate

_J

4

6

(xlOgm/sec)

Fig. 3. Comparison of the dimensionless flat-film height for the standard case (HN = 0.0289cm) as a function of air flow rateffrom correlation I against the time-averaged data of Smith et al. (1984). Points A and B correspond to the limiting cases depicted in Fig. 2. Curve I refers to downffow with f -z f, Curves II, III and IV refer to downflow with f >f,, splitflow and upflow, respectively. Open triangles are data below the feed and closed ones are data above the feed.

4 F,

= 0.0635

F,

= 0.0539 fz

(3Ob)

F,

= 0.0827

(3k)

(H;/H,)f*

(H,/H;)f’

(3ON

where Hg is 0.0289 cm andfis the gas mass flow rate in units of log/s. We will replace u by f in subsequent discussion. This seemingly clumsy unit for us is introduced by Smith et al. (1984) and Zabaras (1985) and for convenience in subsequent comparison we will retain their unit. Equations (28~(30) are only valid for upflow and downflow withf > f,, wheref, is equivalent to uCin eq. (19). Forf
= F,

= 0.

CORRELATION

II

3

R

2

1

(306)

In Figs 3 and 4, we compare the results from correlations I and II to the time-averaged data of the standard case (Hz = 0.0289). It is seen that f, is approximately 3.0 and this can be assigned as the flooding point. The splitflow configuration is only possible if Ye of eq. (24d) lies between 0 and H *, viz. when there are streamlines of opposing directions. Hence, the splitflow curve joins the down- and upflow curves at points corresponding to cases when the interfacial velocity first becomes negative and when the slope of the wall velocity vanishes. These points are labelled A and B, and the corresponding velocity profile is sketched in Fig. 2. Point A is not of interest because the jump at the flood point& usually precludes its existence (see Figs 3

0

0

2

4

6

f Fig. 4. Same as in Fig. 3 for correlation II. Points A and B and curves I-IV are defined in the caption of Fig. 3.

and 4). Point B, however, is most interesting. It occurs for R=R,=3,/2.

(31)

For R > R,, a flat film with opposing streamlines occurs and this is usually considered to be unrealistic. Our results seem to confirm this since R values as high

Nonlinear waves on liquid fiIm surfaces-1

as this portion of the upflow curve are not observed. Moreover, the agreement of the flat-film solution is best at low and high gas rates. In the flooding region near the flooding point, correlation I overshoots the experimental R values, whereas correlation II undershoots them. In Fig. 5, we plot R vs. the dimensionless interfaciai stress for the flat surface, r* = TfH;/Q

/A.

R

1.25

-

1.00

-

0.75

-

2469

(32)

All correlations map into one curve in this case. Maron and Dukler (1984) have noted that the experimental data converge to point B as the gas flow increases, when plotted on this diagram. Their data are also shown in the figure. In Figs 6 and 7, the data taken above the feed from Smith et al.‘s (1984) experiments for various liquid rates are systematically compared to the flat-film solutions of both correlations. Again, the fit is not as satisfactory near flooding but as f increases, the data seem to converge to the flat-film solutions of both correlations. The agreement of local time-averaged film thickness with the flat-film value away from the flooding region has been observed experimentally by Zabaras (1985), who reports a similar plot to Fig. 5. This is surprising since extremely irregular and finite-amplitude waves cover the surface at many locations. To understand this, one must study the dynamic equation instead of the steady one, eq. (21). Insofar as the dynamic Navier-Stokes equation, eqs (3 j(5), is too difficult to analyse, we derive a simpler long-wave equation in the next section for this purpose.

o-5&-c-

’

4.0

3.5

4.5

5.0

5.5

f

Fig. 6. Dimensionless flat-film height vs. air flow rate for various liquid Bow rates. Correlation I is used. Curves III and IV refer to splitflow and upflow,respectively. The experimental data of Smith et al. (1984) are also shown.

4ti

CORRELATlON

0.50

I 2.5

3.0

3.5

1

II

4.0

4.5

5.0

I 5.5

f Fig. 7. Dimensionless flat-film height vs. air flow rates for several liquid flow rates using correlation II is compared to Smith et al.‘s data. R

LONG-WAVE

L OO

I

2

4

1 6

7*

Fig. 5. Dimensionless flat-film height plotted against dimensionless interfacial stress. Both correlations I and II of Figs 3 and 4 map into the same curve in this plot. Points A and B and curves I-IV are defined in the caption of Fig. 3. The timeaveraged data of Smith et al. (1984) and Zabaras (1985) are also plotted for comparison.

APPROXIMATION

In this section, the formal long-wave perturbation formulation of the complete equations, eqs (3)-(9) with the interfacial stress correlation of eq. (13) is carried out. We begin by rendering these equations dimensionless with respect to the flat-film solutions. For convenience, the characteristic length in the x-direction, which cannot be provided by the flat-film solution, is taken to be the yet unspecified wavelength lo of periodic waves. However, the formulation is entirely valid for nonperiodic boundary conditions. In such cases, lo is simply a measure of the variation in the xdirection.

HSUEH-CHIA

2470

The dimensionless variables are

CHANG

W”)

x = X/I,,

Y = Y/H*,

h = H/H*,

u = UICQ/H*h

0 = V/(QIH*k,

P = P/do,

(33)

r = TI(bH*IQ)

(37b)

where the characteristic length scales are IO, the unspecified wavelength divided by 27r, H *, any of the flat-film thicknesses obtained from the earlier sections, the characteristic velocity in the x-direction is Q/H* and the time scale is I,H*/Q. The perturbation parameter is defined as E = He/l,

* Y,

F = g H*3/vQ, (35)

W = G/pgH*2.

We will restrict ourselves to cases where i?, F - O(1) and&W - 0( 1). Hence, the subsequent analysis is only valid for small H* at low liquid flow rates. However, one must note that fi and F should be O(1) relative to E. For long waves and thin base-state thickness, E is extremely small and a Reynolds number in the low hundreds can also be tolerated. We also neglect the gas-phase pressure gradient. Expanding the liquid-phase velocity and pressure in an asymptolic series in E,

av,

(37c)

+--0

dY and the boundary conditions Y =

(34)

which is typically small (- 1O-2-1O-4). Note that the characteristic scales are dependent on the gas flow rate f and the liquid flow rate represented by H,, since the base states, the flat-film solutions, are functions of these parameters. From Figs 3 and 4, it is evident that there is a multiplicity of base-state solutions for a given f (gas flow rate) near the flooding region. Hence, there is an ambiguity of which one should be used for the scaling of eqs (33) and (34). Whenever possible, we will not specify the actual base state. When this is not possible, we adhere to our policy of using the most conservative one. Along this line, the lowest splitflow (curve III in the figures) base state will be used which will be the more stable base state with the lower film thickness. Substituting these dimensionless variables into eqs (3)-(9), the usual dimensionless groups, Reynolds, Froude and Weber numbers, appear as &?= Ql,/H

au, -ix

ug =

0,

Y=h,

vg =

au,=

L40+EU1

+&%,+

v=vo+&V,+E2v2+

(38b) (38~)

PO = 0

where r is the dimensionless interfacial stress scaled with respect to H* instead of H, in the version for flat surfaces [eq. (32)], r(h) = TF H*‘/Q,u.

(39)

From eqs (37b) and (38c), p,, vanishes, reflecting the fact that the pressure gradient is O(E). Hence, eqs (37) and (38) are resealed versions of the flat-film equations, eqs (21) and (22). In the present case, however, u0 and ho are functions of x, y and t. Also, v0 is nonvanishing. Specifically,

y s

u,(x, y, t) = - F y2/2 + (r + Fh) y %(X, Y, ‘) = -

(4w

au,

o

--dy dx

Fob)

where the t and x-dependence arises from h(x, t). The resemblance is because, at zeroth order, the long waves appear as flat surfaces locally. The steady-state flat-film solution can be derived by substituting unity for h in eq. (40a). Hence, u*(y) = - Fy2/2 + (rO + F)y

(41)

where TV =

r(l)

=

7*RZ

(42)

which is related to the other dimensionless flat-film stress r* of eq. (32). Applying the normalization condition of eqs (24), one obtains the following normalization relationships between F and r0 for the various configurations of Fig. 2:

F = (6 - 3s0)/2

upeow

(43a)

F = (-6--3r,)/2

downflow

(43b)

splitflow

(4%)

w

. ...

Pa)

r(h)

aY

8r~+21roF2+24r~F+6F3-6F2 u =

0

(36)

p = plJ +&PI + 2pZ + . . . and collecting terms in the dimensionless equations, one obtains for the following zeroth-order equation which contains no inertial terms,

=0

The first-order equations are d%J, 7-F$=l? ay aP, dy=O

$+u~~+,~‘EIO

(44a) aY > (Mb)

waves on liquid film surfaces-1

Nonlinear

!%+?!L=O

(-1

aY y=o,

u1

for

(a)

= 0.

(440

the zeroth-order

u, and

u1 (x. y,

(44)

dy=O

p+EWh,, Applying

=o

=ul

au,

Y = h,

solutions,

we can

solve

u,,

t) = EM’Fh,,,(hy

- y2/2) + fiF’h, (y3/6 ~ h2y/2)

+ Rr&(y’/6

- h2y/2)

+ J?(r + Fh) F h,(y4/24

- h3 y/6)

+ R(t + Fh) T’h,(y4/24

- h3y/6)

(45)

and UI(X>Y,t)

=

ydu,

-

(46)

axdY

s” where

z’(h) = ;,:. We next interface

invoke

the

condition

at y =

at

1

the

(48)

yields

h,+uohx-uVo+&[u,h,--ul]

= O(E’)

sty=

1.

(49) Since d

h(x.

uih, - v, = 8-x s 0

I,

ui dy

(50)

for i = 0 and 1, we substitute eqs. (40) and (45) into eqs (49) and (50) to obtain the following equation for h which is correct to O(E):

h,+A(h)h,+&

{B(h)h,+C(h)h,,,}

= 0 (51a)

where

A(h) = Th + h* ?/2 + FhZ B(h) = $

J?(r + Fh) (F + r’) h5

(5fb)

C(h) = Fh3EW/3. For

the case with

no interfacial

stress, T = T’ = 0, eq.

(51) reduces to the version derived by Benney (1966) after correcting for the difference in the definition of parameters and variables. Note that given a particular solution h(x, t) to the highly nonlinear partial differential equation (PDE) in eqs (51), the interfacial velocity and the wall stress, two often (40a)

measured and (45),

variables,

u(h) CES 41/10-M

u,(x,

can

Since the interfacial stress, T, is negative by our notation, it is clear from the leading order approximation of eq. (40a) that flow reversal of the liquid indicated by a sign reversal of u(h) can occur. This is most likely at large h, viz. at the peak of the waves, as is consistent with the speculation of a mixing vortex in that region by Dukler (1976). For freely falling films (5 = 0), flow reversal cannot be obtained from the zerothorder expression since uo(h) is always equal to the positive value of FhZ/2 and one needs to include u, of eq. (45), which is negative on the leeward side (h, > 0) of large waves due to the fourth-order term in the expression. Since the contribution u1 (h) in eq. (52a) is multiplied by E, flow separation is less likely in freely falling films except in very large waves far downstream from the feed. That the second-order term in a such as ours can cause flow lubrication formulation separation has bien shown by Chow and Soda (1972) for blood flow in periodically constricted tubes. The leading order approximation to the wall stress is ‘5, -

(47)

kinematic

h, = uh, - v = 0 which

247 1

be obtained

h, t) + EU~(x. h, t)

from

eqs

(52a)

r(h)+

Fh.

(52~)

Since r(h) is linear with respect to h from correlation I, eq. (52~) suggests a strong correlation between the wall stress and the local film thickness, with the maxima of both occurring at the same point in a wave. In Zabaras’s data, the maximum in t, Follows the one for h by a small time lag that is remarkedly constant for all conditions. We believe that this lag is associated with the asymmetry of the interfacial stress-thickness relationship on the windward and leeward sides of the wave. As we have mentioned in an earlier section, symmetry is assumed in both of our correlations. However, it must be noted that for extremely large waves which require resolution of the second-order terms, the built-in asymmetry ofu, due to the h, terms may be sufficient to account for the above time lag. CONSERVATION

PROPERTJES

AND

FLOODING

MECHANJSMS

The long-wave equation of eq. (51) represents a simplified version of the Navier-Stokes equation of eqs (3) and (4) with free-surface boundary conditions, eqs (7b(9). Hence, the evolution of finite-amplitude, nonlinear surface waves can be studied better via this equation. For example, the quantitative effects of varying gas flow rate f on the system dynamics can be obtained and perhaps an explanation of the chaotic surface behaviour near flooding will result. However, the long-wave equation is a nonlinear PDE containing fourth-order derivatives in < and seventh-order algebraic nonlinearities in h. As we will clarify in Part 11 of this series, the fourth-order derivative is already sufficiently destabilizing to most numericalalgorithms, not to mention the high algebraic nonlinearity. Hence, complete numerical solution of eq. (51) remains difficult and in a latter section, we will introduce further simplification to make the problem more accessible to

HSUEH-CHIA

2472

numerical solution. Nevertheless, one can still make several important statements regarding the properties of solutions to eq. (51) without explicitly solving for them. These properties explain the excellent agreement of the time-averaged film thickness to the flat-film value reported in an earlier section. They also lead to a speculated mechanism for flooding. Consider a wave or wave packet travelling on the interface. Typically these waves are bounded on both ends by smooth substrates. Subtracting the small ripples which are not described by the long-wave equation and which tend to reduce the average film height due to nonlinear interactions of the ripples as demonstrated by Kapitza (1964), these bounding substrates are essentially flat films at the base-state thickness (Dukler, 1976; Smith et al., 1984; Zabaras, 1985). In the next section, we will demonstrate that the waves can be assumed to be travelling at a constant celerity A. Since we are isolating a particular wave packet bounded by two substrates, the waves may travel at a different velocity from the substrates. The celerity of the substrates need not be specified in this formulation. Hence, these waves are mass carrying which can overtake or fall behind the substrates. In the moving coordinates t defined with respect to the wave celerity A [see eq. (68)], eq. (51a) becomes h, -Ah,

+ A(h)h, + E$ {B(h)h, +C(h)h,,,}

= 0. (53a)

Note that the wave profile h is not stationary as it propagates but is allowed to evolve as indicated by the h, transient term. Assuming that the substrates are long compared to the wave packet, the appropriate boundary conditions for one particular isolated wave packet solution to eq. (53a) are the ones corresponding to solitary waves, h(C = - cu) = h(< = co) = 1 h, = h,, = hcss =

0

at 6 = + co.

(53b) (53c)

Boundary condition (53b) stipulates that the substrate thickness is the flat-film base-state value and condition (53~) reflects the (relative) smoothness of the substrates. We are not alone in imposing these boundary conditions. Bach and Villadsen (1984) used the same conditions in their numerical work. In fact, all studies on solitary solutions of any equation stipulate the same set of boundary conditions (Jefferey and Kakutami, 1972). Integrating eq. (53a) from -cc to + 00 in the moving coordinate and invoking the boundary conditions of eqs (53b) and (53c), one obtains a timeinvariant property of the wave,

Hence, if the wave profile is flat at any moment during its propagation [the flat profile h = 1 is always a solution to eq. (53a) and boundary conditions (53b) and (53c)], the average film thickness in the moving

CHANG

coordinate is always equal to the base-state value! This unique conservation property is due to the fact that all the spatial derivative terms with respect to 5 in eq. (53a) are exact differentials which, after integration, yield terms at the infinities which cancel or vanish due to boundary conditions (53b) and (53~). Recall that during total upflow and total downflow, a flat film region always exists near the liquid entrance which implies that the time-averaged film height at all downstream locations must correspond to the height in this flat region! (Note that spatial averaging in the moving coordinate is equivalent to time-averaging at a particular location.) This explains why the timeaveraged film heights agree so well with the flat-film value for upflow and downflow in Figs 3-7. It is instructive to note why eq. (54) is invalid in the flooding region. It is not simply because a flat region is never observed throughout the column in this region, which will also be explained subsequently. Since there is a multiplicity of base states in this region as shown in the flat-film section, the film height at the front substrate can be different from the back substrate, viz. Ah=h(<=

-coo)-h(<=

-co)#O.

(55a)

[The unity value in condition (53b) is ambiguous since it is scaled with respect to a particular base state.] With the new boundary conditions, subsequent integration of eq. (53a) yields ;

s

_= hd< = AAh. 0”

t-b)

Hence, there is a constant pick-up or loss of mass as the shock propagates, depending on the sign of Ah. Also, if Ah is negative [h( - co) > h( + co)], the shock travels in the positive direction (downwards). For positive Ah, the opposite is true. Note that these “shocks” are covered with waves, which may have large amplitudes corresponding to “roll” waves. These surface waves travel at celerity A as stipulated by the formulation. Hence, the waves and the underlying shocks may travel in opposite directions! More importantly, since there are multiple base states, individual shocks can travel in totally different directions as governed by their boundary substrates. Hence, the surface is a confusion of shocks and waves interacting in an extremely complex manner which leads to the chaotic surface behaviour. What is missing from our speculated mechanism is what determines the base state a particular substrate takes and whether it will switch to another one subsequently. We suspect that this is related to the initial conditions and subsequent interactions among various shocks and waves. There is experimental evidence corroborating the above mechanism for flooding. Maron and Dukler (1984) report flooding wave traces which repeatedly “switch” among the various distinct base-state values. These portions correspond to the substrate regions of the wave traces, according to our theory. Also, as we will demonstrate in the next section, the waves riding on the shocks in the flooding region travel downwards (positive A), regardless of the base states used for the

2473

Nonlinear waves on liquid film surfaces-I substrates.

[This is amply supported

by the experimenand Zabaras the underlying shocks

tal measurement of Smith et al. (1984) (1985);

see Fig.

8.7 However,

(substrates) may travel upwards for positive Ah. Hence, the combination of the two gives the picture that the substrates are being sucked up by the downward moving “roll” waves. Hewitt et al. (1984) have reported observations of these opposite travelling “substrates” and have formulated a phenomenological model to explain Note

them. that the above mechanism for Aooding is in no way diminished by the fact that we do not have an accurate interfacial stress correlation for this region. The underlying mechanism of multiple substrate thickness originates from the multiplicity of flat-film solutions due to the different possible flat-film flow

configurations in this region [see eqs (24)]. This existence of a multiplicity region is independent of the interfacial stress correlation used. Even the location of the multiplicity region is not significantly affected by the correlation since the three flow configurations, by definition, should coexist only near the flooding point. LINEAR

STABILITY

OF

FLAT-FILM

SOLUTIONS

The linear stability of the fiat-film solution, eq. (41), is especially interesting because it offers a direct results from comparison against the the Orr-Sommerfeld-type stability analysis of Smith (1970). Although the Orr-Sommerfeld approach is valid for all wavelengths, analytical solution of the eigenvalue problem is usually only possible for long waves. Hence, in Smith’s analysis, the long-wave approximation is made after linearization. In this _. section, we linearize after the long-wave approxima-

tion has been made. Since the two processes should be commutative, the results should agree.

Writing h in eq. (51) in the standard deviation form from the flat-film solution, h =

l+exp{wt+ixJ

(56)

where the 2n-periodic The

objective

nature is reflected in the term ix. is to evaluate the complex number o.

Substituting into eq. (51) and linearizing, the dispersion relations

one obtains

Re {w> = E[B -c] Im{c0}

(57a)

= A

where A, B and c are constants

(57b) denoting

A(l),

B( 1) and

C(l), respectively. From eqs (57b) and (41), the phase velocity 2 can be related to the flat-film interfacial velocity, u*(l), a =r0+r’(I)/2+F = u*(1) + F/2 + r’(I)/2

(58)

which is in exact agreement with Smith’s results. Specifically, for the case with no interfacial stress, eq. (2) is obtained in dimensionless form. With interfacial stress, however, F is related to 7. because of the normalization conditions, eqs (43), due to the liquid flow balance of the flat-film base states. Since the normalization condition for splitflow, eq. (43c), is in an implicit form, we will only re-express eq. (58) for down- and upflows. Introducing eqs (43a) and (43b), the wave speed becomes A = & 3 - r,/2 + r’(I)/2

(59)

where the plus and minus signs refer to down- and upflow, respectively. This equation is valid for any given correlation of the interfacial stress 7(h). In this section, however, we wili consider only correlation I of eq. (14) as a representative example. Correlation II will give qualitatively identical results. A detailed examination of the various scalings will show that correlation I of eq. (14) yields the following expression for T,, and 7’(l), 7.

-R’F,-2R3F,

mw

r’(l) = - 2 R’F,

(f-1

=

where R is the dimensionless flat-film solution of eqs (28) and F0 and F1 are defined in eqs (30). Consequently, the wave speed of eq. (59) can be expressed as a function off and H, only a =

+t

l<=

4

Fig. 8. Theoretical linear wave celerity of the standard case compared to the data of Zabaras (1985). The characteristic velocity QfRHN is used to scale Fig. 4.25 of Zabaras’s thesis with Q/H, = 27.2 mm/s and R given by Fig. 3. Note positive celerity after the flooding point.

t 3 + R2Fo/2.

(61)

8, we plot

the wave speed for all flow In Fig. configurations. The solid curve in the multiplicity region corresponds to the most conservative splitflow. If upflow and downflow are used as base states in this region (dotted curves), there will be a sharp singularity to + co at the flooding point. From Fig. 8, it is clear that the waves continue to travel downwards (positive A) immediately after the flooding point at A. At some

distance beyond the flooding point, the celerity de-

2474

HSUEH-CHIA

creases drastically until it becomes negative, indicating that the waves begin to travel upwards at that point. The time-averaged celerity data from Zabaras’s thesis (Fig. 4.25) are also plotted in this figure for comparison. The excellent agreement amply validates the long-wave equation and our approach. It also indicates that, on the average, the finite-amplitude waves do not travel at velocities significantly different from the linear celerity provided by eq. (61). (Scatter about the time-averaged value is small.) This fact will be utilized to simplify eq. (51) further in the next section. Another interesting result lies in the evaluation of Re {o}, the growth term. The general expression is, from eqs (57a) and (52) Re(w)

= E(B(~~+F)(F+~‘(~))-~FFEWJ/LS

(62)

which, after substitution of the normalization condition for up- and downflow, becomes Re{w}

= .s{a (+6-~~)

(+3-33r,/2+r’(l))

- lOFsW}/30.

(63)

Again, this stability result is valid for any given correlation of interfacial stress. We will, as before, use only correlation I. Applying Eqs (60), one obtains Re {w} = $

{i? (R2Fo + 2R3F,

+ 2R3F,

+ 6)(3R2Fo

10Fs W>

_e 6) -

(64)

where down- and upflow are again represented by + and -, respectively. As before, the equivalent expression for splitflow must be obtained numerically. In Fig. 9, the first term of eq. (64)

-

8-

6cl(f)

In

I

i A / //

q(A HN) = (R’F,

f 2R3F,

+ 2R3FI

WEAKLY

NONLINEAR

h=

-: t

Gas

2

f

Rate

(lOgm/sec)

4

6

Dependenceof the amplificationfactor CJonf. H, is taken as the standardvalue, 0.0289cm. The splitflow (curve III) configurationis used in the flooding region. The dotted

Fig.

+ 6)/30

(65)

is plotted along with the splitflow equivalent for the standard case. The quantity q is essentially the growth term of eq. (64), if surface tension is neglected. A positive CJvalue represents exponential growth from the corresponding base state, while negative q indicates a linearly stable base state (neglecting the stabilizing effects of surface tension, of course). The linear stability analysis, however, does not allow us to determine the wave forms that result after the surface moves away from the unstable (positive q) flat-film base states. This must await the numerical study of Part II. Nevertheless, the larger the growth term (positive q), the steeper, larger and more complex one intuitively expects the resulting nonlinear waves to be. Hence, the sharp maximum in q in Fig. 9 strongly suggests that the most nonlinear waves with the largest amplitudes appear at the maximum near the flood point. This conjecture holds true even if surface tension is considered. Also, one must realize that a negative y value in the small region beyond B only indicates that the surface is stable to one-dimensional disturbances. Finite two-dimensional waves can still exist in this region. Finally, the q depicted in Fig. 9 is based on the more conservative splitflow base state (state III) in the multiplicity region. If the other states are used, q would be even larger in this flooding region, as indicated in Fig. 9.

1+&n

t’ = t&.

0

-+ 6) (3RZF,

ANALYSIS

If only relatively small amplitude waves which travel at a constant celerity given by the linear analysis of eq. (61) are pertinent, considerable simplification of the long-wave equation of eq. (51) can be achieved by a weakly nonlinear formulation. The assumption that the waves travel at the constant linear celerity is of course verified by Fig. 8. In this approach, we only consider small-amplitude deviations from the flat-film solution, h = 1, and hence retain the lowest-order nonlinearities. Also, the growth term Re (w} in the linear analysis is of 0 (E- ‘) in the previous time scale. To account for this, we define a new time scale t’with appropriate order of magnitude. Hence,

‘*1 10

CHANG

9.

curves emanating from A and B indicate the higher q values if curves II and IV of Fig. 3 are used as the base states.

(64) (67)

Also, to simplify the analysis further we will introduce the Galilei transformation, which defines a new coordinate moving with a constant wave speed A given by eq. (61), < =x-At.

(68)

Introducing eqs (66k(68) into eq. (51) and retaining only the O(1) terms, one obtains

Nonlinear waves on liquid film surfaces-1 where A’ =$(I)=

_+6+2R2F,

(70)

and q(f, HN) is given in eq. (65) and Fig. 9. In writing eq. (70), we have only used correlation I and only consider the up- ( -) and downflow ( +) configuration. Also, normalization conditions, eqs (43), are used whenever possible. For A’ = 0, eq. (69) is just a linear equation which either grows or decays exponentially according to the stability of h = 1. For A’ # 0, we define

v&!3tl

(71)

4

and eq. (69) can be rewritten as FW v,. _t 4vv< + a q II<<+ E3

v<<<< = 0.

(72)

Since the transformation c -+ - < converts the case with the minus sign (upflow) to the case with the positive sign, we will only retain the positive sign in eq. (72) without any loss of generality. Finally, we use the transformation v = Rgu,

t’ = S/Rq

(73)

to obtain ug + 4 uug + u<<+ bU<<<<= 0

(74)

where the only remaining parameter is eFW a=z

(75)

and the final variable u is related to q by

Since the spatial coordinate x is never disturbed in the various transformations, the boundary conditions for h still apply. In particular, if periodic conditions are imposed u(<. 8) = u(< +2x, 8).

(77)

It is not a coincidence that eq. (74), upon proper scaling of u, can be converted into eq. (1) which has been derived by various workers for the case without interfacial shear. Setting F. and F, to zero in the expression for q in eq. (65) we see that 4 has a value 4 for the case without shear, which is in exact agreement with the previous results. Hence, eq. (74) applies to a system with an arbitrary amount of interfacial shear, and we have hence generalized the previous works. DISCUSSJON

Although the more general long-wave equation, eq. (?+I), is simplified through a weakly nonlinear analysis to eq. (74), limited solutions to the original long-wave equation can be constructed. This has been attempted by Pumir er al. (1983) for the freely falling case. A more general analysis using dynamic singularity theory has

2475

been completed by our group and will be presented in a separate paper. Finally, the validity of periodic boundary conditions, which will be used extensively in Part II of this series, needs to be examined. Unlike the infinite domain problem which is validated in an earlier section, periodic conditions are more for mathematical tractability. It corresponds to a monochromatic wave which retains its wavelength as it grows in amplitude. From the linear Orr-Sommerfeld analysis of the Navier-Stokes equation (Smith, 1970), it is evident that a most unstable wavelength exists. This can also be obtained from the amplification factor of eq. (62), which is second order in E indicating the existence of a maximum with respect to the wave number E. Hence, it is reasonable to assume that at the inception region near the entrance, the waves are periodic with this critical wavelength. In this region, periodic boundary conditions have a solid physical basis. However, as these waves propagate downstream, the heterogeneities which destroy the periodicity begin to grow with the wave amplitudes. The average celerity is still close to the linear value but the distribution begins to spread. This causes wave interactions as waves catch up with another. The average wavelength also increases downstream as the large “roll” waves appear. By invoking the periodicity condition in this region, we are essentially removing the possibility of irregularities and wave interactions. However, there is a major advantage to the periodicity condition other than numerical accessibility. Since eq. (51) only contains a single parameter which includes the wavelength Ia through E and the gas flow ratefthrough q, the study on the effects ofJis identical to the study on IO.In this manner, we have conveniently unified the analysis for finite-amplitude waves of the same wavelength under different operating conditions (f ) and finite-amplitude waves of different wavelengths under the same condition. This also implies that increasing q has the same effect as increasing I0 (or decreasing E). This has some physical support since Zabaras observed that the roll waves that appear very far downstream under freely falling conditions resemble those at the entrance near flooding conditions. However, wave interaction through overtaking cannot be studied with periodic boundary conditions. Hence, numerical study of the infinite domain problem will also be included in Part

NOTATION

wave speed tube diameter gas mass flow rate, log/s dimensionless stress Froude number surface tension gravitational acceleration film height Nusselt film thickness given by eq. (27) Heaviside function wavelength/2rt

HSUEH-CHIA

24%

Q”

P,P R, R” T, t u, u K v W x,x y, Y

Greek E P

v

tube length Iiquid flow rate per unit circumference, H/H, pressure tube radius dimensionless flat-film height Reynolds number time, stress vertical velocity, reduced horizontal velocity Weber number vertical coordinate horizontal coordinate

film height

letters n*/kn

viscosity kinematic viscosity specific gravity

P 0

EFW/~&

7

stress

Subscripts and superscripts * flat-film solution g S

FJ

gas

stress normal Nusselt

REFERENCES

Atherton, R. W., 1972, Studies of the hydrodynamics of a viscous liquid film flowing down an inclined plane. Engineer’s Thesis, Stanford University. Bach, P. and Villadsen, J., 1984, Simuiation of the vertical flow of a thin, wavy film using a finite-elements method. Inc. J. Heat Mass Trans$iir 27, 81>827. Penney, D. J., 1966, Long waves on liquid films. J. Math. Phys. 45, 15(r155. Bharathan, D., 1978, Air-water countercurrent annular flow in vertical tubes. Electric Power Rrs. Inst. Rep. EPRI NP786. Centinbudaklar, A. G. and Jameson, G. J., 1969, The mechanism of flooding in vertical countercurrent two-phase flow. C&m. Engng Sci. 24, 1669-1680. Chow, J. C. F. and Soda, K., 1972. Laminar flow in tubes with constriction. Phys. Fluids 15, 170&1706. Chung. K. S., 1978, Flooding phenomenon in countercurrent two-phase flow. Ph.D. Degree, University of California at Berkeley. Demekhin, E. A., Demekhin, I. A. and Shkadov, V. Ya, 1983, S&tons in viscous fihns flowing down a vertical wall. Izn. Akad. Nauk SSSR Mekh. Zh. Gaza 4, 9-16.

CHANG

Dutfey, R. B., Ackerman, M. C., Piggot, B. D. G. and Fairbaim, S. A., 1978, The effects of countercurrent single and two phase flow on the quenching rate of hot surfaces. Int. J. Multiphase Flow 4, 117-140. Dukler, A. E., 1976, The wavy gas-liquid interface. Chem. Engng Educ. 108-138. Hewitt, G. F., Martin, C. J. and Wilkes, N. S., 1984, Experimental and modelling studies of annular ffow in the region between flow reversal and the pressure drop minimum. Int. Symp. Annular and Dispersed Flows, Piss, Italy. Imura, H., Kusuda, H. and Funatsu, S., 1977, Flooding velocity in a countercurrent annular two-phase flow. Chem, Engng Sci. 32, 79-87. Jefferey, A. and Kakutani, T., f972, Weak nonlinear dispersive waves: a discussion centered around the Korteweg-de Vries equation. SIAM Rev. 14, 583-643. Kapitza, P. L., 1964, Waue Flow of Thin Layers of o Viscous Fluid (collected papers by P. L. Kapitza). Macmillan, New York. Krishna, M. V. G. and Lin, S. P., 1977, Nonlinear stability ofa viscous film with respect to three-dimensional side-band disturbances. Phys. Fluids 20, 1039-1044. Maron, D. M. and Dukler, A. E., 1984, Flooding and upward film flow in vertical tubes--II. Speculations on film flow mechanisms. Inc. J. Multiphase Flow 10, 599-621. Miles, J. W., 1957, On the generation of surface waves by shear flow. J. F&d Mech. 3, 185-204. Nusselt. W., 1916, Die ObserfIachenkondensation des Wasserdamphes. Z. Vu. Dtsch. Ing. 60, 541-552. Pumir, A., Manneville, P. and Pomeau, Y., 1983, On solitary waves running down an inclined plane. J. Fluid Mech. 135, 27-50. Shearer, C. J. and Davidson, J. F.. 1965, The investigation of a standing wave due to gas blowing upwards over a liquid film: its relation to flooding in wetted-wall columns. J. Fluid Mech. 22, 321-335. Shlang, T. and Sivashinsky, G. I., 1982, Irregular flow of a liquid film down a vertical column. J. Phys. 43, 459-466. Sivashinsky, G. 1. and Michelson, D. M., 1980, On irregular wavy flow of a liquid film down a vertical plane. Prog. theor. Phys. 63, 2112-2114. Smith, F. I. P., 1970, Stability of liquid film flow down an inclined plane with oblique air flow. Phys. Fluids 13, 1693-1700. Smith, L., Chopra. A. and Dukler. A. E., 1984, Flooding and upward film flow in tube-I. Experimental studies. Int. J. Multiphase Flow IO, 585-598. A. V., Preobrazhenskii, E. I. and Semenov, D. A., Solevev, 1967, Hydraulic resistance in two-phase flow. Int. them. Engng. 7, 5-63. Tougou, H., 1981, Deformation of supercritically stable waves on a viscous liquid film down an inclined plane wall with the decrease of wave number. J. Phys. Sot. Jap. SO, 1017-1024. Urinstev, A. L., 1978, Branching of solutions in the problem of wavy flow of a viscous liquid

with a free boundary.

Zh.

Prik. Mekh Tekh. Gi.z. 4, 12~%142. Wallis, G. B., 1969, One-Dimensional Two-Phase Flow, Chap. 11. McGraw-Hill, New York. Zabaras, G. J., 1985, Ph.D. Thesis, University of Houston, University Park.