Mass crystals and cage effect in vorticity crystals

Abstract

We study the motion of tiny heavy inertial particles advected by a two dimensional inviscid fluid flow composed of $N$ identical point vortices regularly placed on a ring, and forming a crystal. In the limit of weak particle inertia, we show asymptotically that, in the reference frame of the crystal, inertial particles have $N$ asymptotically stable equilibrium positions located outside the crystal, in agreement with numerical observations by Ravichandran et al. (Sādhanā 42, 2017). In addition to these ”satellite” attracting points, we observe that for $N\geq 3$ the center of the ring, though degenerate, is a stable equilibrium position for inertial particles. This creates a kind of cage effect, where inclusions slowly drift towards the center under the effect of the surrounding vortices. This cage effect is observed to persist even at larger Stokes numbers, in contrast with the satellite attracting points that vanish when the Stokes number is above some critical value.

Vortex crystals are two-dimensional flow structures appearing in inviscid fluids, where point vortices place themselves naturally at the tips of polygonal structures that rotate as a solid body Aref et al. (2002). They were observed experimentally four decades ago in superfluid Helium Yarmchuk et al. (1979), then in pure electron columns (Fine et al., 1995; Durkin and Fajans, 2000a, b), floating disks (Grzybowski et al., 2000; Grzybowski and Whitesides, 2002), and more recently on the poles of Jupiter Adriani et al. (2018). Both experiments and simulations showed that crystals can emerge from random vortex distributions and persist even in the presence of weak viscous effects Schecter et al. (1999); Jin and Dubin (2000); Siegelman et al. (2022). Shortly after the initial disordered state, vortices advect each other in a chaotic manner, and produce vortex layers creating a diffuse background Schecter et al. (1999). Successive vortex mergers lead to intense vortices which eventually converge towards a stable configuration in rigid rotation.

The influence of this self-organization on the transport of inertial particles suspended in the fluid has received little attention so far. However, it can be conjectured that the emergence of such flow structures will certainly affect the distribution of particles, as these objects are known to be very sensitive to the distribution of vorticity. In strong contrast with turbulent flows however, vortex crystals induce steady flows in the reference frame rotating with the crystal. In the case where two point vortices rotate around each other, their centres forming a rigid segment, it has been shown that heavy inertial particles could be trapped by attracting points at rest in the frame of the segment Angilella (2010); Ravichandran et al. (2014). This property has been observed to persist in the case where more than two vortices interact to form a crystal Ravichandran et al. (2017). The goal of the present work is to address this question theoretically, in the case of heavy inertial particles transported at low particle Reynolds number in an inviscid potential flow composed of identical vortices regularly placed on a ring. It will be shown that equilibrium points always exist in the frame of the crystal for inertial particles, in the limit of weak inertia, and that these points are hyperbolic (i.e. the real parts of the eigenvalues of the particle flow gradient are non-zero there). Some of these points will be shown to be asymptotically stable, so that the particle cloud at long time takes the form of a mass crystal at rest in the reference frame of the vortex crystal. These trapping points being located away from the center of vorticity, they will be referred to as satellite points. In addition, in the present case where no vortex is placed at the center of vorticity, we will show that this point, though degenerate (it is a center point with pure imaginary eigenvalues), is asymptotically stable and can attract a significant number of particles. Inclusions trapped there will then be surrounded by the vortices and maintained within a kind of cage. In the laboratory frame, these particles will appear as almost motionless at long times, while vortices rotate around them. Both the satellite attracting points and the cage effect will be studied asymptotically, in the limit of small Stokes numbers. We will then show some numerical simulations, involving up to seven vortices, which confirm theoretical predictions.

Prior studying the motion of inclusions in vortex crystals, we recall some basic considerations about equilibrium positions of inertial particles (Babiano et al., 2000; Haller and Sapsis, 2008; Cartwright et al., 2010). We assume particles are much heavier than the fluid, do not interact, and have a small particle Reynolds number. Under these conditions, the non-dimensional motion equation of an inertial particle with time-dependent position $\mathbf{X}(t)$, transported in a steady flow with velocity field $\mathbf{u(X)}$, assuming a linear drag and neglecting the added mass and Boussinesq-Basset forces, is of the form:

$$\ddot{\mathbf{X}}=\frac{1}{St}(\mathbf{u}(\mathbf{X})-\dot{\mathbf{X}})+\mathbf{F}(\mathbf{X},\dot{\mathbf{X}})$$ (1)

