Detailed hydrodynamic characterization of harmonically excited falling-film flows: A combined experimental and computational study

Alexandros Charogiannis,1 Fabian Denner,2 Berend G. M. van Wachem,2 Serafim Kalliadasis,3 and Christos N. Markides1,* 1Clean Energy Processes (CEP) Laboratory, Department of Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom 2Department of Mechanical Engineering, Imperial College London, London SW7 2AZ, United Kingdom 3Department of Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom (Received 9 February 2016; published 17 January 2017)

We present results from the simultaneous application of planar laser-induced fluorescence (PLIF), particle image velocimetry (PIV) and particle tracking velocimetry (PTV), complemented by direct numerical simulations, aimed at the detailed hydrodynamic characterization of harmonically excited liquid-film flows falling under the action of gravity. The experimental campaign comprises four different aqueous-glycerol solutions corresponding to four Kapitza numbers (Ka = 14, 85, 350, 1800), spanning the Reynolds number range Re = 2.3-320, and with forcing frequencies f w = 7 and 10 Hz. PLIF was employed to generate spatiotemporally resolved film-height measurements, and PIV and PTV to generate two-dimensional velocity-vector maps of the flow field underneath the wavy film interface. The latter allows for instantaneous, highly localized velocity-profile, bulk-velocity, and flow-rate data to be retrieved, based on which the effect of local film topology on the flow field underneath the waves is studied in detail. Temporal sequences of instantaneous and local film height and bulk velocity are generated and combined into bulk flow-rate time series. The time-mean flow rates are then decomposed into steady and unsteady components, the former represented by the product of the mean film height and mean bulk velocity and the latter by the covariance of the film-height and bulk-velocity fluctuations. The steady terms are found to vary linearly with the flow Re, with the best-fit gradients approximated closely by the kinematic viscosities of the three examined liquids. The unsteady terms, typically amounting to 5%-10% of the mean and peaking at approximately 20%, are found to scale linearly with the film-height variance. And, interestingly, the instantaneous flow rate is found to vary linearly with the instantaneous film height. Both experimental and numerical flow-rate data are closely approximated by a simple analytical relationship with only minor deviations. This relationship includes terms such as the wave speed c and mean flow rate Q, which can be obtained by relatively simple and inexpensive methods, thus allowing for spatiotemporally resolved flow-rate predictions to be made without requiring explicit knowledge of the full flow-field information underneath the wavy interface.

A. Problem definition and motivation
Falling liquid films are convectively unstable open-flow hydrodynamic systems that exhibit a rich variety of spatial and temporal structures and a rich spectrum of wave forms and complex wave transitions between different regimes. Not surprisingly, they have been the subject of intense research during the last few decades, spearheaded by the groundbreaking work of the father-son team of the Kapitza family [1,2], and aiming to elucidate the different wave forms and sequence of spatiotemporal transitions. Aside from the fundamental theoretical interest in falling film flows, research on this class of flows is also largely motivated by the broad range of industrial and engineering applications in which they are employed, predominantly owing to their high surface-to-volume ratios and associated high heat and mass transfer capabilities even at modest flow rates. Examples of such applications include wetted-wall absorbers, heat exchangers, film condensers and evaporators, and also film reactors and microcooling schemes, to name but a few.
Extensive theoretical and experimental investigations have indicated that the rates of heat and mass transfer in film flows are strongly coupled to the spatiotemporal hydrodynamic characteristics of these films, and in particular the waviness of the liquid-gas interface [3][4][5]. Yet, despite the numerous and rather detailed experimental studies, only a limited number of studies on the simultaneous and spatiotemporally resolved variations of film thickness and velocity are currently available. This is a direct consequence of the substantial challenges associated with performing these measurements, including the restricted fluid domains under observation and the intermittent nature of the moving wavy interface. Our study is motivated by the strong demand for a detailed, efficient, and systematic simultaneous measurement of spatiotemporally resolved film thickness and velocity in falling films. Such measurements can advance our understanding of the underlying fluid dynamics phenomena, and facilitate the development and validation of advanced and accurate analytical and numerical models. Here, they are complemented by state-of-the-art direct numerical simulations (DNSs) which show good agreement with the experiments but they also guide and motivate experiments to validate some largely unexpected computational findings, thus elucidating the subtle elements of the hydrodynamics of falling liquid films. Therefore, both experiments and DNSs are instrumental in our efforts to further enhance the heat and mass-transfer rates in related flows and associated practical-technological processes.

