Emergence of a hexagonal pattern in shear-thickening suspensions under orbital oscillations

TL;DR

This paper investigates how orbital oscillations induce a hexagonal pattern in shear-thickening suspensions, revealing the intrinsic heterogeneity of the state and the influence of boundary conditions on pattern formation.

Contribution

It demonstrates the emergence of hexagonal density patterns driven by orbital oscillations and provides a linear stability analysis explaining the conditions for pattern formation.

Findings

Localized density waves appear beyond a frequency threshold.

Hexagonal patterns form due to coupling between advection and shear-thickening.

Rigid confinement suppresses density wave formation.

Abstract

Dense particle suspension under shear may lose its uniform state to large local density and stress fluctuations, which challenge the mean-field description of the system. Here, we explore the novel dynamics of a non-Brownian suspension under orbital oscillations, where localized density waves along the flow direction appear beyond an excitation frequency threshold and self-organize into a hexagonal pattern across the system. The spontaneous occurrence of the inhomogeneity pattern arises from a coupling between particle advection and the shear-thickening nature of the suspension. Through linear stability analysis, we show that they overcome the stabilizing effects of particle pressure at sufficient particle volume fraction and oscillation frequency. In addition, the long-standing density waves degenerate into random fluctuations when replacing the free surface with rigid confinement. It indicates that the shear-thickened state is intrinsically heterogeneous, and the boundary conditions are crucial for developing local disturbance.

Equations

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot ({\vec{u}}_p\varphi )=0.$$

$$\varphi {\rho }_p\left(\frac{\partial {\vec{u}}_p}{\partial t}+({\vec{u}}_p\cdot \nabla ){\vec{u}}_p\right)=C_d(\vec{U}-{\vec{u}}_p)+{\eta }_p{\nabla }^2{\vec{u}}_p-\nabla \Pi .$$

$$\frac{\partial U}{\partial t}$$

$$\displaystyle\frac{\partial{U}}{\partial z}\Big{|}_{z=h}=0$$

$$\mathcal{C}+k^2\frac{{\eta }_p}{C_d}<\frac{\varphi {\tilde{U}}^{\prime }}{\sqrt{{\Pi }^{\prime }/{\rho }_p}},$$

Full text

Dense suspensions composed of mixtures of particles and fluid are ubiquitous in natural phenomena and in industrial processes [1, 2]. For a sufficiently large volume fraction of particles, $\Phi$, the suspensions can exhibit a wide range of nonlinear phenomena, including yielding, shear thinning or thickening, and shear jamming  [3, 4, 5], which has stimulated decades of research in the rheology and physics community  [2, 6]. In particular, discontinuous shear thickening (DST), the suspension viscosity, $\eta$, changes abruptly over several orders of magnitude, has recently been understood as a transition from lubricated to frictional particle interactions, when the applied stress overcomes the interparticle repulsion [7, 8, 9]. The constitutive curve of such a suspension, relating the shear stress $\tau$ and the shear rate $\dot{\gamma}$, displays a $\mathsf{S}$-shape in which a negatively sloped region connects the lubricated branch of low-viscosity at low stresses and frictional branch of high-viscosity at high stresses [9, 10, 11]. Because of the characteristic $\mathrm{d}\dot{\gamma}/\mathrm{d}\tau\leq 0$, the steady uniform state of a suspension may become unstable, and large spatiotemporal fluctuations of stress and local densities develop into structures/patterns as consequences [12]. Typical one-dimensional examples are shear-induced dynamic gradient banding [13] and vorticity banding [14, 15]. Periodic density and stress fluctuations are reported in both the Couette and the parallel plate configurations [16, 17]. Yet, little is known about the dynamics of shear thickening suspensions in flow configurations beyond simple shear [18, 19]. Moreover, although boundary confinements are proposed to be essential for DST of non-Brownian suspensions [9], its influences on the instability development in such systems, if any, has not been explored.

In this Letter, we investigate the instability arising in the flow direction in dense Non-Brownian suspensions under swirling excitation. A novel unsteady dynamics where density waves self-organized into a hexagonal pattern is observed, accompanied with large stress fluctuations ( e.g., a cornstarch aqueous suspension under excitation can hold up steel balls floating on its surface constantly without shear-jamming globally). Our result shows that it is closely related to the discontinuous shear thickening transition of the suspension. Moreover, the essential role of boundary conditions on the development of disturbances is discussed.

