Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 5.1.2 - A Numerical Study of Faraday Waves
Full Citation
Permanent Link:
 Material Information
Title: 5.1.2 - A Numerical Study of Faraday Waves Instabilities
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Périnét, N.
Juric, D.
Tuckerman, L.S.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: numerical simulation
pattern formation
Faraday waves
Abstract: We simulate the full dynamics of hexagonal Faraday waves in three dimensions and provide comparisons of the calculated spatio-temporal spectra to experimental results. The Navier–Stokes equations are solved using a finite-difference projection method coupled with a front-tracking method for the interface between the two fluids. We study the long-time behavior of the hexagonal patterns which simulations suggest is a fixed point of a homoclinic orbit. Furthermore, we present recent results on the drift instability in 2D simulations for air/water for which we are able to begin to characterize secondary and possibly tertiary bifurcations.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00119
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 512-Juric-ICMF2010.pdf

Full Text

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

A Numerical Study of Faraday Waves

N. P6rinet*, D. Jurictand L. S. Tuckerman*

PMMH ESPCI, UMR 7636, 10 rue Vauquelin 75231 PARIS CEDEX 5 FRANCE
Keywords: Numerical simulation, instabilities, pattern formation, Faraday waves


We simulate the full dynamics of hexagonal Faraday waves in three dimensions and provide comparisons of
the calculated spatio-temporal spectra to experimental results. The Navier-Stokes equations are solved using a
finite-difference projection method coupled with a front-tracking method for the interface between the two fluids. We
study the long-time behavior of the hexagonal patterns which simulations suggest is a fixed point of a homoclinic
orbit. Furthermore, we present recent results on the drift instability in 2D simulations for air/water for which we are
able to begin to characterize secondary and possibly tertiary bifurcations.

1 Introduction

The Faraday experiment consists of shaking vertically a
container holding two immiscible fluids thereby induc-
ing oscillations of the fluids and the interface between
them. Beyond a certain threshold, the interface can form
many kinds of standing-wave patterns, including crys-
talline patterns and others which are more complex. This
phenomenon was first studied by Faraday (1831) who
noticed that the vibration frequency of the interface was
half that of the forcing. Benjamin & Ursell (1954) car-
ried out the first theoretical linear analysis of the Fara-
day waves, restricted to inviscid fluids. They decom-
posed the fluid motion into normal modes of the con-
tainer and showed that the evolution equation of each
mode reduced to a Mathieu equation whose stability di-
agram is well known.
In the 1990s, new behaviours of the interface were
discovered, such as quasi-crystalline eight-fold patterns
seen by Christiansen, Alstrom & Levinsen (1992). Other
patterns emerge using two-frequency forcing such as
twelve-fold quasi-patterns (Edwards & Fauve (1994)),
triangular patterns (Miller (1993)) and superlattice pat-
terns (Kudrolli, Pier & Gollub (1998)). Spatiotemporal
chaos was studied by Kudrolli & Gollub (1996), who
also surveyed the occurrence of lattice patterns stripes,
squares or hexagons as a function of viscosity and fre-
quency. In addition to patterns or quasi-patterns, very
localized circular waves called oscillons may occur, as
seen by Lioubashevski et al. (1996). The Faraday in-

stability is the first macroscopic system in which such
structures have been observed. These discoveries endow
the Faraday instability with a very great fundamental
interest for understanding the natural formation of pat-
A number of theoretical or semi-numerical analyses
were inspired by these experiments. Kumar & Tucker-
man (1994) extended the linear stability analysis of Ben-
jamin and Ursell to viscous fluids. This analysis was ex-
perimentally confirmed by Bechhoefer et al. (1995) and
used by Kumar (1996) to predict cases in which the re-
sponse would be harmonic rather than subharmonic.
Linear analysis provides no information about the
shape of the patterns which appear; other means are
necessary to understand the occurence of a given pat-
tern or the amplitude of stabilization. Weakly nonlin-
ear approximations have been derived from the Navier-
Stokes equations by Vifials and co-workers, e.g. Zhang
& Vifials (1997) and Chen & Vifials (1999), and by Skel-
don & Guibodoni (2007), focusing on the competition
between different patterns. Vega, Knobloch & Martel
(2001) derived equations governing the interaction be-
tween Faraday waves and the mean flow. There has
been a great deal of analysis of lattices, superlattices
and quasipatterns using equivariant dynamical systems
theory, as well as model equations designed to produce
specific patterns, e.g. Porter, Topaz & Silber (2004). The
approximation of quasipatterns in spatially periodic do-
mains has also been addressed in Rucklidge & Silber

_________ Cubic lixed
3D mesh

Interface: (x, ,t)

Triangular deformable
2D mesh

Figure 1: Spatial discretization of the domain.

