We study the stability of the interface between
(a) two adjacent viscous layers flowing
due to gravity through an inclined or vertical channel that is confined
between two
parallel plane walls, and (b) two superimposed liquid films flowing
down an inclined
or vertical plane wall, in the limit of Stokes flow. In the case of channel
flow, linear
stability analysis predicts that, when the fluids are stably stratified,
the flow is neutrally
stable when the surface tension vanishes and the channel is vertical, and
stable
otherwise. This behaviour contrasts with that of the gravity-driven flow
of two
superimposed films flowing down an inclined plane, where an instability
has been
identified when the viscosity of the fluid next to the plane is less than
that of the top
fluid, even in the absence of fluid inertia. We investigate the nonlinear
stages of the
motion subject to finite-amplitude two-dimensional perturbations by numerical
simulations based on boundary-integral methods. In both cases of channel
and film
flow, the mathematical formulation results in integral equations for the
unknown
interface and free-surface velocity. The properties of the integral equation
for multi-film
flow are investigated with reference to the feasibility of computing a
solution by
the method of successive substitutions, and a deflation strategy that allows
an iterative
procedure is developed. In the case of channel flow, the numerical simulations
show
that disturbances of sufficiently large amplitude may cause permanent deformation
in
which the interface folds or develops elongated fingers. The ratio of the
viscosities and
densities of the two fluids plays an important role in determining the
morphology of
the emerging interfacial patterns. Comparing the numerical results with
the predictions
of a model based on the lubrication approximation shows that the simplified
approach
can only describe a limited range of motions. In the case of film flow
down an inclined
plane, we develop a method for extracting the properties of the normal
modes,
including the ratio of the amplitudes of the free-surface and interfacial
waves and their
relative phase lag, from the results of a numerical simulation for small
deformations.
The numerical procedure employs an adaptation of Prony's method for
fitting a signal
described by a time series to a sum of complex exponentials; in the present
case, the
signal is identified with the cosine or sine Fourier coefficients of the
interface and
free-surface waves. Numerical simulations of the nonlinear motion confirm
that the
deformability of the free surface is necessary for the growth of small-amplitude
perturbations, and show that the morphology of the interfacial patterns
developing
subject to finite-amplitude perturbations is qualitatively similar to that
for channel
flow.