1 Introduction
Small particles embedded in a turbulent flow have interesting behaviours when the particle density,
$\unicode[STIX]{x1D70C}_{p}$
, is different from the density of the surrounding fluid,
$\unicode[STIX]{x1D70C}_{f}$
. For example, within a certain parameter range, heavy particles tend to cluster in regions where the strain rate is higher than the rotation rate (Maxey Reference Maxey1987; Wang & Maxey Reference Wang and Maxey1993a
; Eaton & Fessler Reference Eaton and Fessler1994; Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012), while the opposite is true of lighter particles (Biferale, Scagliarini & Toschi Reference Biferale, Scagliarini and Toschi2010) such as bubbles (Wang & Maxey Reference Wang and Maxey1993b
; Calzavarini et al.
Reference Calzavarini, van den Berg, Toschi and Lohse2008) and oil droplets (Gopalan, Malkiel & Katz Reference Gopalan, Malkiel and Katz2008). This clustering effect (Balkovsky, Falkovich & Fouxon Reference Balkovsky, Falkovich and Fouxon2001; Bec Reference Bec2003) can enhance collision rates (Sundaram & Collins Reference Sundaram and Collins1997; Reade & Collins Reference Reade and Collins2000; Wang, Wexler & Zhou Reference Wang, Wexler and Zhou2000; Falkovich, Fouxon & Stepanov Reference Falkovich, Fouxon and Stepanov2002; Bewley, Saw & Bodenschatz Reference Bewley, Saw and Bodenschatz2013). The rate of clustering can be related to the velocity gradient structure of the surrounding fluid experienced by particles along their trajectories (Maxey Reference Maxey1987; Balkovsky et al.
Reference Balkovsky, Falkovich and Fouxon2001). Other important aspects of multiphase flows in various applications, such as particle rotation and orientation (Pumir & Wilkinson Reference Pumir and Wilkinson2011; Parsa et al.
Reference Parsa, Calzavarini, Toschi and Voth2012; Chevillard & Meneveau Reference Chevillard and Meneveau2013), droplet or bubble deformation (Maffettone & Minale Reference Maffettone and Minale1998; Biferale, Meneveau & Verzicco Reference Biferale, Meneveau and Verzicco2014) and nutrient uptake (Batchelor Reference Batchelor1980; Karp-Boss, Boss & Jumars Reference Karp-Boss, Boss and Jumars1996), similarly depend on the local velocity gradient structure. Much of the recent research on particle evolution in fluid turbulence (Toschi & Bodenschatz Reference Toschi and Bodenschatz2009) is based on direct numerical simulations (DNS). However, the high dimensionality of the Navier–Stokes equations especially for high-Reynolds-number turbulence (Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009; Ireland, Bragg & Collins Reference Ireland, Bragg and Collins2016) complicates basic analysis and the development of physical insights.
Dynamical systems models for the velocity gradient along Lagrangian paths provide an interesting possibility for reducing turbulent dynamics to a low-dimensional representation. Vieillefosse (Reference Vieillefosse1982, Reference Vieillefosse1984) and Cantwell (Reference Cantwell1992) developed and studied the so-called restricted Euler system, which is obtained by taking the spatial gradient of the Navier–Stokes equations and neglecting the viscous and anisotropic pressure Hessian contributions. The model consists of a system of
$3\times 3$
nonlinear coupled ordinary differential equations for velocity gradient tensor elements,
$$\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},\quad i,j=1,2,3,\end{eqnarray}$$
where
$\unicode[STIX]{x1D608}_{ij}=\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}$
is the fluid velocity gradient,
$u_{i}(\boldsymbol{x},t)$
is the velocity field and
$\text{D}/\text{D}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+u_{k}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x_{k}$
represents the Lagrangian time derivative following a fluid element in the flow.
The restricted Euler system can be projected to just two degrees of freedom and was shown to display important features seen in turbulent flows, such as the preferential alignment of the vorticity vector in the direction of the eigenvector associated with the median eigenvalue of the strain rate (Kerr Reference Kerr1985; Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987), negative skewness in longitudinal velocity gradients, as well as the tendency to produce extreme velocity gradient events (Meneveau Reference Meneveau2011), which are clustered along the Vieillefosse tail. Without the neglected unclosed terms, however, the restricted Euler system eventually yields finite-time singularities for almost all initial conditions. The restricted Euler system played an important role in motivating subsequent work on modelling the unclosed terms (Girimaji & Pope Reference Girimaji and Pope1990; Jeong & Girimaji Reference Jeong and Girimaji2003; Chevillard & Meneveau Reference Chevillard and Meneveau2006; Biferale et al. Reference Biferale, Chevillard, Meneveau and Toschi2007; Chevillard et al. Reference Chevillard, Meneveau, Biferale and Toschi2008; Wilczek & Meneveau Reference Wilczek and Meneveau2014; Johnson & Meneveau Reference Johnson and Meneveau2016) and related work on the perceived velocity gradient (Chertkov, Pumir & Shraiman Reference Chertkov, Pumir and Shraiman1999; Pumir, Bodenschatz & Xu Reference Pumir, Bodenschatz and Xu2013) at various scales, which has resulted in models capable of reproducing certain turbulent statistics with good quantitative accuracy, although extension to arbitrarily high Reynolds numbers remains an open challenge (Martins Afonso & Meneveau Reference Martins Afonso and Meneveau2010; Meneveau Reference Meneveau2011).
In this paper, we derive an extension to the restricted Euler system that considers the effect of inertia on the velocity gradient dynamics when following an inertial particle, yet can likewise be projected into a dynamical system with just two degrees of freedom, namely the two tensor invariants
$Q=-(1/2)\text{Tr}(\boldsymbol{A}^{2})$
and
$R=-(1/3)\text{Tr}(\boldsymbol{A}^{3})$
. We explore whether the behaviour in the full
$RQ$
plane observed in DNS (Benzi et al.
Reference Benzi, Biferale, Calzavarini, Lohse and Toschi2009) can be explained by the new model. The model is developed from governing equations in § 2, followed by an exploration of its behaviour compared with DNS in § 3. Model properties such as fixed points and stability are investigated analytically in § 4, before conclusions are drawn in § 5.
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)
$$\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:
$$\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
$$\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:
$$\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:
$$\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.
3 Results
We now pursue a more detailed exploration of inertial effects in Navier–Stokes turbulence according to the restricted Euler model developed in § 2. Figure 3 shows the
$RQ$
phase-space portrait for non-inertial (fluid tracer), heavy and light particles computed numerically from (2.5). Also shown is the stationary joint probability density function (PDF) of
$Q$
and
$R$
computed from DNS at
$Re_{\unicode[STIX]{x1D706}}=185$
(Bec et al.
Reference Bec, Biferale, Lanotte, Scagliarini and Toschi2010). Although statistical stationarity (and hence direct comparison of the joint PDF) cannot be achieved in the system of (2.5) without introducing models for the neglected terms, the qualitative comparison of streamlines with the joint PDF in
$RQ$
space for heavy particles from DNS is informative. In particular, (a) shows the original restricted Euler system (
$\unicode[STIX]{x1D6FC}=0$
), for which trajectories move from left to right along lines of constant
$Q^{3}+(27/4)R^{2}$
, eventually proceeding towards the finite-time singularity in the fourth quadrant (Vieillefosse Reference Vieillefosse1982, Reference Vieillefosse1984; Cantwell Reference Cantwell1992). The sheared teardrop shape in the joint PDF in (b) highlights the dynamical significance of the Vieillefosse tail for the full (Navier–Stokes) dynamics of the velocity gradient tensor (Soria et al.
Reference Soria, Sondergaard, Cantwell, Chong and Perry1994; Blackburn, Mansour & Cantwell Reference Blackburn, Mansour and Cantwell1996; Chong et al.
Reference Chong, Soria, Perry, Chacin, Cantwell and Na1998; Ooi et al.
Reference Ooi, Martin, Soria and Chong1999).
Figure 3. Restricted Euler streamlines (a,c,e) and DNS-computed joint PDF isocontours (b,d,f) for Lagrangian trajectories (a,b), heavy-particle trajectories with
$\unicode[STIX]{x1D6FD}=0$
,
$St=0.3,\unicode[STIX]{x1D6FC}=1.8\sqrt{\unicode[STIX]{x1D708}/\langle \unicode[STIX]{x1D716}\rangle }$
(c,d) and light-particle trajectories with
$\unicode[STIX]{x1D6FD}=3$
,
$St=0.1,\unicode[STIX]{x1D6FC}=-1.2\sqrt{\unicode[STIX]{x1D708}/\langle \unicode[STIX]{x1D716}\rangle }$
(e,f). The time scale
$|\unicode[STIX]{x1D6FC}|$
is used to normalize the axes on the streamline plots, while
$\sqrt{2\unicode[STIX]{x1D708}/\langle \unicode[STIX]{x1D716}\rangle }$
is used to normalize the axes for the DNS results, where
$\langle \unicode[STIX]{x1D716}\rangle$
is the average turbulent dissipation rate from the simulation. The red circles show fixed points of the
$RQ$
dynamics, providing a visual connection between the two normalizations. The DNS data are from a pseudo-spectral simulation performed at
$Re_{\unicode[STIX]{x1D706}}=185$
with a grid resolution of
$512^{3}$
(Bec et al.
Reference Bec, Biferale, Lanotte, Scagliarini and Toschi2010). The PDF isocontours are spaced logarithmically with levels
$10^{z}$
,
$z=1,0,-1,-2,-3,-4$
.
In figure 3(c,d), the inertial restricted Euler phase-space portrait is shown for the case of heavy particles (
$\unicode[STIX]{x1D6FC}>0$
). The finite-time singularity down the Vieillefosse line in the fourth quadrant remains and is strengthened. In addition, some initial conditions proceed to a finite-time singularity down the other branch of the Vieillefosse line in the third quadrant. However, it is a very unstable manifold in the third quadrant, meaning that any noise in the system will prevent particles from proceeding to that singularity. In the first quadrant, the downward ‘flow’ of particles is enhanced while the left-to-right ‘flow’ is suppressed. The DNS results for heavy particles indeed show the tendency down the Vieillefosse tail in the fourth quadrant, as well as reduced probabilities in the upper half (
$Q>0$
).
Finally, the phase-space trajectories for light particles (
$\unicode[STIX]{x1D6FC}<0$
) are shown in figure 3(e,f). The restricted Euler trajectories tend to proceed towards the fixed point in the fourth quadrant. There, a rapid collapse towards the Vieillefosse tail is followed by slower evolution along it towards the fixed point. The restricted Euler dynamics imposes more resistance to (e.g. noise-driven) movement away from the tail than movement along the tail. Trajectories no longer exhibit a finite-time singularity down the Vieillefosse tail. However, some trajectories in the second and third quadrants (e.g.
$R(0)|\unicode[STIX]{x1D6FC}|^{3}<-3.2$
with
$Q(0)|\unicode[STIX]{x1D6FC}|^{2}=0$
) do blow up in finite time with
$Q>0$
along inverted Vieillefosse-like manifolds with
$Q\sim R^{2/3}$
. The joint PDF from DNS data indeed suggests that the Viellefosse tail is still dynamically important for light particles, but that light particles do not tend to reach extreme states as far down the Viellefosse tail compared with neutral and heavy particles, an effect that may be qualitatively linked to the fixed point in the restricted Euler dynamics. The fixed point here takes on a clear physical interpretation, that inertial effects prevent light particles from sampling regions of the flow with more extreme values down the Vieillefosse tail. In general, the lower probabilities in the
$Q<0$
region are offset by higher probabilities in the
$Q>0$
region. Additionally, the upward and left-to-right movement in the first quadrant (towards
$R\gg 0$
) of the inertial restricted Euler streamlines is consistent with the enhanced probabilities observed in the DNS results.
While the qualitative comparisons between streamlines of the inertial restricted Euler system and joint PDFs from DNS are encouraging for both heavy and light particles, quantitative comparison of stationary statistics cannot be accomplished without models for the neglected unclosed terms, as was also the case for the original restricted Euler system (Meneveau Reference Meneveau2011). Besides the pressure Hessian and viscous Laplacian, additional modelling work is probably necessary for the additional terms introduced for inertial trajectories, namely
$\unicode[STIX]{x2202}[(v_{i}-u_{i})\unicode[STIX]{x1D608}_{ij}]/\unicode[STIX]{x2202}x_{k}$
. It is important to note that the finite-time singularities in the original and inertial restricted Euler systems are not physical and exist only because of the absence of the neglected terms.
The
$St$
numbers used in figure 3 are evidently low enough for good qualitative agreement, and increasing accuracy of the linear perturbation solution with further decreasing
$St$
could be further investigated with the careful development of reliable statistical closure schemes for the neglected terms to enable quantitative comparisons. For larger
$St$
(or larger
$|\unicode[STIX]{x1D6FC}|$
), the model predicts that the equilibrium points will move ever closer to the origin, thus increasing the deviations of the joint PDF with those along tracer particles. This can be appreciated by comparing with the DNS of Benzi et al. (Reference Benzi, Biferale, Calzavarini, Lohse and Toschi2009), which were for larger
$St$
than the present results.
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),
$$\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
$$\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
$$\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
$$\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.
5 Conclusions
With only two degrees of freedom, the extension of the restricted Euler system for inertial particle paths yields qualitative agreement with basic trends seen from DNS in the
$RQ$
plane. The trends observed follow directly from first principles, i.e. from the ‘self-stretching’ properties of the nonlinear term in the Navier–Stokes and particle transport equations, whose effects are elucidated here by neglecting all of the ‘non-local’ spatial flux terms. While the restricted Euler system cannot itself offer quantitative predictions in most cases, the qualitative success in representing basic inertial effects suggests that it can be a good starting point for developing more complete models for velocity gradients along inertial particle trajectories for applications such as preferential (fractal) concentration (Bec Reference Bec2003; Bec et al.
Reference Bec, Biferale, Boffetta, Cencini, Musacchio and Toschi2006; Esmaily-Moghadam & Mani Reference Esmaily-Moghadam and Mani2016) of heavy and light anisotropic particles (Parsa et al.
Reference Parsa, Calzavarini, Toschi and Voth2012; Chevillard & Meneveau Reference Chevillard and Meneveau2013) and deformation of liquid droplets (Biferale et al.
Reference Biferale, Meneveau and Verzicco2014) or bubbles.
Acknowledgements
The authors are very grateful to L. Biferale and F. Toschi for making their DNS simulation data available for the plots generated in figure 3. P.J. was supported by a National Science Foundation Graduate Research Fellowship Program under grant no. DGE-1232825. C.M.’s research was made possible by a grant from The Gulf of Mexico Research Initiative.
