Faraday instability on a sphere: numerical simulation

TL;DR

This paper investigates the Faraday instability on a spherical drop using numerical simulations, revealing pattern formations like Platonic solids and oscillations driven by parametric forcing.

Contribution

It introduces a comprehensive 3D simulation approach to study spherical Faraday instability, exploring pattern formation and wave responses across various wavenumbers.

Findings

Observation of spherical harmonic wave excitation from 1 to 6.

Identification of patterns analogous to Platonic solids.

Detection of both subharmonic and harmonic responses at low viscosity.

Abstract

We consider a spherical variant of the Faraday problem, in which a spherical drop is subjected to a time-periodic body force, as well as surface tension. We use a full three-dimensional parallel front-tracking code to calculate the interface motion of the parametrically forced oscillating viscous drop, as well as the velocity field inside and outside the drop. Forcing frequencies are chosen so as to excite spherical harmonic wavenumbers ranging from 1 to 6. We excite gravity waves for wavenumbers 1 and 2 and observe translational and oblate-prolate oscillation, respectively. For wavenumbers 3 to 6, we excite capillary waves and observe patterns analogous to the Platonic solids. For low viscosity, both subharmonic and harmonic responses are accessible. The patterns arising in each case are interpreted in the context of the theory of pattern formation with spherical symmetry.

Tables

Table 1. Table 1: Comparison between theoretical and computed thresholds a c subscript 𝑎 𝑐 a_{c} for various ℓ ℓ \ell values and numerical resolutions.

Spherical mode $\ell$	Resolution	Theoretical $m{s}^{-2}$	Numerical $m{s}^{-2}$	$\frac{\left|\Delta a_c\right|}{a_c}$(%)
6	${128}^3$	0.2604	0.2640	1.38
5	${128}^3$	0.1857	0.1886	1.52
4	${128}^3$	0.09415	0.0919	2.39
6	${256}^3$	0.2604	0.2607	0.08
5	${256}^3$	0.1857	0.1864	0.34
4	${256}^3$	0.09415	0.09424	0.10

Table 2. Table 2: Comparison between decay rates and nonlinear frequencies from simulations and from low-viscosity theory of ( 13 )-( 14 ).

Viscosity	Spherical Mode	Numerical	Low-viscosity		Numerical		Deviation (%)
		resolution	theory
$\begin{array}{c}{\nu }_d\\ \left(m^2{s}^{-1}\right)\end{array}$	$\begin{array}{c}\ell \end{array}$		$\begin{array}{c}\lambda \\ \left(s^{-1}\right)\end{array}$	$\begin{array}{c}{\omega }^{\prime }\\ \left(s^{-1}\right)\end{array}$	$\begin{array}{c}\lambda \\ \left(s^{-1}\right)\end{array}$	$\begin{array}{c}{\omega }^{\prime }\\ \left(s^{-1}\right)\end{array}$	$\frac{\left|\Delta \lambda \right|}{\lambda }$	$\frac{\left|\Delta {\omega }^{\prime }\right|}{{\omega }^{\prime }}$
${10}^{-5}$	5	${128}^3$	0.122	3.711	0.113	3.607	7.53	2.83
${10}^{-5}$	5	${256}^3$			0.114	3.635	6.71	2.06
${10}^{-4}$	5	${128}^3$	1.222	3.507	0.799	3.607	34.6	4.33
${10}^{-4}$	5	${256}^3$			0.805	3.381	34.1	3.56
${10}^{-5}$	4	${128}^3$	0.075	2.661	0.070	2.615	6.66	1.72
${10}^{-5}$	4	${256}^3$			0.073	2.632	3.34	1.09
${10}^{-4}$	4	${128}^3$	0.750	2.555	0.523	2.497	30.2	2.25
${10}^{-4}$	4	${256}^3$			0.529	2.512	29.5	1.68

Table 3. Table 3: Parameters of imposed forcing and of the observed response.

	Temporal		Forcing			Initial	Final
$\ell$	Response		Amplitude	Frequency		Condition	Pattern
Gravitational			$a/a_g$	$\omega /{\omega }_g$
1	harmonic		0.30	1.00		spherical	translational
2	harmonic		0.50	1.37		axisymmetric	prolate-oblate
Capillary			$a/a_c$	$\omega /{\omega }_c$
3	subharmonic		13.53	13.21		axisymmetric	tetrahedral
4	subharmonic		5.07	16.61		axisymmetric	axisymmetric
4	subharmonic		5.07	16.61		cubic	cubic-octahedral
5	subharmonic		22.00	25.02		axisymmetric	$D_4$
6	subharmonic		16.07	30.63		icosahedral	no stable pattern

Equations

$$\zeta (t,\theta ,\varphi )=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{m=\ell }{\zeta }_{\ell }^m(t)Y_{\ell }^m(\theta ,\varphi )$$

$$\rho \left(\frac{\partial \text{u}}{\partial t}+\text{u}\cdot \nabla \text{u}\right)=-\nabla P+\rho \text{G}+\nabla \cdot \mu \left(\nabla \text{u}+\nabla {\text{u}}^T\right)+\text{F},\nabla \cdot \text{u}=0$$

$$\mathbf{G}=-\left(g+a\cos (\omega t)\right)\frac{r}{R}{\text{e}}_r$$

$$\left.\begin{array}[]{c}\rho(\textbf{x},t)=\rho_{1}+\left(\rho_{2}-\rho_{1}\right)I(\textbf{x},t)\\ \mu(\textbf{x},t)=\mu_{1}+\left(\mu_{2}-\mu_{1}\right)I(\textbf{x},t).\end{array}\right.$$

$$\text{F}=\sigma {\kappa }_H\nabla I,{\kappa }_H=\frac{{\text{F}}_L\cdot \text{N}}{\text{N}\cdot \text{N}}$$

