2 Model construction

As illustrated in figure 1, while fluid tracers (position $\boldsymbol{x}(t)$ ) move according to $\text{d}x_{i}/\text{d}t=u_{i}(\boldsymbol{x},t)$ , inertial particle trajectories ( $\boldsymbol{y}(t)$ ) evolve following the particle velocity $\boldsymbol{v}(t)$ according to $\text{d}y_{i}/\text{d}t=v_{i}(t)$ , where, in general, $v_{i}(t)\neq u_{i}(\boldsymbol{y},t)$ . When the particle radius $a\ll \unicode[STIX]{x1D702}=\unicode[STIX]{x1D708}^{3/4}\langle \unicode[STIX]{x1D716}\rangle ^{-1/4}$ ( $\unicode[STIX]{x1D702}$ is the Kolmogorov length scale, where $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the surrounding fluid and $\langle \unicode[STIX]{x1D716}\rangle$ is the average dissipation rate of the fluid flow) and $Re_{a}=a|\boldsymbol{v}-\boldsymbol{u}|/\unicode[STIX]{x1D708}\ll 1$ (particle Reynolds number), the dynamical equation of the inertial particle trajectory (Maxey & Riley Reference Maxey and Riley1983) in the absence of gravitational settling can be simplified to (Maxey Reference Maxey1987; Balkovsky et al. Reference Balkovsky, Falkovich and Fouxon2001)

(2.1) $$\begin{eqnarray}\frac{\text{d}v_{i}}{\text{d}t}=\unicode[STIX]{x1D6FD}\frac{\text{D}u_{i}}{\text{D}t}+\frac{u_{i}-v_{i}}{\unicode[STIX]{x1D70F}_{p}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}=3\unicode[STIX]{x1D70C}_{f}/(2\unicode[STIX]{x1D70C}_{p}+\unicode[STIX]{x1D70C}_{f})$ is the added mass parameter and $\unicode[STIX]{x1D70F}_{p}=a^{2}/3\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FD}$ is the relaxation time for the trajectory of a spherical particle of radius $a$ . For small Stokes number based on the Kolmogorov time scale ( $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D708}^{1/2}\langle \unicode[STIX]{x1D716}\rangle ^{-1/2}$ ), $St=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\ll 1$ , Maxey (Reference Maxey1987) constructed a perturbation solution to linear order in $St$ which yields the following approximation:

(2.2) $$\begin{eqnarray}v_{i}=u_{i}-(1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D70F}_{p}\frac{\text{D}u_{i}}{\text{D}t}.\end{eqnarray}$$

This solution admits an interpretation in terms of a particle velocity field, $v_{i}(\boldsymbol{x},t)$ , such that the velocity of a particle at location $y_{i}(t)$ can be approximated by $v_{i}(t)=v_{i}(\boldsymbol{y}(t),t)$ . In this way, the time derivative of the particle can be interpreted as $\text{d}/\text{d}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+v_{k}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x_{k}$ . While at finite Stokes number the particle velocity field could be multi-valued as two particles can have different velocities at the same point, the linear perturbation expansion for $St\ll 1$ gives a single-valued particle velocity field. Since this velocity field has non-zero divergence, it can describe clustering effects due to particle inertia.

Figure 1.

Sketch of the fluid and inertial particle trajectories. In this paper, we consider the time history of fluid velocity gradients, $\unicode[STIX]{x1D608}_{ij}(t)$ , along these trajectories.

In this paper, we consider the evolution of the fluid velocity gradient, $\unicode[STIX]{x1D608}_{ij}=\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}$ , along the particle trajectory, as sketched in figure 1. Considering a particle velocity field $v_{i}(\boldsymbol{x},t)$ , the evolution equation for the velocity gradient can be related to the Lagrangian evolution by $\text{d}\unicode[STIX]{x1D608}_{ij}/\text{d}t=\text{D}\unicode[STIX]{x1D608}_{ij}/\text{D}t+(v_{k}-u_{k})\unicode[STIX]{x2202}\unicode[STIX]{x1D608}_{ij}/\unicode[STIX]{x2202}x_{k}$ , which upon substitution of the gradient of Navier–Stokes yields