where $St$ is the Stokes number, that is the non-dimensional response time of the inclusion. In the case of spherical particles with mass $m_{p}$ and radius $r_{p}$, the Stokes number reads $St=\Omega_{0}m_{p}/(6\pi\mu\,r_{p})$, where $\mu$ is the dynamic viscosity of the fluid and $1/\Omega_{0}$ is a characteristic flow time which will be defined below.

The term $\mathbf{F}$ contains various terms corresponding to body forces, hydrodynamic forces or fictitious forces (if the reference frame is non-inertial). Equilibrium positions, if any, correspond to

$$\mathbf{\Phi}(\mathbf{X},St)\equiv\mathbf{u}(\mathbf{X})+St\,\mathbf{F}(\mathbf{X},{\mathbf{0}})={\mathbf{0}}.$$ (2)

We assume that the fluid flow has a stagnation point at some position $\mathbf{X}_{eq}^{0}$: $\mathbf{u}(\mathbf{X}_{eq}^{0})=\mathbf{0}$. Then, by virtue of the implicit functions theorem, if the gradient of $\mathbf{\Phi}$ at $(\mathbf{X}_{eq}^{0},0)$ is invertible, that is if the gradient of $\mathbf{u}$ at $\mathbf{X}_{eq}^{0}$ is invertible, there exists a continuously differentiable function $\mathbf{X}_{eq}(St)$, defined for $St$ in the vicinity of 0, and satisfying $\mathbf{\Phi}(\mathbf{X}_{eq}(St),St)=\mathbf{0}$, and $\mathbf{X}_{eq}(0)=\mathbf{X}_{eq}^{0}$. In other words, a single-valued continuous equilibrium branch persists, inertial particles have equilibrium positions that depend continuously on $St$, and which coincides with $\mathbf{X}_{eq}^{0}$ in the limit of small inertia. This shows that equilibrium positions of inertial particles exist in the vicinity of any fluid equilibrium position with non-zero eigenvalues, in the limit of small Stokes numbers. In the following, we will study the position and stability of these points in the case of a vortex crystal.

We consider the case where the flow is a vortex crystal composed of $N$ identical point vortices located at $\mathbf{x}_{n}=a(\cos(n2\pi/N),\sin(n2\pi/N))$, with $n=0,..,N-1$. We will limit ourselves to the case $N\leq 7$ (Thomson’s heptagon), as a larger number of vortices placed on a ring might lead to unstable flows Mertz (1978). These vortices form a rigid body (a regular polygon with radius $a$) that rotates with respect to the laboratory frame, with a constant angular velocity $\Omega_{0}=(N-1)\,{\Gamma}/{(4\pi a^{2})}$, where $\Gamma$ is the strength of vortices. Without loss of generality we assume $\Gamma>0$, so that the crystal rotates in the counter-clockwise direction ($\Omega_{0}>0$). In the frame rotating with vortices, the flow is steady and its velocity field is given by $\mathbf{u}=(\partial\psi/\partial y,-\partial\psi/\partial x)$, where the streamfunction $\psi$ reads:

$$\psi=\frac{1}{2}\Omega_{0}|\mathbf{x}|^{2}-\frac{\Gamma}{4\pi}\sum_{n}\ln(|\mathbf{x}-\mathbf{x}_{n}|^{2})$$ (3)

with $\mathbf{x}=(x,y)$. The corresponding streamlines for $N=$ 2, 3, 4 and 7 are shown in Fig. 1. In addition to the circular streamlines around vortices, the flow is characterized by $N$ anticyclonic cells rotating around elliptic equilibrium positions of fluid points, that will be denoted $\mathbf{X}_{eq}^{0}$ throughout the paper (the superscript 0 indicating that they correspond to particles with $St=0$).

Figure 1: Streamlines of the flow in the frame rotating with the crystal, when $N$ = 2, 3, 4 and 7. Lengths are set non-dimensional by means of the radius $a$ of the circumscribed circle containing vortices. Vortices have a positive sign and are located at $(\cos(n2\pi/N),\sin(n2\pi/N))$, $n=0,..,N-1$. In each case, the $N$ recirculation cells are anticyclonic. Note that, when $N>2$, an additional anticyclonic cell exists around the center-of-vorticity located at $(x,y)=(0,0)$.

