1. Introduction
In the Faraday (Reference Faraday1831) problem, the interface between two superposed fluid layers, subjected to periodic vertical vibration of sufficient amplitude, forms sustained standing wave patterns, called Faraday waves. Many universal dynamical-systems phenomena were first discovered through the Faraday experiment. In particular, a number of experimental investigations have focused on the variety of patterns formed in square (e.g. Douady & Fauve Reference Douady and Fauve1988; Simonelli & Gollub Reference Simonelli and Gollub1989; Douady Reference Douady1990) and cylindrical (e.g. Ciliberto & Gollub Reference Ciliberto and Gollub1985; Das & Hopfinger Reference Das and Hopfinger2008; Batson, Zoueshtiagh & Narayanan Reference Batson, Zoueshtiagh and Narayanan2013; Rajchenbach, Clamond & Leroux Reference Rajchenbach, Clamond and Leroux2013) containers, which have led to important theoretical developments (e.g. Meron Reference Meron1987; Silber & Knobloch Reference Silber and Knobloch1989; Crawford, Knobloch & Riecke Reference Crawford, Knobloch and Riecke1990; Crawford Reference Crawford1991; Gomes & Stewart Reference Gomes and Stewart1994) in the fields of nonlinear physics and pattern formation.
Recently, Périnet, Juric & Tuckerman (Reference Périnet, Juric and Tuckerman2009, Reference Périnet, Juric and Tuckerman2012) performed the first three-dimensional direct numerical simulation of Faraday waves in horizontally periodic domains containing a few basic cells. Here, we describe numerical simulations of Faraday waves in much larger domains which have yielded more complicated patterns. For this purpose, we have used the free-surface code developed by Shin, Chergui & Juric (Reference Shin, Chergui and Juric2014) for simulations of two-phase flows on massively parallel computer architectures. In a square domain (i.e. with square horizontal cross-section), we have obtained patterns with two different length scales, very similar to those reported in experiments by Douady & Fauve (Reference Douady and Fauve1988) and Douady (Reference Douady1990). Patterns of this type were later termed superlattices by Kudrolli, Pier & Gollub (Reference Kudrolli, Pier and Gollub1998) and Arbell & Fineberg (Reference Arbell and Fineberg2002), who observed a large variety of such patterns in their experiments, but by using a temporal forcing with two frequencies. The spectral analysis of our patterns, produced by single-frequency forcing, shows that they result primarily from the superposition of two superlattices constructed from very similar wavelengths. We have also simulated the same conditions in a cylindrical container. The resulting interface contains squares and pentagons.
This paper is organized as follows. First, we briefly present the problem and explain the methods involved in the code. We then present the numerical linear stability analysis which is validated by comparison with the theoretical treatment of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). We demonstrate the robustness of the obtained superlattice-like pattern by testing several boundary conditions. This pattern is described and analysed spectrally. Then, we show the results of simulations with the same conditions but in a cylindrical domain.
2. Problem formulation, governing equations and numerical scheme
The governing equations for an incompressible two-phase flow can be expressed by a single-field formulation:
$$\begin{eqnarray}{\it\rho}\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}\right)=-\boldsymbol{{\rm\nabla}}P+{\it\rho}\boldsymbol{G}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}})+\boldsymbol{F},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$
where
$\boldsymbol{u}$
is the velocity,
$P$
is the pressure,
${\it\rho}$
is the density,
${\it\mu}$
is the dynamic viscosity,
$\boldsymbol{G}$
represents the homogeneous volume forces and
$\boldsymbol{F}$
is the local surface tension force at the interface. Here,
$\boldsymbol{G}$
contains the gravitational term supplemented by a time-dependent force accounting for the vibrations of the domain:
$$\begin{eqnarray}\boldsymbol{G}=-(g+{\it\gamma}\cos ({\it\omega}t))\boldsymbol{e}_{z},\end{eqnarray}$$
where
$g$
is the gravitational acceleration,
$\boldsymbol{e}_{z}$
is the upward vertical unit vector,
${\it\gamma}$
is the amplitude of the inertial forcing and
${\it\omega}=2{\rm\pi}/T$
is its frequency.
The code essentially consists of two modules: one that solves the three-dimensional incompressible Navier–Stokes equations and another that treats the free surface using a parallel Lagrangian tracking method and is only active if the flow is a two-phase flow. Material properties such as density or viscosity are defined in the entire domain:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\it\rho}(\boldsymbol{x},t)={\it\rho}_{1}+({\it\rho}_{2}-{\it\rho}_{1})I(\boldsymbol{x},t)\\ {\it\mu}(\boldsymbol{x},t)={\it\mu}_{1}+({\it\mu}_{2}-{\it\mu}_{1})I(\boldsymbol{x},t).\end{array}\right\}\end{eqnarray}$$
The indicator function,
$I$
, is a numerical Heaviside function, ideally zero in one phase and one in the other phase. It is resolved with a sharp but smooth transition across three to four grid cells and is generated using a vector distance function computed directly from the tracked interface (Shin & Juric Reference Shin and Juric2009a
).
The fluid variables
$\boldsymbol{u}$
and
$P$
are calculated by a projection method (Chorin Reference Chorin1968). The temporal scheme is first order, with implicit time integration used for the viscous terms. For spatial discretization we use the staggered-mesh Marker-in-cell (MAC) method (Harlow & Welch Reference Harlow and Welch1965) on a uniform finite-difference grid with second-order essentially non-oscillatory (ENO) advection (Shu & Osher Reference Shu and Osher1989). The pressure and distance function are located at cell centres while the
$x,y$
and
$z$
components of velocity are located at the faces. All spatial derivatives are approximated by standard second-order centred differences. The treatment of the free surface uses a hybrid front-tracking/level-set technique which defines the interface both by a discontinuous density field on the Eulerian grid and by triangles on the Lagrangian interface mesh.
The surface tension
$\boldsymbol{F}$
is implemented by the hybrid formulation
$$\begin{eqnarray}\boldsymbol{F}={\it\sigma}{\it\kappa}_{H}\boldsymbol{{\rm\nabla}}\boldsymbol{I},\quad {\it\kappa}_{H}=\frac{\boldsymbol{F}_{L}\boldsymbol{\cdot }\boldsymbol{N}}{\boldsymbol{N}\boldsymbol{\cdot }\boldsymbol{N}},\end{eqnarray}$$
where
${\it\sigma}$
is the surface tension coefficient, assumed to be constant, and
${\it\kappa}_{H}$
is twice the mean interface curvature field calculated on the Eulerian grid, with
$$\begin{eqnarray}\boldsymbol{F}_{_{L}}=\int _{{\it\Gamma}(t)}{\it\kappa}_{f}\boldsymbol{n}_{f}{\it\delta}_{f}(\boldsymbol{x}-\boldsymbol{x}_{f})\,\text{d}s,\quad \boldsymbol{N}=\int _{{\it\Gamma}(t)}\boldsymbol{n}_{f}{\it\delta}_{f}(\boldsymbol{x}-\boldsymbol{x}_{f})\,\text{d}s.\end{eqnarray}$$
Here,
$\boldsymbol{x}_{f}$
is a parameterization of the time-dependent interface,
${\it\Gamma}(t)$
, and
${\it\delta}_{f}(\boldsymbol{x}-\boldsymbol{x}_{f})$
is a Dirac distribution that is non-zero only where
$\boldsymbol{x}=\boldsymbol{x}_{f}$
;
$\boldsymbol{n}_{f}$
stands for the unit normal vector to the interface and
$\text{d}s$
is the length of an interface element;
${\it\kappa}_{f}$
is twice the mean interface curvature obtained on the Lagrangian interface. The geometric information, unit normal,
$\boldsymbol{n}_{f}$
, and interface element length,
$\text{d}s$
in
$\boldsymbol{N}$
are computed directly from the Lagrangian interface and then distributed onto the Eulerian grid using the discrete delta function and the immersed boundary method of Peskin (Reference Peskin1977). A detailed description of the procedure for calculating
$\boldsymbol{F}$
,
$\boldsymbol{N}$
and
$I$
can be found in Shin (Reference Shin2007).
The Lagrangian interface is advected by integrating
$\text{d}\boldsymbol{x}_{f}/\text{d}t=\boldsymbol{V}$
with a second-order Runge–Kutta method where the interface velocity,
$\boldsymbol{V}$
, is interpolated from the Eulerian velocity.
The time step
${\rm\Delta}t$
is chosen at each iteration in order to satisfy a criterion based on
$$\begin{eqnarray}\{{\rm\Delta}t_{CFL},{\rm\Delta}t_{int},{\rm\Delta}t_{vis},{\rm\Delta}t_{cap}\},\end{eqnarray}$$
which ensures stability of the calculations. These bounds are defined by
$$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\displaystyle {\rm\Delta}t_{CFL}\equiv \min _{j}\left(\min _{domain}\left(\frac{{\rm\Delta}x_{j}}{u_{j}}\right)\right), & \displaystyle {\rm\Delta}t_{int}\equiv \min _{j}\left(\min _{{\it\Gamma}(t)}\left(\frac{{\rm\Delta}x_{j}}{\Vert \boldsymbol{V}\Vert }\right)\right),\\ \displaystyle {\rm\Delta}t_{vis}\equiv \min \left(\frac{{\it\rho}_{2}}{{\it\mu}_{2}},\frac{{\it\rho}_{1}}{{\it\mu}_{1}}\right)\frac{{\rm\Delta}x_{min}^{2}}{6}, & \displaystyle {\rm\Delta}t_{cap}\equiv \frac{1}{2}\left(\frac{({\it\rho}_{1}+{\it\rho}_{2}){\rm\Delta}x_{min}^{3}}{{\rm\pi}{\it\sigma}}\right)^{1/2},\end{array}\right\}\end{eqnarray}$$
where
${\rm\Delta}x_{min}=\min _{j}({\rm\Delta}x_{j})$
. With the periodic volume force term
$\boldsymbol{G}$
of (2.2), a supplementary upper bound
${\rm\Delta}t_{{\it\omega}}=2{\rm\pi}/(50{\it\omega})$
is required.
Fortran 2003 allows the definition of a set of dynamically allocated derived types and generic procedures associated with the grids, scalar and vector fields, operators as well as the various solvers used in the Navier–Stokes and Lagrangian tracking modules. The parallelization of the code is based on algebraic domain decomposition, where the velocity field is solved by a parallel generalized minimal residual (GMRES) method for the implicit viscous terms and the pressure by a parallel multigrid method motivated by the algorithm of Kwak & Lee (Reference Kwak and Lee2004). Communication across process threads is handled by message passing interface (MPI) procedures.
Finally, the code also contains a module for the definition of immersed solid objects and their interaction with the flow, which we have used to simulate Faraday waves in a cylindrical container. A Navier-slip dynamic contact line model is implemented where hysteresis is accounted for by fixing the contact angle limits to prescribed advancing or receding angles as in Shin & Juric (Reference Shin and Juric2009b ). Further details are available in Shin et al. (Reference Shin, Chergui and Juric2014).
3. Parameters and linear stability
We consider a layer of silicone oil of thickness
$h=14.5$
mm with density
${\it\rho}_{1}=965~\text{kg}~\text{m}^{-3}$
and dynamic viscosity
${\it\mu}_{1}=0.02~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$
. This layer is overlaid with 29.5 mm of air of density
${\it\rho}_{2}=1.205~\text{kg}~\text{m}^{-3}$
and dynamic viscosity
${\it\mu}_{2}=1.82\times 10^{-5}~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$
. The surface tension at the interface between the oil and the air is
${\it\sigma}=0.02~\text{kg}~\text{s}^{-2}$
. This choice of parameters was originally motivated by experiments on Faraday waves in a rotating cylindrical container (F. Clément, G. Pucci & Y. Couder, personal communication), which will be the subject of a future investigation.
Before simulating the patterns, we computed the critical acceleration
${\it\gamma}_{c}$
from direct numerical simulations (DNS) and compared it with that of the Floquet linear stability analysis of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). These simulations are conducted in a small periodic domain of length
${\it\lambda}_{c}$
and width
${\it\lambda}_{c}/2$
, where the critical wavelength is
${\it\lambda}_{c}=2{\rm\pi}/k_{c}=11.3~\text{mm}$
with
$k_{c}$
the critical wavenumber. We compute the growth rate of the total kinetic energy of the system (figure 1
a) after the growth becomes exponential. Near the threshold
${\it\gamma}_{c}$
, the growth rate varies linearly with
${\it\gamma}$
, so that
${\it\gamma}_{c}$
can be deduced by linear interpolation. Table 1 compares the Floquet predictions of
${\it\gamma}_{c}$
with the results from DNS for three vibration frequencies. The discrepancy is under 3 % in all three cases. The ratio of gravitational to capillary forces is measured by the Bond number, whose value here is
$g{\rm\Delta}{\it\rho}/({\it\sigma}k_{c}^{2})=1.5$
, indicating that both forces are at work.
Our large-scale simulations are carried out with a vibration frequency of
${\it\omega}/2{\rm\pi}=30~\text{Hz}$
, i.e. a period of
$T\approx 0.033$
s, and an amplitude of
${\it\gamma}=g=1.36{\it\gamma}_{c}$
. Figure 1(b) tracks the time step bounds (2.7) during the simulation. It can be seen that
${\rm\Delta}t_{cap}$
is the limiting time step throughout the simulation and is of the order of
$5\times 10^{-5}~\text{s}$
. The domain is of size
$L\times L\times L/3$
, where
$L=132~\text{mm}=11.7~{\it\lambda}_{c}$
. This domain is subdivided into
$8\times 8\times 4=256$
subdomains whose traces can be seen in figure 2(a), each of which contains
$64^{3}$
gridpoints, leading to an overall resolution of
$512\times 512\times 256$
gridpoints on a Cartesian mesh. The critical wavelength is thus resolved by
$44=512{\it\lambda}_{c}/L$
grid cells. Each subdomain is assigned to a process thread. Our initial condition has zero velocity and random noise on the interface of order
$10^{-2}~\text{mm}$
. The pattern emerges after approximately 40 forcing periods, increases exponentially, and saturates at about 90 forcing periods, requiring about 15 days of computation time on 256 threads of an IBM x3750-M4. The cylindrical container has been modelled by an impermeable solid object, as shown in figure 2(b).
Figure 1. (a) Evolution of the total kinetic energy for
${\it\omega}/2{\rm\pi}=30~\text{Hz}$
and various accelerations in a small domain. The fact that the peaks for a given acceleration lie along a line shows that the energy growth is very close to exponential. (b) Evolution of the various bounds on the time step during the course of the simulation. The smallest,
${\rm\Delta}t_{cap}$
, is used in the simulations.
Table 1. Comparison of Floquet and DNS instability thresholds for various frequencies.
4. Results and discussion
4.1. Visualizations
The top and bottom of the domain are taken to be rigid; we impose no-slip boundary conditions on the velocity field. Three types of boundary conditions have been tested on the lateral walls: no-slip, free-slip and periodic (figures 3
a–c respectively). Although our parameters were originally chosen for a different purpose, we find that the patterns that we observe are strikingly similar to those observed in earlier experiments by Douady & Fauve (Reference Douady and Fauve1988), figure 1(b), and Douady (Reference Douady1990), figures 10(b), 11(b) and 15(b), with different parameters: water, with
${\it\rho}_{1}=1000~\text{kg}~\text{m}^{-3}$
,
${\it\mu}_{1}=9\times 10^{-4}~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$
and
${\it\sigma}=0.03~\text{kg}~\text{s}^{-2}$
, in a container of dimensions
$80~\text{mm}\times 80~\text{mm}\times 5~\text{mm}$
oscillated with
$f=37~\text{Hz}$
or
$f=70~\text{Hz}$
.
For all three boundary conditions, we observe the superlattice-like patterns shown in figure 3. The interface in the
$x$
–
$y$
plane consists of four large square sub-blocks, each composed of smaller squares whose sides have approximate length
${\it\lambda}_{c}$
. In each case, the overall pattern has symmetry
$D_{2}$
, i.e. it is invariant under reflections through the two diagonals. Each of the four sub-blocks is in phase opposition with its two adjacent neighbours. In figure 3, the patterns in (a) and (c) have periodicity length
$L$
while (b) has a period of
$2L$
. Although one would expect the periodic case (c) to be the simplest to analyse, we observe that the blocks in the free-slip case (b) are more homogeneous than those in cases (a) and (c); we will therefore focus on this case.
Figure 2. The calculation domain of size
$L\times L\times L/3$
(
$L=132~\text{mm}$
), divided into
$8\times 8\times 4$
subdomains. (a) Interface profile. (b) An immersed solid cylinder of diameter
$D=L$
and height
$h=L/3$
surrounding the flow. The flow outside of the cylinder is not computed. For both (a) and (b) the Cartesian resolution per subdomain is
$64^{3}$
, which gives a global resolution of
$512\times 512\times 256$
.
Figure 3. Snapshots of the top view of the interface coloured by the vertical velocity for (a) no-slip, (b) free-slip and (c) periodic boundary conditions.
Figure 4 displays the temporal evolution of the interface at saturation. (This evolution is seen in greater detail in supplementary movie 1 to this article available at http://dx.doi.org/10.1017/jfm.2015.213.) After one period of forcing oscillation
$T$
, troughs are replaced by crests (see figure 4
a,c) as a consequence of the subharmonic behaviour of the interface. This subharmonic behaviour is also displayed by the large sub-blocks, since the appearance of each block is transformed into that of its two adjacent neighbours after time
$T$
and then returns to its initial shape after time
$2T$
has passed. The pattern is thus invariant under the combined operations of rotation by
${\rm\pi}/2$
and time translation by
$T$
, a spatio-temporal symmetry.
Figure 4. Interface profiles at times separated by intervals of
$T/2$
, where
$T$
is the forcing period and
$2T$
is the subharmonic response period. In (a) and (c), the small squares consist of wells surrounded by ridges with peaks at each corner, while in (b), low central peaks are surrounded by four higher peaks. In (a), a weak diagonal ‘bridge’ connects the large squares on the bottom left and top right; in (c), the bridge links the top left and bottom right. (See supplementary movie.)
Figure 5 shows the interface and the velocity on a vertical slice at the same three instants. The contrasting behaviour of the velocity in the left and right halves of the slice is displayed. The interface of figure 5(a) has a maximum (minimum) at the left (right) boundary, and vice versa for (c), illustrating again that the field is not periodic in a domain of length
$L$
.
Figure 5. Velocity profiles on the plane
$y=3L/4$
, corresponding to each of the three instants shown in figure 4. The vectors are coloured according to their vertical component.
4.2. Fourier spectra
For free-slip boundary conditions, the pattern is periodic when it is doubled horizontally by reflection at the boundaries. The doubled pattern is invariant under the symmetry group
$D_{4}$
of the square, consisting of reflections through the
$x$
and
$y$
coordinate axes and rotations of multiples of
${\rm\pi}/2$
. We decompose the doubled interface height profile
${\it\zeta}(x,y,t)$
into spatial Fourier modes, and then into spatio-temporal modes as follows:
$$\begin{eqnarray}{\it\zeta}(x,y,t)=\mathop{\sum }_{l,m=-\infty }^{\infty }\hat{{\it\zeta}}_{l,m}(t)\text{e}^{\text{i}\boldsymbol{k}_{l,m}\boldsymbol{\cdot }\boldsymbol{ x}}=\mathop{\sum }_{l,m=-\infty }^{\infty }\left[\mathop{\sum }_{n=-\infty }^{\infty }\hat{\hat{{\it\zeta}}}_{l,m,n}\text{e}^{\text{i}n{\it\Omega}t/2}\right]\text{e}^{\text{i}\boldsymbol{k}_{l,m}\boldsymbol{\cdot }\boldsymbol{ x}},\end{eqnarray}$$
where
$\boldsymbol{k}_{l,m}\equiv (({\rm\pi}l)/L,({\rm\pi}m)/L)$
,
$\boldsymbol{x}=(x,y)$
and
$k_{l,m}=|\boldsymbol{k}_{l,m}|$
. The square symmetry of the doubled interface implies that the eight modes in an octet
$\hat{{\it\zeta}}_{\pm l,\pm m}$
and
$\hat{{\it\zeta}}_{\pm m,\pm l}$
all have the same amplitude. Figure 6 shows two such octets, with representative elements
$\hat{{\it\zeta}}_{23,1}$
and
$\hat{{\it\zeta}}_{23,3}$
, located on circles which we will call
$\mathscr{C}_{1}$
and
$\mathscr{C}_{3}$
respectively, i.e. with
$l=23$
and
$m=1$
or
$m=3$
.
Figure 6. Spatial spectrum of a snapshot of the interface in the doubled domain. Circles
$\mathscr{C}_{1}$
(dark blue) and
$\mathscr{C}_{3}$
(dark brown) have slightly different radii. The circle
$\mathscr{C}_{1}$
contains the fundamental mode
$\hat{{\it\zeta}}_{23,1}$
and its images under reflection and rotation (dark blue dots surrounded by circles). The circle
$\mathscr{C}_{3}$
contains the highest third-harmonic mode
$\hat{{\it\zeta}}_{23,3}$
and its images (dark brown dots surrounded by squares). Only these modes and, to a lesser extent, the weaker second harmonic
$\hat{{\it\zeta}}_{46,0}$
(orange) have significant amplitude. (a) Restriction to
$[-25,25]^{2}$
highlights the square symmetry of the pattern. (b) Restriction to the first quadrant
$[0,50]^{2}$
includes
$\hat{{\it\zeta}}_{46,0}$
.
The circle
$\mathscr{C}_{1}$
has radius
$k_{c}$
, the critical wavenumber. The dominant wavevector octet of our pattern, represented by
$\boldsymbol{k}_{23,1}$
, is in accordance with the experimental observations of Douady & Fauve (Reference Douady and Fauve1988), who scanned over the forcing frequency in order to vary the selected wavenumbers and patterns. They noted that the range of existence was greatest for pairs
$(l,m)$
with a small ratio
$m/l$
and was maximum for
$m=1$
, as in our case, i.e. for wavevectors that are almost parallel to the walls.
The modes on
$\mathscr{C}_{3}$
can be seen as the result of cubic interactions between modes of
$\mathscr{C}_{1}$
:
$$\begin{eqnarray}\boldsymbol{k}_{23,3}=\boldsymbol{k}_{23,1}+\boldsymbol{k}_{23,1}+\boldsymbol{k}_{-23,1}.\end{eqnarray}$$
Modes on
$\mathscr{C}_{3}$
are fed continually and damped slowly, if at all, since
$|k_{23,3}-k_{c}|/k_{c}=0.75\,\%$
. Generalizing (4.2) to any
$m\ll l$
shows that such cubic interactions lead to
$\boldsymbol{k}_{l,3m}$
, with
$|k_{l,3m}-k_{c}|/k_{c}\sim 4(m/l)^{2}$
. Due to its proximity to the critical circle, this mode is weakly damped and should often be present when Faraday waves are excited in a square domain that is much larger than the critical wavelength. A second harmonic of
$\mathscr{C}_{1}$
, arising from the quadratic interaction
$$\begin{eqnarray}\boldsymbol{k}_{46,0}=\boldsymbol{k}_{23,1}+\boldsymbol{k}_{23,-1},\end{eqnarray}$$
is seen in figure 6(b) along with its image under
${\rm\pi}/2$
rotation. These modes occur in a quartet rather than an octet because of their reflection symmetry in the coordinate axes and are of smaller amplitude than the modes on
$\mathscr{C}_{1}$
or
$\mathscr{C}_{3}$
. The experimental techniques of Douady & Fauve (Reference Douady and Fauve1988) did not give them access to amplitudes of Fourier modes other than the dominant one.
According to the terminology of Simonelli & Gollub (Reference Simonelli and Gollub1989), the four components
$\hat{{\it\zeta}}_{\pm l,\pm m}$
comprise one pure mode, while
$\hat{{\it\zeta}}_{\pm m,\pm l}$
comprise another pure mode. Douady & Fauve (Reference Douady and Fauve1988) observed mixed modes (which they called symmetric) for small
$m/l$
, and pure modes (which they called dissymmetric) for larger
$m/l$
. Pure modes may be combined to form two types of mixed modes, the in-phase mode in which all of the components have the same sign (as in our figure 6) and the out-of-phase mode in which the components of the two pure modes are of opposite sign. Silber & Knobloch (Reference Silber and Knobloch1989) showed that these two types of mixed modes are equivalent if
$l$
and
$m$
have opposite parities. Crawford (Reference Crawford1991) showed that this was also true for
$l$
and
$m$
with the same parity (as in our case) by invoking the concept of hidden symmetries. In a square container with Neumann boundary conditions, the doubled pattern formed by reflecting at each boundary, as we have done above, satisfies periodic boundary conditions, whose inherited (hidden) translation invariances are incorporated into the normal form.
For each of
$\mathscr{C}_{1}$
or
$\mathscr{C}_{3}$
, the resulting pattern has a superlattice structure, in which squares of different sizes coexist, as shown in figure 7(a,b). Although shown in domains of length
$L$
, the periodicity length is
$2L$
. The smallest squares have length scale
$2{\rm\pi}/k_{23,m}$
for
$m=1,3$
. The patterns formed by modes on
$\mathscr{C}_{3}$
have an additional intermediate length scale
$2L/3$
. Close examination of figure 7(b) shows that adjacent medium squares are almost, but not exactly, identical. Figure 7(c) shows a superposition of the two patterns of figure 7(a,b), weighted by the amplitudes of their respective Fourier components. During most of the time, the interface resembles this superposition; compare with figure 4. When the amplitudes of modes on
$\mathscr{C}_{1}$
and
$\mathscr{C}_{3}$
approach zero, which occurs twice during one subharmonic period (see figure 8
a and figure 4
b), it is necessary to take modes like
$\hat{{\it\zeta}}_{46,0}$
into account.
Figure 7. Patterns constructed from Fourier modes in (a)
$\mathscr{C}_{1}$
, (b)
$\mathscr{C}_{3}$
and (c) a superposition of the two. Lighter (darker) zones correspond to the peaks (troughs) of the interface. Superlattices, consisting of smaller and larger squares, are observed in all of the plots; the smallest boxes are approximately of size
${\it\lambda}_{c}$
. (b) Medium squares of size
$2L/3$
are present (long-dashed box). (c) Reconstructed pattern, made by combining patterns in (a) and (b) weighted by the corresponding Fourier coefficients
$\hat{{\it\zeta}}_{23,1}$
and
$\hat{{\it\zeta}}_{23,3}$
. Compare with figure 4(c).
Figure 8. (a) Temporal evolution of the three main spatial modes. Modes
$\hat{{\it\zeta}}_{23,1}$
(blue curve and circles) and
$\hat{{\it\zeta}}_{23,3}$
(red curve and
$+$
symbols) are subharmonic: their period of oscillation is twice that of the forcing. Mode
$\hat{{\it\zeta}}_{46,0}$
(black curve and
$\times$
symbols) is harmonic. (b) Temporal spectrum of the main modes. The modes
$\hat{{\it\zeta}}_{23,1}$
,
$\hat{{\it\zeta}}_{23,3}$
contain only odd temporal harmonics while
$\hat{{\it\zeta}}_{46,0}$
has mainly even harmonics. Each spatial mode has a finite
$n=0$
component and hence a non-zero temporal mean.
The temporal evolution of
$\hat{{\it\zeta}}_{23,1}$
,
$\hat{{\it\zeta}}_{23,3}$
and
$\hat{{\it\zeta}}_{46,0}$
is displayed in figure 8(a) and the spatio-temporal spectrum is shown in figure 8(b). The mode amplitudes
$\hat{{\it\zeta}}_{23,1}$
and
$\hat{{\it\zeta}}_{23,3}$
are related by a multiplicative constant and contain only subharmonic temporal components (i.e.
$n$
odd in (4.1)). The temporal component
$n=1$
corresponding to
${\it\Omega}={\it\omega}/2$
dominates the temporal evolution. Both
$\hat{{\it\zeta}}_{23,1}$
and
$\hat{{\it\zeta}}_{23,3}$
evolve in phase opposition, crossing zero at the same time and yielding patterns such as that in figure 4(b). In contrast,
$\hat{{\it\zeta}}_{46,0}$
is harmonic and is dominated by two temporal components
$n=0$
and
$n=2$
, as expected from the quadratic resonance of
$\hat{{\it\zeta}}_{23,1}$
with itself or with
$\hat{{\it\zeta}}_{23,3}$
. Each spatial mode has a finite
$n=0$
mean component; the offset of
$\hat{{\it\zeta}}_{46,0}$
is especially noticeable.
4.3. Simulation in a cylindrical container
We performed simulations under the same conditions, but in a cylindrical container of diameter
$D=L$
and height
$h=L/3$
(see figure 2
b). The lateral boundary conditions are free-slip and the advancing and receding contact angles are fixed at
$100^{\circ }$
and
$60^{\circ }$
respectively. Figure 9 shows snapshots of the interface every
$T/2$
. (This evolution is seen in greater detail in supplementary movie 2 to this article.) The pattern has the spatio-temporal symmetry of invariance under combined time evolution by
$T/2$
and rotation by
${\rm\pi}/2$
, while each instantaneous pattern is reflection-symmetric about both diagonal axes. At the centre, we observe squares whose characteristic wavelength remains close to
${\it\lambda}_{c}$
of table 1 for
$f=30~\text{Hz}$
. Closer to the boundary, five-sided cells can be seen as the pattern organizes itself to fit in the circular domain. Aside from this, the pattern is relatively insensitive to the circular shape of the domain because
$L/{\it\lambda}_{c}=11.7\gg 1$
, in contrast to small-radius experiments and theory. In most of these (e.g. Ciliberto & Gollub Reference Ciliberto and Gollub1985; Crawford et al.
Reference Crawford, Knobloch and Riecke1990; Das & Hopfinger Reference Das and Hopfinger2008; Batson et al.
Reference Batson, Zoueshtiagh and Narayanan2013), the interface takes the form of Bessel functions, including axisymmetric modes, while Rajchenbach et al. (Reference Rajchenbach, Clamond and Leroux2013) report intriguing patterns of the form of stars and pentagons.
Figure 9. Interface profile in the cylinder at time intervals of
$T/2$
. The central part of the domain is still occupied by squares, but five-sided cells are present near the boundaries, allowing the pattern to fit inside a circle. (See supplementary movie.)
5. Conclusion
We have presented the first three-dimensional numerical results on the Faraday instability in domains much larger than the critical wavelength in both square and cylindrical containers. In the square domain, the interface displays patterns that are very similar to those found in the work of Douady & Fauve (Reference Douady and Fauve1988) and Douady (Reference Douady1990). By means of a spatial spectral decomposition, we have described these patterns quantitatively and have found that they have a complex superlattice-like structure due to resonances between modes of very similar wavelengths. We conjecture that this pattern arises from two effects. (i) In a square container that is large compared with the critical wavelength, Faraday wave patterns tend to be mixed modes whose wavevectors are almost parallel to the boundaries of the box, as stated by Douady & Fauve (Reference Douady and Fauve1988). (ii) In this case, cubic nonlinearities generate another wavelength which is very close to the dominant one. The combination of the two effects leads to the superlattice pattern that we observe in our simulations, and we infer that such patterns should be observed in a wide variety of systems. The present work demonstrates that it is now possible to numerically compute all the fields of Faraday waves in domains sufficiently large to produce complex patterns such as superlattices, quasicrystalline patterns or oscillons. The hydrodynamic mechanisms responsible for the formation of these patterns can then be explored numerically.
Acknowledgements
We thank S. Douady, S. Fauve, and G. Pucci for helpful discussions. This work was performed using high-performance computing resources provided by the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Équipement National de Calcul Intensif). N.P. was supported by a grant from CONICYT–FONDECYT/postdoctoral research project 3140522. L.S.T. acknowledges support from the Agence Nationale de la Recherche (ANR) for the TRANSFLOW project. This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and future planning (NRF-2014R1A2A1A11051346).
Supplementary movies
Supplementary movies are available at http://dx.doi.org/10.1017/jfm.2015.213.