B. Background theory
A thin liquid film flowing down an inclined plane is generally unstable to long-wavelength perturbations (much longer than the film height) above a critical Reynolds number Re, which instigate the formation and growth of long-wavelength waves on the free surface of the film [6]. The main instability mechanism between two fluids in shearing motion is the phase shift between the interface displacement and vorticity caused by inertia. The motion induced by the out-of-phase vorticity is manifested by the advection of mass into the liquid beneath the interface disturbance, acting to amplify its growth [7][8][9]. Surface tension, viscosity, and the cross-stream component of gravity have a stabilizing effect on the flow, while inertia and the streamwise component of the flow tend to destabilize it. In most practical applications, free-surface instabilities driven by natural noise (i.e., broadband noise that is naturally present in the system due to an interaction with the environment) emerge from the short-lived flat-film region with a frequency content centered narrowly around the fastest growing mode (f m ). Via a secondary modulation instability, these waves eventually evolve into a two-dimensional (2D) wave train comprising highly asymmetric solitary waves with long flat tails and steep fronts preceded by capillary ripples, and eventually break down into complex three-dimensional (3D) wave patterns [10][11][12][13].
In contrast to naturally developing instabilities, imposing a monochromatic perturbation at the inlet, e.g., by forcing the inlet flow rate, triggers, in the first instance, a rapid evolution of wave sequences which adopt the forcing frequency. The interfacial features characteristic of forced film flows at short distances from the inlet, which are of primary interest to this work, have been investigated analytically, numerically, and experimentally in many studies [14][15][16][17][18][19][20][21][22][23], from which we have the following key observations: (i) at low forcing frequencies (considerably below f m ), small waves appear near the inlet and grow quickly into solitary waves traveling over a thin substrate layer; and (ii) solitary wave growth saturates within a few wavelengths from the inlet and their shapes fluctuate modestly through interactions with their neighbors. Our experiments focus on this forcing frequency band and wave regime. At longer distances from the inlet, solitary waves interact continuously through inelastic collisions [20,[24][25][26][27]; faster moving waves absorb smaller and slower-moving ones, reorganizing as their mass grows, and eventually form bound states of two or more solitary waves traveling at the same speed [28][29][30]. Such bound states are also formed on films flowing down vertical fibers [31,32], non-Newtonian and shear-thinning films [33], but also weakly nonlinear prototypes, such as the Kuramoto-Sivashinsky and generalized Kuramoto-Sivashinsky equations [34][35][36]. For inlet forcing nearly matching the excitation frequency beyond which wave growth is inhibited (f c ), wave growth near the inlet is weak and sideband and subharmonic instabilities develop quickly. Large humps and short-separation pulses appear away from the inlet and the forcing frequency content is lost.

C. Experimental methods
Optical diagnostics allow for noninvasive visualization of one-dimensional (1D), 2D, and 3D spatial distributions of seeded or naturally occurring tracers and particles, and have thus been extensively employed in a broad range of contemporary investigations in the fields of fluid mechanics, heat transfer, combustion, and reacting flows, as well as material science, among many others. Novel optical methods have been developed and applied to liquid-film flows, while existing ones from different research fields have been appropriately adapted and implemented to the study of the instantaneous motion and spatial extent of the liquid domains under observation [37][38][39][40]. For the sake of completeness and clarity, we review below the optical diagnostic applications to film flows. But, it should be noted that this review is by no means exhaustive; instead, it is merely indicative of the considerable number of available studies and breadth of experimental practices developed towards elucidating the fundamental physics and associated underlying mechanisms of the wide class of two-phase flows.
Approaches utilizing brightness or intensity-based molecular tagging are often adopted when pursuing film-height measurements, with small quantities of a chemical fluorescent dyestuff dissolved in optically transparent liquids. Upon illumination by an excitation source, the fluorescence emitted by the dye is collected by a photodetector and translated to measurements of the local film height, following an appropriate calibration procedure. Liu and co-workers [14][15][16][17], for example, investigated the evolution of interfacial instabilities in aqueous-glycerol solution films by calibrating the fluorescence intensity emitted by the dye-doped liquid bulk. Using rhodamine 6G and the same measurement approach, Alekseenko and co-workers conducted experiments over a broad range of unsteady gasliquid flows, such as downwards annular-flows [41][42][43] and gas shear-driven flows [44]. Focusing on the dynamics of the larger waves (referred to some times as "disturbance waves"), their aim was to uncover the underlying mechanism that potentially links these waves to liquid entrainment phenomena.
Planar laser-induced fluorescence (PLIF) imaging, which involves illumination and visualization of a planar (2D) section of the flow, has been preferred by some researchers to intensity-based fluorescence imaging, as it allows for spatiotemporally resolved velocity (along with the use of appropriate velocity tracers) or temperature (adiabatic, heated, or cooled flows) information to be recovered simultaneously with the film-height data. Rhodamine-B PLIF, for example, was employed alongside infrared imaging in a study of spatiotemporal film-height and heat transfer coefficient (HTC) variations in heated films flowing down an inclined heated flat foil [45][46][47]. Spanning a wide range of flow conditions (Re, forcing frequencies and heat fluxes), these experiments showed that despite the fact that a steady-flow analysis does describe well the correlation between the mean film height and mean HTC reasonably well, significantly higher (by up to a factor of 2) local and instantaneous HTC enhancements ensue for thinner films relative to the predicted flat-film heat-transfer performance. It should be noted, however, that this is not the only study where the employment of rhodamine-B was favored for thermographic purposes [48]. This is due to its excellent temperature sensitivity even at low (near-room) temperatures, as well as a high fluorescence efficiency which promotes the employment of rhodamine-B PLIF in a plethora of film-flow investigations, such as downwards cocurrent gas-liquid annular flows [49] and horizontal liquid-liquid flows [50,51], allowing for a clear identification of liquid-gas and liquid-liquid interfaces. Therefore, rhodamine-B PLIF has been selected here as it is considered the most suitable tracer-diagnostic combination towards meeting our objectives, while allowing for the simultaneous employment of a suitable, optical velocimetry technique.
Optical velocimetry techniques can broadly be grouped depending on whether particle seeding or molecular tagging is employed. The photochromic dye activation (PDA) technique is a prime example of the latter; thin fluorescing traces are generated perpendicular to the direction of the flow and imaged at different time delays using a high-speed visualization setup yielding direct velocity-profile measurements. The PDA technique has been implemented in a plethora of two-phase flow visualization experiments, including countercurrent annular flows [52], slug and stratified flows [53], free-falling films down a vertical tube [54], laminar films flowing down an incline [55], and falling films with countercurrent shear-gas flow [56].
Laser Doppler velocimetry (LDV), on the other hand, which belongs to the particle-seeding category of optical velocimetry techniques, has also been frequently employed. This technique relies on local measurements of the Doppler shift between the laser light incident on particles seeded within the probe volume and the light scattered from them. It was adopted by Mudawar and Houpt [4] in mass and momentum transport investigations in thick aqueous propylene-glycol films, and by Dietze and co-workers [57,58] in backflow and flow separation studies in inclined and vertically falling laminar films. In the latter case, negative velocities in the laboratory reference frame ensued in the capillary wave region, confirming the presence of local backflow.
In contrast to the aforementioned techniques, particle image velocimetry (PIV) and particle tracking velocimetry (PTV) allow for 2D and 3D flow-field characterization. Sequential illumination and imaging (within a small time delay) of the scattered or fluorescence signal generated by particles seeded in the flow allows for tracking of the motion of particle groups (PIV) or individual particles (PTV), typically using a multipass cross-correlation approach during image post processing. A comprehensive review of contemporary particle-based velocimetry techniques can be found in the book by Raffel et al. [59]. It is worth mentioning that PTV typically offers superior resolution at the cost of increased processing time, while suffering less from bias errors in the presence of strong velocity gradients within the examined fluid domain. This effect becomes particularly pronounced in the near-wall region of film flows, with PTV or micro-PIV being considered superior options to PIV.
Schubring et al. [60] employed micro-PIV to assess the validity of the universal velocity profile (UVP) in horizontal, gas-shear-driven liquid flows, while Adomeit and Renz [61] investigated the hydrodynamics of vertically falling laminar films using simultaneous fluorescence imaging and micro-PIV. In the latter study, strong deviations from the steady-flow velocity profile predictions were observed depending on film topology, despite the presence of parabolic profiles throughout the investigated Re range (Re = 27-200). Zadrazil and Markides [62] recently studied downwards cocurrent gas-liquid annular flows, also by employment of PLIF, PIV and PTV. In this case, PTVderived velocity profiles revealed multiple recirculation zones inside the disturbance waves (in a reference frame moving with the wave speed), that the authors hypothesized might be linked to liquid entrainment as well as mass and momentum transfer from the near-wall zone towards the liquid-gas interface. Also using micro-PIV, Ashwood et al. [63] provided mean velocity-profile data in cocurrent upward flows in a rectangular duct, and suggested that modifications to the UVP must be implemented to accurately track their results. Finally, by applying micro-PIV, Dietze et al. [58] linked changes in the surface curvature and ensuing adverse pressure gradients to backflow in the capillary trough closest to the solitary wave crest, and flow deceleration in the second capillary trough. The abundance of information that can be extracted from the application of PIV and PTV in resolving the flow field of interest within a 2D domain led us to the selection of this technique, implemented simultaneously alongside PLIF, for the measurement of film thickness and velocity.

D. Aims and objectives
One of our main aims is to provide a comprehensive account of the link between the interface height and film thickness and its fluctuations, and the liquid mass-transfer rate in the film (in the flow direction) and its fluctuations. Towards this end, a systematic study has been undertaken comprising simultaneously conducted local and instantaneous film-height and flow-rate measurements complemented by advanced numerical simulations, which builds upon two previous contributions. The first was dedicated to establishing an experimental methodology capable of providing simultaneous and spatiotemporally resolved film height and velocity data in harmonically excited film flows (Ref. [64]). Along with the development of the combined (PLIF-PIV-PTV) optical technique, an in-depth error analysis was carried out and presented in this earlier work. The second contribution explored various aspects of solitary wave dynamics by bringing together experimental results and computational studies: both low-dimensional modeling and DNSs. Direct comparisons between experimental and computational results were pursued with respect to topological characteristics and velocity fields for selected flow conditions, while the roles of flow recirculation, separation, and reversal were investigated and elucidated in our very recent study by Denner et al. [65].
Our experimental campaign comprises four Kapitza-number liquids (Ka = 14, 85, 350, and 1800), Reynolds numbers spanning the range Re = 2.3-320, and forcing frequencies f w = 7 and 10 Hz. Of these two important dimensionless variables, Re represents a quantitative description of the relative significance of inertial to viscous forces and is defined in terms of a characteristic flow velocity U , a characteristic length scale l, and the fluid kinematic viscosity ν f : while Ka is a ratio of surface tension to inertia, and acts as an indicator of the hydrodynamic wave regime: where γ f denotes surface tension, ρ f density, β inclination angle, and g gravitational acceleration.
It is noteworthy that unlike Re which is a flow-control parameter, Ka depends only the physical properties of the liquid, gravity, and inclination angle and hence it is fixed for a given liquid. Following a brief description of the experimental arrangement (Sec. II) and the experimental and numerical methodologies (Secs. III and IV) employed in this work, the topological characteristics of the films under investigation are examined in Sec. V A, followed by results on the steady and unsteady components of the local and instantaneous flow rate in Sec. V D. The latter are linked to the local and instantaneous film-height characteristics based on both experimental and numerical data from the same flow conditions (Sec. V E). Concluding remarks and a summary of our results are offered in Sec. VI.

A. Flow loop
The experimental setup and methodologies developed to retrieve film height and velocity data have been described in detail in a previous study (Charogiannis et al. [64]), and briefly in Ref. [65], and therefore only a synopsis is presented here for clarity and completeness. A schematic of the experimental configuration, itself an evolution of a previously developed setup [46,47], is presented in Fig. 1 collecting and resupplying the liquid to the test section. A settling chamber comprised of meshes and a honeycomb is installed at the test-section inlet, serving a dual purpose: (i) to break down any secondary flow and dissipate any large-scale eddies, thus producing a well-characterized uniform flow; and (ii) to introduce this over the test section with uniform thickness. The liquid is then collected in a tank, passed through a heat exchanger to remove any heat imparted upon it by the pump, and directed to the flow pulsation arrangement (V3, V4, and oscillator valve). There, the flow is split forming a steady (Q) and a pulsating supply (Q ).
Flow pulsation is instigated using an oscillator valve, which consists of a metallic spindle with a hole drilled transversely through it, rotating within the valve body. A stepper motor and dedicated electronic circuitry are used to drive the spindle at the desired rotational speed. Using this setup, accurate control of the wave frequency and wave amplitude can be achieved, selectively triggering the rapid growth of regular waves and fully developed wave regimes within a short distance from the flow inlet. The mean flow rate and pulsation intensity are calculated from the output signal of an ultrasonic flowmeter that measures the instantaneous flow rate into the settling chamber. In more detail, the pulsation intensity A is obtained by dividing the measured forcing flow-rate amplitude by the mean flow rate. This parameter could be controlled to within a few % points, and was varied in this work between A = 5% and 10%. It should finally be noted that a series of separate experimental runs have shown that the wave characteristics at the measurement location, such as the wave-crest height, trough height, and mean film height are not influenced when varying the pulsation intensity in the range A = 5%-15%.
As noted earlier, our experimental campaign employs two forcing frequency and forcing amplitude combinations: (i) f w = 7 Hz and A = 10%, and (ii) f w = 10 Hz and A = 5%. Our excitation setup allows for disturbances of up to ≈11 Hz to be imposed by forcing the inlet flow rate periodically; the higher frequency f w = 10 Hz is, therefore, selected based on the capabilities of our excitation mechanism. It should also be noted that this frequency is considerably lower than the cutoff frequency f c , above which the film remains flat; based on the experiments of Liu et al. [17] pertaining to flows with β = 4 • , Ka = 2340 and for Re 25, film flows are always unstable for f w 10 Hz (see Ref. [13]), while in the case of water film flows falling down a vertical plane, the simulations of Nosoko et al. [19] show that film flows are unstable even at forcing frequencies exceeding 100 Hz.  Finally, while conducting a separate set of experiments during which the excitation frequency was set below f w = 7 Hz, we observed the generation of parasitic pulses when f w 4 Hz, which imposed a lower bound for the selection of our excitation frequency. Other studies have also observed that below a threshold frequency, the regular wave pattern generated by forcing the flow at the inlet is destroyed by the emergence of parasitic waves (see, for example, Refs. [5,17,21]). As will become apparent later, the disparity between the two frequencies employed here is sufficient in order to generate film flows with sufficiently different film-height and flow-field statistics.

B. Test section
The main test section comprises a 0.7-mm-thick, 400 mm × 285 mm rectangular soda-lime glassplate installed on an aluminium support frame and inclined at β = 20 • to the horizontal. The inclination angle was not altered throughout the course of this experimental campaign. The settling chamber outlet is equipped with a knife edge, the height of which can be adjusted by a pair of micrometer stages in order to regulate the flow contraction and prevent the generation of a hydraulic jump or air backflow. The Reynolds number Re can be evaluated at the flow inlet using the mean flow-rate measurement per unit width of the channel (Q s ), obtained from an ultrasonic flowmeter: The liquid phase (aqueous-glycerol solutions; see Sec. II C) was seeded with rhodamine-B dye at a concentration of ∼0.5 g/L, while glass hollow spheres (11.7-μm mean diameter) were seeded at a concentration of ∼0.18 g/L for tracking the fluid motion. The Stokes number St (ratio of particle motion characteristic time scale to fluid motion characteristic time scale) is estimated to be of the order of 10 −3 for the majority of the examined flows; hence, it can be assumed that the particles follow the fluid motion faithfully.
The dye-doped liquid flowing over the glass plate was illuminated and excited by a laser light sheet from the unwetted (under)side to avoid spatially and temporally nonuniform beam stirring and lensing from the film's free surface (Fig. 2), with imaging also carried out from underneath the plate to limit image distortions. A frequency-doubled double-cavity Nd:YAG laser (532 nm) was used at a rate of 100 Hz per laser head, allowing a maximum PLIF temporal resolution of 5 ms. Sheet optics were employed in forming a thin (∼100 μm) laser sheet extending along the TABLE I. Experimentally determined ρ f , ν f , and γ f of the four aqueous-glycerol solutions employed in our experiments, presented alongside the examined liquid flow-rate ranges and resulting Re and Ka numbers.  (32 mm) were employed for imaging. The PIV camera was equipped with a short-pass filter (with a cutoff wavelength at 550 nm), thus making the visible spectrum below the cutoff wavelength available for calibration, and a Scheimpflug adaptor as it was positioned at a slight angle (13 • -17 • ) relative to the excitation plane. A long-pass optical filter with a 540-nm cutoff was installed on the PLIF camera, thereby only allowing the red-shifted fluorescence emission to be collected. Both cameras and lasers were operated using the LaVision DaVis software, and synchronized by a LaVision High-Speed Controller (HSC). In order to vary the Kapitza number Ka, while avoiding any restrictions being imposed on fluid selection (see, for example Refs. [49,66]), an alternative approach to the employment of a typical correction box was deemed necessary. A Perspex box with the surface that would otherwise adhere to the glass plate removed was positioned at the imaging plane using a micrometer stage and filled with the tested fluid through a hole. The imaging planes of the two cameras were then mapped using a calibration graticule immersed inside the box and aligned with the excitation plane. Optical distortion corrections for the two-camera setup were performed using a pinhole model, with a resulting apparent magnification varying between 28.0 μm/pixel and 29.7 μm/pixel, and a root mean square (RMS) positional error between 0.5 and 0.9 pixels for the four Ka flow cases considered in this paper.

C. Fluid properties
The density ρ f , viscosity ν f , and surface tension γ f of the four aqueous-glycerol solutions were determined experimentally, and are presented in Table I. The glycerol concentrations, starting from the solution with the highest content, correspond to 82% (S1), 65% (S2), 45% (S3), and 21% (S4) by weight. Along with the liquid properties and their associated measurement errors, the examined liquid flow-rate ranges and resulting Re and Ka numbers, as well as the errors associated with their calculation, have also been included. It should also be noted that samples were collected throughout the duration of optical measurements, which confirmed that the surface tension of the employed aqueous solutions did not vary as a consequence of the fluid exposure to the laboratory environment and the possible deposition of dust particles acting as surfactants. The observed surface tension values were, however, lower than those typically found in the literature, most probably owing to the presence of the PIV particles and the rhodamine dye, known to reduce the surface tension of water-rich solutions. All measurements were conducted at 25 • C ± 0.5 • C.

A. PLIF and PTV data acquisition and processing
Sample PLIF and PIV frames from all four examined Ka liquids (following perspective distortion corrections and the removal of any reflections) are presented in Fig. 3. The identification of the two interfaces (liquid-gas and liquid-solid) from raw PLIF images involves a series of challenges. Near the liquid-solid interface, reflections from the substrate blur the PLIF signal locally, while at the liquid-gas interface the fluorescence emitted by the illuminated liquid volume is reflected about the interface and back to the camera. The reflection intensities between the wave crests and troughs are 014002-8 DETAILED HYDRODYNAMIC CHARACTERIZATION OF . . . typically higher compared to other regions of the flow due to lensing effects stemming from locally high film curvature. In fact, with increasing film curvature the intensity of out-of-plane reflections increases to the point where accurate identification of the interface becomes locally difficult.
To systematically identify the location of the two interfaces, different methods where devised and tested using in-house MATLAB processing routines. In particular, the location of the liquid-solid interface was obtained using an edge-detection algorithm (similarly to Ref. [67]) whereby the signal level corresponding to the onset of the liquid film was calibrated against the graticule. For the liquid-gas interface, the intercept between a linear fit to the maximum signal gradient and a linear fit to the reflection intensity profile on a per image-column basis was employed as the best estimate of this liquid boundary. A detailed account of our interface-identification methodology, including the approach used for the removal of reflections, can be found in Ref. [64].
Both CMOS cameras were operated in dual-frame mode, so that for every PIV frame, a corresponding PLIF frame was collected and processed in the manner discussed above. The interframe separation between two successive acquisitions was varied in the range 0.45-1.50 ms, depending on flow conditions (Re and Ka), so that the corresponding particle displacement in the near-interface regions was in the range 8-15 pixels. Unlike the PLIF images where reflection intensities are substantially lower compared to the original emission, reflections in the particle images are comparable to the original signal level, thus rendering the identification of secondary scattering sources highly nontrivial. In order to address this, following full processing of the PLIF frames and retrieval of the liquid-gas interface, binarized images (signal level equal to unity ascribed to the illuminated liquid domain and zero to the rest of the image) were used to mask out secondary scattering regions in the raw PIV images (regions outside the illuminated liquid bulk). The final, masked particle images were used to generate 2D velocity-vector maps by employment of a four-pass cross-correlation approach using LaVision DaVis. For the first and second passes, a (32 × 32)-pixel interrogation window with 50% overlap was selected, while for the remaining two passes the PIV window size was reduced to 16 × 16 pixels. The resulting PIV vector-to-vector spatial resolution corresponds to 8 pixels (amounting to 222.4-236.6 μm). Individual particles were subsequently tracked (PTV calculation) based on the employment of the PIV results as reference estimators of the velocity field.
Along with a diminished spatial resolution compared to PTV, PIV also suffers from bias errors in the presence of velocity gradients (shear), as the particle-group motion is inevitably averaged within the interrogation window. As a consequence, significant systematic errors arise in the nearwall regions of film flows, which are further pronounced due to the fact that local particle-image displacements are small. A detailed account of the relevant shortcomings of PIV can be found in Refs. [68,69]. Unlike PIV, PTV data contribute significantly lower velocities in these regions, rendering the latter the preferred measurement method; it should, however, be noted that the presence of particle-image reflections within the masked liquid domain constitutes an error source affecting both PIV and PTV velocity field calculations. For PTV, this effect manifests as a bias towards higher velocities, especially when time or space averaging over larger domains, while for PIV it is indistinguishable from the aforementioned error. Such uncertainties that affect near-wall velocity measurements and originate from optical effects are often (and, inevitably) encountered in the literature (see, for example, Refs. [68,69]), with an importance that depends on the flow facility and optical setup. In our present experiments, the near-wall region, which is expected to be more strongly affected by these limitations, extended up to about 5%-7% of the total spatial extent of the film, and consequently, the impact on bulk velocity and flow-rate measurements is minimal.

B. Wave and phase-locked flow-field averaging
To examine the velocity field underneath the wavy interface in detail, and to recover flow-rate data with high spatial resolution, we follow the standard practice of averaging instantaneous PTV velocity vector maps. The primary challenge associated with this approach is to identify which images and which location in these images correspond to the same spatial region of the flow (i.e., position underneath the periodic wave). The desired topology is selected first by identifying a suitable (instantaneous) reference film-height (interfacial wave) profile, which is then cross correlated with all available film-height traces from the images in the same data set [ Fig. 4(a)]. Trace pairs satisfying a maximum displacement condition (80 pixels, or less) are then repositioned using the crosscorrelation results, and finally averaged. Standard deviations calculated along the reference signal serve as deviation estimators of the implemented methodology [ Fig. 4(b)], peaking in this case at approximately 60 μm (3.8% of the mean film height and 2.6% of the maximum) in the region between the wave crest and trough. As expected, this particular region displays the highest deviations as it typically contains the strongest gradients. In the other regions, the standard deviations seldom exceed 15 μm (below 1% of the mean height). The same translation values are finally used to rearrange and average out the corresponding instantaneous PTV velocity maps, producing averaged velocity fields such as the one displayed in Fig. 4(c). It is essential here to highlight a study by Alekseenko et al. [70], whereby the evolution of waves down on a vertical plate was investigated by simultaneous application of the shadow method (not to be confused with the shadowgraph method that relates localized refractive-index variations in transparent media to localized temperature, pressure, or composition variations), for determining the film height, and stroboscopic particle visualization, for velocity measurements. This particular technique allowed, for the first time, the retrieval of instantaneous velocity profiles along the wavy interface, and constitutes the basis of the current study as well.

C. Validation experiments
Aside from the measurement runs carried out in harmonically excited films, experiments were also conducted in flat, unperturbed, viscous films (Ka = 14-20) to assess the validity and accuracy of our combined optical methodology. First, the mean (time-averaged) heights of these films were compared to predictions from the 1D, laminar, fully developed, and steady solution of the Navier-Stokes (NS) equations under the assumption of negligible inertia, also known as the Nusselt solution [71]. The analytical Nusselt-solution expressions for the film height h N as a function of the total flow rate Q and the flow Re, the liquid bulk velocity U N , and the axial-velocity profile U x as a function of the distance from the wall y, are given in Eqs. (4)-(7), respectively: In these expressions, β, s, and ν f denote inclination angle, film span, and liquid kinematic viscosity, respectively.

014002-11
The mean absolute film-height deviation from the Nusselt results was, for a total of six different flat films, 23 μm (corresponding to a mean relative deviation of 0.9%). Direct comparisons to micrometer-stage measurements where also pursued, with the mean-absolute and RMS deviations amounting to no more than 10 and 30 μm, respectively, and corresponding approximately to 1% and 3% of the thinnest film examined in this set of experiments (mean film heights between 1 and 3 mm). Thinner films were, however, sometimes encountered during the course of the experimental campaign, and consequently a worst-case measurement uncertainty of 2% is believed to be a more representative estimate of the error associated with average film-height measurements (time-averaged and wave and phase-lock averaged).
Furthermore, an estimate of the random film-height measurement error, originating from camera noise, can be obtained by calculating standard deviations over a wide range of flat-film runs. The latter were found to be of the order of the spatial resolution, around 30 μm. The mean uncertainty associated with the instantaneous film-height measurement is therefore estimated at a maximum 40 μm, which for the range of examined film flows corresponds to errors between 1% and 10% depending on local film topology (greater where the film is thinner). As the mean film height of all flows examined in our experiments amounts to approximately 1.5 mm, a mean relative error of 2.5% is quoted for the instantaneous film-height measurements.
In an effort to validate the PTV-derived interfacial and bulk velocities, comparisons between measurements carried out in flat films and their complementary Nusselt calculations were performed; in both cases, the average relative deviations were of the order of 3.2%. In our previous discussion, we noted that both PTV and PIV near-wall measurements are affected by local reflections, which give rise to locally overestimated velocities. With respect to the flows examined for validation purposes, this phenomenon was apparent in regions extending up to approximately 100-200 μm away from the wall or, on a relative basis, a spatial extent of approximately 5%-7% of the film height. Despite the fact that the ensuing effect on bulk-velocity measurements is quite modest, an effort to alleviate this artifact in the analysis of wavy flows was made (see Sec. V C 2).
A final test aimed at assessing the effectiveness of the combined optical technique in recovering spatiotemporally resolved local and instantaneous streamwise flow-rate results (product of the PLIF-derived film height and PTV-derived bulk-velocity data) and at obtaining an estimate of the relevant measurement error, was conducted by comparing the results with independent flowmeter measurements. Assuming that any 3D out-of-plane flow components are negligible, the time-averaged flow rate measured at the excitation plane is representative of the time-averaged flow rate at all locations along the film span. The mean relative deviation between time-averaged flow-rate measurements carried out by combined PLIF-PTV and the ultrasonic flowmeter for a total of six flat and nine wavy flows included in this test was 1.6%, with no systematic trends identified. As this value is smaller than the flowmeter uncertainty quoted by the manufacturer (3%), the latter is used as the estimate of our mean flow-rate measurement error. Figure 5(a) shows local and instantaneous flow-rate data (from PLIF-PTV) normalized by mean flowmeter measurements Q fl , plotted against local and instantaneous film heights normalized to their respective Nusselt heights h N , for flat films corresponding to Ka = 14 and Re = 5.8, 4.9, 4.2, and 2.6. Effectively, the deviations discussed earlier (mean flow rate and mean film height recovered by PLIF and PTV) are presented alongside any random uncertainties associated with these measurements. The latter have been isolated in Fig. 5(b), where the same (PLIF-and PTV-derived) data are normalized to their respective mean values. The flow-rate standard deviation for the displayed results is 1.3%, an estimate of the relevant random error.
The larger scatter observed in Fig. 5(a) originates from the use of flowmeter measurements Q fl to normalize the PLIF-and PTV-derived flow-rate data, and the use of equivalent Nusselt heights h N to normalize the PLIF-derived film-height data (as opposed to the actual mean flow rate). Thus, the uncertainties associated with the measurement of Q fl , as well as those linked to the calculation of h N (mainly from the flow-rate, viscosity, and density measurements), all contribute to the scatter observed in this plot. In turn, the degree of scatter reflects the variations in the flow-rate and film-height values, and highlights the absence of significant bias in the measurements. It should finally be noted that the region with the highest probability (dark-red region) arises due to the significant overlap of two out of the four data sets employed in this assessment.

IV. NUMERICAL METHODOLOGY
The numerical methodology applied in this study is based on the general-purpose finite-volume framework for the simulation of interfacial flows of Denner and van Wachem [72,73]. The underlying numerical framework, whose application extends beyond film flows to general interfacial flows, resolves the dynamic behavior of the entire two-phase system consisting of the liquid film, the liquid-gas interface, and the gas phase.

A. Numerical framework
The incompressible flow of the immiscible two-phase system is governed by the continuity equation and the momentum equations where t represents time, x ≡ (x,y,z) denotes a Cartesian coordinate system with x the streamwise direction, y is cross-stream direction, and z is the transverse direction [we adopt the Einstein notation in Eqs. (8) and (9)], u is the velocity, p is the pressure, ρ is the density, μ is the dynamic viscosity, and f σ is the volumetric surface force acting at the liquid-gas interface.
The governing equations of the flow are solved as a single linear system of equations. The momentum equations [Eq. (9)] are discretized using a second-order backward Euler scheme for the temporal derivative of the velocity, and central differencing for the convection, diffusion, and pressure terms (see Ref. [72] for details). The continuity equation is discretized using a specifically designed momentum-weighted interpolation method proposed by Denner and van Wachem [72].

B. Two-phase methodology
The volume-of-fluid (VOF) method [74] is adopted to capture the interface between the immiscible gas and liquid phases. The local volume fraction of both phases in each mesh cell is represented by the color function φ, defined as φ = 0 in the gas phase and φ = 1 in liquid phase, with the interface located in mesh cells with a color function value of 0 < φ < 1. The local density ρ and viscosity μ are calculated based on the local color function φ as respectively, where subscript g denotes properties of the gas phase and f denotes properties of the liquid film, and the advection of the color function φ is governed by the linear hyperbolic equation based on the underlying flow with velocity u. Equation (12) is discretized using the compressive VOF methodology proposed by Denner and van Wachem [75], based on the CICSAM scheme [76]. The surface force per unit volume is described by the continuum surface force (CSF) model [77], assuming a constant surface tension coefficient σ and neglecting mass transfer between the bulk phases, as where κ is the interface curvature, which translates the surface force resulting from surface tension into a volume force which can be discretized in a finite volume framework. A balanced-force implementation of the surface force is assured by discretizing Eq. (13) on the same computational stencil as the pressure gradient [72]. Aliasing errors in the evaluation of the interface normal vectorm = ∇φ/|∇φ| and the interface curvature κ = − ∇ ·m are diminished by convoluting the volume fraction φ using of a cosine convolution kernel with a support of ε = 3 x, where x is the mesh spacing; details are given in Refs. [65,72,78]. No convolution is applied to smooth the surface force [79], and numerical artifacts that manifest in parasitic or spurious flow features near or at the interface are addressed by an artificial viscosity model [65,80].

C. Numerical setup
The 3D computational domain has dimensions 0.256 m + 50h N × 4h N × 0.1h N and is represented by a Cartesian mesh with a resolution of 10 cells per Nusselt film height h N . Flow-rate and film-height data are taken at downstream position x = 0.256 m at every time step. The distance x = 0.256 m corresponds to the furthermost downstream location that was optically accessible in our experimental apparatus. Prior to conducting our experiments, we verified by both PLIF imaging and computations that for film flows excited at 7 and 10 Hz, wave growth has already seized at the measurement location. The numerical time step applied in the presented simulations satisfies a Courant number of Co = t u |u|/ x 0.25, as well as the capillary time-step constraint proposed by Denner and van Wachem [73]. A no-slip condition is enforced on the substrate and a free-slip boundary condition is applied at the top (gas-side) liquid boundary. A monochromatic forcing is imposed at the domain inlet by periodically varying the flow rate with frequency f and amplitude A around a mean value, chosen to match the experiments for any given flow case. At the domain inlet, a semiparabolic velocity profile is prescribed for the liquid phase: and a spatially invariant velocity is prescribed for the gas phase The film height at the inlet is constant h(x = 0) = h N and the domain outlet is modeled as an open boundary following Nosoko and Miyara [19] and Denner et al. [65], which assures that the flow can leave the domain with minimal reflections. Initially, at time t = 0, the velocity field is fully developed and the film is flat. Results from the complementary DNS effort described here are included in Sec. V E, while a more detailed description of the numerical setup can be found in Ref. [81].

A. Topological features
Here, we outline and discuss the main interfacial features observed across the range of examined flow conditions as a function of Ka, Re, and f w , both qualitatively and quantitatively. Using the wave and phase-locked averaging methodology introduced earlier, the shape of the interface around the crests of solitary waves is shown in Fig. 6 for two excitation frequencies (7 and 10 Hz), and flows with Re = 12, 21, 25 for a fluid with Ka = 85, Re = 37, 54, 76 for Ka = 350, and Re = 95, 195, 317 for Ka = 1800.
As was noted in the Introduction, solitary waves consist of a main hump with a steep front preceded by capillary ripples, and a long, gently sloping tail. Increasing the Re for a fixed Ka fluid results in an increase in the wave-crest height (maximum film height) and mean film height. The wave-crest height ceases to grow following the onset of flow recirculation with respect to a frame of reference moving with the wave speed, and once the interface velocity increases beyond the wave speed at sufficiently high inertia [65]. The wave-trough height (minimum film height) initially decreases as the Re of the film flow is increased, and then quickly settles to a near-constant value which is independent of f w . Therefore, the normalized wave amplitude, which is defined as the difference between the wave crest and trough height divided by the mean film height, also increases with increasing Re. Forcing at 7 Hz rather than 10 Hz (for the same Re) results in a decrease in the number of capillary ripples and an increase in the wave amplitude for the films with Ka = 85 and 350 films, while the Ka = 1800 films appear unaffected. In fact, wave topologies seem near-identical irrespective of the imposed f w , an observation that is also supported by the statistical results depicted in Fig. 7. In addition, it can be observed that in the flows with Ka = 85 and 350, the film ahead of the capillary-wave region is always thinner and flatter for the films forced at f w = 7 Hz compared to those forced at f w = 10 Hz.
In calculating statistical quantities such as the mean film height h and standard deviation of the film height σ h , 50 local and instantaneous measurements were selected uniformly, at random, and without replacement from each PLIF frame. The (absolute) mean film height h was found to increase monotonically with the flow Re, following a similar trend to that suggested by Nusselt [71] and reported in many previous studies [37,55,[82][83][84][85].
It should, however, be noted that when considering the (relative) mean film height normalized by the Nusselt height h N , a variety of experimental studies have reported different results and trends. For example, Moran et al. [55] observed h > h N for laminar falling films down a 45 • incline (Re = 11-220), as did Takahama and Kato [86] for the range Re below Re ≈ 100-200, and h h N for higher Re. Finally, Karapantsios et al. [85] have reported that h > h N for Re = 509-13 090. Evidence stemming from both experimental and computational studies [87,88] suggests that large-amplitude waves preceded by capillary ripples reduce the wall shear stress, on average, in comparison to an undeformed interface (Re = 10-100); consequently, fluid particles travel faster, on average, in a wavy film flow of the same average flow rate as a flat film. In order for the average flow-rate to be conserved, the mean film height must then reduce (h h N ). Given the large number of experimental studies performed at higher Re, however, a reversal of this relation may be expected past the onset of transition to turbulence. In our case (Fig. 8), we observe h ≈ h N for Ka = 85 and f w = 10 Hz, h h N for Ka = 85 and f w = 7 Hz, h < h N for Ka = 350, and h h N for Ka = 1800. Given our estimated error of 2% for the mean film-height measurement, however, none of the examined film flows exceed h N beyond experimental error.
The film-height standard deviation σ h , which represents the variation of film heights about their mean value, can be employed as a measure of the (absolute) film waviness. Results are plotted in Fig. 7(a) for the three Ka and two f w flows examined in this experimental campaign, as a function of Re. The following remarks can be made: (1) Among the Ka = 85 and 350 films, those forced at f w = 7 Hz display significantly higher σ h values compared to those excited at f w = 10 Hz. This trend is not, however, evident in the Ka = 1800 films.
(2) For the Ka = 85 films, and the Ka = 350 films below Re ≈ 70, σ h rises quickly and monotonically with Re. In contrast, for the Ka = 1800 films, σ h remains approximately constant over the range Re = 90-320. The data of Moran et al. [55] on unforced films falling down a 45 • incline (Ka = 18, Re = 11-220) agree with our low Ka result, both in terms of trends and σ h magnitudes. In their experiments, σ h saturates beyond Re ≈ 120, while the experiments of Drosos et al. [82] indicate that σ h becomes independent of Re past Re ≈ 200 (vertically falling, naturally developing films, Ka = 3500, 2240, and 1450).

014002-17
(3) The ratio σ h /h, also referred to as the (relative) film roughness or the coefficient of variation, decreases throughout the range of examined Re for the Ka = 1800 data set [ Fig. 7(b)]. The highest σ h /h values amount to approximately ∼50% in the range Re = 70-80.

B. Axial velocity profiles
By implementing the experimental methodologies and processing steps described above, the flow field underneath the wavy liquid-gas interface can be scrutinized over a broad range of harmonically excited films. Direct comparisons can then be made to analytical estimations based on Nusselt theory by considering corresponding Nusselt velocity profiles based on the measured local film heights. A similar analysis has been attempted by other researchers (see, for example, Refs. [4,54,55,61]), although to our knowledge mostly concentrating on unforced film flows (see Ref. [70] for harmonically excited films), and axial-velocity profile comparisons to analytically derived ones based on the (quasisteady, varying based on the local film thickness) Nusselt solution to the NS equations. Using the PLIF data of the local film height h and Eq. (7), experimentally measured velocity profiles are compared to analytically derived counterparts at Locations 1-15 along the wave in Fig. 9. The results, which are provided for demonstration purposes, were generated by means of the wave and phase-locked averaging methodology and relate to a flow with Ka = 85, Re = 21, and f w = 10 Hz. The following observations can now be made: (1) The analytical predictions overestimate the experimental results significantly, close to the wave crest (Location 1); e.g., the interfacial velocities corresponding to Profiles 1-3 are overestimated by more than 50%.
(2) Many of the profiles obtained at Locations 6-15 exhibit nonparabolic shapes, such as Profile 7 (wave front), 8 (wave trough), and 12 (capillary wave region). The deviations between the Nusselt and experimental profiles are, however, considerably smaller compared to ones occurring near the wave crest. These observations were discussed extensively alongside direct comparisons to DNS in our previous study [65], and the relevant findings will therefore not be repeated here. It should be noted, however, that despite the occurrence of locally nonparabolic velocity profiles, the flow rate along the streamwise direction as predicted by the steady flow analysis at the same locations is expected to match closely the experimental values, (the local flow-rate corresponds to the integral of the local velocity profile along the local film height, multiplied by that film height).

Space-resolved flow-rate measurements
The preceding comparison between the experimentally obtained local axial-velocity profiles and corresponding analytical estimates based on the (quasisteady) Nusselt solution can be expanded to allow for local flow-rate comparisons over the range of different flows. In more detail, each local, wave and phase-locked averaged axial-velocity profile along an averaged PTV map, such as the ones previewed in Fig. 9, was first integrated using the trapezoidal rule and multiplied by the local film height and the known film span (28.5 cm). In addition, using the experimental local film-height values, corresponding Nusselt bulk velocities were then calculated by using the Nusselt relationship in Eq. (6), yielding analytical flow-rate predictions that are used here as a basis for comparison with the experimental results. Sample results are shown in Fig. 10, for Ka = 85; Re = 21, Ka = 350; Re = 54, and Ka = 1800; Re = 196, all for forcing at f w = 7 Hz. It should be noted that despite the fact that different Ka and Re flows were selected, the mean flow rates Q N of all three flows are very similar (to within 5%).
Referring to Fig. 10(a), deviations between experiments and theory for the Ka = 85 film flow are small ahead of the solitary wave and at the wave tail. The mass-carrying capacity of the wave is, however, significantly overestimated owing to the overpredicted bulk velocities in the vicinity of the wave crest (see Fig. 9). The peak deviation, which is observed at the crest height, amounts to around 70%, while in the capillary wave region, as well as the flat-film region ahead of it, the analytical results are typically 20%-25% lower than the experimental ones. It should be noted that these trends were examined in greater detail in a previous contribution [64], which also included comparisons to the data provided by Adomeit and Renz [61], and Moran et al. [55].
For flows with higher Ka and Re numbers of 350 and 54, respectively [ Fig. 10(b)], the aforementioned trends persist and the deviations near the wave crest grow to beyond 100%, while those in the region ahead of the main wave front increase to 25%-30%. Furthermore, for the Ka = 1800, Re = 196 film flows [ Fig. 10(c)], the analytical results overestimate the experimental obtained ones in nearly all locations along the solitary wave, with the maximum deviation growing to nearly 180%. It is very interesting and rather surprising that different locations along the wave with the same film height appear to carry the same mass along the streamwise direction of the flow, even if the local velocity profiles at those locations have different velocity profiles. For example, the 014002-19 film near x = 0 in Fig. 9 (wave tail) is as thick as it is between the wave crest and trough (Locations 6 and 7); the velocity profiles at those locations are, however, highly nonparabolic, unlike those at the wave tail [64,65].

Time-varying flow-rate measurements
Time-varying flow-rate data were generated by averaging the measured film heights and bulk velocities over a 100-pixel (2.8-mm) region of the flow on a per image basis. Thus, every LIF-PTV image pair contributed a single, local, instantaneous flow-rate measurement; upon averaging all individual flow rates over an entire data set corresponding to a fixed flow condition, a mean flow rate can be obtained and compared to the flow meter reading for validation purposes. These were found to deviate by only 2.1% on average over 50 studied flows.
The methodology adopted in retrieving velocity data for the purposes of this study is demonstrated in Fig. 11. Velocity profiles from four distinct locations [P1, P2, P3, and P4 in Fig. 11(a)] along a solitary wave are plotted in Fig. 11(b) (P1 and P2) and 11(c) (P3 and P4). The spatial separation between each of these four locations corresponds to 100 pixels, matching the spatial extent used (c) Velocity profiles generated for Locations P3 and P4 (wave front and first capillary wave), presented alongside the space-averaged profile between P3 and P4, and the near-wall correction generated for that profile.
for averaging. Between P1 and P2, the film is relatively flat and therefore the local space-averaged velocity profile, plotted alongside P1 and P2, matches both P1 and P2 very closely [ Fig. 11(b)]. In contrast, the topology between P3 and P4 is more complex; both P3 and P4 are nonparabolic, and the local space-averaged profile does not resemble the shapes of P3 and P4 [ Fig. 11(c)]. In order, therefore, to obtain a mean flow rate over an extended and irregular domain such as P3-P4, the latter is broken down to four separate subdomains yielding independent flow-rate measurements which are then averaged. We believe that this allows for a good compromise between accuracy and spatial resolution.
The averaged profiles were subjected to further corrections at the near-wall regions where they erroneously exhibit higher velocities (see Sec. III A for a discussion of the origins of these errors). The profile gradient displays a maximum at the onset of the biased data region, which was identified by tracing the gradient locally. The profile data between this position and the wall were then removed and substituted by a parabolic fit [blue dotted line in Fig. 11(b) and 11(c)]. It should be noted that such regions typically do not extend beyond 250 μm away from the wall at worst, and consequently affect a very small portion of the profiles derived at Locations P1 and P2. In terms of the other two locations, and in particular P4 that lies very near the wave trough, the relative fraction of the fitted regions is larger, however, the smoothness and cohesion of the resulting composite profiles suggest that the correction procedure can still be implemented effectively.
Sample local and instantaneous flow-rate time traces generated by the employment of the above approach are presented alongside complementary analytical results in Fig. 12. The flow conditions correspond to Ka = 85, f w = 7 Hz, and Re = 11, 14, 20, and 25. As the Re increases, this time by increasing the flow rate, the deviations between experiments and analytical results grow, and peak in close vicinity to the wave crests, in agreement with the previous assessment. It is therefore evident from the results presented thus far that the examined analytical approach fails to reproduce the experimentally obtained data, as the velocities underneath the waves fall significantly short of the values expected from theoretical analyses relying on the (quasisteady) Nusselt solution.

D. Mass-transfer characterization
The impact of unsteadiness on the mass-carrying capacity of harmonically excited falling-film flows can be examined more rigorously by a "Reynolds decomposition" of the time-varying flow rate per unit span s into mean Q s and fluctuating Q s components, as follows: The time-varying flow rate Q s is equal to the product of the time-varying bulk velocity U b and the time-varying film height h, which can also be decomposed in the same fashion. Thus, Eq. (16) can be expanded into a function that includes the means and fluctuations of the bulk velocity U b and U b and film thickness h and h : Upon averaging Eq. (17), the following expression is obtained for the mean flow rate: The first term in Eq. (18), hereby referred to as the "steady term," is the product of the timeaveraged bulk velocity U b and the time-averaged film height h, and corresponds to the flow rate of an equivalent steady film flow (without waves) that has a bulk velocity and thickness equal to the averaged bulk velocity and thickness of the actual flow under investigation. The second term, designated the "unsteady term," is the covariance of the two fluctuating terms U b and h , and represents the (generally nonlinear) coupling between the local and instantaneous film thickness and the local and instantaneous velocity. For a steady and flat film, the unsteady term vanishes (Q s = 0), and the time-averaged bulk velocity U b corresponds to the Nusselt velocity U N . The introduction of waviness entails a nonzero unsteady term and a consequent deviation between the Nusselt velocity U N and the time-averaged bulk velocity U b obtained in the experiment (Fig. 12). The Nusselt velocity U N overestimates U b due to significant velocity deviations underneath the waves, which grow with increasing unsteadiness (stronger flow-rate fluctuations).
As a consequence of these observations, an alternative approach to the direct comparative study of experiments and theory conducted so far has been pursued. Returning to the flow-rate decomposition of Eq. (18), the PLIF and PTV results allow us to calculate time-averaged flow-rates Q s , as well as time-averaged bulk velocity U b and film-height h values; unsteady terms can then be obtained directly from Eq. (18). For the same experimentally derived time-varying flow-rate trace, film heights and bulk velocities can be calculated at each measurement point using the Nusselt expressions as before, which if processed in the same manner as the experimental results, yield equivalent steady and unsteady terms: where the subscript "QN" designates Nusselt results obtained from time-varying flow rates, as opposed to the subscript "N" which denotes quantities calculated using the Nusselt relations and the mean flow rate Q. The bulk velocities calculated using the relevant Nusselt relation and the instantaneous flow rates are based on the assumption of a parabolic profile everywhere in the flow. Thus, comparing the analytical with the experimental results allows us to assess the suitability of the particular assumption with respect to the presently described investigation. This approach was implemented over a total of 58 flow conditions, 18 originating from the Ka = 85 data set, 22 from the Ka = 350 data set, and 18 from the Ka = 1800 data set.
Based on the results depicted in Figs. 13-15, the following remarks can be made: (1) The analytical predictions consistently match the experimentally observed steady terms, with a mean absolute deviation below 1%, as evidenced by the data in Fig. 13. Consequently, the assumption of parabolic profiles everywhere in the flow has very limited effect on the calculation of the product U b h.
(2) The steady terms scale linearly with Re for each of the three Ka data sets; a relatively unsurprising outcome given that both U b and h are included in the film Re definition.
(3) The gradients of the three best fits to the steady term data approximate closely the respective kinematic viscosities of the three aqueous-glycerol solutions. When the fits are forced through zero, the gradients obtained underestimate the kinematic viscosities of the Ka = 85, 350, and 1800 data sets by 8.1%, 5.4%, and 8.9%, respectively. Without forcing the fits through zero, the deviations amount to −3.4%, −4.2%, and −6.8% for Ka = 85, 350, and 1800, respectively.
(4) This result suggests that knowledge of the mean flow rate per unit span Q s and the kinematic viscosity also allows one to obtain first the steady term (Fig. 13) and then the unsteady term [from Eq. (19)].
(5) The mean deviation between the experimentally obtained and analytically derived unsteady terms over all examined flow conditions corresponds to approximately 8%. (6) The variation in the unsteady terms for the Ka = 1800 flow case is limited compared to the other two liquids; this is inherently linked to the similar values of σ h observed throughout the examined Re range in Fig. 7.
(7) For all three liquids, the unsteady terms scale linearly with the variance of the film height σ 2 h , which is effectively a measure of the (absolute) film waviness (see Fig. 14). Thus, the magnitude of the unsteady term can be linked directly with the waviness of the gas-liquid interface, without knowledge of the flow field underneath the interface. (8) The gradients of these fits, which have units of s −1 , correspond to 128 s −1 for Ka = 85, 204 s −1 for Ka = 350, and 312 s −1 for Ka = 1800. The relevant time scales (inverse of these gradients) calculated with 95% confidence bounds amount to 7.7 ± 0.1 ms for Ka = 85, 4.9 ± 0.2 ms for Ka = 350, and 3.2 ± 0.1 ms for Ka = 1800. Considering the balance between gravity, which drives the flow, and viscosity, which resists it, characteristic length l v and time t v scales ("viscous-gravity scales") can be derived by simple dimensional analysis [13]. The latter can be calculated according to t v = [v/(g sinβ)] 1/3 , and amount to 11.2, 8.0, and 5.3 ms for the three Ka number liquids, respectively. Interestingly, the time scales obtained from the experimental data correspond to 0.6-0.7 of the corresponding t v values. (9) The magnitude of the unsteady terms in the flows excited at f w = 7 Hz are generally higher than those calculated for films excited at f w = 10 Hz, for Ka = 85 and 350, whereas they become independent of the imposed excitation frequency f w for the liquid with Ka = 1800. In both cases, this finding is in agreement with the trend observed between σ h and h in Fig. 7.  In the first case, the relative unsteady terms amount to between 5% and 12%, while the latter, to between 2% and 8%.
(11) Increasing Re for the Ka = 350 flows results in a consistent reduction in the relative magnitudes of the unsteady terms for films excited at 7 Hz, and a somewhat less smooth reduction for films excited at 10 Hz. Both data sets, however, display considerably higher relative unsteady terms compared to their lower Ka counterparts, ranging from 10% to 20% for f w = 7 Hz and from 7% to 12% for f w = 10 Hz. (12) Finally, the unsteady terms do not appear to show a dependency with f w for the Ka = 1800 flows, and decrease in a regular fashion with increasing Re, inside the range 3%-12%.

E. Correlation between the local and instantaneous film height and flow rate
The primary objective here is to examine the link between the local and instantaneous film height and the local and instantaneous flow rate in harmonically excited, fully developed thin-film flows. A prerequisite for meeting this objective is, similarly to the preceding analysis, the availability of simultaneously acquired and spatiotemporally resolved film-height and velocity data from within the wavy films. Using such data, we have demonstrated earlier but also in our previous contributions (Refs. [64,65]) that axial-velocity profiles of nonparabolic shape are encountered ahead of the solitary wave crests and that, depending on the wave topology, the mass transferred along the main wave hump is severely overpredicted or underpredicted by the Nusselt solution based on the local and instantaneous film height. We also showed that the covariance of the bulk velocity fluctuation and the film-height fluctuation vary linearly with the film-height variance, and that both the covariance term and the product U b h can be evaluated with knowledge of the mean flow rate and kinematic viscosity alone.
In the upcoming analysis, DNSs results are used to complement the experimental data, as outlined in Table II. The flow cases considered comprise three aqueous-glycerol solutions (Ka = 85, 350, and 1800), four Re per liquid composition spanning the ranges Re = 11-25, 37-76, and 94-276, respectively, and two forcing frequencies (f w = 7 and 10 Hz). DNS data are provided for three Re per liquid solution for flows with Ka = 85 and 350.
Mean flow-rate data generated by combined PLIF-PTV (normalized by the inlet flow rate Q N measured using an ultrasonic flowmeter) are provided in Fig. 16 for quality assessment purposes. The dotted lines indicate the ±3% error quoted by the flowmeter manufacturer. The results indicate that the majority of the mean flow rates obtained using the combined optical technique lie within the error bounds of the ultrasonic flowmeter measurement; the maximum deviation corresponds to 6%, and the mean absolute and RMS deviations to 2.5% and 3.1%. In an effort to quantify the mass fraction transported by the waves in unforced, laminar falling films, Moran et al. [55] generated instantaneous flow-rate and film-height time series, and plotted the results against each other after normalizing by the mean values. In a similar study, Karimi and Kawaji [54] produced the same graph without normalization for turbulent films falling down a vertical tube. In both cases, the authors noted that thicker film regions (waves) transport more mass than thinner ones (substrate or base film), but less than predicted by the Nusselt flat-film solution in the same regions (thick-film regions). (Recall that a Nusselt flow observes a cubic relationship between the flow rate and film height, Q N ≈ h 3 N .) In this study, we have reproduced such plots for some of the data shown in Fig. 12 (Ka = 85, Re = 11 and 25, and f w = 7 Hz), and have also included the respective Nusselt predictions of the instantaneous flow rate based on the measured film heights in Fig. 17. The flow-rate data are normalized by the mean flow rate measured at the downstream location where the flow was imaged, which by continuity is equal to the mean flow rate at the inlet Q N , while the film-height data are normalized by the local mean film height h. The color scheme associated with the data points corresponds to the probability density (i.e., the relative likelihood of observing a particular pair of normalized flow-rateQ and film-heightĥ values based on the bivariate distribution of the local and instantaneous flow rate and film height).
The deviations highlighted earlier between experimental and analytical results, also in agreement with the results of Moran et al. [55], can once again be observed in this figure. It is interesting to note, however, that based on the experimental data on display, the local and instantaneous flow rate varies linearly with the corresponding local and instantaneous film height, which is surprising. As a consequence of this observation, it can be ascertained that irrespective of the flow field underneath the wavy interface, and by extension the shape of the local axial-velocity profile, there exists a unique flow rate for every film height. In other words, two regions of the flow with the same film height, for example, one behind the wave crest and one in front of it where nonparabolic velocity profiles ensue, will carry the same mass in the streamwise direction. It should be noted that the behavior was not clearly observed or explicitly noted in the aforementioned literature sources; instead, the data provided in those sources displayed a significant scatter. The origin of this scatter will be discussed later on.
Based on the results generated over the entire range of film-flow conditions listed in Table II, we can report that these linear trends were identified consistently in all experiments and in the accompanying numerical simulations. Moreover, the slopes of first-order polynomial (i.e., linear: Q = Aĥ + B) fits to the data from all the examined film flows, such as those shown in Fig. 17, were found to vary across flows, as did the y-axis intercepts.
The slopes (A) and intercepts (B) can be obtained a priori. First, let us consider the continuity equation between two locations along a fully developed wavy film flow, x and x + δx, in a frame of reference moving with the wave speed c: where h (x) is the film height at x, and h (x+δx) the film height at x + δx. Integrating this equation over the film height, and noting that h (x+δx) = h x + δh and Q (x+δx) = Q x + δQ results in Integrating Eq. (21) in the streamwise direction between a location downstream of the inlet where Q = Q = Q N and h = h (we showed earlier that h = h N ), i.e., where the local and instantaneous flow rate and film thickness are equal to the mean flow rate and mean film thickness, and any other location farther downstream (x) gives which is an equation that links the local and instantaneous flow rate per unit span Q(x) (the subscript "s" has been dropped) to the local and instantaneous film height h(x), via the wave speed c and the mean flow rate Q, which also corresponds to the inlet flow rate Q N . The value of Q that will be used in the following expressions corresponds to the flowmeter measurement rather than the mean flow rate calculated using the PLIF-PTV data.
which can be used to describe the data plotted in Fig. 17, with A = c/Q N and B = 1 − ch/Q N . Alternatively, by decomposing h(x) into mean and fluctuating components h(x) = h + h (x), Eq. (22) gives which is an expression linking the normalized flow rate Q(x)/Q N to the film-height fluctuation h (x). Sample plots of the local and instantaneous flow rate per unit span Q(x) against the local and instantaneous film height h(x) from experiments and DNSs are provided in Fig. 18, for Ka = 350, Re = 45, and f w = 10 Hz, and Re = 65 and f w = 7 Hz. Unlike in Fig 17, where experimental results of the local and instantaneous flow rate are plotted against the local and instantaneous film height and compared to analytical results based on the Nusselt solution and the measured film heights, the agreement between the experiments and DNSs is excellent; in the first case, the Nusselt predictions of the local and instantaneous flow rates are linked to the corresponding film heights via a third-order polynomial function, whereas both experiments and DNSs clearly indicate a linear relationship. Furthermore, according to Eq. (22), the plotted data should display a first-order polynomial relationship with a slope described by A = c and an intercept B = Q N − ch. At this point, one can proceed to assessing the validity of these relations, but first, the wave speeds obtained in the experiments and the DNSs are compared in Fig. 19(a) (c in m/s), and in Fig. 19(b) c in units of U N , the Nusselt velocity calculated from Eq. (6).
In the experiments, wave speeds were recovered by first wave and phase-locked averaging a frame with the wave front located near its center (as in Fig. 6). The same procedure was then applied to the previous and next frames, and the displacement of the wave front was calculated by cross correlating the region bounded by the wave crest and trough in each signal pair. The average wave-front displacement was finally divided by the inter-frame separation time dt, yielding a wave speed c. The mean absolute deviation in c calculated from the experimental data and the data from DNS is 3.5% over the 12 flows examined in the present work, while the RMS amounts to 4.4%. As we noted in Ref. [65] for Ka = 85 and 350, the wave speed c decreases with increasing inertia (increasing Re); the same trend is observed here for the data from the films with Ka = 1800 as well.
Returning to Eq. (22), linear fits were generated for both experimental and numerical data in plots of h(x) against Q(x), and compared to the predictions A = c and B = Q N − ch. The wave speeds used in the analytically calculated coefficients were obtained directly from the experimental or the DNSs data, while Q N , which is common to both, was measured in the experiments and provided as an input to the DNSs scheme. The results for A and B normalized by their respective analytical values are presented in Fig. 20. Key observations are summarized below: (i) The mean absolute deviation between analytical and best-fit data for experimentally derived values of A is 4%, while the mean deviation is −3.8%, suggesting that the best-fit slopes typically underestimate the experimentally derived wave speeds. For the DNS data, these deviations amount to 1.8% and 0.1%.
(ii) The mean absolute deviation between analytical and best-fit data for experimentally derived values of B is 8%, and the mean deviation is −6%. For the DNSs, these deviations amount to 3.4% and −0.2%.
A similar analysis was also performed based on Eq. (24). The slope and intercept results are plotted in Fig 21. Based on this description, B = Q/Q N 1, that is the mean flow rate obtained from the experiments and DNS normalized by Q N . As was noted earlier, Q N is an input to the DNSs and therefore any deviations from B = 1 are insignificant (below 1%), and probably originate from numerical errors. The experimentally derived values of B were, effectively, examined in Fig. 16 and displayed a mean absolute deviation of 2.5%. With respect to the slope results, the best-fit values recovered in the experimental and DNS studies were normalized by c/Q N , and displayed the same deviations as the A = c data examined earlier (around 4% and 2%, respectively). The agreement between experiments, simulations, and the proposed correlation [Eq. (24)] is excellent. Furthermore, it is also very interesting to note that the description based on Eq. (24) does not require explicit knowledge of the film height, but rather only the fluctuation about some mean value. Consequently, accurate flow-rate data can be determined for any location or instant along the film, even if the real mean film height is not known, that is, even if the location of the liquid-solid interface is not identified. In fact, from an experimental point of view, the application of this simple linear relationship presupposes knowledge of the mean flow rate and the wave speed only. The former can be easily obtained using a flowmeter installed before the settling chamber (like in this study), while the latter can be measured by direct imaging (without the need for an expensive and complex PLIF-PIV setup).
Finally, it is essential to discuss briefly the results of Moran et al. [55] and of Karimi and Kawaji [54], with regard to their efforts to characterize the mass transported by the waves in their laminar and turbulent film flows. We initiated this analysis by noting that the aforementioned researchers generated normalized (by the mean values) flow-rate and film-height time series, and plotted the results against each other, only to observe a strong scatter rather than the linear behavior identified in our study. It is possible that the former stems from the fact that in both experiments, the films were not harmonically excited and were, consequently, still developing at the measurement location. This implies that in our case the wave shape and speed display only insignificant variations over time (in essence, we repeatedly observe the same solitary wave), whereas in their case the same flow condition will span a wide range of wave heights and wave speeds. Even if individual waves or even wave packets moving with the same speed obey a linear trend, the latter would only be clearly discernible provided that an analysis similar to ours is carried out for each wave structure separately. 014002-31

VI. CONCLUSIONS
A measurement technique featuring combined, simultaneous planar laser-induced fluorescence imaging (PLIF), particle image and tracking velocimetry (PIV-PTV) was employed alongside direct numerical simulations (DNSs), towards the detailed characterization of the hydrodynamic characteristics of harmonically excited, falling liquid-film flows. In our experimental campaign, we adopted four different aqueous-glycerol solutions, corresponding to four Ka (14,85,350, and 1800), in flows with Re covering the range 2.3-320 and inlet forcing frequencies f w = 7 and 10 Hz. PLIF was employed to generate spatiotemporally resolved film-height measurements, and PIV to generate 2D velocity maps of the flow field underneath the wavy interface. PTV was implemented as a final processing step to enhance the spatial resolution of the velocity measurements and moderate the effect of bias errors in regions of the flow characterized by high strain. Wave and phase-locked flow field averaging based on the PTV and PLIF results were introduced as a processing methodology aimed at the generation of highly localized velocity profile, bulk velocity, and flow-rate data along the wave topology. Based on these measurements, the effect of waviness and local film topology on the flow field underneath the waves was studied in detail.
A series of validation and quality assessment experiments were conducted to estimate the errors associated with our optical technique. First, film-height data from flat (unforced) films (Ka = 14) were compared to micrometer stage measurements, and the 1D, steady, fully developed solution to the NS equations, referred to as the Nusselt solution. From these, a mean relative error of 2.5% can be quoted for the local and instantaneous film height, while a measurement error of 2% is believed to be a representative estimate of the error associated with average film-height measurements (time averaged and wave-locked average). Relative deviations were calculated between PTV-derived interfacial and bulk velocities and Nusselt predictions, with mean values amounting to 3.2% for both velocities. Finally, flow-rate comparisons were conducted between PLIF and PTV-derived data and flowmeter readings. The mean relative deviation over a total of six flat and nine wavy flows included in this test was 1.6%, with no systematic trends identified. This value is smaller than the flowmeter uncertainty quoted by the manufacturer (3%), so the latter can be used as the estimate of our mean flow-rate measurement error.
Our experiments revealed that the mean film height h was lower than the Nusselt height h N for the majority of examined films, while the film-height standard deviation σ h , effectively a measure of film waviness, was found to be higher for films excited at f w = 7 Hz compared to 10 Hz for Ka = 85 and 350. The Ka = 1800 results suggested that σ h is decoupled from both Re and f w . The coefficient of variation σ h /h peaked at 47% in the range Re = 70-80.
Nonparabolic velocity profiles were encountered in regions bounded by the wave crests and troughs, as well as underneath the capillary waves. Comparisons to velocity profiles generated based on the steady-flow solution indicated that the latter overestimate the real mass-carrying capacity of the waves, often by more than 100%.
Time-varying flow-rate traces were generated to study systematically the impact of unsteadiness on the mass-carrying capacity of harmonically excited film flows. The local and instantaneous flow rate was decomposed into steady and unsteady terms, the former represented by the product of the mean film height and bulk velocity U b h, and the latter by the covariance of the film height and bulk-velocity fluctuations U b h . The most significant results of this analysis are briefly summarized below: (i) The steady terms vary linearly with the flow Re, and the gradients of the three fits to the steady-term data (one per Ka) are approximately equal to the kinematic viscosities of the respective aqueous-glycerol solutions.
(ii) For all three liquids, the unsteady terms scale linearly with the variance of the film height σ 2 h , which is effectively a measure of absolute film waviness. The gradients of these fits, which have units of s −1 , amount to 128 s −1 (7.8 ms) for Ka = 85, 204 s −1 (4.9 ms) for Ka = 350, and 312 s −1 (3.2 ms) for Ka = 1800. The experimentally derived time scales are of the same order of magnitude as the characteristic time scales t v calculated based on the balance between gravity and viscosity, varying in the range (0.6-0.7) × t v .
Based on local, time-varying flow-rate and film-height traces, plots of the instantaneous film height against the instantaneous flow rate were generated using data from both experiments and numerical simulations. The two quantities were found to vary linearly against each other, with the slopes and y intercepts being functions of the flow and inlet conditions (Re, Ka, and f w ), a result which is particularly interesting. This result effectively suggests that the investigated solitary waves are fully developed, which together with mass conservation arguments, suggests that regions with the same film height are associated with the same mass-transfer capacity in the streamwise direction of the flow, irrespective of the flow field underneath the wavy interface.
From an analysis of the coefficients of linear fits (of the form Q = Ah + B) to experimental and numerical data from all 24 examined film flows in this work, it was found that the data obey the relationship Q = (Q N − ch) + ch with small deviations. An alternative formulation of the same analytical relationship, following a decomposition of h into steady (h) and fluctuating (h ) components, was also put forward: Q/Q N = 1 + c/Q N h . The advantage of this over the previous expression is that it can be employed in order to yield local and instantaneous flow-rate data from film-height measurements, based on knowledge only of the film-height fluctuation h , rather than the real absolute film height; the former is much easier to obtain experimentally than the latter.