These anticyclonic cells will play a key role in the trapping of inertial particles: the Coriolis force acting on these particles will be directed towards the interior of the cell. Note also that, as soon as $N>2$, another anticyclonic recirculation cell exists around the center-of-vorticity $(0,0)$, which is also an elliptic equilibrium position of fluid points. The equilibrium positions of fluid points are fixed points of the system $\dot{\mathbf{X}}=\mathbf{u}(\mathbf{X})$. Their stability can be studied by means of the eigenvalues $\mu$ of the gradient of $\mathbf{u}$, which satisfy

$$\mu^{2}=(\frac{\partial u}{\partial x})^{2}+(\frac{\partial u}{\partial y})^{2}-2\frac{\partial u}{\partial y}\equiv I(\mathbf{X}).$$ (4)

Here we chose $1/\Omega_{0}$ as time unit, and made use of the fact that the flow is incompressible and irrotational in the laboratory frame, that is $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ and $\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}=-2$. Using both expressions we could re-write the gradient of $\mathbf{u}$ in terms of derivatives of the $x$-velocity only, and obtain its characteristic polynomial (4). Elliptic equilibrium points $\mathbf{X}_{eq}^{0}$ are characterized by pure imaginary eigenvalues, so that $I_{0}\equiv I(\mathbf{X}_{eq}^{0})<0$.

In particular, the determinant of $\nabla\mathbf{u}$ at $\mathbf{X}_{eq}^{0}$ is equal to $|I_{0}|$ and is non-zero. The implicit functions theorem mentioned above implies the existence of the branch $\mathbf{X}_{eq}(St)$, for inertial particles transported in this flow in the limit of small Stokes number $St$. For a given Stokes number, the stability of the branch $\mathbf{X}_{eq}(St)$ emerging from $\mathbf{X}_{eq}^{0}$ can be investigated by linearising the particle motion equation (1), in the vicinity of $(\mathbf{X}_{eq}(St),\dot{\mathbf{X}}=\mathbf{0})$, with $\mathbf{F}=\mathbf{X}-2\hat{\mathbf{z}}\times\dot{\mathbf{X}}$ corresponding to the non-dimensional centrifugal and Coriolis forces ($\hat{\mathbf{z}}$ is the direct unit vector perpendicular to the $(x,y)$ plane). The system is then recast as a dynamical system with 4 degrees-of-freedom $(\mathbf{X},\dot{\mathbf{X}})$. Exploiting the fact that $St\ll 1$, one obtains that the eigenvalues of the dynamical system (1) at the equilibrium point have the following real parts:

$$Re(\lambda_{1}) = Re(\lambda_{3})=-(1+I_{0})\,St,$$ (5)

$$Re(\lambda_{2}) = Re(\lambda_{4})=-\frac{1}{St}+(1+I_{0})\,St$$ (6)

plus terms of order $O(St^{2})$. Hence, if $1+I_{0}\not=0$, and even if the fluid stagnation point is degenerate (non-hyperbolic), the neighbouring equilibrium point of inertial particles is not, since all real parts are non-zero in the limit of small Stokes numbers. It is hyperbolic and its stability can be derived by examining the sign of these real parts. In the limit where $St\ll 1$, the real part of $\lambda_{2}$ and $\lambda_{4}$ is always strictly negative. If $1+I_{0}>0$, the real part of $\lambda_{1}$ and $\lambda_{3}$ is also strictly negative and the equilibrium point $\mathbf{X}_{eq}(St)$ is asymptotically stable and attracts particles released in its basin of attraction. For the four cases of figure 1, the position of the fluid elliptic stagnation points in the anticyclonic cells could easily be found (exactly for $N$=2 and 4), and one gets: $1+I_{0}=1/4$ ($N=2$), $1+I_{0}\simeq 0.121$ ($N=3$), $1+I_{0}=11/4-\sqrt{7}$ ($N=4$), and $1+I_{0}\simeq 0.335$ ($N=7$). These values being positive, we conclude that inertial particles must be trapped in the vicinity of these points, provided their Stokes number is sufficiently small. Since the elliptic points in the anticyclonic cells are all located at $|\mathbf{X}|>1$, that is outside the ring formed by the vortices, the attracting points $\mathbf{X}_{eq}(St)$ will be called satellite attracting points in the following.