$${\text{F}}_{{}_L}={\int }_{\Gamma (t)}{\kappa }_f{\text{n}}_f{\delta }_f\left(\text{x}-{\text{x}}_f\right)\mathrm{d}a,\text{N}={\int }_{\Gamma (t)}{\text{n}}_f{\delta }_f\left(\text{x}-{\text{x}}_f\right)\mathrm{d}a$$

$$\left\{\Delta t_{\mathrm{CFL}},\Delta t_{\mathrm{int}},\Delta t_{\mathrm{vis}},\Delta t_{\mathrm{cap}}\right\}$$

$$\begin{array}[]{ll}\displaystyle{\Delta t_{\rm CFL}\equiv\min_{j}\left(\min_{\rm domain}\left(\frac{\Delta x_{j}}{u_{j}}\right)\right)}&\displaystyle{\Delta t_{\rm int}\equiv\min_{j}\left(\min_{\Gamma(t)}\left(\frac{\Delta x_{j}}{\|{\bf V}\|}\right)\right)}\\ \displaystyle{\Delta t_{\rm vis}\equiv\min\left(\frac{\rho_{2}}{\mu_{2}},\frac{\rho_{1}}{\mu_{1}}\right)\frac{\Delta x_{\min}^{2}}{6}}&\displaystyle{\Delta t_{\rm cap}\equiv\frac{1}{2}\bigg{(}\frac{(\rho_{1}+\rho_{2})\Delta x_{\min}^{3}}{\pi\sigma}\bigg{)}^{1/2}}\end{array}$$

$$\zeta (\theta ,\varphi ,t)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{m=\ell }{\zeta }_{\ell }^m(t)Y_{\ell }^m(\theta ,\varphi )$$

$${\zeta }_{\ell }(t)\equiv {\left[\sum _{m=0}^{\ell }|{\zeta }_{\ell }^m(t){|}^2\right]}^{1/2}$$

$${\zeta }_0(t)={\zeta }_0^0={\int }_0^{2\pi }d\varphi {\int }_0^{\pi }d\theta \sin \theta \zeta (\theta ,\varphi ,t)Y_0^0(\theta ,\varphi )$$

$${\omega }^2=\frac{\sigma }{\rho }\frac{\ell (\ell -1)(\ell +2)}{R^3}$$

$${\zeta }_{\ell }(t)\propto e^{-\lambda t}|\cos ({\omega }^{\prime }t)|$$

$$\lambda =\frac{{\nu }_d(\ell -1)(2\ell +1)}{R^2}{\omega }^{\prime }={\left({\omega }^2-{\lambda }^2\right)}^{1/2}.$$

$$\text{mbox}all\ell$$

$$\ell =3$$

$$\ell =4$$

$$D_4$$

$$\ell =6$$

$${\omega }_g^2\equiv \frac{g}{R}$$

$$a_g\equiv R{\omega }_g^2=g$$

$$\zeta -R\propto \sqrt{7}{Y}_4^0+\sqrt{5}{Y}_4^4+c.c.$$

$$\zeta -R\propto \sqrt{3}{Y}_5^0+\sqrt{5}{Y}_5^4+c.c.$$

$$\zeta -R\propto \sqrt{11}{Y}_6^0+\sqrt{14}{Y}_6^5+c.c.$$

Full text

published in Journal of Fluid Mechanics 870, 433-459 (2019)

1 Introduction

The Faraday (1831) instability, in which the interface between two superposed fluid layers subjected to periodic vertical vibration of sufficient amplitude forms sustained standing wave patterns, has been instrumental in the study of pattern formation, leading to the discovery and analysis of hydrodynamic quasipatterns (Edwards & Fauve, 1994; Rucklidge & Skeldon, 2015), superlattices (Kudrolli et al., 1998; Silber & Proctor, 1998; Arbell & Fineberg, 2002), supersquares (Douady, 1990; Kahouadji et al., 2015) and other exotic patterns (Périnet et al., 2012; Rajchenbach et al., 2011).

Here, we consider a spherical analogue to the Faraday instability, a fluid drop subjected to a time-periodic radial body force. The study of such problems therefore relies particularly on numerical simulation.

Numerical studies of axisymmetric oscillating viscous drops have been carried out by Lundgren & Mansour (1988) via a boundary-integral method, by Patzek et al. (1991) and Meradji et al. (2001) via a Galerkin/finite-element method, and by Basaran (1992) via a marker-in-cell initial-value problem. Tsamopoulos & Brown (1983) used a Poincare-Lindstedt expansion method, to calculate the shapes of axisymmetric inviscid drops subjected to moderate amplitude oscillations for the three lowest capillary modes (see (1)). They were unable to confirm mathematically the existence of an asymptotic finite-amplitude motion; this can be verified only by means of numerical calculations. Indeed, the problem of axisymmetric ellipsoidal drop oscillation and decay has now come to be seen as a routine validation test case for numerical codes and interface methods for multiphase flows. However, the full non-linear problem is non-axisymmetric and requires three-dimensional numerical simulation with interface algorithms that ensure volume conservation and precise calculation of capillary forces as well as the ability to integrate highly spatially resolved systems over long physical times.

Recent advances have led to powerful general purpose codes such as Gerris (Popinet, 1993) and BLUE (Shin et al., 2017). It is with the multiphase code BLUE, which is based on recently developed hybrid front-tracking/level-set interface methods implemented on parallel computer architectures, that we conduct the current study. We have previously used this code to study large-scale square patterns of Faraday waves (Kahouadji et al., 2015). Here, we use BLUE to carry out the first numerical investigation of the Faraday problem on a sphere.

One of the most appealing aspects of the Faraday instability is that the pattern length scale is not set by the geometry, but by the imposed forcing frequency. (It is this feature which has allowed the generation of quasipatterns and superlattices, since multiple length scales can be excited simultaneously over the entire domain by superposing different frequencies.) This means that for a drop with a fixed radius $R$, patterns can be created with any wavenumber $\ell$, where the interface is described via its spherical harmonic decomposition