Investigation of the full nonlinear viscous problem,
however, requires numerical simulations, of which there
have been very few up to now, specifically those of Chen
& Wu (2000) and Chen (2002), Murakami & Chikano
(2001), Valha, Lewis & Kubie (2002), Ubal, Giavedoni
& Saita (2003) and O'Connor (2008). With the excep-
tion of O'Connor (2008), all previous simulations have
been two-dimensional. Here we use the method devel-
oped and validated in Perinet, Juric & Tuckerman (2009)
to report on the results of further studies of Faraday
waves, specifically the long time behavior of the hexag-
onal pattern. In addition we investigate, using 2D simu-
lations, the drift instability for the air/water system and
its bifurcation diagram.
In the next section of this article, we present the hy-
drodynamic equations that govern the Faraday instabil-
ity and briefly describe the computational method. We
then compare our results with the experiments of Kityk
et al. (2005, 2009) and present the three-dimensional ve-
locity field for the hexagonal pattern. We discuss the
long-term stability of the hexagonal pattern and investi-
gation of the homoclinic orbit. Finally we present results
from a 2D study of the drift instability in the air/water
system and discuss the emerging bifurcation diagram.

2 Equations of motion

The mathematical model of the Faraday experiment con-
sists of two incompressible and immiscible viscous flu-
ids in a three-dimensional domain x (x, y, z) E
h'-' x [0, h], bounded at z = 0 and z = by flat walls at
which adherence conditions are applied. The two fluids,
each uniform and of densities pi, P2 and viscosities H1
and P2, initially form two superposed horizontal layers
with an interface between them. This two-dimensional
interface is defined by x' = (x, y, ((x, y, t)). Within
the parameter range we wish to simulate, the height C

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


i /I
I .. I I

I I / I


Figure 2: Lattice formed by the spatial modes compris-
ing a hexagonal pattern. The principal modes,
with wavenumbers kc, 2k, and v 3k, in-
volved in later quantitative investigations are
indicated by hollow black circles. The arrows
illustrate resonance mechanisms leading to
harmonic contributions to higher wavenum-

remains a single-valued function of (X, y, t).
The gravitational acceleration, g, is augmented by a
temporally periodic inertial acceleration

G = (a cos(wt) g)e,

where a is the amplitude of the forcing and c is its fre-
The Navier-Stokes equations for incompressible,
Newtonian fluids are:


Vp + pG+V .(Vu + Vu) + s (2a)


Here p is the pressure and s is the capillary force (per
unit volume) and is defined below. Equations (2a) and
(2b) are valid for the entire domain, including the inter-
face, in spite of the fact that the density and viscosity
change discontinuously and the surface tension, s, acts
only at the interface. In this single-fluid formulation, the
density and viscosity fields are defined in terms of the
densities and viscosities of the two fluids

P = P + (P2 P)H
P = /1 + (/12 1 )H

with the aid of a Heaviside function,