In addition to these satellite points, another equilibrium point exists for particles in these flows. Indeed, in the vicinity of the fluid stagnation point located at the center-of-vorticity $(0,0)$, one can check that the non-dimensional streamfunction reads $\psi\simeq\frac{1}{2}r^{2}+O(r^{N})$, leading to $I_{0}=-1$ when $N>2$. Eq. (5) therefore implies that $Re(\lambda_{1})=Re(\lambda_{3})=O(St^{2})$, so that the linear analysis does not allow to conclude about the stability of the equilibrium branch near the center. Particles rotate around the center, while they can drift towards the interior or the exterior of the cell under the effect of the non-linear terms $O(r^{N})$ appearing in the streamfunction. To study the effect of these terms, we consider an inertial particle injected very near $(0,0)$, i.e. $|\mathbf{X}(0)|=\varepsilon\ll 1$ and re-scale the positions by $\varepsilon$, setting $x_{*}=x/\varepsilon$, $y_{*}=y/\varepsilon$. This leads to

$$\psi^{*}=\frac{\psi}{\varepsilon^{2}}=\psi_{0}^{*}+\varepsilon^{N-2}\psi_{N-2}^{*}+O(\varepsilon^{2N-2})$$ (7)

where $\psi_{0}^{*}=(x_{*}^{2}+y_{*}^{2})/2$ and $\psi_{N-2}^{*}(x_{*},y_{*})$ is a polynomial of degree $N$: $\psi_{N-2}^{*}=x_{*}^{3}-3x_{*}y_{*}^{2}$ (if $N=3$), $\psi_{N-2}^{*}=2/3\,{x_{*}}^{4}+2/3\,{y_{*}}^{4}-4\,{x_{*}}^{2}{y_{*}}^{2}$(if $N=4$), etc.

Figure 2: Clouds of particles in the laboratory frame, advected by 3 vortices (first column), 4 vortices (middle column) and 7 vortices (third column), after two turns (upper row) and 12 turns (middle row) of the crystal. White circles are vortex centers, and black dots are particles. The Stokes number of particles is $St=0.062$ for $N=3$ and 4 vortices, and $St=0.021$ for the case $N=7$. The last row shows the initial positions of the particles trapped by satellite points (blue dots) and by the center point (red dots). White squares indicate the position of the degenerate fluid stagnation points $\mathbf{X}_{eq}^{0}$ in anticyclonic cells.

Then we calculate the change of the streamfunction $\Delta\psi^{*}=\psi^{*}(\mathbf{X}(T))-\psi^{*}(\mathbf{X}(0))$ during one period $T=2\pi$ of the rotation of the crystal. Negative values of this quantity will indicate that the particle drifts towards the ring center (as $\psi^{*}$ increases with the distance to $(0,0)$, when $\varepsilon$ is sufficiently small). To achieve this we consider a fixed value of $\varepsilon$ and assume that $\varepsilon^{2N-2}\ll St\ll\varepsilon$. On perturbing Eq. (1) we get, for the re-scaled variables:

$$\dot{\mathbf{X}}_{*}\simeq\mathbf{u}_{*}+St\,(-(\nabla\mathbf{u}_{*})\mathbf{u}_{*}+\mathbf{X}_{*}-2\hat{\mathbf{z}}\times\mathbf{u}_{*})+O(St^{2})$$ (8)

where $\mathbf{u}_{*}=\nabla\psi^{*}\times\hat{\mathbf{z}}$. We then obtain

$$\Delta\psi^{*}=St\,\int_{0}^{T}\nabla\psi_{*}(\mathbf{X}(t))\cdot\dot{\mathbf{X}}_{*}dt$$

$$=St\,\int_{0}^{T}[\nabla\psi^{*}\cdot\mathbf{X}_{*}-2|\nabla\psi^{*}|^{2}-\nabla\psi^{*}\cdot(\nabla\mathbf{u}_{*})\mathbf{u}_{*}]dt.$$ (9)

The first term in the integral (9) is the contribution of the centrifugal force, the second one is the contribution of the Coriolis force, and the last one is the contribution of the drag due to the flow induced by vortices. As expected, the second term is always negative because the Coriolis force acts towards the right-hand-side of $\dot{\mathbf{X}}$ here, whereas the streamfunction $\psi^{*}$ increases towards the left-hand-side of streamlines. These integrals are calculated by approximating the particle trajectory by the fluid point trajectory within integrands: $\mathbf{X}_{*}(t)={\mathbf{X}}_{f}(t)+O(St)$ with $\dot{\mathbf{X}}_{f}\simeq\mathbf{u}_{*0}+\varepsilon^{N-2}\mathbf{u}^{*}_{N-2}$. The solution is approximated in the form ${\mathbf{X}}_{f}=(\cos t,-\sin t)+\varepsilon{\mathbf{X}}_{1}(t)+\varepsilon^{2}{\mathbf{X}}_{2}(t)+etc$, where the first term corresponds to the motion of a fluid point when vortices are at infinity, and ${\mathbf{X}}_{i}(t)$ are to be found by injecting this expansion into the motion equation of fluid points. In the case $N=3$, we get