The suspension consists of deionized water mixed with cornstarch particles, of an average size $d=15\$\mu\mathrm{m}, confined within an open cylindrical container. The diameter and height of the container is $32\ $\mathrm{c}\mathrm{m} and $5\$\mathrm{c}\mathrm{m} respectively. The system is subjected to a horizontal circular vibration, i.e., a circular movement of the entire platform. The oscillation frequency can be tuned from $0.5$ to $8.33~{}$\mathrm{H}\mathrm{z} with a fixed amplitude $A=5~{}$\mathrm{m}\mathrm{m} indicating the radius of the circular motion. The suspension is illuminated by a LED panel from below. A high speed camera (Microtron EoSens 1.1cxp2) is fixed in the laboratory frame of reference. The local volume fraction, $\phi$, is calculated according to the transmission ratio of the back light [[20](#bib.bib20)]. Surface topography of suspension is measured on demand via the deflection of laser sheets. The system studied here is defined by three major experimental parameters: the oscillation frequency $f$, the global particle packing fraction $\Phi$, and the thickness of the dense suspension system $h$. We first demonstrate the primary observations in a reference system specified at $\Phi=0.42$ and $h=5$\mathrm{m}\mathrm{m}. The influence of various parameters is discussed afterwards.

The suspension is uniform at low oscillation frequencies, as shown in Fig. 1(a). However, when $f$ is increased beyond $3.67\$\mathrm{H}\mathrm{z}, high density regions appear suddenly, followed by low density tails. Individual dense area move in a circular path, with little drifts, at the frequency $f$ precisely (see [[20](#bib.bib20)] for multimedia view). These density patches self-organize into a hexagonal pattern with further increase of $f$, as shown in Fig.[1](#S0.F1)(b). In such a state, the surface becomes uneven too, introducing a peak-to-valley height difference $\Delta h\approx 1$\mathrm{m}\mathrm{m}. We measure the velocity of tracer particles near the suspension surface, $u_{p}$, and local packing fraction, $\phi$ 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 Fig.2, i.e., particles in denser regions tend to move faster. Nevertheless, the particle velocity is smaller than that of the density pattern. This indicates 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 Fig.1(c)], a signature of supersonic waves.

Although the suspension can be viewed as a single phase in the uniform state, the particle phase has to be treated separately to understand the spontaneous occurrence of the observed inhomogeneity. The continuity equation of particles reads

[TABLE]

When the inertia of particles is neglected, Eq. 1 can be written in the form of advection-diffusion equation of $\phi$ [21]. Shock waves are thus expected in certain circumstances. However, it has been shown that the inertia of particles are critical for instability development [22, 23]. Therefore, we consider the momentum equation of particles in addition.

The general analytical form of the stress tensor of the particle phase is still a difficult task to tackle after decades of efforts [24, 6]. Here, we only employ dominating terms for the relatively dense suspensions:

[TABLE]

On the right hand side of Eq. 2, the first term represents the hydrodynamic drag on particles, proportional to the relative velocity between the particle $\vec{u}_{p}$ and the average flow of the mixture $\vec{U}$, and the Richardson-Zaki approximation for $C_{d}$ [25, 26] is used [20]. 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, $\Pi$. Both $\eta_{p}$ and $\Pi$ will be evaluated with their suspension counterparts, as the stresses of particle phase dominates the suspension dynamics for large $\phi$ [27]. The instability analysis of Eqs. 1-2 will be performed with respect to the uniform state. Therefore, the average velocity $\vec{U}$ of the suspension in the uniform state is computed first.

The motion of the uniform state of the suspension is a solution to Stokes problem in two-dimension, i.e., the motion is induced by the oscillating bottom plate 111We ignore the side wall boundary, as changing the diameter of the container from 32 $\mathrm{c}\mathrm{m}$ to 16 $\mathrm{c}\mathrm{m}$ does not alter the phenomena.. Since the swirling excitation is the superposition of two perpendicular harmonic oscillations with a phase difference of $\pi/2$, the flow velocity $\vec{U}=(U_{x},\,U_{y})$ can be transformed into a complex function $U(z,t)=U_{x}+iU_{y}$ which satisfies

[TABLE]

The effective kinematic viscosity of the suspension $\nu=\nu_{0}(\phi_{J}-\phi)^{-2}$ diverges as $\phi$ approaching the jamming packing density $\phi_{J}$ [6, 29, 30]. A phenomenological constitutive model developed recently [31] 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}-e^{{-\tau^{*}}/{\tau}}(\phi_{0}-\phi_{m})$. For the sample studied here, we found $\nu_{0}=0.95\times 10^{-6}~{}$\mathrm{m}^{2}\mathrm{/}\mathrm{s}, $\tau^{*}=3.7~{}$\mathrm{P}\mathrm{a}, $\phi_{m}=0.45$ and $\phi_{0}=0.58$ according to rheology measurements [20].

The solution of Eq. 3 can be described by a dimensionless number $l=\sqrt{{2\nu}/{\omega}}/{h}$ [32], 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 swirling motion of the bottom plate with little phase lag and magnitude decay with height $z$. However, due to the rate-dependence of $\nu$, the shear stress, $\tau=\rho\nu|\partial U/\partial z|$, and $l$ are mutually determined for a given $\phi$, where $\rho$ is the suspension density. To simplify the calculation, $\tilde{U}(l,\phi)=|U(z=h,t)|$ and $\tilde{\tau}(l,\phi)=\tau(z=0)$ are thus used to characterize the mainstream flow of the suspension. We assume that the initial development of the instability occurs in the mainstream direction, which results into the observed motion of density patterns. This assumption is significant, as Eqs 1-2 are hence reduced to one-dimension and a linear stability analysis can be readily performed. The uniform flow ($\phi=\Phi$ and ${u_{p}}=\tilde{U}$) becomes unstable against density perturbations of a wavenumber $k$, provided [20, 33]

[TABLE]

where $\tilde{U}^{\prime}$ and $\Pi^{\prime}$ are the derivatives of $\tilde{U}$ and $\Pi$ with respect to $\phi$ at $\phi=\Phi$. $\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 Eq. 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$). The onset of instability is thus insensitive to the value of $\eta_{p}$, which, however, shapes the developed density waves (see discussion on Fig. 4 below). On the right hand side of Eq. 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 faster flow locally, as confirmed in Fig. 2. The particle pressure term $\Pi^{\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 $\Pi=\tilde{\tau}(l,\phi)$ is used [34]. With all terms defined, we calculate the onset frequency of the disturbance growth, $\omega_{c}$, and find that $\mathcal{C}=1.3$ matches the experimental observation [Fig. 3]. In particular, the theoretical $\omega_{c}(\Phi)$ agrees with the its experimental counterparts quantitatively for the density-matched suspensions, where buoyancy is absent. Otherwise, the theoretical $\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 diffusion of particles in the vorticity direction and the implemented equality of $\Pi=\tilde{\tau}$. Both underestimate the denominator of the right hand side of Eq. 4.

The model predicts a minimum packing fraction around $\Phi=0.389$ [Fig. 3], below which the uniform state is always stable. This result is closely related to the DST of the suspension and reveals the underlying nature of the instability development. The solution for $\tilde{U}$ in Eq. 3 increases rapidly at intermediate $\phi$, thus corresponding to a regime of large $\tilde{U}^{\prime}$. An increase of $l$ would shift this regime toward lower densities. The pressure term $\Pi$ follows a similar trend. It can be shown that the variation of $\phi\tilde{U}^{\prime}/{\sqrt{\Pi^{\prime}/\rho_{p}}}$, the right hand side of Eq. 4, is dominated by $\sqrt{l/(\phi_{J}-\phi)}$ for $\omega\gtrsim 10~{}$\mathrm{r}\mathrm{a}\mathrm{d}~{}\mathrm{s}^{-1} [[20](#bib.bib20)]. The effective viscosity of the suspension $\nu$ in general grows with $\omega$, i.e., $\nu\sim\omega^{\alpha}$ and $\alpha>0$. However, without DST $\alpha<1$ and $l\sim\omega^{(\alpha-1)/2}$ decreases with $\omega$. In contrast, when DST is triggered, $\alpha>1$ and $l$ rises, leading to the proliferation of $\tilde{U}^{\prime}$ at the given $\Phi$ and the inequality in Eq. [4](#S0.E4). This distinct behavior is illustrated in Fig. [3](#S0.F3) inset. A minimum $\Phi$ is expected for the occurrence of DST  [[17](#bib.bib17), [35](#bib.bib35)], thus so is it for the observed density waves here. The relation with DST is additionally confirmed with experiments of a aqueous solution of polydispersion silica beads (of an average diameter of 20 $\mu$\mathrm{m}) in the same setup. It was reported that dissolving electrolytes within the solvent reduces the magnitude of repulsive forces between silica grains due to surface charge, and DST disappears accordingly [36]. The density waves in the suspension of silica beads here are significantly weakened in the same manner [20].

Once the supersonic disturbance arises and grows in the flow direction, secondary instabilities may further develop [21, 37], and two-dimensional structures form [38]. In our experiments, the density waves self-organize into a hexagonal pattern as $f$ increases (cf. Fig. 1 and 4 Inset). A thoroughgoing analysis of the formation of the observed pattern is beyond the scope of the linear theory above. Instead, we argue that it can be understood from the symmetry perspective. Consider that at a moment the uniform state becomes unstable and break up into alternating density bands perpendicular to the flow direction that are separated by a characteristic wavelength $\lambda_{c}$. In the following moments, however, the flow direction changes constantly in the horizontal plane, and those bands can not preserve the alignment with the flow without breaking the translational symmetry in the lateral (vorticity) direction. Therefore, the band structure itself is unstable and must be reduced to localized density patches. Diffusion dominates, when those high density regions are misaligned. The most stable structure is thus that maximizes the duration of the alignment with the flow. In other words, the favorite structure retains the highest rotational symmetry and a discrete translational symmetry of $\lambda_{c}$ in two-dimension, i.e., the observed hexagonal pattern. A subtle inference along this line of argument is that the maximum local density fluctuates 6 times during one oscillation period, if the observer travels along with the density wave. Such fluctuations suggest the out-of-alignment and realignment events between the flow and the hexagonal pattern, and this is indeed confirmed in experiments [20].

$k_{c}={2\pi}/\lambda_{c}$ of the observed pattern is measured from the two-dimensional Fourier spectrum of the experimental images, typified by Fig. 4 inset. The measured $k_{c}$ is compared with linearly selected wave number $k_{m}$, corresponding to the peak growth rate Fourier mode in the model. The magnitude of $\eta_{p}$, bounded by $\rho\nu(\phi)=\rho\nu_{0}(\phi_{J}-\phi)^{-2}$, is critical here [20]. From the perspective of linear analysis, $\eta_{p}\sim\rho\nu(\Phi)$ is implemented, and the resultant $k_{m}$ is of $\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-1}$. Nevertheless, as seen in Fig. 1c, the local $\phi$ can be as high as 0.449 for the developed density waves, and the corresponding $\eta_{p}(\phi)$ would shift $k_{m}$ towards longer wavelength. The evaluation of $k_{m}$ with such a high $\eta_{p}=1800~{}$\mathrm{P}\mathrm{a}\cdot\mathrm{s}$$ is in reasonable agreement with $k_{c}$ [Fig. 4]. It again indicates that the description of the fully-developed density pattern is beyond purely linear predictions [39, 21].

Our last remark is about the role of boundary confinement in the growth of density disturbances. Recent advances find that the steady shear-thickening state is a precursor to shear jamming [34, 40, 10, 30], where the dilation of the particle phase under shear is (partially) frustrated by the confining stress [41, 34]. The boundary confinement is thus considered essential. The open system studied so far is confined by the air-suspension interface tension, $\Gamma$, and the curvature of the free surface, $k_{c}^{2}\Delta h$. For the developed density waves, this confining pressure necessarily balances the onset stress of DST, $\Gamma k_{c}^{2}\Delta h\sim\tau^{*}$. Reducing $\Gamma$, e.g., by covering the suspension with a layer of silicone oil, leads to an increase of $k_{c}$ [20]. On the other hand, confining the suspension with an acrylic plate, either comoving or fixed in the lab frame of reference, completely suppresses the long-standing density waves, and the suspension appear uniform for $\omega>\omega_{c}$, though transient fluctuations of $\phi$ and thrust on the plate are present instead, as reported for DST in other configurations [42, 43, 44, 17]. This dramatic contrast indicates that the boundary confinement alters the underlying growth of disturbances. We speculate that only density disturbances below a critical may size persist under rigid confinements. This may be readily seen in experiments under the comoving top-plate confinement. In this configuration, a suspension stays uniform for $\omega>\omega_{c}$ as described above. However, soon after removing the confinement and restarting the oscillation at $\omega<\omega_{c}$, the growth of inhomogeneous density profile can be seen for a finite duration (see the details of the experiment in [20]). It implies that finite-size clusters already exist in the comoving-confinement configuration during shear [45], which experiences a transient growth even for $\omega<\omega_{c}$ after the removal of the rigid confinement.

The significance of boundary confinement highlighted above may not be simply captured by a cutoff of $\Pi$ (indicating the boundary stiffness). As evidenced by recent experiments with spatial resolution [44, 43, 46, 14, 16], the shear-thickening state seems heterogeneous. The increase of the global stress response on shear consists in presence of individual high stress zones. The growth of localized disturbances and the resultant stress response are highly nonlinear with size [47, 48]. In particular, when its dimension is comparable to the gap between boundaries, a sudden rise of stress response is expected [44, 43]. The constitutive relation, calibrated with bulk rheology measurements, only capture these intermittent microscopic events on average and may not be sufficient for a quantitative description of the growth of local disturbances.

This work is partially supported by the National Natural Science Foundation of China (Grant No. 12172277).