2.1 Pattern formation on a sphere

Our previous paper (Ebo-Adou & Tuckerman Reference Ebo-Adou and Tuckerman2016) 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 \leqslant m\leqslant \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 Reference Busse1975; Busse & Riahi Reference Busse and Riahi1982; Ihrig & Golubitsky Reference Ihrig and Golubitsky1984; Riahi Reference Riahi1984; Golubitsky, Stewart & Schaeffer Reference Golubitsky, Stewart and Schaeffer1988; Chossat, Lauterbach & Melbourne Reference Chossat, Lauterbach and Melbourne1991; Matthews Reference Matthews2003) 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$ . Here $O(2)$ consists of the symmetries of a circle and patterns which are axisymmetric (for some orientations) 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\unicode[STIX]{x03C0}/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 \geqslant 2$ (Chossat et al. Reference Chossat, Lauterbach and Melbourne1991). 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 Reference Matthews2003) or which extremizes a functional (Busse Reference Busse1975). 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 Reference Busse1975; Busse & Riahi Reference Busse and Riahi1982; Riahi Reference Riahi1984) 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.

Second, the theory applies close to the threshold. Our simulations are carried out far from threshold, so that instabilities can grow on a reasonable time scale 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:

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}P+\unicode[STIX]{x1D70C}\boldsymbol{G}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})+\boldsymbol{F},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

where $\boldsymbol{u}$ is the velocity, $P$ is the pressure, $\unicode[STIX]{x1D70C}$ is the density, $\unicode[STIX]{x1D707}$ is the dynamic viscosity and $\boldsymbol{F}$ is the local surface tension force at the interface. Here, $\boldsymbol{G}$ is an imposed time-dependent radial acceleration:

(2.2) $$\begin{eqnarray}\boldsymbol{G}=-(g+a\cos (\unicode[STIX]{x1D714}t))\frac{r}{R}\boldsymbol{e}_{r},\end{eqnarray}$$

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

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

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}(\boldsymbol{x},t)=\unicode[STIX]{x1D70C}_{1}+(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})I(\boldsymbol{x},t)\\ \unicode[STIX]{x1D707}(\boldsymbol{x},t)=\unicode[STIX]{x1D707}_{1}+(\unicode[STIX]{x1D707}_{2}-\unicode[STIX]{x1D707}_{1})I(\boldsymbol{x},t).\end{array}\right\}\end{eqnarray}$$

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

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 the Level-Set distance function field on the Eulerian grid as well as by triangles on the Lagrangian interface mesh.

The surface tension $\boldsymbol{F}$ is implemented by the hybrid/compact formulation (Shin Reference Shin2007)

(2.4a,b ) $$\begin{eqnarray}\boldsymbol{F}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}_{H}\unicode[STIX]{x1D735}I,\quad \unicode[STIX]{x1D705}_{H}=\frac{\boldsymbol{F}_{L}\boldsymbol{\cdot }\boldsymbol{N}}{\boldsymbol{N}\boldsymbol{\cdot }\boldsymbol{N}},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}$ is the surface tension coefficient and $\unicode[STIX]{x1D705}_{H}$ is twice the mean interface curvature field calculated on the Eulerian grid, with

(2.5a,b ) $$\begin{eqnarray}\boldsymbol{F}_{_{L}}=\int _{\unicode[STIX]{x1D6E4}(t)}\unicode[STIX]{x1D705}_{f}\boldsymbol{n}_{f}\unicode[STIX]{x1D6FF}_{f}(\boldsymbol{x}-\boldsymbol{x}_{f})\,\text{d}a,\quad \boldsymbol{N}=\int _{\unicode[STIX]{x1D6E4}(t)}\boldsymbol{n}_{f}\unicode[STIX]{x1D6FF}_{f}(\boldsymbol{x}-\boldsymbol{x}_{f})\,\text{d}a.\end{eqnarray}$$