$$\Delta\psi^{*}=-36\,\pi\,St\,\varepsilon^{2}$$ (10)

plus terms of order $St^{2}$. We conclude that particles released in the vicinity of $\mathbf{X}=\mathbf{0}$, i.e. within the anticyclonic cell around the center, will slowly drift towards the center. The set of vortices maintain the particles trapped in a kind of ”cage”. We note also that the drift of particles towards the center is very slow: back to the original scales (i.e. when lengths are reduced by the radius of the vortex crystal $a$), Eq. (10) reads $\Delta\psi=-36\,\pi\,St\,|\mathbf{X}(0)|^{4}$. Therefore, the change of streamfunction at each period is very small since $|\mathbf{X}(0)|\ll 1$.

These theoretical results suggest that particles, in a system of $N$ vortices located on a ring, can be trapped by $N$ satellite attracting points located outside the ring, and by an additional attracting point located at the center-of-vorticity. To check these results we have performed numerical simulations (in the laboratory frame with coordinates $(X,Y)$) involving $N$ point vortices moving under their mutual influence and initially placed on a ring with unit radius. In addition, 10000 inertial particles have been released within the flow domain, and their positions have been computed by solving Eq. (1). It should be noted that when $N=3$ the critical Stokes number, below which equilibrium points do not exist, can be obtained analytically from Eq. (2). After a few manipulations using polar coordinates we show that the non-dimensional radius $r=|\mathbf{X}_{eq}|$ of equilibrium positions satisfies either $r=0$ or $r^{6}-{6}r^{4}/(1+St^{2})+{9}r^{2}/(1+St^{2})-1=0$. This is a cubic equation in $r^{2}$, the discriminant of which is zero when the Stokes number is equal to the critical value:

$$St_{c}=1/2\,\sqrt{2\,{\frac{\sqrt{48\,A-{2}^{4/3}\,{19}^{2/3}A+280}}{\sqrt{A}}}-4-2\,A}$$

where $A=({24+{2}^{4/3}\,{19}^{2/3}})^{1/2}$, leading to $St_{c}\approx 0.206$ in the three-vortex system, in agreement with Fig. 3 of Ref. Ravichandran et al. (2017) (up to a $2\pi$ factor due to a different definition of the Stokes number). For $N=4$ or more, no analytical expression can be found, and one has to use a numerical algorithm to approach $St_{c}$, leading to $St_{c}\approx 0.15$ ($N=4$) and $St_{c}\approx 0.05$ ($N=7$) (see also Fig. 3 of Ref. Ravichandran et al. (2017)). The first column of Fig. 2 shows the resulting clouds of particles, when $N=3$ and $St=0.3St_{c}$. It confirms the existence of the three satellite attracting points, as well as trapping near the center-of-vorticity. The second and third columns of Fig. 2 correspond to $N=4$ and $N=7$ respectively, with Stokes number chosen arbitrarily, such that $St<St_{c}$. As expected, some particles are trapped by satellite points, or by the center point. Trapping by this point is very slow and long-term simulations (not shown here) would show tiny spots of particles near the satellite points, with a large cloud of particles slowly decaying near the center. The last row of Fig. 2 shows the initial positions of trapped particles, which give a view of the basin of attraction of the satellite points (blue dots), together with the basin of the center point (red dots) for $N=$ 3, 4 and 7. The size of the basin of the central point significantly increases with $N$, as the size of the anticyclonic cell grows with $N$ and most trapped particles come from the interior of this cell. We note that these basins are not fractal, as expected since the dynamics is non-chaotic, but they have a spiral structure near vortices. Therefore, blue and red basins might be highly intricated near vortices.

These analyses show that inertial particles might eventually form a mass crystal moving as a rigid-body with the vorticity crystal. In addition, another attracting point exists near the center of the vortex crystal, in high contrast with the 2-vortex system where this point does not exist Angilella (2010); Ravichandran et al. (2014). These observations could reveal interesting features of particle dynamics in either laboratory of geophysical flows, where vortex crystals have been shown to emerge. Further analyses about particles lighter than the fluid, or with a density of the order of the fluid density, should also reveal interesting particle dynamics. Finally, the case where a vortex is present at the center-of-vorticity, which is known to increase the stability of the crystal, is also of interest for particle transport and will be the next step of this work.