[TABLE]

For fixed $R$, the length scale associated with $\ell$ is $R/\ell$. More importantly, $\ell$ is associated with a set of allowed patterns; each value of $\ell$ leads to a qualitatively different situation. We will explore the motion and shape of an oscillating drop for values of $\ell$ up to 6.

2 Methods

2.1 Pattern formation on a sphere

Our previous paper (Ebo-Adou & Tuckerman, 2016) concerns the linear stability, via Floquet analysis, of the spherical Faraday problem. As is the case for all spherically symmetric problems, the equation governing the linear stability does not depend on the order $m$ of the spherical harmonic. Therefore, a bifurcation from the spherically symmetric state involves $2\ell+1$ linearly independent solutions ($Y_{\ell}^{m}$ with $-\ell\leq m\leq\ell$) with the same growth rate. The combination of these modes, i.e. the pattern, that can result from such a bifurcation is determined by the nonlinear terms. A pattern with a given $\ell$ cannot be associated with a unique combination of modes $m$, since rotation of a spherical harmonic changes $m$ (but not $\ell$).

Symmetry groups provide a classification of patterns which does not depend on orientation. A number of researchers (Busse, 1975; Busse & Riahi, 1982; Riahi, 1984; Ihrig & Golubitsky, 1984; Golubitsky et al., 1988; Chossat et al., 1991; Matthews, 2003) have studied the patterns which are allowed and those which are preferred for various values of $\ell$. Patterns that appear at bifurcations can be associated with a subgroup of the group $O(3)$ of symmetries of the sphere. The subgroups of interest are $O(2)$, $D_{k}$, and the exceptional subgroups $T$, $O$, and $I$. $O(2)$ consists of the symmetries of a circle; patterns which are axisymmetric (for some orientation) are in this category. $D_{k}$ describes the symmetries of a $k$-gon, and thus patterns which are invariant under reflection and rotation by $2\pi/k$ (about a fixed axis). The three exceptional subgroups are associated with the five Platonic solids: $T$ (tetrahedron), $O$ (octahedron or cube), and $I$ (icosahedron or dodecahedron).

The Platonic solids are regular polyhedra. Although a drop does not have angular vertices and flat faces, a polyhedron with the same symmetry properties can be constructed from a drop in a $T$, $O$, or $I$ configuration by assigning a local maximum on the surface to a vertex and a local minimum to a face. The dual of a polyhedron is obtained by inverting its vertices with its faces, an operation which preserves symmetry. Similarly, the dual of a drop can be formed by inverting maxima and minima. An oscillating drop provides an ideal opportunity to observe duality: since a location which is the site of a maximum contains a minimum after half of an oscillation period, the interface alternates between a pattern and its dual. For the Platonic solids, the dual of an octahedron is a cube, that of an icosahedron is a dodecahedron and the dual of a tetrahedron is another tetrahedron.

The patterns are highly dependent on the value of $\ell$ considered. Axisymmetric (sometimes called zonal) solutions are never stable if $\ell\geq 2$ (Chossat et al., 1991). For odd $\ell$, solutions exist which are stable at onset. In contrast, for even $\ell$, all solutions are unstable near the bifurcation point. In this case, the preferred solution can be considered to be that with the smallest number of unstable eigenvalues (Matthews, 2003) or which extremizes a functional (Busse, 1975). Unstable solution branches produced at a transcritical bifurcation can be stabilized, for example at a saddle-node bifurcation some distance from the threshold.

The theory of pattern selection with $O(3)$ symmetry differs from the framework of our simulations in some important ways. First, the theory applies to steady bifurcations rather than oscillatory solutions resulting from time-periodic forcing. The archetypical application is Rayleigh-Bénard convection in a sphere (Busse, 1975; Busse & Riahi, 1982; Riahi, 1984) and the steady symmetry-breaking bifurcations it undergoes. However, many of the conclusions can be generalized to the time-periodic context, by considering the discrete-time dynamical system derived by sampling the continuous-time system at a single phase of the forcing period.

Secondly, the theory applies close to the threshold. Our simulations are carried out far from threshold, so that instabilities can grow on a reasonable timescale and stabilize at an amplitude that can be clearly seen. Our patterns contain modes generated by nonlinear interactions and their existence or stability may result from secondary bifurcations.

Despite these differences, we will see that there is much common ground between the patterns we observe and those predicted by theory.

2.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:

[TABLE]

where u is the velocity, $P$ is the pressure, $\rho$ is the density, $\mu$ is the dynamic viscosity and F is the local surface tension force at the interface. Here, ${\bf G}$ is an imposed time-dependent radial acceleration:

[TABLE]

where $g$ is a constant acceleration, referred to for simplicity as gravitational, which is set to zero when we carry out capillary simulations. $R$ is the radius of the drop, $\textbf{e}_{r}$ is the radial unit vector, and $a$ and $\omega$ are the amplitude and frequency of the oscillatory forcing.

Material properties such as density or viscosity are defined in the entire domain:

[TABLE]

The indicator function $I$ in (4) is the Heaviside function, whose value is zero in one phase and one in the other phase. In our discrete numerical implementation $I$ is approximated by $I_{\rm num}$, which represents a smooth transition across 3 to 4 grid cells, as pioneered in the immersed boundary method of Peskin (1977). In our method $I_{\rm num}$ is generated using a vector distance function computed directly from the tracked interface (Shin & Juric, 2009).

The fluid variables u and $P$ are calculated by a projection method (Chorin, 1968). 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, 1965) on a uniform finite-difference grid with second-order essentially non-oscillatory (ENO) advection (Shu & Osher, 1989). The pressure and distance function are located at cell centers while the $x$, $y$ and $z$ components of velocity are located at the faces. All spatial derivatives are approximated by standard second-order centered differences. The treatment of the free surface uses a hybrid Front-Tracking/Level-Set technique which defines the interface both by the Level-Set distance function field on the Eulerian grid as well as by triangles on the Lagrangian interface mesh.

The surface tension F is implemented by the hybrid/compact formulation (Shin, 2007)

[TABLE]

where $\sigma$ is the surface tension coefficient and $\kappa_{H}$ is twice the mean interface curvature field calculated on the Eulerian grid, with

[TABLE]

Here, $\textbf{x}_{f}$ is a parameterization of the time-dependent interface, $\Gamma(t)$, and $\delta_{f}\left(\textbf{x}-\textbf{x}_{f}\right)$ is a Dirac distribution that is non-zero only where $\textbf{x}=\textbf{x}_{f}$; $\textbf{n}_{f}$ stands for the unit normal vector to the interface and $\mathrm{d}a$ is the area of an interface element; $\kappa_{f}$ is twice the mean interface curvature obtained on the Lagrangian interface. The geometric information, unit normal, $\textbf{n}_{f}$, and interface element length, $\mathrm{d}a$ in 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 (1977). A detailed description of the procedure for calculating F, N and $I_{\rm num}$ can be found in Shin & Juric (2007, 2009), where in particular we demonstrate that this method of calculating the surface tension force reduces any parasitic currents in the standard static drop test case to a level of $O(10^{-7})$ for fluids with properties similar to those we use here.

The Lagrangian interface is advected by integrating $\mathrm{d}\textbf{x}_{f}/\mathrm{d}t=\textbf{V}$ with a second-order Runge-Kutta method where the interface velocity, V, is interpolated from the Eulerian velocity.

The parallelization of the code is based on algebraic domain decomposition, where the velocity field is solved by a parallel generalized minimum residual (GMRES) method for the implicit viscous terms and the pressure by a parallel multigrid method motivated by the algorithm of Kwak & Lee (2004). Communication across process threads is handled by message passing interface (MPI) procedures.

The code 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 spherical container. In order to simulate a fluid within a solid sphere we take the simple approach of defining all grid cells whose centers lie within the solid region as solid. Then Dirichlet (no-slip) boundary conditions for the velocity and Neumann conditions for the pressure are applied to those cell faces as in the projection method. On a Cartesian grid this necessarily creates a stair-stepped solid/fluid boundary; however it is found that the method works well in practice since the discrete momentum flux is conserved and is simpler than and equivalent to other approaches which impose a near-wall force to ensure a no-slip condition at the solid (Tryggvason et al., 2011). See Shin et al. (2017) for further details.

2.3 Physical parameters of the fluids

Our numerical code can treat inner and outer spherical domains of any size containing fluids of any density $\rho$, viscosity $\nu$, and surface tension $\sigma$, leading potentially to a large number of non-dimensional parameters, in addition to those describing the forcing amplitude $a$ and frequency $\omega$. We have chosen to limit the parameter space as follows.

The density and viscosity ratios that we have chosen are typical of oil droplets in air. The inner and outer densities are $\rho_{d}=940$ kg/m3 and $\rho_{\rm out}=1.205$ kg/m3, leading to the density ratio $\rho_{d}/\rho_{\rm out}=780$, while the inner and outer kinematic viscosities are $\nu_{d}=10^{-5}$ m2/s and $\nu_{\rm out}=1.5\times 10^{-5}$ m2/s, leading to the viscosity ratio of $\nu_{d}/\nu_{\rm out}=0.66$. The viscosity of both fluids is sufficiently low that the stability diagram approaches that of the Mathieu equation. (One of the conclusions of Ebo-Adou & Tuckerman (2016) is that the Mathieu equation describes the inviscid Faraday instability even in a spherical geometry.) More specifically, $\nu_{d}/(R^{2}\omega)$ is between $3\times 10^{-4}$ and $7\times 10^{-3}$ for the capillary cases and between $6\times 10^{-4}$ and $9\times 10^{-4}$ for the gravitational cases. Usually, Faraday waves are subharmonic, i.e. their period is twice the forcing period $T$. However, in this low-viscosity regime, we can easily excite harmonic waves as well (Kumar, 1996), whose period is the same as the forcing period.

For the capillary cases, we set the surface tension to be $\sigma=0.02$ kg/s2, the constant radial acceleration to be $g=0$, and the drop radius to be $R=0.06$ m, while for the gravitational cases, we set the constant radial acceleration to be $g=-1$ m/s2, the surface tension to be $\sigma=0$, and the drop radius to be $R=0.05$ m. We will use these parameters to non-dimensionalize the forcing frequencies and amplitudes as in equation (16) and list these in section 2.7 and table 3.

2.4 Numerical parameters

For the capillary cases, the domain is a sphere whose radius is $2R$ (see section 2.2 and figures 9 and 9). The Eulerian mesh is uniform and Cartesian, with a resolution of $N_{x}^{3}$, which is generally $128^{3}$. The Lagrangian triangular grid used to represent the interface is constructed in such a way that the sides of the triangles remain close to the length of the diagonals of the cubic Eulerian mesh. The grid spacing $\Delta x=4R/N_{x}$ for the capillary cases should be compared to the circumference $2\pi R$ of the drop and to the approximate wavelengths $2\pi R/\ell$. The radius and circumference are spanned by $R/\Delta x=N_{x}/4$ and $2\pi R/\Delta x=N_{x}\pi/2$, which are 32 and 201, respectively for the $128^{3}$ grid. The number of gridpoints per wavelength $N_{x}\pi/(2\ell)$ goes from 201 for $\ell=1$ to 33 points for $\ell=6$ for this grid. The interface is represented by approximately $4\pi R^{2}/(\Delta x)^{2}=N_{x}^{2}\pi/4$ points, i.e. 12,868 points for $N_{x}=128$. For the gravitational cases, the radius of the bounding sphere is $2.4R$, leading to slightly different but similar numbers for the resolution, e.g. $R/\Delta x=26.7$ rather than 32. We have also carried out simulations with grids of $256^{3}$ (to confirm our results) and of $64^{3}$ (for visualisation of our patterns). Parallelization is achieved through domain decomposition, in which each subdomain is assigned to its own process thread. Here, we use $512$ subdomains or processes, each with $16^{3}$ or $32^{3}$ gridpoints for the global resolutions of $128^{3}$ and $256^{3}$, respectively.

The time step $\Delta t$ is chosen at each iteration in order to satisfy a criterion based on

[TABLE]

which ensures stability of the calculations. These bounds are defined by:

[TABLE]

where $\Delta x_{\min}=\min_{j}\left(\Delta x_{j}\right)$. In our simulations, this minimum is realized by $\Delta t_{\rm CFL}\approx\Delta t_{\rm int}$ for our gravitational waves and by $\Delta t_{\rm cap}\approx\Delta t_{\rm vis}$ for our capillary waves.

2.5 Spherical harmonic transform

In order to analyze the shape of the drop quantitatively, we will compute and present the time-dependent spectral coefficients, obtained from the decomposition

[TABLE]

where $\zeta$ is the distance from the domain center. We compute both $|\zeta_{\ell}^{m}(t)|$ and

[TABLE]

We note that, according to the normalization convention we use,

[TABLE]

i.e. $\zeta_{0}(t)$ is $\sqrt{4\pi}$ times the spherically averaged distance of the interface from the origin.

We calculate the spherical harmonic transform as follows. The three-dimensional Lagrangian interface consists of triangles composed of points $\left(x_{i},y_{i},z_{i}\right)$ on the fixed Cartesian grid. We transform each of these points to spherical coordinates $\left(\zeta_{i},\theta_{i},\phi_{i}\right)$ and interpolate the function $\zeta(\theta_{i},\phi_{i})=\zeta_{i}$ onto a regular spherical grid $(\theta_{j},\phi_{k})$, discretizing the latitudinal interval $[0,\pi]$ and longitudinal interval $[0,2\pi]$ with $N_{\theta}=80$ and $N_{\phi}=80$ points. Using the discrete data $\zeta(\theta_{j},\phi_{k})$ and truncating the series (9) at four times the dominant value of $\ell$, we compute the spherical harmonic coefficients $\zeta_{\ell}^{m}$ by the method of least squares (Politis, 2013, 2016). This procedure is carried out for each sampled temporal snapshot of the interface. We note that the spherical harmonic coefficients are weighted averages of the surface height over the entire drop surface, and hence more accurate than the individual surface height values.

The spectrum in $m$ (but not in $\ell$) depends on the orientation. We can rotate a pattern by localizing the coordinates $(\theta,\phi)$ of one of its features, such as a maximum, rotating in longitude $\phi$ to place it in the $(x,z)$ plane and then rotating in colatitude $\theta$ to place this feature at the north pole. This allows us to interpret our spectra by using the explicit representations given by Busse (1975) of various patterns in terms of $Y_{\ell}^{m}$. Conversely, if the shape of the pattern is known or constant, the spectrum in $m$ can be used to track changes in its orientation.

2.6 Validation

The front-tracking approach inherently conserves mass to high accuracy compared to other numerical interface methods (Shin & Juric, 2009). We confirm mass conservation by showing the fractional deviation of the volume. Figure 1 shows the time evolution of the bounding envelope of $\Delta V/V$ for a simulation with resolution $128^{3}$. The fluid parameters are those of section 2.3 and table 3 for the $\ell=6$ case. (See figure 12 for an illustration of the extraction of the bounding envelope.) We see that $|\Delta V(t)|/V$ remains less than about $3\times 10^{-6}$ throughout the simulation.

In our previous investigation (Ebo-Adou & Tuckerman, 2016), we extended the method of Kumar & Tuckerman (1994) for computing the Faraday threshold via Floquet theory to a spherical geometry. In figure 3 and table 1 we compare these theoretical results with thresholds obtained by interpolating growth rates from numerical simulations, as shown in figure 3. The parameters are as previously stated for capillary waves in section 2.3, except that we increase the viscosity to $\nu_{d}=10^{-4}m^{2}s^{-1}$. For $\ell=6$, the error in the threshold is about 1% for a resolution of $128^{3}$ and only about 0.1% for $256^{3}$.

Figure 4 compares the time evolution of the envelope of the spherical harmonic coefficients $|\zeta_{\ell}^{m}|$ for capillary wave simulations using resolutions $128^{3}$ and $256^{3}$. The fluid parameters are again those of section 2.3 and table 3 for the $\ell=6$ case. Starting from an icosahedral initial condition, the solution evolves similarly for the two spatial resolutions. The main difference is that near-zero modes grow more slowly for the higher resolution, which can be understood as a manifestation of the lower level of noise introduced in the simulation.

We also validated our numerical code by testing the free decay of a perturbed sphere in the unforced case (i.e. $a=0$, no imposed oscillatory forcing) under the influence of surface tension. In the absence of viscosity, Rayleigh (1879) showed that a drop initialized with a perturbation of spherical wavenumber $\ell$ oscillates with frequency $\omega$ where

[TABLE]

For small kinematic viscosity $\nu_{d}$, Lamb (1932) showed that the oscillation amplitude decays to the spherical rest state according to

[TABLE]

where the decay rate and the oscillation frequency are given by

[TABLE]