(2.3) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D608}_{ij}}{\text{d}t}=-\unicode[STIX]{x1D608}_{ik}\unicode[STIX]{x1D608}_{kj}-\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}x_{i}\unicode[STIX]{x2202}x_{j}}-\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}\unicode[STIX]{x1D608}_{ij}-\frac{\unicode[STIX]{x2202}T_{ijk}}{\unicode[STIX]{x2202}x_{k}},\end{eqnarray}$$

where $p$ is the pressure divided by density and $T_{ijk}$ represents spatial fluxes of velocity gradient due to viscosity and inertial effects according to $T_{ijk}=-\unicode[STIX]{x1D708}\unicode[STIX]{x2202}\unicode[STIX]{x1D608}_{ij}/\unicode[STIX]{x2202}x_{k}+\unicode[STIX]{x1D608}_{ij}(1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D70F}_{p}\text{D}u_{k}/\text{D}t$ . A key step is to evaluate the divergence of the particle velocity field using (2.2) for a divergence-free fluid velocity field (Balkovsky et al. Reference Balkovsky, Falkovich and Fouxon2001), i.e. $\unicode[STIX]{x2202}v_{k}/\unicode[STIX]{x2202}x_{k}=-(1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D70F}_{p}\unicode[STIX]{x1D608}_{k\ell }\unicode[STIX]{x1D608}_{\ell k}$ . The final steps in deriving the new inertial restricted Euler system are, similarly to the classical restricted Euler model, (i) to replace the pressure Hessian $\unicode[STIX]{x2202}_{i}\unicode[STIX]{x2202}_{j}p$ by its isotropic part $\unicode[STIX]{x1D6FB}^{2}p\,(\unicode[STIX]{x1D6FF}_{ij}/3)$ and to invoke the pressure Poisson equation $\unicode[STIX]{x1D6FB}^{2}p=-\unicode[STIX]{x1D608}_{k\ell }\unicode[STIX]{x1D608}_{\ell k}$ and (ii) to neglect any spatial fluxes, i.e. set $T_{ijk}=0$ , where we make the strong assumption of neglecting fluxes due to viscosity and inertia effects.

The resulting system reads as follows:

(2.4) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D608}_{ij}}{\text{d}t}=-\unicode[STIX]{x1D608}_{ik}\unicode[STIX]{x1D608}_{kj}+\frac{1}{3}\unicode[STIX]{x1D608}_{k\ell }\unicode[STIX]{x1D608}_{\ell k}\unicode[STIX]{x1D6FF}_{ij}+(1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D70F}_{p}\unicode[STIX]{x1D608}_{k\ell }\unicode[STIX]{x1D608}_{\ell k}\unicode[STIX]{x1D608}_{ij},\end{eqnarray}$$

thus extending the restricted Euler system of equations to include inertial trajectory effects. The original restricted Euler equation is recovered by considering particles with equal density to the surrounding fluid, $\unicode[STIX]{x1D70C}_{p}=\unicode[STIX]{x1D70C}_{f}$ , hence $\unicode[STIX]{x1D6FD}=1$ . Equation (2.4) thus shows, within the limitations of the restricted Euler assumptions, how the inertia of a particle impacts the rate of change for each of the velocity gradient components it experiences along its trajectory. It is important to note that the effect of viscosity essential to studying inertial particles, Stokes drag, is represented in (2.4), while the less crucial effect of spatial diffusion of fluid velocity gradient is neglected. These simplifications enable us to focus on the terms that can be represented exactly in a low-dimensional dynamical system.

The inertial restricted Euler dynamics given by (2.4) can be projected into the two-dimensional space of tensor invariants $Q$ and $R$ , and yields the following two-dimensional dynamical system:

(2.5a,b ) $$\begin{eqnarray}\frac{\text{d}Q}{\text{d}t}=-3R-\frac{2}{3}\unicode[STIX]{x1D6FC}Q^{2},\quad \frac{\text{d}R}{\text{d}t}=\frac{2}{3}Q^{2}-\unicode[STIX]{x1D6FC}QR,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=6(1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D70F}_{p}$ is the time scale representing inertial effects. The second invariant, $Q=(\unicode[STIX]{x1D6FA}_{ij}\unicode[STIX]{x1D6FA}_{ij}-S_{ij}S_{ij})/2$ , represents the relative balance between local rotation, $\unicode[STIX]{x1D6FA}_{ij}=(\unicode[STIX]{x1D608}_{ij}-\unicode[STIX]{x1D608}_{ji})/2=-\unicode[STIX]{x1D716}_{ijk}\unicode[STIX]{x1D714}_{k}/2$ , and straining, $S_{ij}=(\unicode[STIX]{x1D608}_{ij}+\unicode[STIX]{x1D608}_{ji})/2$ . The third invariant, $R=-S_{ij}S_{jk}S_{ki}/3-\unicode[STIX]{x1D714}_{i}S_{ij}\unicode[STIX]{x1D714}_{j}/4$ , represents the balance between strain production and enstrophy production (Tsinober Reference Tsinober2001).

For particles that are heavier than the surrounding fluid, $0<\unicode[STIX]{x1D6FD}<1$ and $\unicode[STIX]{x1D6FC}>0$ . For particles lighter than the surrounding fluid, $1<\unicode[STIX]{x1D6FD}<3$ and $\unicode[STIX]{x1D6FC}<0$ . For heavy particles ( $\unicode[STIX]{x1D6FC}>0$ ), the inertial term in the evolution equation for $Q$ tends to oppose rotation-dominant states ( $Q>0$ ) and reinforce strain-dominant states ( $Q<0$ ). The exact opposite is true for light particles, where the inertial term opposes highly straining states and favours highly rotating states. In this way, heavy particles cluster in straining regions ( $Q<0$ ) and lighter particles cluster in rotating regions ( $Q>0$ ), qualitatively mimicking well-known preferential concentration trends. In homogeneous turbulence, $\langle Q\rangle =0$ and $\langle R\rangle =0$ , where angle brackets denote ensemble averaging (Betchov Reference Betchov1956), but when averaging over inertial trajectory ensembles, one observes that $\langle Q\rangle <0$ for heavy particles and $\langle Q\rangle >0$ for light particles (Ireland et al. Reference Ireland, Bragg and Collins2016). The qualitative features of $RQ$ space are sketched in figure 2, including this qualitative effect of inertia on $\langle Q\rangle$ , which is valid for any random velocity field, but in (2.5) is combined with turbulence-like dynamics in a manner consistent with first principles.

Figure 2.

Sketch outlining the features of the $RQ$ invariant space, including representative local flow topology cubes. The Vieillefosse tail (dashed line) represents the boundary between real and complex eigenvalues of the velocity gradient tensor.

4 Analysis

Due to its inherent simplicity, many of the features of the inertial restricted Euler system can be investigated analytically. A salient feature of the original restricted Euler equation ( $\unicode[STIX]{x1D6FC}=0$ ) is the invariant $Q^{3}+(27/4)R^{2}$ (Vieillefosse Reference Vieillefosse1982; Cantwell Reference Cantwell1992). For the extended system given by (2.5),

(4.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left(Q^{3}+\frac{27}{4}R^{2}\right)=-2\unicode[STIX]{x1D6FC}Q\left(Q^{3}+\frac{27}{4}R^{2}\right),\end{eqnarray}$$

so that for the particular choice $Q^{3}+(27/4)R^{2}=0$ , this remains an invariant of the dynamics. Strikingly, this means that the Vieillefosse tail, $Q_{v}(R)=-(27/4)^{1/3}R^{2/3}$ , is an invariant manifold for all values of $\unicode[STIX]{x1D6FC}$ , an important observation clearly supported by the DNS evidence in figure 3.

It is straightforward to show that (2.5) has two fixed points, one at the origin and another at $R_{0}=-2/\unicode[STIX]{x1D6FC}^{3}$ , $Q_{0}=-3/\unicode[STIX]{x1D6FC}^{2}$ , which lies on the Vieillefosse tail, i.e. $Q_{0}^{3}+(27/4)R_{0}^{2}=0$ , as clearly seen in figure 3. Linear stability analysis of this fixed point reveals eigenvalues of $\unicode[STIX]{x1D706}_{1}=6/\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D706}_{2}=1/\unicode[STIX]{x1D6FC}$ with (unnormalized) eigenvectors $\boldsymbol{e}^{(1)}=(1,-(3/2)\unicode[STIX]{x1D6FC})^{T}$ and $\boldsymbol{e}^{(2)}=(1,\unicode[STIX]{x1D6FC})^{T}$ . The fixed point is unstable for heavy particles and stable for light particles, as seen in figure 3. The slope of the Vieillefosse tail at the fixed point is $\text{d}Q_{v}/\text{d}R|_{Q_{0}}=\unicode[STIX]{x1D6FC}$ , so that the eigenvector associated with the more weakly stable/unstable eigenvector points along the Vieillefosse manifold.

Along the Vieillefosse manifold, the dynamics is given by $\text{d}R/\text{d}t=(3/2^{1/3})R^{4/3}+(3/2^{2/3})\unicode[STIX]{x1D6FC}R^{5/3}$ , which for $\unicode[STIX]{x1D6FC}\neq 0$ can be written as

(4.2) $$\begin{eqnarray}\frac{\text{d}R}{\text{d}t}=\frac{6}{\unicode[STIX]{x1D6FC}^{4}}\left[\left(\frac{R}{R_{0}}\right)^{4/3}-\left(\frac{R}{R_{0}}\right)^{5/3}\right].\end{eqnarray}$$

This shows the reinforcement of the original finite-time singularity behaviour in the fourth quadrant for heavy particles as well as a similar path to singularity in the third quadrant, for $R<R_{0}<0$ . It also shows that there is no longer a finite-time singularity along the Vieillefosse manifold for light particles due to the stable fixed point. This stable fixed point for light particles highlights the role of inertia in counteracting the tendency towards more heavily strain-dominated regions down the Vieillefosse tail. This tendency, meanwhile, is strengthened for heavy particles.

The linear stability of the Vieillefosse manifold is examined by considering the trajectory $Q(R)=Q_{v}(R)+\unicode[STIX]{x1D716}(R)$ . Using $\text{d}\ln \unicode[STIX]{x1D716}/\text{d}t=\text{d}\ln \unicode[STIX]{x1D716}/\text{d}R\,\text{d}R/\text{d}t$ , the linearized behaviour of $\unicode[STIX]{x1D716}$ can be shown to be

(4.3) $$\begin{eqnarray}\frac{\text{d}\ln \unicode[STIX]{x1D716}}{\text{d}t}=2^{1/3}\unicode[STIX]{x1D6FC}R^{1/3}\left(R^{1/3}-\left(\frac{16}{\unicode[STIX]{x1D6FC}^{3}}\right)^{1/3}\right).\end{eqnarray}$$

When $\text{d}\ln \unicode[STIX]{x1D716}/\text{d}t>0$ , the Vieillefosse line is an unstable manifold. When $\text{d}\ln \unicode[STIX]{x1D716}/\text{d}t<0$ , it is a stable manifold. The stability of the manifold changes sign twice, once at the origin, and also at the point

(4.4) $$\begin{eqnarray}(R_{s},Q_{s})=\left(\frac{16}{\unicode[STIX]{x1D6FC}^{3}},-\frac{12}{\unicode[STIX]{x1D6FC}^{2}}\right).\end{eqnarray}$$

For $\unicode[STIX]{x1D6FC}>0$ (heavy particles), the Vieillefosse manifold is stable for $0<R<R_{s}$ , while it is unstable if $R<0$ or $R>R_{s}$ . Meanwhile, for $\unicode[STIX]{x1D6FC}<0$ (light particles), it is stable for $R<R_{s}$ or $R>0$ and unstable for $R_{s}<R<0$ . The point of neutral stability, $R_{s}$ , has an opposite sign to the fixed point, $R_{0}$ , and is also eight times larger in magnitude. Thus, it is unlikely to be of much relevance to the stationary statistics at low $St$ number, as shown in figure 3.