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
t LIMSICNRS, BP 133, Bat 508 F91403 ORSAY CEDEX FRANCE
perinet@pmmh.espci.fr
Keywords: Numerical simulation, instabilities, pattern formation, Faraday waves
Abstract
We simulate the full dynamics of hexagonal Faraday waves in three dimensions and provide comparisons of
the calculated spatiotemporal spectra to experimental results. The NavierStokes equations are solved using a
finitedifference projection method coupled with a fronttracking method for the interface between the two fluids. We
study the longtime 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 standingwave 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 quasicrystalline eightfold patterns
seen by Christiansen, Alstrom & Levinsen (1992). Other
patterns emerge using twofrequency forcing such as
twelvefold quasipatterns (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 quasipatterns, 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
terns.
A number of theoretical or seminumerical 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 coworkers, 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
(2009).
_________ 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 twodimensional. 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 threedimensional ve
locity field for the hexagonal pattern. We discuss the
longterm 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 threedimensional 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 twodimensional
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
'k~k
i /I
I .. I I
I 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
bers.
remains a singlevalued 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
quency.
The NavierStokes equations for incompressible,
Newtonian fluids are:
Du
DL
Vp + pG+V .(Vu + Vu) + s (2a)
V.u0
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 singlefluid 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 threedimensional 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
threedimensional 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
fronttracking/immersedboundary 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 noslip 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 threedimensional finitedifference mesh.
Within the domain, the two distinct immiscible flu
ids are separated by a twodimensional 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 fronttracking/immersedboundary method (Peskin
(1977); Tryggvason et al. (2001)). Because we have as
sumed that ((x, y, t) is singlevalued, 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
S(n~m)
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
height.
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 m3, 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 twodimensional 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
height.
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:
6
((x, t) () + A(k, t) eikex
j=1
6
+ A(2k6,t) Z ei2cek_ x
j1
6
+ A( 3k6,t)Ze1 i3k'
j1
+ 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
1,...6.
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
]P
0.2
0.4
0.2
S0.4
0.6
.05'O 6.I 6.15 6.2 6.25 6.3
tinle(s)
10 x
0
IIf
0
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 911. 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
0
ii/ '
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
w(mm s')
80
60
40
20
0
MvN
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 roundtrips 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
by:
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
I M
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
constraint:
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
y
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
tions:
6
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 R27/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
imaginaries,
, ,
. *
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
It i
0 0.5 1 1.5 2
x/l
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
y/c
10
20
10 0
k/k
y
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 paritybreaking 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 NavierStokes 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).
4
00 o o0
0 0& 0 O o
. 0 0 0
So OO
0 0 0 0
$000, o
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.015
0 I 
0 1
8 10 12 14 16 1S 20
ilmc(sl
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
1.2
1
0 8
0.6
04
0.125 0.13 0.135 0.14 0.145 0.15 0.155 0.16
alg
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
0.005
i SO
0.01 0.02 0.03 0.04
x (m)
Figure 18: Timeaveraged 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/2periodic.
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
percritical.
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 spatiotemporal mean the
volumetric flux of either fluid (water or air), expressed
as:
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)
(14)
So.o1
>1
0.2
02
h =
(q(x,z)) t t+
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.015
i IJ ^ ^f   ' ~ m
0.005
x (m)
x (m)
Figure 19: Timeaveraged 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: Timeaveraged 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 threedimensional
simulations of Faraday waves using the method devel
oped in Perinet, Juric & Tuckerman (2009). The incom
pressible NavierStokes equations for two fluid layers
of different densities and viscosities are solved using a
finitedifference method. The interface motion and sur
face tension are treated using a fronttracking/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 spatiotemporal 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 longlived 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.
0.015
0.Z7
I~ol~ ::
0.005
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.
References
BECHHOEFER, J., EGO, V., MANNEVILLE, S. &
JOHNSON, B. An experimental study of the onset of
parametrically pumped surface waves in viscous fluids.
J. Fluid. Mech. 288, 325350, 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, 505515, 1954.
BRACKBILL, J.U., KOTHE, D.B. & ZEMACH, C.
A continuum method for modeling surface tension. J.
Comp. Phys. 100, 335354, 1992.
CHEN, P. & VIIALS, J. Amplitude equation and pattern
selection in Faraday waves. Phys. Rev. E 60, 559570,
1999.
CHEN, P. & WU, K.A. Subcritical bifurcations and
nonlinear balloons in Faraday waves. Phys. Rev. Lett. 85,
38133816,2000.
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, 745762, 1968.
CHRISTIANSEN, B., ALSTROM, P. & LEVINSEN,
M. T. Ordered capillarywave 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, 309315, 1989.
EDWARDS, W. S. & FAUVE, S. Patterns and quasi
patterns in the Faraday experiment. J. Fluid Mech. 278,
123148, 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, 299340, 1831.
FAUVE, S., DOUADY, S.& THUAL, O. Drift instabili
ties of cellular patterns. J. Phys. I 1, 311322, 1991.
GODA, K. A multistep technique with implicit differ
ence schemes for calculating two or threedimensional
cavity flows. J. Comput. Phys. 30, 7695, 1979.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
GUERMOND, J. L., MINEV, P. & SHEN, J. An
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, 201225, 1981.
KITYK, A. V., EMBS, J., MEKHONOSHIN, V. V. &
WAGNER, C. Spatiotemporal characterization of inter
facial Faraday waves by means of a light absorption
technique. Phys. Rev. E 72, 036209, 2005.
KITYK, A. V., EMBS, J., MEKHONOSHIN, V. V. &
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.
KNOBLOCH, E., MARTEL, C. & VEGA, J. M. Cou
pled mean flowamplitude equations for nearly inviscid
parametrically driven surface waves. New York Acad.
Sci. 974 201219, 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, 133154, 1996.
KUDROLLI, A., PIER, B. & GOLLUB, J.P. Superlattice
patterns in surface waves. Physica D 123, 99111, 1998.
KUMAR, K. Linear theory of Faraday instability in vis
cous fluids. Proc. R. Soc. Lond. A 452, 11131126,
1996.
KUMAR, K. & TUCKERMAN, L. S. Parametric in
stability of the interface between two fluids. J. Fluid.
Mech. 279, 4968, 1994.
LIOUBASHEVSKI, O., ARBELL, H. & FINEBERG,
J. Dissipative solitary states in driven surface waves.
Phys. Rev. Lett. 76, 39593962, 1996.
Martin, Martel & Vega (2002)
MARTIN, E., MARTEL, C. & VEGA, J. Drift instability
of standing Faraday waves. J. Fluid Mech. 467, 5779,
2002.
MULLER, H.W. Periodic triangular patterns in the Fara
day experiment. Phys. Rev. Lett. 71, 32873290, 1993.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
MURAKAMI, Y. & CHIKANO, K. Twodimensional 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, 6574, ately large aspect ratio. Physical D 154, 313336, 2001.
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 301330, 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, 1249, 1988.
PERINET, N., JURIC, D. & TUCKERMAN, L.S. Numer
ical simulation of Faraday waves. J. Fluid Mech. 635,
126, 2009.
PESKIN, C. S. Numerical analysis of blood flow in the
heart. J. Comput. Phys. 25, 220252, 1977.
PORTER, J., TOPAZ, C.M. & SILBER, M. Pat
tern control via multifrequency 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, 298347, 2009.
SHU, C. W. & OSHER, S. Efficient implementation of
essentially nonoscillatory shock capturing schemes, II.
J. Comput. Phys. 83, 3278, 1989.
SKELDON, A. C. & GUIDOBONI, G. Pattern selection
for Faraday waves in an incompressible viscous fluid.
SIAM J. ofAppl. Math. 67, 10641100, 2007.
SUSSMAN, M., FATEMI, E., SMEREKA, P. & OSHER,
S. An improved level set method for incompressible
twophase flows. Comput. Fluids 27, 663680, 1998.
TEMAM, R. Une m6thode d'approximation de la solu
tion des equations de NavierStokes. Bull. Soc. Math.
France 96, 115152, 1968.
TRYGGVASON, G., BUNNER, B., ESMAEELI, A., JU
RIC, D., ALRAWAHI, N., TAUBER, W., HAN, J.,
NAS, S. & JAN, Y.J. A fronttracking method for the
computations of multiphase flow. J. Comput. Phys. 169,
708759,2001.
UBAL, S., GIAVEDONI, M. D. & SAITA, F. A. A nu
merical analysis of the influence of the liquid depth on
twodimensional Faraday waves. Phys. Fluids 15, 3099
3113, 2003.
VALHA, J., LEWIS J. S. & KUBIE, J. A numerical
study of the behaviour of a gasliquid interface subjected
to periodic vertical motion. Int. J. Num. Meth. Fluids 40,
697721, 2002.