Figure 5 shows the evolution of $\Delta\zeta\equiv\zeta_{max}-\zeta_{min}$ obtained via numerical simulations with a resolution of $256^{3}$ of capillary waves with the fluid parameters of section 2.3 initialized with axisymmetric perturbations $Y_{5}^{0}$ for two different viscosities, along with curves obtained from the low-viscosity theory of (13)-(14). Table 2 gives the values of the decay rates and nonlinear frequencies from figure 5, along with additional cases in which the spherical mode and the spatial resolution have been varied. As could be expected, the deviation is lowest when the resolution is higher ($256^{3}$) so that the simulation is more accurate, and when the viscosity is lower ($\nu_{d}=10^{-5}ms^{-2}$) so that the low-viscosity theory is more applicable. This best-case deviation is about 7% for $\lambda$ and 2% for $\omega^{\prime}$ for $\ell=5$ and only about 3% for $\lambda$ and 1% for $\omega^{\prime}$ for $\ell=4$.

2.7 Survey of cases studied

Table 3 lists the simulations which we will describe in this paper. As an initial condition, we perturb the spherical interface by one of the following combinations of spherical harmonics:

[TABLE]

These formulas for functions with a given symmetry are given by Busse (1975) and Riahi (1984) for patterns with a single $\ell$, aligned along the $z$ axis. When the patterns are rotated to a different orientation, the value of $\ell$ is conserved, but the combinations of $m$ values change. Patterns with other symmetries are also possible, but are not used or not achieved here. For $\ell=1$, the initial condition for the interface is a sphere perturbed only by its representation on a triangular mesh. In all cases, the initial velocity is zero. The quantitative results we will present in the next section use a resolution of $128^{3}$, confirmed by simulations with $256^{3}$. Grids with $64^{3}$ are used to plot the visualisations of the drop.

The physical parameters are as given in section 2.3. We study either gravity or capillary waves, i.e. a Bond number $\rho gR^{2}/\sigma$ of either infinity or zero. Figure 6 locates the parameter values that we have used for our simulations within the instability tongues for the gravitational or capillary cases, where

[TABLE]

See Ebo-Adou & Tuckerman (2016) for more details on this non-dimensionalization. Because our viscosity is low, the tongues are very close to the inviscid ones; we plot the inviscid tongues for simplicity.

Frequencies were chosen to induce instabilities from $\ell=1$ to $\ell=6$, as predicted from linear Floquet theory. In each case, the value of $\ell$ from the full three-dimensional nonlinear numerical simulations agreed with the theoretical value.

3 Results

3.1 Case $\ell=1$

We begin by presenting the $\ell=1$ case in the purely gravitational regime, i.e. in the presence of a constant radial force $g{\bf e}_{r}$ included in the time-periodic force (3) and without surface tension. Our simulations of this case exhibit a periodic subharmonic translational motion of the sphere about its original position, as shown in figure 7 and supplementary movie 1. This explains why the $\ell=1$ case is prohibited for capillary waves: translational motion does not involve deformation of the interface, and so surface tension cannot act as a restoring force. (Note the factor of ($\ell-1$) in formulas (12) and (14).) Although the drop arrives quite close to the bounding sphere in this simulation, it does not touch it.

As stated by Busse (1975), there is no competition between different patterns in this case. All solutions are obtained by rotation of a single axisymmetric pattern, which is expected to be stable. (The three spherical harmonics $Y_{1}^{1},Y_{1}^{0},Y_{1}^{-1}$ are also related by rotation.)

3.2 Case $\ell=2$

To produce a case in which $\ell=2$ is the dominant mode, we have simulated harmonic oscillations in the gravitational regime. A sequence of figures from our numerical simulations showing the shape and velocity field at the interface is shown in figure 9 and figure 9, which also shows the velocity field in the outer fluid; see supplementary movie 2. The interface remains axisymmetric, with two principal axes of the same length, and a third of a different length. Figure 9(a) shows a prolate shape, like a rugby ball (the third axis is longer than the other two), while the shape in figure 9(c) is oblate, like a disk (the third axis is shorter). During the prolate phase, the velocity is directed in the polar direction, while during the oblate phase it is directed in the equatorial direction. At maximum deformation, the velocity changes direction, so that the magnitude of the velocity is lowest when the surface is most deformed and highest when it is least deformed. Because the interface is opaque, velocity vectors directed inwards cannot be seen.

Figure 12 shows the amplitude $|\zeta_{\rm max}-\zeta_{\rm min}|$ as a function of time. The prolate and oblate configurations are represented by higher and lower maximum values of $|\zeta_{\rm max}-\zeta_{\rm min}|$, respectively. Oblate-to-prolate oscillations have been studied extensively (Trinh & Wang, 1982; Tsamopoulos & Brown, 1983; Patzek et al., 1991). These authors observe that the drop spends a longer time in the prolate than in the oblate configuration, which agrees with our observation that the drop spends about 53% of each period in the prolate configuration. This is explained by the combination of two facts: the restoring force is governed by the pressure at the curved surface, and the poles cover a smaller surface than the equator. The prolate rugby-ball-like form, with high curvature at the poles, is therefore subjected to a smaller restoring force than the oblate disk-like form, with high curvature at the equator. Hence it takes more time for the drop to return from the prolate to the spherical shape than it takes to return from the oblate shape.

In figures 12 and 12, we present the spherical harmonic coefficients of the interface height $\zeta(\theta,\phi,t)$, calculated from equations (9) and (10) as explained in section 2.5. Figure 12 shows the time evolution of $\zeta_{\ell}$, defined in (10). Its most visible feature is its rapid oscillation. In order to examine the dynamics over timescales much larger than the forcing period $T$, we extract the envelope of each $\zeta_{\ell}$, as shown by the bold curves in figure 12. The variation of $\zeta_{0}$ cannot be seen on this scale, since it is $\sqrt{4\pi}$ times the mean radius, which should be nearly constant due to incompressibility. Although $\ell=2$ is the dominant non-zero spherical wavenumber for this case, other even $\ell$ (multiples of 2) are also present, generated by nonlinear interactions. Standard bifurcation theory predicts that the growth rate of the harmonic $n\ell$ should be approximately $n$ times that of $\ell$, shown here by the fact that $\zeta_{4}$ grows about twice as quickly as $\zeta_{2}$.

Although we began our simulations by perturbing the sphere with an axisymmetric initial condition, proportional to $Y_{2}^{0}$, the drop orientation quickly tilts away from the $z$ axis, acquiring components with $m\neq 0$ while remaining axisymmetric about its own axis. The symmetry axis continues to rotate and to oscillate periodically, manifested by the behavior of components $m=0$, 1 and 2, shown in figure 12. We have also obtained prolate-oblate oscillations in simulations (not presented here) of subharmonic capillary waves with the same axisymmetric initial condition.

Busse (1975) and Chossat et al. (1991) show that the only allowed solution for $\ell=2$ is axisymmetric and unstable at onset. In our simulations, which are far from onset, we observe an axisymmetric pattern whose orientation continually changes.

3.3 Case $\ell=3$

For $\ell=3$, there exist solutions with three possible symmetries: axisymmetric, $D_{6}$, and tetrahedral. Either the $D_{6}$ or the tetrahedral solution can be stable at onset (Busse, 1975; Chossat et al., 1991). Starting from an axisymmetric initial condition, our solution for this case rapidly develops tetrahedral symmetry, as can be seen in figure 15 and supplementary movie 3, which show subharmonic capillary oscillations. Over a half oscillation period, the tetrahedron in figure 15(a) reverses its orientation, as shown in figure 15(c), since this polyhedron is self-dual. The initial condition is axisymmetric, but the same behavior is seen starting from a sphere slightly perturbed by numerical noise.

The pattern contains components other than $\ell=3$ and its harmonics, as illustrated in figure 15. In particular, an important $\ell=4$ component is present, of approximately the same magnitude as the first harmonic $\ell=6$. Figure 15 shows the rapid progression from the axisymmetric initial condition to a tetrahedral one, as component $m=0$ falls and the other components rise. There follows a long phase during which the orientation of the tetrahedal pattern varies. Eventually, by $t/T\approx 2000$, the tetrahedron aligns with the $z$ axis, with only $m=2$ remaining; recall from equation (15b) that $Y_{3}^{2}+{\rm c.c.}$ corresponds to a tetrahedron aligned with the $z$-axis (Busse, 1975).

3.4 Case $\ell=4$

When $\ell=4$, bifurcation theory predicts two possible solutions at onset (Busse, 1975; Chossat et al., 1991), with axisymmetric or cubic symmetry, produced by a transcritical bifurcation. Neither branch is stable at onset but the cubic solution is the preferred one (Busse, 1975; Matthews, 2003). We have been able to produce patterns of both kinds by starting with initial conditions which are axisymmetric or cubic

[TABLE]

Visualisations of our numerical simulations of subharmonic capillary oscillations with cubic symmetry are shown in figure 18 and supplementary movie 4a. The pattern oscillates between resembling a cube, with six square faces, of which three meet at each vertex, and its dual, the octahedron with eight equilateral triangles, four of which meet at each corner.

Figure 18 shows that, in addition to $\ell=4$ and its harmonics, the $\ell$ wavenumber spectrum contains a large $\ell=6$ component, as well as modes resulting from interactions between the $\ell=4$ and $\ell=6$ families. The pattern is aligned with the numerical domain, with two faces or minima (for the cube) or two vertices or maxima (for the octahedron) located at the north and south poles. With this alignment, the cubic solution is expected to be a sum of $Y_{4}^{0}$ and $Y_{4}^{\pm 4}$. Indeed figure 18 shows that $\zeta_{4}^{0}$ and $\zeta_{4}^{4}$ are dominant and constant. We have also obtained a cubic pattern in simulations (not presented here) of harmonic gravitational waves with a perturbed spherical initial condition and a resolution of $256^{3}$.

From an axisymmetric initial condition, the solution remains axisymmetric over the time of our simulation, but develops a tilting and rotating axis of symmetry, as was the case for $\ell=2$. Visualisations of this state are shown in figure 21 and supplementary movie 4b. Figure 21 shows the evolution of the amplitudes of the various $\ell$ modes. We observe that $\zeta_{2}$ is higher (a feature we have observed whenever a pattern is axisymmetric) and $\zeta_{6}$ is lower than was the case for the cubic solution. The oscillatory behavior of the various $m$ modes in figure 21 are manifestations of the fact that the axis of symmetry rotates, as in figure 12.

3.5 Case $\ell=5$

The theoretical analysis of the $\ell=5$ case is the most complicated of those surveyed here (Busse & Riahi, 1982; Riahi, 1984), since it leads to eight allowed solutions with different symmetries (Chossat et al., 1991). We find two of these solutions in our simulations, one as a short-lived unstable equilibrium and the other as a long-term asymptotic state.

Starting from an initial condition which is an axisymmetric perturbation of the sphere, the solution quickly equilibrates to axisymmetric oscillations which are shown in figure 23. The solution is then replaced by another solution with $D_{4}$ symmetry, shown in figure 23. See supplementary movie 5. In the $D_{4}$ solution, the upper hemisphere contains four patches of outward (inward) flow, and the lower hemispere contains four similar patches, located at longitudes which are halfway between those of the patches of the upper hemisphere. This pattern corresponds to the preferred solution called $\mathbf{F}$ by Riahi (1984):

[TABLE]

The spectra in $\ell$ and $m$ are shown in figures 25 and 25. Because the pattern remains aligned with the $z$ axis, the transition from the axisymmetric pattern to (18) can be tracked by following the different $m$ modes, as is done in figure 25. The initial plateau in $\zeta_{5}^{0}$ indicates that the axisymmetric phase of the oscillations shown in figure 23 comprises a solution, albeit unstable, of our system. (This is accompanied in figure 25 by a brief slowing of the decay of the $\ell=2$ mode which is present in all of our nonlinear axisymmetric solutions.) The subsequent decrease in $\zeta_{5}^{0}$ is accompanied by an increase in $\zeta_{5}^{4}$, which together comprise the stable solution with $D_{4}$ symmetry.

3.6 Case $\ell=6$

We compute subharmonic capillary oscillations for $\ell=6$. For $\ell=6$, there exist four possible solutions: axisymmetric, six-fold symmetric ($D_{6}$) octahedral, and icosahedral. Of these, the icosahedral solution is preferred near threshold, via a maximization argument (Busse, 1975) or a stability argument (Matthews, 2003). With this in mind, we use as an initial condition an icosahedral perturbation of the sphere (Busse, 1975):

[TABLE]

Figure 27 shows the icosahedral/dodecahedral oscillations during the first phase of the simulation, with a clear five-fold symmetry; see also supplementary movie 6. Figure 27(a) resembles a dodecahedron, i.e. pentagons which meet in sets of three at the vertices, while 27(c) resembles an icosahedron, i.e. triangles which meet in sets of five at the vertices. Further evolution leads to a second phase of the solution, shown in figure 27, which oscillates between a cube and its dual, an octahedron. Indeed, simulations from a spherical initial condition (perturbed by its discrete representation on a triangular mesh) also led directly to the cubic/octahedral solution.

Figure 28 shows the time evolution of the spectrum in $\ell$, while figure 29 shows the time evolution of the most important $m$ modes associated with $\ell=6$. The icosahedral pattern, consisting of modes $m=0$ and $m=5$, quickly tilts away from the $z$ axis and by $t/T\approx 300$ is replaced by a cubic pattern, whose dominant mode is $m=4$. The subsequent increase in $m=3$ is associated with the tilting of the cubic pattern away from the $z$ axis. Eventually, however, the solution shifts to the pattern shown in the rightmost insets of figure 29, which is neither icosahedral nor cubic.

4 Discussion

Using the parallel front-tracking code BLUE, we have been able to simulate the spherical Faraday problem and thereby to produce patterns with all spherical harmonic wavenumbers between $\ell=1$ and $\ell=6$. We have simulated both gravitational and capillary waves, in both the harmonic regime and the more usual subharmonic regime.

Our simulations agree in most cases with two types of theory. First, the spherical harmonic wavenumber $\ell$ obtained in each case agrees with the results of Floquet analysis, presented in Ebo-Adou & Tuckerman (2016) and in figure 6. Second, the interface shapes we observe are readily interpreted using the theory of pattern formation on the sphere (Busse, 1975; Busse & Riahi, 1982; Riahi, 1984; Ihrig & Golubitsky, 1984; Golubitsky et al., 1988; Chossat et al., 1991; Matthews, 2003). For $\ell=1$ and $\ell=2$, only one type of pattern is possible, and that is the one we observe. For our vibrating drop, the $\ell=1$ “pattern” is manifested as a back-and-forth motion of the spherical drop, while the $\ell=2$ pattern alternates between an oblate and a prolate spheroid. For $\ell=3$, we observe a tetrahedral pattern, one of the two solutions predicted to be stable (out of the three which can exist), while for $\ell=4$, depending on the initial conditions, we observe both of the possible solutions, a cubic/octahedral or an axisymmetric pattern. For $\ell=5$, we observe an unstable axisymmetric and a stable $D_{4}$ pattern. These are two of the eight possible solutions; $D_{4}$ is one of those which can be stable at onset. For $\ell=6$, an initially icosahedral solution makes a transition to a cubic/octahedral pattern, which is succeeded by a solution with neither symmetry. It is surprising to find such good agreement with theory given the differences with our configuration mentioned in the introduction, i.e. the fact that our patterns are oscillatory rather than steady and that our parameters are far above threshold. A complete study of each of these cases, varying the forcing amplitude from threshold to higher values, would be desirable. Another crucial issue, both empirical and mathematical, is the possible difference between the harmonic and subharmonic regimes in each of the cases.

An important avenue of exploration that has arisen in our study is that of the long-term dynamics. The $\ell=2$ (prolate-oblate) and $\ell=4$ (axisymmetric) show the drop tumbling into and out of alignment with the coordinate system, as illustrated in figures 12 and 21. The $\ell=3$ (tetrahedral) and $\ell=6$ (cubic/octahedral) oscillations in figures 15 and 29 also show very long phases (2000 or more forcing periods $T$) during which the orientation continues to change. Another case in which very long-term Faraday-wave dynamics has been found is that of hexagonal waves in a minimal domain (Périnet et al., 2012). To the best of our knowledge, there exists no explanation of such long-term dynamics in terms of pattern formation, fluid dynamics, or any other kind of theory.

We have demonstrated the feasibility of well-resolved numerical simulation of the spherical version of the Faraday instability over extremely long times. We believe that the Faraday problem serves as a rigorous proving ground for numerical interface techniques. The wide variety of drop shapes and their detailed patterns simulated here demonstrate the necessity and advantages of using high fidelity numerical techniques developed for accurately computing two-phase flows and free-surfaces particularly where precise volume conservation, interface advection and calculation of surface tension forces are fundamental. We hope that this work on the spherical Faraday instability can open up new possibilities in the study of pattern formation and that the numerical code can serve as a useful tool in the exploration of this rich dynamical system. In addition, we believe that the results obtained here serve both to validate the techniques implemented in BLUE and also to provide encouragement in applying the code to two-phase flow scenarios in highly non-linear regimes.

Acknowledgements

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 Equipement National de Calcul Intensif) through grants A0042B06721 and A0042A01119. L.S.T. acknowledges support from the Agence Nationale de la Recherche (ANR) through the TRANSFLOW project. This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (No. 2017R1D1A1B03028518).