Here, $\boldsymbol{x}_{f}$ is a parameterization of the time-dependent interface, $\unicode[STIX]{x1D6E4}(t)$ , and $\unicode[STIX]{x1D6FF}_{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}a$ is the area of an interface element; $\unicode[STIX]{x1D705}_{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}a$ , 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_{num}$ can be found in Shin & Juric (Reference Shin and Juric2007, Reference Shin and Juric2009), 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 $\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 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 (Reference Kwak and Lee2004). 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 centres 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, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011). See Shin et al. (Reference Shin, Chergui and Juric2017) 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 $\unicode[STIX]{x1D70C}$ , viscosity $\unicode[STIX]{x1D708}$ , and surface tension $\unicode[STIX]{x1D70E}$ , leading potentially to a large number of non-dimensional parameters, in addition to those describing the forcing amplitude $a$ and frequency $\unicode[STIX]{x1D714}$ . 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 $\unicode[STIX]{x1D70C}_{d}=940~\text{kg}~\text{m}^{-3}$ and $\unicode[STIX]{x1D70C}_{out}=1.205~\text{kg}~\text{m}^{-3}$ , leading to the density ratio $\unicode[STIX]{x1D70C}_{d}/\unicode[STIX]{x1D70C}_{out}=780$ , while the inner and outer kinematic viscosities are $\unicode[STIX]{x1D708}_{d}=10^{-5}~\text{m}^{2}~\text{s}^{-1}$ and $\unicode[STIX]{x1D708}_{out}=1.5\times 10^{-5}~\text{m}^{2}~\text{s}^{-1}$ , leading to the viscosity ratio $\unicode[STIX]{x1D708}_{d}/\unicode[STIX]{x1D708}_{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 (Reference Ebo-Adou and Tuckerman2016) is that the Mathieu equation describes the inviscid Faraday instability even in a spherical geometry.) More specifically, $\unicode[STIX]{x1D708}_{d}/(R^{2}\unicode[STIX]{x1D714})$ 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 Reference Kumar1996), whose period is the same as the forcing period.

For the capillary cases, we set the surface tension to be $\unicode[STIX]{x1D70E}=0.02~\text{kg}~\text{s}^{-2}$ , 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~\text{m}~\text{s}^{-2}$ , the surface tension to be $\unicode[STIX]{x1D70E}=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 (2.15), and we list these in § 2.7 and table 3.

2.4 Numerical parameters

For the capillary cases, the domain is a sphere whose radius is $2R$ (see § 2.2 and figures 8 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 $\unicode[STIX]{x0394}x=4R/N_{x}$ for the capillary cases should be compared to the circumference $2\unicode[STIX]{x03C0}R$ of the drop and to the approximate wavelengths $2\unicode[STIX]{x03C0}R/\ell$ . The radius and circumference are spanned by $R/\unicode[STIX]{x0394}x=N_{x}/4$ and $2\unicode[STIX]{x03C0}R/\unicode[STIX]{x0394}x=N_{x}\unicode[STIX]{x03C0}/2$ , which are 32 and 201, respectively, for the $128^{3}$ grid. The number of gridpoints per wavelength $N_{x}\unicode[STIX]{x03C0}/(2\ell )$ goes from 201 for $\ell =1$ to 33 points for $\ell =6$ for this grid. The interface is represented by approximately $4\unicode[STIX]{x03C0}R^{2}/(\unicode[STIX]{x0394}x)^{2}=N_{x}^{2}\unicode[STIX]{x03C0}/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 (for example,  $R/\unicode[STIX]{x0394}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 visualization 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 $\unicode[STIX]{x0394}t$ is chosen at each iteration in order to satisfy a criterion based on

(2.6) $$\begin{eqnarray}\{\unicode[STIX]{x0394}t_{CFL},\unicode[STIX]{x0394}t_{int},\unicode[STIX]{x0394}t_{vis},\unicode[STIX]{x0394}t_{cap}\},\end{eqnarray}$$

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

(2.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x0394}t_{CFL}\equiv \min _{j}\left(\min _{domain}\left(\frac{\unicode[STIX]{x0394}x_{j}}{u_{j}}\right)\right),\quad \unicode[STIX]{x0394}t_{int}\equiv \min _{j}\left(\min _{\unicode[STIX]{x1D6E4}(t)}\left(\frac{\unicode[STIX]{x0394}x_{j}}{\Vert \boldsymbol{V}\Vert }\right)\right)\\ \displaystyle \unicode[STIX]{x0394}t_{vis}\equiv \min \left(\frac{\unicode[STIX]{x1D70C}_{2}}{\unicode[STIX]{x1D707}_{2}},\frac{\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x1D707}_{1}}\right)\frac{\unicode[STIX]{x0394}x_{min}^{2}}{6},\quad \unicode[STIX]{x0394}t_{cap}\equiv \frac{1}{2}\left(\frac{(\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2})\unicode[STIX]{x0394}x_{min}^{3}}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}}\right)^{1/2},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x0394}x_{min}=\min _{j}(\unicode[STIX]{x0394}x_{j})$ . In our simulations, this minimum is realized by $\unicode[STIX]{x0394}t_{CFL}\approx \unicode[STIX]{x0394}t_{int}$ for our gravitational waves and by $\unicode[STIX]{x0394}t_{cap}\approx \unicode[STIX]{x0394}t_{vis}$ for our capillary waves.

Figure 1.

Time evolution of the normalized deviation of $\unicode[STIX]{x0394}V(t)/V$ for a simulation with $\ell =6$ . The value of $|\unicode[STIX]{x0394}V|/V$ remains less than approximately $3\times 10^{-6}$ throughout the simulation.

2.5 Spherical harmonic transform

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

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=\mathop{\sum }_{\ell =0}^{\infty }\mathop{\sum }_{m=-\ell }^{m=\ell }\unicode[STIX]{x1D701}_{\ell }^{m}(t)Y_{\ell }^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}),\end{eqnarray}$$

where $\unicode[STIX]{x1D701}$ is the distance from the domain centre. We compute both $|\unicode[STIX]{x1D701}_{\ell }^{m}(t)|$ and

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D701}_{\ell }(t)\equiv \left[\mathop{\sum }_{m=0}^{\ell }|\unicode[STIX]{x1D701}_{\ell }^{m}(t)|^{2}\right]^{1/2}.\end{eqnarray}$$

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

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D701}_{0}(t)=\unicode[STIX]{x1D701}_{0}^{0}=\int _{0}^{2\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D719}\int _{0}^{\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)Y_{0}^{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\end{eqnarray}$$

i.e.  $\unicode[STIX]{x1D701}_{0}(t)$ is $\sqrt{4\unicode[STIX]{x03C0}}$ 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 $(x_{i},y_{i},z_{i})$ on the fixed Cartesian grid. We transform each of these points to spherical coordinates $(\unicode[STIX]{x1D701}_{i},\unicode[STIX]{x1D703}_{i},\unicode[STIX]{x1D719}_{i})$ and interpolate the function $\unicode[STIX]{x1D701}(\unicode[STIX]{x1D703}_{i},\unicode[STIX]{x1D719}_{i})=\unicode[STIX]{x1D701}_{i}$ onto a regular spherical grid $(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D719}_{k})$ , discretizing the latitudinal interval $[0,\unicode[STIX]{x03C0}]$ and longitudinal interval $[0,2\unicode[STIX]{x03C0}]$ with $N_{\unicode[STIX]{x1D703}}=80$ and $N_{\unicode[STIX]{x1D719}}=80$ points. Using the discrete data $\unicode[STIX]{x1D701}(\unicode[STIX]{x1D703}_{j},\unicode[STIX]{x1D719}_{k})$ and truncating the series (2.8) at four times the dominant value of $\ell$ , we compute the spherical harmonic coefficients $\unicode[STIX]{x1D701}_{\ell }^{m}$ by the method of least squares (Politis Reference Politis2013, Reference Politis2016). 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 $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ of one of its features, such as a maximum, rotating in longitude $\unicode[STIX]{x1D719}$ to place it in the $(x,z)$ plane and then rotating in colatitude $\unicode[STIX]{x1D703}$ to place this feature at the north pole. This allows us to interpret our spectra by using the explicit representations given by Busse (Reference Busse1975) 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.

Figure 2.

Comparison between the theoretical and computed threshold $a_{c}$ Curves show theoretical thresholds calculated from Floquet theory (Ebo-Adou & Tuckerman Reference Ebo-Adou and Tuckerman2016) for $\ell =4$ , 5, 6, while circles show results from numerical simulations with resolution $256^{3}$ at the indicated frequencies calculated by the method shown in figure 3.

Figure 3.

Calculation of thresholds. (a) Time series of interface height oscillate and grow or decay exponentially, depending on the imposed acceleration amplitude, as shown by the linear time dependence of the logarithms of the peaks. (b) The slopes of the logarithms of the time series maxima constitute the growth or decay rates, which are interpolated to determine the acceleration amplitude threshold for $\ell =6$ , shown as the rightmost circle in figure 2.

Table 1.

Comparison between theoretical and computed thresholds $a_{c}$ for various $\ell$ values and numerical resolutions.

2.6 Validation

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

In our previous investigation (Ebo-Adou & Tuckerman Reference Ebo-Adou and Tuckerman2016), we extended the method of Kumar & Tuckerman (Reference Kumar and Tuckerman1994) for computing the Faraday threshold via Floquet theory to a spherical geometry. In figure 2 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 § 2.3, except that we increase the viscosity to $\unicode[STIX]{x1D708}_{d}=10^{-4}~\text{m}^{2}~\text{s}^{-1}$ . For $\ell =6$ , the error in the threshold is approximately 1 % for a resolution of $128^{3}$ and only approximately 0.1 % for a resolution of $256^{3}$ .

Figure 4.

Comparison of evolution of spherical harmonic spectra $Y_{6}^{m}$ for resolutions (a) $128^{3}$ and (b) $256^{3}$ .

Figure 4 compares the time evolution of the envelope of the spherical harmonic coefficients $|\unicode[STIX]{x1D701}_{\ell }^{m}|$ for capillary wave simulations using resolutions $128^{3}$ and $256^{3}$ . The fluid parameters are again those of § 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.

Figure 5.

Free decay of a shape-deformed drop to the spherical shape. The evolution of $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}\equiv \unicode[STIX]{x1D701}_{max}-\unicode[STIX]{x1D701}_{min}$ with time is shown for $\ell =5$ for a resolution of $256^{3}$ . The blue curve corresponds to the output from numerical simulations which is fitted by a nonlinear least-square method (dashed red). The green curve corresponds to the small-viscosity-limit prediction for viscosity (a) $\unicode[STIX]{x1D708}_{d}=10^{-5}~\text{m}^{2}~\text{s}^{-1}$ and (b) $\unicode[STIX]{x1D708}_{d}=10^{-4}~\text{m}^{2}~\text{s}^{-1}$ .

Table 2.

Comparison between decay rates and nonlinear frequencies from simulations and from the low-viscosity theory of (2.12)–(2.13).

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 (Reference Rayleigh1879) showed that a drop initialized with a perturbation of spherical wavenumber $\ell$ oscillates with frequency $\unicode[STIX]{x1D714}$ , where

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{\ell (\ell -1)(\ell +2)}{R^{3}}.\end{eqnarray}$$

For small kinematic viscosity $\unicode[STIX]{x1D708}_{d}$ , Lamb (Reference Lamb1932) showed that the oscillation amplitude decays to the spherical rest state according to

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D701}_{\ell }(t)\propto \text{e}^{-\unicode[STIX]{x1D706}t}|\text{cos}(\unicode[STIX]{x1D714}^{\prime }t)|,\end{eqnarray}$$

where the decay rate and the oscillation frequency are given by