H (x- x') 0 ifz < ((r, y, t)
1 if z > ((x, y, t)

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

S(n2. m)

0.39 2.17 3.95 5.73 7.51

0 (nm)

0.21 0.87 1.53 2.19 2.85

Figure 3: Example of hexagonal patterns: height of the
interface as a function of the horizontal co-
ordinates for a 38.0 ms 2. Resolution:
58 x 100 x 180 in x, y, z directions respec-
tively. Snapshot taken when height of the in-
terface peaks is maximal. Note that the ver-
tical and horizontal scales are different. Each
horizontal direction in the figure is twice that
of the calculation domain. Color indicates in-
terface height.

where we recall that x = (x, y, z) is a point anywhere in
the three-dimensional volume and x' (x, y, ((x, y, t))
is the vertical projection of x onto the interface. The
capillary force is

s = or n 6 (x x') dS (5)
s (t)

where a is the surface tension coefficient, assumed to be
constant, n is the unit normal to the interface (directed
into the upper fluid) and K its curvature. 6 (x x') is a
three-dimensional Dirac distribution that is nonzero only
where x x'. S'(t) is the surface defined by the instan-
taneous position of the interface.
The interface is represented implicitly by H and ad-
vected by material motion of the fluid

The computational method is based on a MAC pro-
jection method Harlow & Welch (1965) coupled to a
front-tracking/immersed-boundary method for the fluid
interface. The calculation domain is a rectangular paral-
lelepiped, horizontally periodic in x and y and bounded
in z by flat walls for which we impose no-slip bound-
ary conditions. The entire domain is discretized by a

Figure 4: Snapshot of hexagonal pattern taken t = 0.3 x
2T after instant of maximal interface height.

uniform fixed three-dimensional finite-difference mesh.
Within the domain, the two distinct immiscible flu-
ids are separated by a two-dimensional interface which
is discretized by a second mesh as sketched in fig-
ure 1. This moving and deformable mesh is com-
posed of triangular elements whose motion is treated
by a front-tracking/immersed-boundary method (Peskin
(1977); Tryggvason et al. (2001)). Because we have as-
sumed that ((x, y, t) is single-valued, the nodes of the
mesh can be fixed in x and y and only their vertical dis-
placements need to be calculated
After setting appropriate initial and boundary con-
ditions, the computational solution algorithm for each
timestep is composed of three main phases. First the in-
terface is advected and the density and viscosity fields
updated according to the new interface position. The
capillary force, s, is then calculated. Finally the velocity
and pressure are found by means of a standard projection
method. Full details of the method including calculation
of the capillary force are described in Perinet, Juric &
Tuckerman (2009)

3 Results: nonlinear analysis

In the full nonlinear evolution of the interface for a >
ac, the amplitude of the interface height grows in time
until nonlinear terms in equations (1)-(6) become im-
portant. After that, the mode whose linear growth rate
is maximal gives rise, via nonlinear resonances, to a se-
ries of other discrete modes, selected according to the
magnitudes and orientations of their wavevector, k. This
selection is responsible for the formation of patterns that
will be the object of our validations. We seek to compare

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

S(n1. m)

0.10 1.33 2.55 3.77 5.00


0.10 0.64 1.18 1.72 2.26

Figure 5: Snapshot of hexagonal pattern taken t
0.48 x 2T after instant of maximal interface

our calculations with the experimental results of Kityk
et al. (2005, 2009) where quantitative data concerning
the Fourier spectrum ((k, t) are available.
We run our numerical simulations with the same ex-
perimental parameters as Kityk et al. (2005): w/27
12 Hz (T 0.0833 s), pi 1346 kg m-3, p
7.2 mPas for the lower fluid and P2 949 kgm 3,
P2 20 mPa s for the upper fluid. The surface tension
at the interface is a 35 mNm 1, the total height of
the vessel is 1.0 cm, and the mean height of the interface,
the initial fill height of the heavy fluid, is () 1.6 mm
(with some uncertainty; see below). The Floquet anal-
ysis for these parameters yields a critical wavelength of
Ac 27/kc 13.2 mm and a critical acceleration of
a, = 25.8 ms 2.
Rather than starting from a sinusoidal interface, we
chose to add two-dimensional white noise of small
amplitude to (C) to define the initial interface height,
((x, y, t 0 ), in order to excite every mode allowed
by the box's horizontal dimensions and number of cells.
It is thus possible to check that the correct critical mode
(that whose growth rate is maximal) emerges from the
linear dynamics. In order to reproduce the experimental
results in a computational domain of a minimal size, the
dimensions in x and y of the box must correspond to the
periodicity and symmetries of the expected pattern.

3.1 Hexagonal patterns

Although the emergent pattern is square at lower am-
plitudes of the forcing acceleration, as a increases the
modes can reorganize. The symmetries change and, in
the experiments of Kityk et al. (2005), the initial square

Figure 6: Snapshot of hexagonal pattern taken t
0.68 x 2T after instant of maximal interface

pattern becomes hexagonal. Though kc remains con-
stant, the horizontal dimensions of the minimal com-
putational box necessary to support the periodic pattern
must change too. These dimensions become 47/k, in y
and 47/(/3k,) in x. The wavevector lattice for hexag-
onal patterns is shown in figure 2. The principal modes
of three amplitudes: kc, 2k6, and /3k,. When a pat-
ter is hexagonal, a mode will have the same amplitude
and temporal behaviour as each of its images through
rotations by any integer multiple of 7/3. The interface
height is thus:
((x, t) () + A(k, t) eikex
+ A(2k6,t) Z ei2cek_ x
+ A( 3k6,t)Ze1 i3k'
+ higher order terms (7)

where ej ex cos(wj/3) + ey sin(7j/3) and e' =
ex cos(7w/6 + 7j/3) + ey sin(7w/6 + 7j/3). for j
Our simulations are carried out at acceleration a
38.0 ms 2 and mean height (C) 1.6 mm. In figures
3 to 5, we show visualizations of the patterns at four in-
stances in time. The 7/3 rotational symmetry confirms
that the rectangular numerical grid does not forbid the
formation of hexagonal patterns, which are not aligned
with this grid. The patterns reproduce several prominent
features from the visual observations of hexagons in the

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010




.05'O 6.I 6.15 6.2 6.25 6.3

10 x

10 x

Figure 7: Temporal evolution of the amplitudes of the
spatial modes with wavenumbers ke, 2k, and
/3k,. Solid curves represent experimental
results (Kityk & Wagner, private communica-
tion) at a 38.5 ms 2. Dashed curves rep-
resent the simulation for a 38.0 m s 2 at
resolution (in x, y, z directions) of 58 x 100 x
180. Arrows indicate times corresponding to
figures 3 to 5.

experiments. For example, one can observe in figures 3
and 4 the up and down hexagons shown in the experi-
mental snapshots (figure 10 of Kityk et al. (2005)). The
pattern in figure 5, when the surface elevation is mini-
mal, is dominated by wavenumbers higher than k, as is
also the case in figure 10 of Kityk et al. (2005). This is
reflected by the disappearance of A (ke) and the resulting
dominance of A(2k,) and A(/3kc) at the correspond-
ing instant in the spectral timeseries of figure 7.
The spectra from experiments and simulations are
represented in figures 7 and 8. Given the experimental
uncertainties concerning the hexagons (Kityk & Wagner,
private communication), the agreement is remarkable.
The principal mode is well reproduced while the other
two modes show rough agreement. Every wavevector is
a superposition of harmonic and subharmonic temporal
modes, so that each has temporal period 2T.
In addition to the interface height, our simulations
also produce the entire velocity field, which is the focus
of figures 9-11. These figures show the velocity fields
on horizontal planes at three instants spanning the oscil-
lation period of a hexagonal pattern, as well as the ver-
tical velocity on the interface. Figures 9, 10 and 11 cor-
respond approximately to the visualizations of figures 3,
4 and 5, where the structures are more visible since the
interface has been repeated periodically in the horizontal
directions for clarity.
Figure 9 is taken at t = 0.07 x (2T), just after the in-
terface reaches its maximum height (at t = 0), when the
peaks are beginning to descend. Consequently, the fluid
converges horizontally towards the interface peaks, then
descends dramatically below them. The fluid then di-

Figure 8: Temporal Fourier transform of the amplitudes
in figure 7. Circles represent experimental re-
sults (Kityk & Wagner, private communica-
tion) for a z 38.5 ms 2. Crosses represent
numerical results for a 38.0 m s 2 at reso-
lution 58 x 100 x 180.

verges horizontally outwards near the bottom and moves
upwards in the large regions between the peaks. The
motion shown in figure 10, at t = 0.41 x (2T), is quite
different from that in figure 9. The peaks of figure 9 have
collapsed into wide flat craters. The fluid diverges out-
wards horizontally above the peaks, then descends into
the craters and diverges outwards horizontally just be-
low them. Figure 11, at t 0.73 x (2T), shows that the
rims of the wide flat craters seen in figure 10 have in turn
collapsed inwards, forming circular waves which invade
the craters, whose remnants are visible as dimples. The
velocity field of figure 11 shows fluid converging hori-
zontally below these dimples. These are erupting at ve-
locities which are the largest in the cycle, and will even-
tually reconstitute the high peaks seen in figure 9.

4 Behaviour of hexagons at long time

The hexagonal patterns which appear in our simulation
of the experiment of Kityk et al. seem to be only tran-
sient. The hexagons vanish after a few seconds and are
replaced by patterns with other symmetries. The exper-
imental observations confirm these outcomes (Kityk &
Wagner, private communication) but the durations of the
hexagonal regime are not equal in both cases.
This symmetry breaking is closely linked to the am-
plitude AC. We indeed notice in Figure 12 that the am-
plitude saturates but only for a short while, after which
AC decreases and the hexagonal symmetry is broken
(see Figures 13 et 14). AC then stays at a lower level for
a few seconds but this state is periodically punctuated
by jumps in amplitude correlated with the partial reap-
pearance of the hexagonal invariance. The pattern ob-
served during such amplitude peaks is shown in Figure
15 where it can be noticed that the rotational invariance


ii/ '

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

w(mm s')


Figure 9: Velocity field at time t = 0.07 x (2T) after
the instant of maximum height. Interface is
colored according to the vertical velocity w.
Arrows show velocity field at z = 0.53 mm
and z = 6.08 mm. (Total height is 10 mm,
average interface height is 1.6 mm.) For clar-
ity, velocity vectors are plotted only at every
fourth gridpoint in each direction. Note that
the vertical and horizontal scales are different.
One computational domain is shown.

of the state is 27/3 instead of 7/3. This implies that
the central symmetry is also broken. We perform a spec-
tral decomposition by which we can track the modes of
wavenumber ke, but oriented in different directions, as
well as those they generate through the nonlinear inter-
actions, shown in Figure 16.
The modes associated with the same wavenumber
have different behaviours. Some of the spatial modes
are subjected to ,igiiiik.ii damping that accompanies
the decay of A(. Their angular distribution is responsi-
ble for the changes in the symmetry of the pattern. These
successive round-trips around the hexagonal regime sug-
gest that this state is the fixed point of a homoclinic orbit.

4.1 Study of the homoclinic orbit

The first step in studying the homoclinic orbit is to cal-
culate the fixed point in phase space, e.g.the amplitude
of each spatiotemporal mode of a perfect hexagonal pat-
tern. To achieve this, hexagonal symmetry of the inter-
face is imposed at regular time intervals. Assuming that
the breaking of hexagonal symmetry is due to an unsta-
ble direction, forcing the hexagonal regime will avoid
the unstable eigenvector and project the solution to the
stable eigenspace at the vicinity of the fixed point which
should lead the dynamics to converge towards this point.

Figure 10: Velocity field at time t = 0.41 x (2T) af-
ter the instant of maximum height. Vectors
shown at z = 0.083 mm and z = 6.25 mm.
Vector and color scales differ from those of
figure 9.

The group D6 of the hexagonal symmetry is generated

7/3 rotational invariance, associated with the oper-
ator R,/3,

reflection symmetry through an axis joining oppo-
site vertices of the hexagon whose direction is fixed
by the ratio between the size of the box and the crit-
ical wavelength and which corresponds to the oper-
ator S. The direction of this axis varies according
to the mode that is considered.

A condition of translational invariance is added to the
elements of D6 in the framework of pattern formation.
The translational invariance is insured by the use of pe-
riodic functions q(x) or numerically by a horizontally
periodic domain. These functions are then spatially ex-
panded into horizontal Fourier modes

q(x)= im .eikmx

27l 27m
where kim = -e, + e is the set of horizontal
xu ly
wavevectors permitted by the discretization of the do-
main with ,, and 1, denoting the dimensions of the pe-
riodic domain along x and y. Each spatial mode q4, is
complex, but imposing q_- = q7* produces a q(x)
which is real.
Let L be an element of D6. The vector field q is in-
variant under the action of L. We suppose here that the

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

- 0.

Figure 11: Velocity field at time t = 0.73 x (2T) af-
ter the instant of maximum height. Arrows
show velocity field at z = 1.58 mm. Vector
and color scales differ from those of figures
9 and 10.

point of horizontal coordinates (0,0) is a center of sym-
metry if
q (L(x)) L (q(x)) Vx (9)
Due to the linearity of both Fourier and L operators, (8)
and (9) are combined to yield the condition of invariance
through L for each Fourier mode,
meikim.L(x) = L(qim)eik/'.x (10)

The solution of (10) with L R,/3 leads to the first

if R ~/ (kim) or Ra2/3 (kim) does not belong to
the Fourier mesh, then q~i 0.

0 5 10 15 20
time (s)

Figure 12: Evolution of A((t) in the nonlinear satu-
rated regime. Resolution: 180 points in z,
50 points per wavelength in x and y. The
hexagonal state is reached when the height
amplitude is maximum.

-10 0
k /k

10 20

Figure 13: Example of patterns occurring after the
hexagons have vanished. Top: contour lines
of ((x, y, t). bottom: spatial spectrum of C.

This condition makes possible the reorganization of in-
dices in a way that is more practical for further calcula-
q(x) q seik x(11)
r s=1
where the index r represents the number of the hexagon
considered and R,/3 (krs) = krmod(s+1,6). Hence, it is
deduced from (9) applied to L R,/3 and (11) that

R7/S (q*s) -= rmod(s+ 1,6)

Thus the additional constraint:

for all r if R,/3 (krs) and R2a/3 (krs) belong to the
Fourier mesh, then R,/3 (q 6) R27/ (q) 5)
-qr4 R-27/3(q3) R- (/3q(2) 41.

The solution of (10) with L Ss (the reflection
through an axis oriented along krs) imposes additional
constraints that depend on the component of the vector
qs to which it is applied:

the vertical component of each q$ is real,

the horizontal components of every q$ are pure

, ,
. *

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

It i

0 0.5 1 1.5 2

-10 0 10 20

Figure 14: Further example of patterns occurring after
the hexagons have vanished. Top: contour
lines of ((x, y, t). bottom: spatial spectrum
of .

*for a given pair (R, S), if kRS provides hexagons,
e.g.satisfies the conditions for rotational invariance,
then for all pairs (r, s), there exists (r', s') such that
,, S=RS (qTS,).

The last condition is non trivial and subtle. This result
states first that the behavior of any 4q, is conditioned by
a reference 4,R fixed for example as the critical mode,
and second that for a fixed r, the set of k,, can form an
hexagon which is not symmetric through SRs. Hence,
the regular rectangular mesh we use implies that there
exists r' such that the hexagon formed by k,,,, is sym-
metric to the former hexagon. The variables that are to
be symmetrized in our case are the interface height (
which acts as a vertical vector and the velocity field u
whose direction depends on its position.

Once the fixed point is found, the hexagonal symme-
try constraints defined above are relaxed. Subsequently,
the dynamics should move away from the fixed point
along the unstable direction. Therefore, the component
of the unstable eigenvector will be easily computed and
we expect to be able to determine the cause of the peri-
odic symmetry breaking from the hexagonal regime.

0 0.5 1


-10 0

1.5 2

10 20

Figure 15: Triangular pattern observed during a jump in
amplitude. Top: contour lines of ((x, y, t).
bottom: spatial spectrum of C.

5 Drift instability

The Faraday problem Faraday (1831) has been exten-
sively studied since it provides a simple and very gen-
eral experiment for the understanding of pattern forma-
tion. It has also served as a model for other theoret-
ically predicted phenomena, in particular, the genera-
tion of a mean flow. Called the drift instability, this
mean flow was first observed in Douady, Fauve & Thual
(1989) for an annular configuration. Theoretical investi-
gations have since then been conducted and it has been
shown to be associated with the breaking of horizontal
reflection symmetry (Fauve, Douady & Thual (1991);
Martin, Martel & Vega (2002)). Knobloch, Martel &
Vega (2002) showed, using weakly nonlinear theory, that
a drift instability occurred through a parity-breaking bi-
furcation and that the appearance of this drift was mainly
caused by the excitation of a visocus mean flow near the
boundaries of the fluids.

Here we study the drift instability by simulation of
the complete Navier-Stokes equations in two dimensions
Perinet, Juric & Tuckerman (2009) in a periodic rectan-
gular domain (or equivalently, a circular domain whose
radius of curvature tends to infinity).


00 o o0
0 0& 0 O o
. 0 0 0
0 0 0 0
$000, o

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010


0 I -

0 1

8 10 12 14 16 1S 20

Figure 16: Evolution of the main horizontal modes of
the interface: the three modes ke, 2key,
(vke, + ?2 e )/2 (solid lines) and their
images through rotations of r/3 (circles)
and 27/3 (crosses).

1 4x 10


0 8



0.125 0.13 0.135 0.14 0.145 0.15 0.155 0.16

Figure 17: Bifurcation diagram for the mean flow insta-
bility. Mean volumetric flux of air (circles)
and water (crosses).

5.1 Physical conditions

Our calculations are based on the physical conditions
provided by Vega, Knobloch & Martel (2001) for a
single layer of fluid. Our numerical approach, how-
ever, is designed to simulate two superposed fluid lay-
ers on a horizontally periodic domain. Thus we take
the lower fluid to be water as in Vega, Knobloch &
Martel (2001) and the upper fluid is taken to be air.
The densities of the fluids are pi 1032 kg m3 and
P2 1.205 kgm 3, respectively, and their viscosities
are/ 1.023 x 103 Pas and 2 1.821 x 105 Pas.
The surface tension a is 7.2 x 10 2 Nm 1, the height
of the water layer is hi = 8 mm and the height of the
air layer is h2 10 mm, such that the total height
of the box h is 18 mm. The frequency of vibration is


i SO

0.01 0.02 0.03 0.04
x (m)

Figure 18: Time-averaged variables for a = 0.127-.,
which is below the secondary bifurcation.
Contours: streamfunction, vectors: velocity
field, line: interface. The profile is symmet-
ric and the mean variables are A/2-periodic.

F = 10.6 Hz for which the critical wavelength Ac is
5.03 cm. Since weakly nonlinear theory is only valid
for supercritical bifurcations, we impose the length of
the box to be 4.65 cm, for which the bifurcation is su-

5.2 Secondary instability (drift instability)

The primary threshold for the onset of the classic Fara-
day instability has been compared to the theoretical
value Kumar & Tuckerman (1994). The numerical
threshold is approximately 0.1243g, whereas the theo-
retical value is 0.126g.
We define the temporal Eulerian mean of a field vari-
able q as

q(x,z,t) dt (13)

where T denotes a harmonic period or the period of vi-
bration of the apparatus. The time interval must be an
integer number of subharmonic periods 2T, the oscilla-
tion period of the fluid. To quantify the drift instability,
we have chosen to study the spatio-temporal mean the
volumetric flux of either fluid (water or air), expressed

V1 (u),- Tj T A = u(x,z,t) ex dzdxdt
8T A it Jo Jo

2 (u), +,t J t1 /( ,t)u(x, z, t) ex dz dx dt
8T A t ((, (14t)





h =-

(q(x,z)) t t+

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010


-i IJ ^ ^f - -- '- ~ m


x (m)

x (m)

Figure 19: Time-averaged variables for a = 0.145g.
Same conventions as in Figure 18. The
streamfunction contours show that the re-
flection symmetry and the A/2 translational
invariance are weakly broken. The stream-
lines are open, conveying the appearance of
a mean flow.

where the spatial average over x is performed in order to
minimize the numerical discretization errors. Since both
fluids are incompressible, the fluxes are in fact indepen-
dent of x.

Simulations have been run for a set of accelera-
tion values above the secondary threshold except a
0.127./ (no mean flux observable at this value of a, see
Figure 18). The corresponding mean fluxes have been
calculated and are shown in Figure 17. This figure dis-
plays three outstanding features. First, the fluxes are
monotonically increasing functions of the acceleration
and each has the shape of a supercritical pitchfork bifur-
cation. Second, the ratio of (1 to D2 remains constant.
Thirdly, for a > 0.15g, the slope of D(a) increases dras-
tically, suggesting the presence of a tertiary instability
for which additional calculations are necessary.

In order to illustrate the characteristics of the mean
flow we plot in Figures 18 to 20 the mean streamfunc-
tion for various accelerations along with the interface
position and velocity field. Figures 18 to 20 exhibit an
inversion of the colors for a = 0.145g, caused by a re-
versal of the mean horizontal velocity. This confirms
that the instability arises from a pitchfork bifurcation, in
which no direction is favored. Initially the fluids are mo-
tionless and the interface is perfectly sinusoidal, so the
sign of the mean flux is only governed by the numerical
approximation whose effects are random.

Figure 20: Time-averaged variables for a = 0.16g.
Same conventions as in Figure 18. The fea-
tures are exacerbated here.

6 Conclusion

We have carried out full nonlinear three-dimensional
simulations of Faraday waves using the method devel-
oped in Perinet, Juric & Tuckerman (2009). The incom-
pressible Navier-Stokes equations for two fluid layers
of different densities and viscosities are solved using a
finite-difference method. The interface motion and sur-
face tension are treated using a front-tracking/immersed-
boundary technique.
Our computations reproduce the hexagonal patterns
observed by Kityk et al. (2005,2009). Quantitative com-
parisons were made between experiment and simulation
of the spatio-temporal spectra. Our numerical results lie
well within the experimental uncertainty. Our direct nu-
merical simulations provide velocity fields and pressure
throughout the entire domain of calculation. Thus, we
have been able to ascertain precisely the fluid motion for
the Faraday waves, both above and below the interface
between the two fluids.
The hexagonal patterns are long-lived transients and
show intriguing dynamical behavior. We are now study-
ing the hexagonal patterns from the point of view of sta-
bility and lifetimes. To test the Ilp1, llhesi that the re-
currence of these patterns is due to the presence of a ho-
moclinic cycle, we are imposing hexagonal symmetry in
an attempt to calculate and to "freeze" an otherwise un-
stable steady state. Relaxing this constraint would then
lead to the unstable eigenvector. In our other current
study, we calculate the temporal mean of the flow result-
ing from the Faraday instability in the air/water system.
We show that a bifurcation involving the breaking of re-
flection symmetry leads to a horizontal mean flow in the
upper and lower layers.



I~ol-~ ::-----------


Our future studies of Faraday waves will include
a more detailed investigation of the dynamics of the
hexagonal patterns, and the simulation and interpreta-
tion of oscillons.


JOHNSON, B. An experimental study of the onset of
parametrically pumped surface waves in viscous fluids.
J. Fluid. Mech. 288, 325-350, 1995.

BENJAMIN, T. B. & URSELL, F. The stability of the
plane free surface of a liquid in vertical periodic motion.
Proc. R. Soc. London, Ser. A 225, 505-515, 1954.

A continuum method for modeling surface tension. J.
Comp. Phys. 100, 335-354, 1992.

CHEN, P. & VIIALS, J. Amplitude equation and pattern
selection in Faraday waves. Phys. Rev. E 60, 559-570,

CHEN, P. & WU, K.-A. Subcritical bifurcations and
nonlinear balloons in Faraday waves. Phys. Rev. Lett. 85,

CHEN, P. Nonlinear wave dynamics in Faraday instabil-
ities. Phys. Rev. E 65, 036308, 2002.

CHORIN, A. J. Numerical simulation of the Navier-
Stokes equations. Math. Comput. 22, 745-762, 1968.

M. T. Ordered capillary-wave states: quasicrystals,
hexagons, and radial waves. Phys. Rev. Lett. 68, 2157-
2160, 1992.

DOUADY, S., FAUVE, S. & THUAL, O. Oscilla-
tory phase modulation of parametrically forced surface
waves. Europhys. Lett. 10, 309-315, 1989.

EDWARDS, W. S. & FAUVE, S. Patterns and quasi-
patterns in the Faraday experiment. J. Fluid Mech. 278,
123-148, 1994.

FARADAY, M. On a peculiar class of acoustical figures;
and on certain forms assumed by groups of particles
upon vibrating elastic surfaces. Philos. Trans. R. Soc.
London 121, 299-340, 1831.

FAUVE, S., DOUADY, S.& THUAL, O. Drift instabili-
ties of cellular patterns. J. Phys. I 1, 311-322, 1991.

GODA, K. A multistep technique with implicit differ-
ence schemes for calculating two- or three-dimensional
cavity flows. J. Comput. Phys. 30, 76-95, 1979.

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

overview of projection methods for incompressible
flows. Comput. Methods Appl. Mech. Engrg. 195, 6011-
6045, 2006.

HARLOW, F. H. & WELCH, J. E. Numerical calculation
of time dependent viscous incompressible flow of fluid
with free surface. Phys. Fluids. 8, 2182, 1965.

HIRT, C. W. & NICHOLS, B. D. Volume of fluid (VOF)
method for the dynamics of free boundaries. J. Comput.
Phys. 39, 201-225, 1981.

WAGNER, C. Spatiotemporal characterization of inter-
facial Faraday waves by means of a light absorption
technique. Phys. Rev. E 72, 036209, 2005.

WAGNER, C. Erratum: Spatiotemporal characterization
of interfacial Faraday waves by means of a light ab-
sorption technique [PRE 72, 036209 (2005)]. submitted
Phys. Rev. E, 2009.

pled mean flow-amplitude equations for nearly inviscid
parametrically driven surface waves. New York Acad.
Sci. 974 201-219, 2002.

KUDROLLI, A. & GOLLUB, J.P. Patterns and spa-
tiotemporal chaos in parametrically forced surface
waves: a systematic survey at large aspect ratio, Phys-
ica D 97, 133-154, 1996.

KUDROLLI, A., PIER, B. & GOLLUB, J.P. Superlattice
patterns in surface waves. Physica D 123, 99-111, 1998.

KUMAR, K. Linear theory of Faraday instability in vis-
cous fluids. Proc. R. Soc. Lond. A 452, 1113-1126,

KUMAR, K. & TUCKERMAN, L. S. Parametric in-
stability of the interface between two fluids. J. Fluid.
Mech. 279, 49-68, 1994.

J. Dissipative solitary states in driven surface waves.
Phys. Rev. Lett. 76, 3959-3962, 1996.
Martin, Martel & Vega (2002)

MARTIN, E., MARTEL, C. & VEGA, J. Drift instability
of standing Faraday waves. J. Fluid Mech. 467, 57-79,

MULLER, H.W. Periodic triangular patterns in the Fara-
day experiment. Phys. Rev. Lett. 71, 3287-3290, 1993.

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

MURAKAMI, Y. & CHIKANO, K. Two-dimensional di- VEGA, J. M., KNOBLOCH, E., MARTEL, C. Nearly
rect numerical simulation of parametrically excited sur- inviscid Faraday waves in annular containers of moder-
face waves in viscous fluid. Phys. Fluids 13, 65-74, ately large aspect ratio. Physical D 154, 313-336, 2001.
ZHANG, W. & VIIALS, J. Pattern formation in weakly
O'CONNOR, N. L. The complex spatiotemporal dy- damped parametric surface waves. J. Fluid Mech. 336,
namics of a shallow fluid layer. Master's Thesis, Virginia 301-330, 1997.
Polytechnic Inst. and State University, 2008.

OSHER, S. & SETHIAN, J. Fronts propagating with cur-
vature dependent speed: algorithms based on Hamilton-
Jacobi formulations. J. Comput. Phys. 79, 12-49, 1988.

ical simulation of Faraday waves. J. Fluid Mech. 635,
1-26, 2009.

PESKIN, C. S. Numerical analysis of blood flow in the
heart. J. Comput. Phys. 25, 220-252, 1977.

tern control via multi-frequency parametric forcing.
Phys. Rev. Lett. 93, 034502, 2004.

RUCKLIDGE, A.M. & SILBER, M. Design of paramet-
rically forced patterns and quasipatterns, SIAM J. Ap-
plied Dynamical Systems 8, 298-347, 2009.

SHU, C. W. & OSHER, S. Efficient implementation of
essentially non-oscillatory shock capturing schemes, II.
J. Comput. Phys. 83, 32-78, 1989.

SKELDON, A. C. & GUIDOBONI, G. Pattern selection
for Faraday waves in an incompressible viscous fluid.
SIAM J. ofAppl. Math. 67, 1064-1100, 2007.

S. An improved level set method for incompressible
two-phase flows. Comput. Fluids 27, 663-680, 1998.

TEMAM, R. Une m6thode d'approximation de la solu-
tion des equations de Navier-Stokes. Bull. Soc. Math.
France 96, 115-152, 1968.

NAS, S. & JAN, Y.-J. A front-tracking method for the
computations of multiphase flow. J. Comput. Phys. 169,

merical analysis of the influence of the liquid depth on
two-dimensional Faraday waves. Phys. Fluids 15, 3099-
3113, 2003.

VALHA, J., LEWIS J. S. & KUBIE, J. A numerical
study of the behaviour of a gas-liquid interface subjected
to periodic vertical motion. Int. J. Num. Meth. Fluids 40,
697-721, 2002.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs