3.1. Experimental results

The suspension is uniform at low oscillation frequencies. However, when $f$ is increased beyond $3.67\ {\rm Hz}$, high-density regions appear suddenly, followed by low-density tails. These density patches self-organize into a hexagonal pattern with a further increase of $f$, as shown in figure 1(b). Individual dense regions move in a circular path at the same frequency $f$ precisely, with little changes in location and size over time (figure 1c–f; refer to supplementary materials for multimedia view). The motion of the dense regions is defined by the position of the highest density locally. The typical diameter of the circular path, denoted as $d_w$, measures 13.3 mm, which is larger than the orbit of the container. In such a state, the surface becomes uneven too, introducing a peak-to-valley height difference $\Delta h\approx 1\ {\rm mm}$. Dense/loose regions correspond to peaks/valleys (figure 1g). As the suspension is bounded by the air–liquid interface, the observation indicates a non-uniform normal stress distribution accompanying the density inhomogeneity, a hallmark of non-Newtonian behaviour.

As reported by Oyarte Gálvez et al. (Reference Oyarte Gálvez, de Beer, van der Meer and Pons2017), the interaction force profile shows a hysteresis when corn starch particles are pressed into contact. However, the inhomogeneity transition here is reversible, which implies that the dense regions in our experiments are not permanent aggregates of contacting particles. To reveal the microscopic nature of the observed inhomogeneity, we measure the velocity of tracer particles near the suspension surface, $u_p$, and local packing fraction, $\phi$, simultaneously at a point on the path of a high-density region. Both $u_p$ and $\phi$ vary in phase with the same period as the oscillation, as shown in figure 2, i.e. particles in denser regions tend to move faster. However, the particle velocity remains lower than that of the density pattern, indicating that the motion of the observed density pattern is the propagation of density waves. Furthermore, the gradient of $\phi$ in the wavefront is significantly steeper than that at the rear (see figure 1g), a signature of shock waves.

Figure 2.

The local volume fraction $\phi$ (a) and particle velocity $u_p$ (b) varies with the phase angle of the oscillation. Both $\phi$ and $u_p$ rise when a density wave passes through the measuring area. Here $u_d={\rm \pi} d_w f$ is the propagation velocity of the density wave concerned. Date here is collected at $\varPhi =0.41$, $f=5.83\ {\rm Hz}$ and $h_{0}=3\ {\rm mm}$. The shallower suspension used here, thus more subtle surface deformation, improves the accuracy of the measurement of $\phi$ and $u_p$.

3.2. Instability onset

The non-uniform state cannot be modelled using a single-phase description of the suspension. In the exploration of flow instabilities in particulate suspensions, two-fluid models are normally used (Chacko et al. Reference Chacko, Mari, Cates and Fielding2018; Batchelor Reference Batchelor1988). The continuity equation of particles reads

(3.1)\begin{equation} \frac{\partial \phi}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot}({\boldsymbol{u}}_p\phi)=0. \end{equation}

When the inertia of particles is neglected, (3.1) can be written in the form of an advection–diffusion equation of $\phi$ (Anderson, Sundaresan & Jackson Reference Anderson, Sundaresan and Jackson1995). Shock waves are thus expected in certain circumstances. However, it has been demonstrated that the inertia of particles (the advective term) is critical for instability development (Batchelor Reference Batchelor1988; Johri & Glasser Reference Johri and Glasser2002). Therefore, we consider the momentum equation of particles in addition.

The general analytical formulation of the stress tensor of the particle phase remains a difficult task to tackle after decades of efforts (Jackson Reference Jackson2000; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). Here, we exclusively consider the dominant terms for relatively dense suspensions,