(2.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{\unicode[STIX]{x1D708}_{d}(\ell -1)(2\ell +1)}{R^{2}},\quad \unicode[STIX]{x1D714}^{\prime }=(\unicode[STIX]{x1D714}^{2}-\unicode[STIX]{x1D706}^{2})^{1/2}.\end{eqnarray}$$

Figure 5 shows the evolution of $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}\equiv \unicode[STIX]{x1D701}_{max}-\unicode[STIX]{x1D701}_{min}$ obtained via numerical simulations with a resolution of $256^{3}$ of capillary waves with the fluid parameters of § 2.3 initialized with axisymmetric perturbations $Y_{5}^{0}$ for two different viscosities, along with curves obtained from the low-viscosity theory of (2.12)–(2.13). 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 ( $\unicode[STIX]{x1D708}_{d}=10^{-5}~\text{m}~\text{s}^{-2}$ ) so that the low-viscosity theory is more applicable. This best-case deviation is approximately 7 % for $\unicode[STIX]{x1D706}$ and 2 % for $\unicode[STIX]{x1D714}^{\prime }$ for $\ell =5$ , and only approximately 3 % for $\unicode[STIX]{x1D706}$ and 1 % for $\unicode[STIX]{x1D714}^{\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:

(2.14a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Axisymmetric all }\ell \quad \unicode[STIX]{x1D701}-R\propto Y_{\ell }^{0} & \displaystyle\end{eqnarray}$$

(2.14b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Tetrahedral }\ell =3\quad \unicode[STIX]{x1D701}-R\propto Y_{3}^{2}+\text{c.c.} & \displaystyle\end{eqnarray}$$

(2.14c ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Cubic }\ell =4\quad \unicode[STIX]{x1D701}-R\propto \sqrt{7}Y_{4}^{0}+\sqrt{5}Y_{4}^{4}+\text{c.c.} & \displaystyle\end{eqnarray}$$

(2.14d ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{4}~\ell =5\quad \unicode[STIX]{x1D701}-R\propto \sqrt{3}Y_{5}^{0}+\sqrt{5}Y_{5}^{4}+\text{c.c.} & \displaystyle\end{eqnarray}$$

(2.14e ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Icosahedral }\ell =6\quad \unicode[STIX]{x1D701}-R\propto \sqrt{11}Y_{6}^{0}+\sqrt{14}Y_{6}^{5}+\text{c.c.} & \displaystyle\end{eqnarray}$$

These formulae for functions with a given symmetry are given by Busse (Reference Busse1975) and Riahi (Reference Riahi1984) 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 visualizations of the drop.

Table 3.

Parameters of imposed forcing and of the observed response.

Figure 6.

Instability tongues resulting from oscillatory forcing with amplitude $a$ and angular frequency $\unicode[STIX]{x1D714}$ for an inviscid drop with the parameters given in table 3. Diamond shapes designate the parameter values of our simulations. Solid curves bound subharmonic tongues and dashed curves bound harmonic tongues. (a) Tongues corresponding to gravitational instability with spherical wavenumbers $\ell =1$ , 2, 3, 4 originate at $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{g}=2\sqrt{\ell }/n$ , with $n=1$ for subharmonic tongues and $n=2$ for harmonic tongues. (b) Tongues corresponding to capillary instability with spherical wavenumbers $\ell =2$ , 3, 4, 5, 6 originate at $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}=2\sqrt{\ell (\ell -1)(\ell +2)}/n$ .

The physical parameters are as given in § 2.3. We study either gravity or capillary waves, i.e. a Bond number $\unicode[STIX]{x1D70C}gR^{2}/\unicode[STIX]{x1D70E}$ 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

(2.15a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{g}^{2}\equiv \frac{g}{R},\quad \unicode[STIX]{x1D714}_{c}^{2}\equiv \frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}R^{3}}, & \displaystyle\end{eqnarray}$$

(2.15c,d ) $$\begin{eqnarray}\displaystyle & \displaystyle a_{g}\equiv R\unicode[STIX]{x1D714}_{g}^{2}=g,\quad a_{c}\equiv R\unicode[STIX]{x1D714}_{c}^{2}=\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}R^{2}}. & \displaystyle\end{eqnarray}$$

See Ebo-Adou & Tuckerman (Reference Ebo-Adou and Tuckerman2016) 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.