(3.2)\begin{equation} \phi \rho_{p} \left( \frac{\partial{\boldsymbol{u}_p}}{\partial t}+(\boldsymbol{u}_p \boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_p \right) =C_d({\boldsymbol{U}}-\boldsymbol{u}_p)+\eta_p\nabla^{2}\boldsymbol{u}_p-\boldsymbol{\nabla} \varPi. \end{equation}

On the right-hand side of (3.2), the first term represents the hydrodynamic drag exerted on particles, proportional to the relative velocity between the particle $\boldsymbol {u}_p$ and the average flow of the mixture ${\boldsymbol {U}}$, and the Richardson–Zaki approximation for $C_d$ (Richardson Reference Richardson1954; Buscall et al. Reference Buscall, Goodwin, Ottewill and Tadros1982) is used. Note that the mixture flow ${\boldsymbol {U}}$ is defined as the weighted average of the velocity of the two phases. The second term is the viscous force between particles, where $\eta _p$ is the dynamic viscous coefficient. The third term describes the gradient of particle pressure, $\varPi$. Both $\eta _p$ and $\varPi$ will be evaluated by their suspension counterparts, as the stresses of particle phase dominate the suspension dynamics at high $\phi$ (Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014). The stability analysis of (3.1)–(3.2) will be performed with respect to the uniform state. Therefore, the mixture velocity $\boldsymbol {U}$ of the suspension in the uniform state is computed first.

Consider a suspension with a volume fraction $\varPhi$ that behaves as a uniform fluid with a kinematic viscosity $\nu$ and is subjected to bottom oscillation. Since the orbital motion is the superposition of two perpendicular harmonic oscillations of angular frequency $\omega =2{\rm \pi} f$ with a phase difference of ${\rm \pi} /2$, the flow velocity $\boldsymbol {U}=(U_x,\,U_y)$ is written as a complex function $U(z,t) = U_x + i U_y$ which satisfies

(3.3a)$$\begin{gather} \frac{\partial U}{\partial t} = \nu \frac{\partial^{2}{U}}{\partial z^{2}}, \end{gather}$$

(3.3b)$$\begin{gather}\left.\frac{\partial{U}}{\partial z}\right\vert_{z=h} = 0 \quad\mathrm{and}\quad U(z=0,t) = A\omega {\rm e}^{{\rm i}\omega t}. \end{gather}$$

Equation (3.3b) indicates the no-slip condition on the bottom and the zero-shear condition on the free surface. For shear-thickening fluid, the effective kinematic viscosity of the suspension, $\nu = \nu _{0}(\phi _{J}-\phi )^{-2}$, is rate-dependent (Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018; Singh et al. Reference Singh, Mari, Denn and Morris2018). A recently developed phenomenological constitutive model (Wyart & Cates Reference Wyart and Cates2014), based on the mean-field description of shear-thickening process, suggests that $\phi _J$ experiences a crossover from its frictionless value $\phi _{0}$ to $\phi _{m}$ for frictional contacts, when the typical stress $\tau$ in the suspension exceeds a characteristic value $\tau ^*$, i.e. $\phi _{J}=\phi _{0}-{\rm e}^{{-\tau ^{*}}/{\tau }}(\phi _0-\phi _m)$. For the sample used here, we found $\nu _0=0.95 \times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$, $\tau ^{*}=3.7\ {\rm Pa}$, $\phi _m=0.45$ and $\phi _0=0.58$ according to rheology measurements (see Appendix C).

Analogous to the Stokes problem in two dimensions, the solution of (3.3) can be described by a dimensionless number $l=\sqrt {{2\nu }/{\omega }}/{h}$ (Yih Reference Yih1968), the square of which is merely the inverse Reynolds number representing the significance of viscosity relative to inertia. For $l\sim 1$, $U$ follows the orbital motion of the bottom plate with little phase lag and magnitude decay along with $z$. Due to the rate-dependence of $\nu$, the shear stress $\tau =\rho \nu \vert \partial U/\partial z\vert$ and $l$ are interrelated for a given $\phi =\varPhi$ (details are provided in Appendix B), where $\rho$ is the suspension density. To simplify the calculation, $\tilde {U} = \vert U(z=h,t)\vert$ and $\tilde {\tau } = \tau (z=0)$ are used to characterize the mainstream flow of the suspension in the uniform state. We assume that the initial development of the instability occurs in the mainstream direction, which results in the observed motion of density patterns. This assumption is significant, as (3.1)–(3.2) are hence reduced to one dimension (the main flow direction), and particle migration along the gradient and the vorticity directions are neglected. A linear stability analysis can be readily performed for the reduced equations. We leave the calculation details in Appendix B and present the main result here. The uniform flow ($\phi =\varPhi$ and ${u_p}=\tilde {U}$) becomes unstable against density perturbations of a wavenumber $k$, provided (Zhao & Pöschel Reference Zhao and Pöschel2021)

(3.4)\begin{equation} \mathcal{C}+k^2 \frac{\eta_p}{C_d}<\frac{\varPhi \tilde{U}^{\prime}}{\sqrt{\varPi^{\prime}/\rho_p}}, \end{equation}

where $\tilde {U}^{\prime }$ and $\varPi ^{\prime }$ are the derivatives of $\tilde {U}$ and $\varPi$ with respect to $\phi$ at $\phi =\varPhi$. In theory, $\mathcal {C}= 1$ is a constant whose value is to be adjusted by comparing with experiments, accounting for the neglected features in the model.

It is clear in (3.4) that the viscosity of the particle phase, $\eta _p$, stabilizes the short-wave disturbances. Therefore, the instability first occurs at the long wave limit ($k\sim 0$). Denoting the container diameter as $L={32}\ {\rm cm}$, the second term on the left-hand side, $\sim \eta _0/C_d L^2 \sim d^2/L^2 \sim O(10^{-5})$, is negligible. The onset of instability is thus insensitive to the value of $\eta _p$, which, however, shapes the developed density waves (see § 3.3). On the right-hand side of (3.4), $\tilde {U}^\prime$ introduces a shock-wave-like kinematic instability promoting the growth of high $\phi$ regions, i.e. a higher $\phi$ leads to a faster flow locally, as confirmed in figure 2. The particle pressure term $\varPi ^{\prime }$, on the other hand, provokes the particle migration out of the local high $\phi$ region. The competition between these two terms sets the onset of instability. Here, $\varPi = \tilde {\tau }(l,\phi )$ is used (Brown & Jaeger Reference Brown and Jaeger2012). With all terms defined, we calculate the onset frequency of the disturbance growth, $\omega _c$, and find that $\mathcal {C}=1.15$ aligns with the experimental observation (figure 3). In particular, the theoretical $\omega _c(\varPhi )$ agrees with its experimental counterparts quantitatively for the density-matched suspensions, where buoyancy is absent. Otherwise, the onset frequency $\omega _c$ is obscured by the suspending threshold of the relatively dense suspension without density-match in practice. The resultant higher value of $\mathcal {C}$ than the theory may be explained by the neglect of the particle migration in the vorticity direction and the implemented equality of $\varPi =\tilde {\tau }$. Both underestimate the denominator of the right-hand side of (3.4).

Figure 3.

State diagram. The onset frequency $\omega _c$ of density waves is denoted by open (solid) circles for the (density-matched) aqueous cornstarch suspension. The contour lines of $\varPhi \tilde {U}^{\prime }/{\sqrt {\varPi ^{\prime }/\rho _p}}$ (the right-hand side of (3.4)) are plotted for comparison. Inset: the dimensionless number $l$ varies with $\omega$ for $\varPhi =0.35$ (no DST) and $\varPhi =0.42$ (DST occurs), respectively.

The model predicts a minimum packing fraction around $\varPhi =0.392$ (figure 3), below which (3.4) cannot be satisfied, and the uniform state is always stable. This result is closely related to the shear-thickening nature of the suspension and reveals the underlying mechanism of the instability development. The solution for $\tilde {U}$ in (3.3) increases rapidly at intermediate $\phi$, thus corresponding to a regime of large $\tilde {U}^\prime$. An increase of $l$ would shift this regime towards lower densities. The pressure term $\varPi$ follows a similar trend. It can be shown that the variation of $\varPhi \tilde {U}^{\prime }/{\sqrt {\varPi ^{\prime }/\rho _p}}$ is dominated by $\sqrt { l/(\phi _J-\varPhi )}$ for $\omega \gtrsim 10\ {\rm rad}\ {\rm s}^{{-1}}$ (see Appendix B). The effective viscosity of the suspension $\nu$ in general grows with $\omega$, i.e. $\nu \sim \omega ^\alpha$ and $\alpha >0$. When approaching DST, $\alpha$ becomes considerably larger than 1, and $l\sim \omega ^{(\alpha -1)/2}$ displays a dramatic increase. Figure 3(inset) illustrates this distinct behaviour. A minimum $\varPhi$ is associated with the occurrence of DST (Fall et al. Reference Fall, Bertrand, Hautemayou, Mezière, Moucheront, Lemaître and Ovarlez2015; Hermes et al. Reference Hermes, Guy, Poon, Poy, Cates and Wyart2016), thus the same holds for the observed density waves here. The relation with shear-thickening is additionally confirmed via experiments of an aqueous solution of polydisperse silica beads (of an average diameter of $20\ \mathrm {\mu } {\rm m}$) in the same set-up. It was known that dissolving electrolytes within the solvent reduces the magnitude of repulsive forces between silica grains due to surface charge, and DST disappears accordingly (Clavaud et al. Reference Clavaud, Bérut, Metzger and Forterre2017). We confirm that the density waves in the suspension of silica beads under orbital oscillations are significantly weakened in the same manner (supplementary movie 3).

3.3. Development of density waves

Once the disturbance arises and grows in the flow direction, secondary instabilities may further develop (Anderson et al. Reference Anderson, Sundaresan and Jackson1995; Duru & Guazzelli Reference Duru and Guazzelli2002), and two-dimensional structures form (Glasser, Sundaresan & Kevrekidis Reference Glasser, Sundaresan and Kevrekidis1998). In our experiments, the density waves self-organize into a hexagonal pattern as $f$ increases (cf. figures 1 and 4 inset). A comprehensive theoretical analysis of the formation of the observed pattern is beyond the scope of this work. Instead, we argue that it can be understood from the symmetry perspective. Consider that the uniform state becomes unstable at one moment and breaks into alternating density bands perpendicular to the flow direction and separated by a characteristic wavelength $\lambda _c$. In the following moments, however, the flow direction changes constantly in the horizontal plane. Those bands cannot preserve the alignment with the flow without breaking the translational symmetry in the lateral (vorticity) direction. Therefore, the band structure is unstable and must be reduced to localized density patches. Particle migration dominates when those high-density regions are misaligned. The most stable structure thus maximizes the duration of the alignment with the flow. In other words, the favourite structure retains the highest rotational symmetry and a discrete translational symmetry of $\lambda _c$ in two dimensions, i.e. the hexagonal pattern, as observed. A subtle inference along this line of argument is that the maximum local density fluctuates six times during one oscillation period if the observer travels along with the density wave. Superharmonic oscillations of the local maxima of $\phi$ are indeed observed in experiments. As shown in figure 4, the distance between neighbouring density waves, $\lambda _c$, is larger than the diameter of their circular motion, i.e. $d_w/\lambda _c\approx 0.5$, ensuring no intersection between their trajectories. Therefore, such superharmonic fluctuations suggest the out-of-alignment and realignment events between the flow and the hexagonal pattern.

Figure 4.

The length scale of the observed density pattern $\lambda _c$ varies with $\omega$ in the aqueous-starch suspension (open circles) and the density-matched suspension (solid circles). Experimental parameters here are $\varPhi =0.42$, $h=5\ {\rm mm}$. The dashed line is the wave number calculated by the linearized model with a proliferated viscosity, $\lambda _m$. See the main text for discussions. Inset: two-dimensional Fourier spectrum of figure 1(b).

The wavelength of the observed pattern, $\lambda _c$, is measured from the two-dimensional Fourier spectrum of the experimental images and averaged over one oscillation cycle, typified by figure 4(inset). Data from the density-matched solution and the aqueous suspension are plotted in figure 4. Overall, $\lambda _c$ decreases with $\omega$, demonstrating a shift of dominance from long waves to short waves. Just above $\omega _c$, the relatively large error bars indicate that the structure is not well-ordered yet. For the intermediate $\omega >\omega _c$, $\lambda _c$ of the density-matched suspension is larger than the pure aqueous solvent, which indicates that the vertical migration of particles could modify the inhomogeneous pattern. With increasing $\omega$, $\lambda _c$ decreases towards an asymptotic value where the aqueous suspension and the density-matched suspension match. The measured $\lambda _c$ is compared with linearly selected wavenumber $\lambda _m$, corresponding to the peak growth rate Fourier mode in the linearized model. The magnitude of $\eta _p$, bounded by $\rho \nu (\phi )=\rho \nu _{0}(\phi _J-\phi )^{-2}$, is critical for evaluating $\lambda _m$. From the perspective of linear analysis, $\eta _p\sim \rho \nu (\varPhi )$ for the uniform state is implemented, and the resultant $\lambda _m$ is of the order of micrometres. Nevertheless, as seen in figure 1(g), the local density, $\phi$, can be as high as 0.449 for the developed density waves, leading to viscosity proliferation locally. The corresponding $\eta _p(\phi )$ would shift $\lambda _m$ towards longer wavelength. The evaluation of $\lambda _m$ with $\eta _p = 2800\ {\rm Pa}\ {\rm s}$ is in reasonable agreement with $\lambda _c$ (figure 4). It again indicates that the description of the fully developed density pattern is beyond purely linear predictions (Anderson et al. Reference Anderson, Sundaresan and Jackson1995; Duru et al. Reference Duru, Nicolas, Hinch and Guazzelli2002).